生信：TCGA学习（R、RStudio安装与下载、常用语法与常用快捷键）

2024-11-12 09:58:14

前置环境

macOS系统，已安装homebrew且会相关命令。

近期在整理草稿区，所以放出该贴。

R语言、RStudio、R包安装

R语言安装

brew install r

RStudio安装

官网地址：https://posit.co/download/rstudio-desktop/

R包下载

注意R语言环境自带的install语法没有内置对于已安装包的检查（这点真的……太乐了
所以建议任何包的下载都包含在require里面。
下载方式有如下几种，这里以dplyr为例：

# r代码实现官网下载
# 打开RStudio运行如下命令即可：
if (!require("dplyr")) {
  install.packages("dplyr")
}
# Bioconductor下载，注意需要先用r代码下载BiocManager包本身
if (!require("BiocManager")) {
  install.packages("BiocManager")
}
BiocManager::install("dplyr")
# Github下载，注意需要先用r代码下载devtools包本身
if (!require("devtools")) {
  install.packages("devtools")
}
devtools::install_github("dplyr")
# 手动安装R包

其中最推荐官网下载，其次为Bioconductor。
为了下载个R包下载devtools本身太不值了，可以看出来R语言并没有兼容github。
而手动安装R包存在需要同时安装该R包的依赖包的情况，过于麻烦（建议R语言出个好点的包管理工具和环境管理工具）

R语言常用语法和RStudio常用快捷键

R语言常用语法

# 引入R包
library("dplyr")

RStudio常用快捷键

ctrl+enter执行命令
ctrl+shift+c注释单行/多行

新版TCGA相关页面操作

TCGA在24年2月中旬更新过一次，现在页面操作建议先看https://www.bilibili.com/video/BV1b34y1g7RM的p3到p5
注意的是这个网站没有做国际化配置，使用浏览器自带的翻译可能会导致页面报错。

Cohort Builder

我们需要关注的是首先需要选取Cohort Builder这个标签页。将鼠标悬停到该标签上，可以看到提示：

Build and define your custom cohorts using a variety of clinical and biospecimen features.

翻译过来就是

使用各种临床和生物样本功能建立和定义您的定制队列。

也就是选一些筛查条件，选中了Cohort Builder这个标签页下，会出现几个可供筛选的表单：

  • Program（计划）：表示研究计划的名称或标签。
  • Project（项目）：表示一个具体的研究项目，通常是在特定的研究计划下进行的。
  • Disease Type（疾病类型）：表示样本与哪种类型的疾病相关。选项如癌症、心血管疾病、神经系统疾病等。
  • Primary Diagnosis（主要诊断）：表示样本的主要疾病诊断，即医生根据病理学、临床表现和其他检查结果所做的初步诊断。
  • Primary Site（主要部位）：表示样本中主要病变或疾病发生的部位或组织。选项如肺、乳腺、结肠等。
  • Tissue or Organ of Origin（起源组织或器官）：表示样本中疾病起源的具体组织或器官。

一般我们只需要关注前两个即可，即Program和Project。Program我们选择TCGA。Project以TCGA-LUSC为例。然后点击Repository。

Repository

将鼠标悬停到该标签上，可以看到提示：

Browse and download the files associated with your cohort for more sophisticated analysis.

翻译过来就是

浏览和下载与您的队列相关的文件，以进行更复杂的分析。

也就是进一步筛查，选中了Repository这个标签页下，会出现几个可供筛选的表单：

  • Data Category(数据类别)：biospecimen生物胺、clinical临床、copy number variation拷贝数变化、dna methylationDNA甲基化、proteome profiling蛋白质组谱分析、sequencing reads测序读数、simple nucleotide variation简单核苷酸变异、somatic structural variation体细胞结构变异、structural variation结构变化、transcriptome profiling转录组图谱分析。
  • Data Type（数据类型）：Gene Expression Quantification基因表达量化、Isoform Expression Quantification异构体表达的量化、miRNA Expression QuantificationMiRNA表达的量化、Splice Junction Quantification拼接接头的量化

一般Data Category选中transcriptome profiling转录组图谱分析， Data Type选中Gene Expression Quantification基因表达量化。

Cart

类似购物车的概念，将要下载的所有文件都放在了这里面，注意每次下载新的数据时要把老的数据全部清除掉。

一般下载两个文件：

  • Cart文件，存储基因表达的文件压缩包，下载所需比较久，两三百MB要下几分钟，推荐用我之后的方式下载。
  • Metadata文件，存储文件名称和样本名的对应关系的json文件。
    在这里插入图片描述

在这里插入图片描述

TCGA数据处理流程

查看癌症/肿瘤名称

https://blog.****.net/u010608296/article/details/112740418

在TCGA下载数据集

见如上新版TCGA页面操作

文件处理步骤

我们需要对下载下来的Cart文件和Metadata文件进行处理，得到列名为样本名，行名为基因名的文件。

自己写的脚本（R语言）

R语言代码如下：

# 安装jsonlite包，用于处理json相关
if (!require("jsonlite")) {
  install.packages("jsonlite")
}
library(jsonlite)

# 安装 tidyverse 包，它包含了一系列用于数据科学的 R 包，如dplyr、ggplot2 等工具
if (!require("tidyverse")) {
  install.packages("tidyverse")
}
library(tidyverse)


# 设置工作目录到 "F:/35P/3TCGA"或者"F:\\35P\\3TCGA"，即meta.data所在的路径
setwd("/Users/shanshan/生信/test3")

# 读入 meta.data 文件
json <- fromJSON("metadata.cart.2024-05-12.json")

# 如果需要查看 json 对象的内容，可以取消注释以下行
# View(json)

# 获取样本名称及文件名称
# 假设 json 对象中的 associated_entities 字段是一个列表，每个元素都是一个数据框
sample_id <- sapply(json$associated_entities, function(x) x[,1])

# 如果需要查看 sample_id 的前 10 个元素，可以取消注释以下行
# sample_id[1:10]

# 创建一个数据框 file_sample，包含样本 ID 和文件名
file_sample <- data.frame(sample_id, file_name = json$file_name)

# 查看 file_sample 数据框的内容
View(file_sample)

# 获取 counts 文件所在的位置
# 注意：pattern 参数的值应该用双引号包围
count_file <- list.files('gdc_download_20240512_165220.387861/', pattern = "*.tsv", recursive = TRUE)

# 如果需要查看 count_file 的前 10 个元素，可以取消注释以下行
# count_file[1:10]

# 获取每个文件的名称
# 使用 strsplit 函数按照 '/' 分割每个文件路径
count_file_name <- strsplit(count_file, split = '/')
# 使用 sapply 函数提取每个分割结果的最后一个元素作为文件名称
# 注意：这里假设文件名是路径的最后一个部分
count_file_name <- sapply(count_file_name, function(x) x[length(x)])

# 如果需要查看 count_file_name 的前 10 个元素，可以取消注释以下行
# count_file_name[1:10]

# 构建一个空的数据框
# 注意：这里的 matrix 应该是 matrix，ncol 应该是 ncol
matrix <- data.frame(matrix(nrow = 60660, ncol = 0))

# 逐个读取及合并
for (i in 1:length(count_file)) {
  # 拼接文件路径
  # 注意：这里的 paste 应该是 paste0，'//' 应该是 '/'
  path <- paste0('gdc_download_20240512_165220.387861/', count_file[i])
  
  # 读取 Counts 文件
  # 注意：这里的 read.delim 的 fill 参数应该是 fill，header 参数应该是 header
  data <- read.delim(path, fill = TRUE, header = FALSE, row.names = 1)
  
  # 设置列名
  # 注意：这里的 data[2, ] 应该是 data[2, ]，假设第二行包含列名
  colnames(data) <- data[2, ]
  
  # 移除前 6 行
  data <- data[-c(1:6), ]
  
  # 选择特定的列
  # 注意：这里的 data[3]、data[6]、data[7] 应该是 data[, 3]、data[, 6]、data[, 7]
  # data <- data[, 3]  # 选择 unstranded counts
  data <- data[6]  # 选择 tpm_unstranded
  # data <- data[, 7]  # 选择 fpkm_unstranded
  
  # 设置列名为对应的样本 ID
  # 注意：这里的 file_samples 应该是 file_sample，file_samples$file_name 应该是 file_sample$file_name
  colnames(data) <- file_sample$sample_id[which(file_sample$file_name == count_file_name[i])]
  
  # 将数据绑定到 matrix 数据框
  # 注意：这里的 cbind 应该是 cbind，matrix 应该是 matrix
  matrix <- cbind(matrix, data)
}

# 转化为 gene_symbol
# 拼接文件路径
# 注意：这里的 paste0 应该是 paste0，'//' 应该是 '/'
path <- paste0('gdc_download_20240512_165220.387861/', count_file[1])

# 读取数据并转换为矩阵
# 注意：这里的 as.matrix 和 read.delim 应该是正确的函数名
data <- as.matrix(read.delim(path, fill = TRUE, header = FALSE, row.names = 1))

# 提取基因名称
gene_name <- data[-c(1:6), 1]

# 如果需要查看 gene_name 的前 10 个元素，可以取消注释以下行
gene_name[1:10]

# 将基因名称绑定到矩阵
# 注意：这里的 matrix 应该是之前创建的 matrix 变量
matrix <- cbind(gene_name, matrix)

# 提取基因类型
gene_type <- data[-c(1:6), 2]

# 如果需要查看 gene_type 的前 10 个元素，可以取消注释以下行
# gene_type[1:10]

# 将基因类型绑定到矩阵
matrix <- cbind(gene_type, matrix)

# 将 gene_name 列去除重复的基因，保留基因最大表达量结果
# 注意：这里的 aggregate 函数应该使用正确的公式和数据框
matrix <- aggregate(. ~ gene_name, data = matrix, max)

# 查看 gene_name 的分布
table(gene_name)


# 保留 protein_coding 类型的基因
# 注意：这里的 subset 函数应该使用正确的数据框和条件
matrix0 <- subset(matrix, gene_type == "protein_coding")

# 查看 gene_type 的分布
table(gene_type)

# 将 gene_name 列设置为行名，并转换为导出格式
# 注意：这里的 matrix 应该是之前创建的 matrix 变量
rownames(matrix) <- matrix[, 1]
matrix <- matrix[, -c(1, 2)]

# 创建一个新的数据框，包含 ID 和之前的矩阵
# 注意：这里的 data.frame 创建应该是正确的
matrix1 <- data.frame(ID = rownames(matrix), matrix)

# 替换列名中的非法字符
# 注意：这里的 gsub 函数应该使用正确的模式和替换字符串
colnames(matrix1) <- gsub('[.]', '_', colnames(matrix1))

# 导出数据
# 注意：这里的 write.table 函数应该使用正确的文件名和参数
write.table(matrix1, 'TCGA_LUSC_TPM.txt', sep = "\t", quote = FALSE, row.names = FALSE)

# 如果需要导出 count 数据，可以取消注释以下行
# write.table(matrix1, 'TCGA_LUSC_count.txt', sep = "\t", quote = FALSE, row.names = FALSE)

医学概念

TPM和FPKM

TPM和FPKM是两种常用的基因表达量计算方法。它们用于衡量基因在不同样本中的表达水平，并在转录组分析和基因表达研究中广泛使用。

TPM（Transcripts Per Million）：TPM是一种归一化的基因表达量表示方法，它考虑了样本中基因表达的总量，并将其转换为每百万个转录本的表达量。TPM的计算方法可根据基因的Read Counts和基因的长度来进行，以考虑到测序深度和基因长度对表达量的影响。TPM值反映了基因在样本中相对的表达水平，且可比较不同样本中基因的表达量变化[2]。

FPKM（Fragments Per Kilobase of transcript per Million mapped reads）：FPKM是基因表达量的另一种归一化表示方法，它也考虑了样本的测序深度和基因长度对表达量的影响。FPKM是将基因的Read Counts除以基因长度（以千碱基对为单位），然后乘以每百万个测序reads的比例缩放到一定大小的值。FPKM常用于RNA-Seq数据分析中，并且对于单个基因的表达量而言，FPKM值可以表示其相对的表达水平[1]。

总之，TPM和FPKM都是用于衡量基因表达水平的归一化方法，可以解决测序深度和基因长度对表达量的影响。TPM计算基于转录本的表达水平，而FPKM则根据测序reads的映射情况计算。

旧版TCGA相关

·因为现在大多数教程用的都是官网1.0版本的页面，所以建议进入1.0的官网地址：https://portal.gdc.cancer.gov/v1/repository

页面操作

Cases用来筛选。Primary Site意味着疾病或者病人种类。在这里以胰腺癌pancreas为例下载。
如果只有一个选项则默认勾选，即使前面的勾选框没有勾选。

Data Category中可以筛选是转录组还是其他。在这勾选Transcriptome Profiling。
Data Type中可以筛选是基因表达水平还是其他。在这勾选Gene Expression Quantification。
Experimental Strategy可以筛选是RNA序列还是其他。在这勾选Experimental Strategy。
Workflow Type可以筛选数据处理的类型。在这勾选Counts。

然后勾选页面上的Add All Files to Cart。

编号

在这里插入图片描述

主要看Sample样品编号。01-09意味着癌症组织。10-19意味着正常组织。20-29意味着正常人的组织。

在这里插入图片描述
Counts意味着无矫正。FPKM意味着用FPKM进行矫正。FPKM-UQ意味着用两种方法：FPKM和UQ进行矫正。
一般下载Counts即可。

这五个都是要下载下来的。其中

  • Clinical意味着临床数据。下载JSON格式的。
  • Biospecimen意味着生物多样性数据。下载JSON格式的。
  • File Metadata用于TCGA编号文件和病人的对应关系。
  • Manifest意味着样本信息，即要下载的TCGA数据文件的详细信息，包括文件名、大小、类型、数据格式等。
  • Cart意味着基因文件。下载是需要最久的。
上一篇：HTML查缺补漏-一、熟记