零基础入门转录组分析——数据处置处罚（TCGA数据库） ...

零基础入门转录组分析——数据处置处罚（TCGA数据库）


  

---

---

TCGA应该是肿瘤数据最权势巨子的泉源之一，但是从TCGA上下载数据集相对来说比较麻烦，因此出现了许多针对TCGA数据进行二次开发的衍生网站，XENA.UCSC就是很直观强盛的一个在线网站，里面收录了众多数据库的数据集，其中就包罗了TCGA数据集，并且是整合好的，可以直接用于分析。
并且XENA.UCSC这个网站不但能下载数据，还能进行在线分析，例如：KM分析，表达量分析等，详细环境可以参考好用的TCGA分析工具：UCSC Xena，但是在这篇教程中仅介绍如何从UCSC上下载所需要的TCGA数据集，并且下载后在R中对数据集进一步处置处罚成后续分析所要的形式。

---

---

1. 本项目以非小细胞肺癌（non-small cell lung cancer，NSCLC）中的肺腺癌（lung adenocarcinoma，LUAD）作为展示
3. 选用的数据库是TCGA。
5. 物种：人类（Homo sapiens）
7. 实验分组：疾病组，对照组。
9. 我这里使用的R版本是4.2.2
11. 要用到的R包：tidyverse, rtracklayer, dplyr

复制代码

1. 数据集获取

起首辈入XENA.UCSC官网（如下图所示），点击箭头指向的位置进一步运行Xena。

进一步点击DATA SETS进入Xena中存储数据集的页面

如下图所示为Xena中存储数据集的页面，其中（1）就是收录的各种疾病的数据集，（2）是筛选条件，由于该网站不但收录了TCGA数据库的，同样还收录了其他数据库的数据集，因此可以根据本身的需要选择想要的数据库及数据集。我们想要下载是TCGA数据集，因此勾选GDC Hub 和TCGA Hub（至于为什么会选择这两个，在下面会有个人的一点理解和表明）

如下图所示，为勾选GDC Hub 和TCGA Hub的筛选环境，可以从图中看到两个都是TCGA数据集，那么区别是什么呢？

直接以本教程要展示的疾病——肺腺癌（lung adenocarcinoma，LUAD）为例，分别找到对应的GDC-TCGA数据集和TCGA数据集，如下所示，分别点进去可以看到详细信息。

如下图所示为GDC-TCGA肺腺癌数据集，可以看到整个界面干净简便，并且图中标注出了4个红字部门就是我们后续要用到的数据，分别是原始Count，FPKM，Phenotype，以及survival data。

4个数据表明如下：

---

• counts——转录组的原始表达矩阵，里面对应的基因表达量又被称作raw_count，行名为基因symbol，列名为样本名（也是病人的id，可以将每一列看作一个病人）
  
• fpkm——raw_count经过转换后的表达矩阵，其盘算公式可参考一文相识Count、FPKM、RPKM、TPM | 相互间的转化 | 收藏教程
  
• phenotype——病人的临床信息，包含分组信息，肿瘤分期，分级，年龄，性别等
  
• survival——病人的生存信息，里面通常会有4列信息，两列的病人的id，另外两列：OS——生存状态（0表现存活，1表现死亡），OS.time——生存时间
  

---

如下图所示为TCGA肺腺癌数据集，可以看到整个界面信息比较多，但也包罗了基因表达量及患者的临床信息和生存信息。

那么这两者区别是什么？可参考别人的教程：UCSC Xena数据库中GDC TCGA数据和TCGA数据的区别和聊UCSC xena的数据下载问题，但是表明的都不是很清楚，并且官方表明也是艰涩难明，对于后续分析该选择哪个并没有做出很好的解答。
总的来说就是TCGA这个数据集的发行年份更早，而GDC则是发行年份更晚一些，做科研的小同伴也都知道一般都要追前沿，那我们这也不例外，选择一个最新的GDC数据库。
同时于偶然间发现了一件事：GDC的数据处置处罚要比TCGA更加方便！如下图所示，分别点进去GDC中的Counts和TCGA中的IlluminaHiSeq。

结果如下所示是GDC的数据，其中（1）是基因的ensemble_ID，这个玩意可以直接转换成基因symbol，（2）表现该数据集GDC数据库是进行了log2（X+1）处置处罚的，因此后续处置处罚是要反转一下才气得到想要的count数据，（3）是下载的链接，只要点击即可下载数据。

TCGA的数据如下图所示，其中（1）对应的不知道是什么，推测应该是探针的ID（没深入研究），也有大概是发行年份较早有些基因缺失。但是这个玩意要想转换成基因symbol，应该需要先找到对应的文件，才气转换成对应基因symbol，同时细心的小同伴也可以发现TCGA这个数据库他只有20531个基因，而前面GDC数据库则是60489个基因，（2）同样表现该数据集是进行了log2（X+1）处置处罚的，因此后续处置处罚也是要反转一下，（3）是下载的链接。

总的来说已目前的理解来看GDC数据库中的数据集更加方便简便一些，可以直接通过ensemble_ID对应基因，因此选用一个更加方便的方法——GDC数据库，前面我们也说了图中（3）是下载链接；“照猫画虎，依葫芦画瓢”，以同样的方式下载FPKM，Phenotype，以及survival data，下载完后注意后缀是不是XX.tsv.gz或XX.tsv。如下图所示是下载好的数据集们，接下来要在R中让他们展示自我并做进一步处置处罚。

2. 数据处置处罚（Rstudio）

1. rm(list = ls()) # 删除工作空间中所有的对象
2. setwd('/XX/XX/XX') # 设置工作路径
3. if(!dir.exists('./00_rawdata')){
4.   dir.create('./00_rawdata')
5. } # 判断该工作路径下是否存在名为00_rawdata的文件夹，如果不存在则创建，如果存在则pass
6. setwd('./00_rawdata/') # 设置路径到刚才新建的00_rawdata下

复制代码

上传前面下载的4个tsv数据到Rstudio中（注意：最好上传到刚才创建的00_rawdata文件夹下），如下图所示，绿色的框就是要处置处罚的原始文件。

起首加载要用的包，全能的tidyverse，dplyr，rtracklayer如果没有安装这3个包，可以通过install.packages(‘XXX’)指令安装，XXX就是包的名字，例如：install.packages('tidyverse')

1. library(tidyverse)
2. library(dplyr)
3. library(rtracklayer)

复制代码

1. ### 读入下载的数据
2. tcga_count <- read_tsv(file = './TCGA-LUAD.htseq_counts.tsv.gz') # count数据
3. tcga_fpkm <- read_tsv(file = "./TCGA-LUAD.htseq_fpkm.tsv.gz") # fpkm数据
4. tcga_survival <- read_tsv(file = "./TCGA-LUAD.survival.tsv") # 患者生存状态
5. tcga_phenotype <- read_tsv(file = "./TCGA-LUAD.GDC_phenotype.tsv.gz") # 患者临床信息
6. ### count数据处理
7. tcga_count <- as.data.frame(tcga_count) # 将导入的文件转成数据框格式
8. tcga_count <- tcga_count[!duplicated(tcga_count$Ensembl_ID), ]
9. tcga_count <- column_to_rownames(tcga_count, var = 'Ensembl_ID') # 将文件第一列转为行名
10. tcga_count <- 2^tcga_count-1 # xena下载的数据经过了log2+1转化，需要将其还原

复制代码

点击右上角的tcga_count可以查看数据集，这里要注意一下：要看一下count数据集中的原始count是不是都是整数，如果不是整数，可以看一下是否经过反转（2^tcga_count-1）。

接下来要导入一个比较紧张的文件——gencode.v22.annotation.gtf
从bing上搜刮gencode.v22.annotation，检索到的第一个就是基因解释文件

进入ENCODE的gencode.v22.annotation页面，直接点击download可下载。

下载后是个.gz的压缩包，解压后同样上传到Rstudio的00_rawdata中，运行如下代码导入：

1. gene_annotation <- import('./gencode.v22.annotation.gtf.gz')
2. gene_annotation <- as.data.frame(gene_annotation)#将文件转换为数据框格式
3. gene_annotation <- gene_annotation [gene_annotation$type == 'gene', ] # 筛选为gene的类型
4. gene_annotation <- dplyr::select(gene_annotation, gene_id, gene_name, seqnames, start, end, width, strand, gene_type)

复制代码

看看最后处置处罚好的gene_annotation长啥样，点击Rstudio右上角的gene_annotation 如下图所示：

• gene_id——基因的ensemble_ID
  
• gene_name——基因名
  
• seqnames——基因位于哪条染色体
  
• start——该基因于染色体上的起始位置
  
• end——该基因于染色体上的停止位置
  
• width——该基因的长度（注意：依附这个就可以盘算FPKM）
  
• strand——链的方向，+表现正链，-表现负链（关于链的方向的描述信息可以自行查找资料，在这里跟教程无关不做过多介绍）
  
• gene_type——基因的类型，有的是miRNA，有的是lincRNA…
  

这里我们用到的是前两列gene_id和gene_name

1. gene_symbol <- gene_annotation[,(1:2)] # 筛选gene_annotation文件中的第一列和第二列
2. colnames(gene_symbol) <- c("ID", "symbol") # 将gene_symbol文件的列名改成ID和symbol

复制代码

gene_symbol如下图所示，第一列是ensemble_ID，第二列是基因symbol

预备好前面的count数据和基因解释文件后，运行下面代码，每句代码后都有解释信息，可以挨个查看，比对处置处罚前后的变化

1. data_tcga <- tcga_count
2. data_tcga$ID <- rownames(data_tcga) # 给data_tcga添加一列，列名为ID
3. data_tcga$ID <- as.character(data_tcga$ID) # 将data_tcga文件中ID列从数值类型数据转化为字符串类型数据
4. gene_symbol$ID <- as.character(gene_symbol$ID) # 将gene_symbol文件中ID列从数值类型数据转化为字符串类型数据
6. data_tcga <- inner_join(data_tcga, gene_symbol, by = "ID") # 将data_tcga文件和gene_symbol文件根据ID列进行合并（这样就能获得ID对应的基因名，但是都排在最后）
7. data_tcga <- dplyr::select(data_tcga, -ID) # 删除ID列
8. data_tcga <- dplyr::select(data_tcga, symbol, everything()) # 重新排列,将最后一列的基因名放到第一列（如果不加everything只会选择symbol一列）
9. data_tcga <- mutate(data_tcga, rowMean=rowMeans(data_tcga[grep('TCGA',names(data_tcga))])) # 添加一列为表达量的平均值
10. data_tcga <- arrange(data_tcga, desc(rowMean)) # 把表达量的平均值从大到小排序
11. data_tcga <- distinct(data_tcga, symbol, .keep_all = T) # distinct函数被用于去重，.keep_all参数表示去重后返回数据框的所有列向量
12. data_tcga <- dplyr::select(data_tcga, -rowMean) # 去除rowMean这一列
13. data_tcga <- tibble::column_to_rownames(data_tcga, var = "symbol")   ## 把第一列变成行名并删除
14. dim(data_tcga)

复制代码

上面的代码说白了就是一个去重加基因解释，结果如下图所示，行为基因symbol，列为样本名，但是同样每个样本后有个01A和11A这个是与癌症相关的，01-09为肿瘤，10-19为正常对照。

接下来根据样本最后的这个数字开端筛选癌症和对还是本

1. ### 筛选癌症和对照样本。01-09为肿瘤，10-19为正常对照
2. mete <- data.frame(colnames(data_tcga))
3. for (i in 1:length(mete[, 1])) {
4.   num <- as.numeric(as.character(substring(mete[i, 1], 14, 15))) # 提取每行第一列中第14，15字符
5.   if(num %in% seq(1, 9)){mete[i, 2] = "T"} # 判断：如果提取的数字在0-9之内就在每行第二列加上T表示肿瘤样本
6.   if(num %in% seq(10, 29)){mete[i, 2] = "N"} # 判断：如果提取的数字在10-29之内就在每行第二列加上N表示正常对照
7. }
8. colnames(mete) <- c("id", "group")
9. table(mete$group)
10. mete$group <- as.factor(mete$group) # 将mete中group列转为因子模式
11. mete <- subset(mete, mete$group == "T") # subset函数，从数据框中选择出group组为T的行
13. exp_tumor <- data_tcga[, which(colnames(data_tcga) %in% mete$id)] # 从大样本中选出组为T的样本
14. exp_tumor <- as.data.frame(exp_tumor)
15. exp_tumor <- exp_tumor[, colnames(exp_tumor) %in% tcga_survival$sample] # 进一步筛选具有生存信息的样本
17. exp_control <- data_tcga[, which(!colnames(data_tcga) %in% mete$id)] # 反向从大样本中选出组为N的样本
18. exp_control <- as.data.frame(exp_control)
19. exp_control <- exp_control[, colnames(exp_control) %in% tcga_survival$sample] # 进一步筛选具有生存信息的样本
21. data_final <- cbind(exp_control, exp_tumor) # 将肿瘤样本文件和正常样本文件合并
22. data_final <- na.omit(data_final) # 去除含有NA的行

复制代码

固然已经根据样本ID可以区分肿瘤样本和正常样本，但是为了确保分组可靠，接下来依据前面导入的临床信息（tcga_phenotype）进一步确认分组

1. table(tcga_phenotype$sample_type.samples) # 查看样本类型，一共四种：FFPE Scrolls，Primary Tumor，Recurrent Tumor，Solid Tissue Normal，这里面我们选择Primary Tumor——原发肿瘤和Solid Tissue Normal——正常组织
2. table(tcga_phenotype$primary_site) # 查看癌症发生的部位
4. clinical_data <- subset(tcga_phenotype, sample_type.samples == 'Primary Tumor' | sample_type.samples == 'Solid Tissue Normal')
5. clinical_data <- clinical_data[clinical_data$submitter_id.samples %in% colnames(data_final), ]
6. clinical_data <- clinical_data[order(clinical_data$sample_type.samples, decreasing = T), ]
8. Group <- data.frame(sample = clinical_data$submitter_id.samples,
9.                     group = clinical_data$sample_type.samples)
10. Group$group <- ifelse(Group$group == 'Solid Tissue Normal', 'control', 'disease')
11. data_final <- subset(data_final, select = Group$sample) # 确保前面整理的count数据与整理的分组信息的样本是一致的
13. table(Group$group)   ##control:59   disease:511

复制代码

如下图所示是整理好的分组信息，第一列是样本的ID，需要和前面整理好的count数据集中的样本ID对应，第二列为样本所属的分组

目前我们已经得到了整理好的count数据以及样本的分组信息，目前还差样本对应的生存信息和fpkm表达矩阵，样本的生存信息比较利益置处罚，如下图所示为导入的样本生存信息，第一列和第三列都是样本ID，第二列OS是患者的生存状态，0表现存活，1表现死亡，第四列OS.time是患者的生存时间。

此时只需要去撤除第三行并且让生存信息表中患者的ID和前面整理好的分组信息表中的患者ID做个匹配即可，最后将整理好的count，group，以及生存信息表生存成csv格式即可。

1. tcga_survival <- tcga_survival[, -3]
2. tcga_survival <- tcga_survival[tcga_survival$sample %in% Group$sample, ]
4. write.csv(tcga_survival, './data_survival.csv')
5. write.csv(clinical_data, './data_clinical.csv')
6. write.csv(data_final, file = './data_count.csv')
7. write.csv(Group, file = './data_group.csv')

复制代码

最后一步就是处置处罚fpkm数据，这是因为在后续分析中除了差异表达分析会用到count数据，别的许多环境用的都是fpkm数据。处置处罚的代码如下，实在团体思路和前面处置处罚count时差不多，最后将结果生存为.csv形式即可。

1. ### fpkm数据整理
2. tcga_fpkm <- as.data.frame(tcga_fpkm)
3. tcga_fpkm <- tcga_fpkm[!duplicated(tcga_fpkm$Ensembl_ID), ]
4. tcga_fpkm <- column_to_rownames(tcga_fpkm, var = 'Ensembl_ID')
5. tcga_fpkm <- 2^tcga_fpkm-1
7. dat_fpkm <- tcga_fpkm
8. dat_fpkm$ID <- rownames(dat_fpkm)
9. dat_fpkm$ID <- as.character(dat_fpkm$ID)
11. gene_symbol$ID <- as.character(gene_symbol$ID)
13. dat_fpkm<-dat_fpkm %>%
14.   inner_join(gene_symbol,by='ID')%>%
15.   select(-ID)%>%     ## 去除多余信息
16.   select(symbol,everything())%>%     ## 重新排列
17.   mutate(rowMean=rowMeans(.[grep('TCGA',names(.))]))%>%    ## 求出平均数
18.   arrange(desc(rowMean))%>%       ## 把表达量的平均值从大到小排序
19.   distinct(symbol,.keep_all = T)%>%      ## symbol留下第一个
20.   select(-rowMean)%>%     ## 反向选择去除rowMean这一列
21.   tibble::column_to_rownames(colnames(.)[1])   ## 把第一列变成行名并删除
22. dim(dat_fpkm)
23. dat_fpkm <- dat_fpkm[rownames(data_final), ]
24. dat_fpkm <- dat_fpkm[,colnames(data_final)]
25. dat_fpkm <- na.omit(dat_fpkm)
26. dat_fpkm <- log2(dat_fpkm + 1)
28. write.csv(dat_fpkm, file = './data_fpkm.csv')

复制代码

到目前为止我们一共得到了5个文件，count，临床信息，分组信息，生存信息，以及最后的fpkm。
注：做到这一步，可以说是生信分析已经完成一半了，后续全部的分析都会基于前面处置处罚的这几个文件，因此这最初的第一步也是最紧张的一步，只有基础打好，后面分析才会可靠一些。

---

结语：
以上就是TCGA数据下载及处置处罚的全部过程，如果有什么需要增补或不懂的地方，大家可以私聊我大概在下方评论。
如果觉得本教程对你有所帮助，点赞关注不迷路！！！
配套资源链接：配套资源（代码+原始数据+处置处罚好的数据）

---

下期预告： 零基础入门转录组分析——数据处置处罚（GEO数据库）

---

• 目录部门跳转链接：零基础入门生信数据分析——导读
  

免责声明：如果侵犯了您的权益，请联系站长，我们会及时删除侵权内容，谢谢合作！更多信息从访问主页：qidao123.com:ToB企服之家，中国第一个企服评测及商务社交产业平台。