从TCGA数据库下载并整合清洗高通量肿瘤表达谱-临床性状数据

TCGA是The Cancer Genome Atlas的简称，由美国NIH管理，是目前出现频率最高，最权威的肿瘤公共数据库。其中收录了人体多个部位肿瘤的高通量测序/基因芯片mRNA、miRNA表达谱数据、基因组变异、甲基化等数据，以及相应样本的人口学特征及临床数据（包括暴露情况、预后、肿瘤分期、治疗手段）。为肿瘤研究的数据挖掘和课题前期分析提供了丰富的资源。

目前课题非常常见的手（tao）段（lu）就是使用公共数据库挖掘肿瘤高通量转录组数据，发现对抑癌/促癌具有潜在价值的基因或基因模块，随后结合功能实验进行验证和深入机制的探讨。所以对没有条件自己送临床样本/动物实验测序的人来说，TCGA是妥妥的宝藏~

就从TCGA RNA-seq表达谱数据下载开始吧。写下这个备忘录，以免以后又出现代码不用就忘记的尴尬。

TCGA官方网站为：https://portal.gdc.cancer.gov

1.点击Repository框。

TCGA官网长这样~

2. 接下来就进入数据类型选择的页面。

选择下载转录组二代测序数据（芯片已经基本成为过去式了），在左侧的File中需要选择Data Category > Transcriptome Profiling, Data Type > Gene Expression Quantification, Experimental Strategy > RNA-seq.

Workflow type >: 关于检测gene expression的RNA-seq的数据类型，TCGA提供已经将reads比对好并注释后的Counts （未经标准化的aligned reads数）, FPKM（已标准化fragments per kilobase per million）, FPKM-UQ（FPKM-upper quartile）下载。

如何选择：（1）如果需要做差异表达分析，则需下载HTSeq-Counts，数据整合后使用R的EdgeR或DESeq2包分析，它们均需输入未经标准化的原始counts matrix。（2）如需做WGCNA模组分析，推荐下载HTSeq-Counts，经DESeq2包标准化后输入分析。（3）如需做CIBERSORT数字免疫细胞分析，需下载HTSeq-FPKM作为输入数据。（4）...

这里以FPKM数据下载为例，但Counts数据的流程也完全一样！

File选项一览

3. 选择需要分析的Case即肿瘤类型和临床样本特征。

将左上角的File调到Case，出现Primary site, Program, Gender, Age等信息需要勾选。需根据课题需要进行相关人口学特征的勾选。

注：为避免后续数据处理中ID转换麻烦，Program需选择TCGA（具有统一的ID特征，可批量处理），project和disease type根据需要的疾病类型勾选。如需做生存分析，Vital status不要勾选某一特定生存状态。

本次以TCGA项目中特定人种 (white) 的肺鳞癌 (LUSC, Lung squamous cell carcinoma) 为例。勾选条件如下：

Case选项一览

4. 添加数据到购物车(Cart)。

此时样本类型和特征都已经选好。点击页面中上的Add all files to cart，然后点击右上角Cart进入购物车下载。Cart上显示已经纳入了392个样本等待下载。

购物车安排上

5. 下载样本的测序数据，临床数据和样本ID信息。

点击Download > Manifest，该文件随后将作为下载器的输入信息文件下载测序数据。（不推荐通过下载Cart直接获取测序数据，因为数据库网络速度慢，遇到断网或无响应就不能断点续传了）。然后同时需点击下载Metadata以及Clinical > json。

从网站下载数据文件

注：第一次使用TCGA，需要点击下载右上角的GDC Data Transfer Tool，根据自己的平台系统选择合适的zip包下载。解压后放入用于储存TCGA下载数据的新建文件夹中。

gdc下载器的下载页面 此时文件夹中的文件已经准备好了！

6.测序FPKM表达矩阵下载。

注：本例使用苹果macbook OS系统操作。我没有windows电脑，所以windows的操作我搜一搜，再补上！

打开terminal终端，把目录转换到储存上述文件所在的文件夹（根据自己的文件夹路径进行更改），并以下载的manifest文件作为输入，使用gdc编译器下载数据：

cd /Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial
./gdc-client download -m gdc_manifest_20200308_044238.txt

接下来就是等待...

Tips：国内最好在白天下载，米国人睡了，传输速度会快一丢丢...但总体来说还是非常感人。而且建议挂学校的学术梯子，有条件挂米国学校的梯子是最好了。

仅作参考，392个样本，国内下午，我用国内普通网络下载了3个小时；挂了Yale的梯子下载同一批数据，只用了15分钟。:)))

terminal数据下载-下载完成展示

7.数据全部下载完成：

数据下载完成后terminal报告success. 文件夹中会出现392个名字神似乱码（然而不是）的新文件夹，每个新文件夹中都包含有.FPKM.txt.gz（一种压缩包形式），这是我们需要的表达矩阵。我们随后使用R将每个单样本的FPKM表达矩阵进行提取，整合，并与临床数据进行匹配。

下载完成后的文件夹及所需文件

接下来的部分使用Rstudio完成：

需要R包：dplyr, R.utils, tidyr, jsonlite, rtracklayer，第一次使用需要先使用install.packages()下载安装。

install.packages('dplyr')
install.packages('R.utils')
install.packages('tidyr')
install.packages('jsonlite')
#rtracklayer is based on bioconductor, so we need to change the way of installation.
source("https://bioconductor.org/biocLite.R")
biocLite("rtracklayer")

1. 设置自定义工作路径，加载数据，解压，并将解压后.txt表达矩阵统一放入新文件夹：

TIPS：gunzip函数是会默认删除压缩包文件，仅保留解压文件的。数据文件下载非常不容易，所以衷心建议各位在批量解压前先整体备份一盘...如果因为匹配pattern或者各种问题导致代码没有循环完成，操作就比较麻烦了！建议找出问题后直接用备份文件重来即可~

library('R.utils')
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial')
#to make sure this directory will be in the first place when performing the for loop.
dir.create('00000_extracted_FPKM')
#The char length of the name of each downloaded expr matrix directory is 36, which is different from other files downloaded from the TCGA website.
dir_FPKM <- dir()[nchar(dir()) == 36]
#decompress the separate files and integrate them into a expression matrix.
for (data in dir_FPKM) {
  filename <- list.files(data, pattern = ".*FPKM")
  R.utils::gunzip(paste0(data,"/",filename))
  file_extracted <- list.files(data, pattern = ".*FPKM.txt$")
  file.copy(paste0(data,"/",file_extracted),"00000_extracted_FPKM")
}

2. 根据每个单独表达矩阵文件的gene_id，将所有样本合并为一个总表达矩阵。

每个样本的FPKM.txt由两列组成，第一列为各基因的ENSEMBL GENE ID，第二列为FPKM表达值，不含标题行。

library('dplyr')
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial/00000_extracted_FPKM')
file_extracted_FPKM <- list.files()
first_file <- read.table(file_extracted_FPKM[1], header = F, sep = '\t', stringsAsFactors = F)
names(first_file) <- c('gene_id', substr(file_extracted_FPKM[1], 1, nchar(file_extracted_FPKM[1])-9))
for (extracted_FPKM in file_extracted_FPKM[2:length(file_extracted_FPKM)]) {
  file_appended <- read.table(extracted_FPKM, header = F, sep = '\t', stringsAsFactors = F)
  names(file_appended) <- c('gene_id', substr(extracted_FPKM, 1, nchar(file_extracted_FPKM[1])-9))  #delete the suffix .FPKM.txt which is made of 9 chars and just maintain the original file names of each sample for following ID conversion.
  first_file <- dplyr::inner_join(first_file, file_appended, by='gene_id')
}
expr_FPKM <- first_file

3. 暂停点：此时已获得最原始的总体FPKM表达矩阵，第一列为各基因的ENSEMBL GENE ID，第一行为样本名称（为txt文件中的file_names），可将其输出备用：

#for original data backup:
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial/00000_extracted_FPKM')
write.csv(expr_FPKM, 'expr_FPKM.csv', quote = F, row.names = F)

若需将原始表达矩阵重新载入：

注：由于由36个字符组成的file_names很多以数字开头，而R不允许这样的命名方式，所以载入后部分样本名称前会自动加一个字符X。所以需要去除。这里合并操作：

#to load the expr_FKPM data exported before and remove the char 'X' automatically added in front of colnames starting with number:
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial/00000_extracted_FPKM')
expr_FPKM <- read.csv('expr_FPKM.csv', header = T, stringsAsFactors = F)
for_replacement_indices <- grep(names(expr_FPKM), pattern = '^X')
names(expr_FPKM)[for_replacement_indices] <- substr(names(expr_FPKM)[for_replacement_indices], 2, 37)

4. 从网站下载的metadata文件中获得各样本的TCGA_IDs，并将合并表达矩阵的样本名转换为完整的TCGA ID (eg: TCGA-33-4538-01A-01R-1201-07)。

#get the TCGA IDs of each sample.
library('jsonlite')
setwd('/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial')
metadata <- jsonlite::fromJSON(txt = "metadata.cart.2020-03-08.json")
metadata_id <- dplyr::select(.data = metadata, c(file_name, associated_entities))
file_name <- c()
TCGA_ID <- c()
sub_TCGA_ID_func <- function(x){
  x$entity_submitter_id
}
ID_conversion_table <- data.frame('file_names'=substr(metadata_id$file_name, 1, 36), 
                                  'TCGA_IDs'=as.character(lapply(metadata_id$associated_entities, FUN = sub_TCGA_ID_func)), 
                                  stringsAsFactors = F)
rownames(ID_conversion_table) <- ID_conversion_table$file_names
colnames(expr_FPKM) <- append('gene_id', ID_conversion_table[colnames(expr_FPKM)[2:length(expr_FPKM)],'TCGA_IDs'])

5. 进一步处理表达矩阵，将第一列的ENSEMBL GENE ID (eg. ENSG00000273842.1)转化为HUGO GENE SYMBOL (eg. DVL1)，便于认识和处理。并去除表达矩阵中可能存在的重复基因。最终生成的表达矩阵前三列为Gene symbol, gene id和bio_type.

TIPS：ENSEMBL数据库下载注释GTF文件速度真的很一言难尽。作为参考，挂了Yale的梯子，45mb左右的文件下载了2个小时。建议不要着急，要微笑:)))

注：ENSEMBL ID与GENE SYMBOL转换时不含其中的小数点和之后的数字，因此需要将其先去除后再使用GTF注释文件转换。

library('tidyr')
library('rtracklayer')
#trim the gene ID by decimals and translate them into gene symbols
expr_FPKM <- tidyr::separate(expr_FPKM, gene_id, into = c('gene_id', 'junk'), sep='\\.')
expr_FPKM <- expr_FPKM[,-2]
#download the annotation file containing ENS IDs, gene symbols, and biotypes.
download.file('ftp://ftp.ensembl.org/pub/release-99/gtf/homo_sapiens/Homo_sapiens.GRCh38.99.chr.gtf.gz',
              'Homo_sapiens.GRCh38.99.chr.gtf.gz')
R.utils::gunzip('Homo_sapiens.GRCh38.99.chr.gtf.gz')
gtf1 <- rtracklayer::import('Homo_sapiens.GRCh38.99.chr.gtf')
gtf_df <- as.data.frame(gtf1)
#Only select the protein coding genes.
mRNA_exprSet <- gtf_df %>%
  dplyr::filter(type=="gene",gene_biotype=="protein_coding") %>%
  dplyr::select(c(gene_name,gene_id,gene_biotype)) %>%
  dplyr::inner_join(expr_FPKM,by ="gene_id")
mRNA_exprSet <- mRNA_exprSet[!duplicated(mRNA_exprSet$gene_name),]

6. 载入临床特征数据，并进行其中缺失数据的清洗。本例中需要保存的临床数据包括病人的性别（gender），肿瘤分期(tumor stage)，生存时间（overall survival），以及生存状态（原vital status，数据清洗后改为censoring status）。

TIPS：之前下载的392个数据中既包含所有病人的癌组织也包含部分病人的癌旁组织，因此文件数量多于实际病人数。此时获得的临床特征数据以病人为单位（n=349），小于下载总样本数，等于总样本中癌症样本数。

注：overall survival的计算方法：临床数据中包含days_to_death以及days_to_last_follow_up两行数据，仅vital status为Dead的样本会有days_to_death数据，而Alive的样本记录有days_to_last_follow_up，以及有极少数样本两项数据均缺失。因此在overall_survival中，死亡样本以days_to_death为结局发生天数，而存活样本以days_to_last_follow_up作为删失发生天数。overall_survival与censoring_status配合，在生存分析中会用到。

clinical_data_original <- jsonlite::fromJSON('clinical.cart.2020-03-08.json')
filter_follow_up <- function(x){
  x$days_to_last_follow_up
}
filter_tumor_stage <- function(x){
  x$tumor_stage
}
days_to_last_follow_up <- as.numeric(lapply(clinical_data_original$diagnoses, filter_follow_up))
tumor_stage <- as.character(lapply(clinical_data_original$diagnoses, filter_tumor_stage))
clinical_trait <- clinical_data_original$demographic[,c('submitter_id', 'gender', 'days_to_death', 'vital_status')] %>%
  cbind(days_to_last_follow_up, tumor_stage)
clinical_trait$tumor_stage <- as.character(clinical_trait$tumor_stage)
clinical_trait <- tidyr::separate(clinical_trait, col= 'submitter_id',
                                  into = c('submitter_id', 'junk'), sep='_', remove = T)[,-2]
clinical_trait <- clinical_trait[!duplicated(clinical_trait$submitter_id),]
clinical_trait[is.na(clinical_trait$days_to_death), 'days_to_death'] <- clinical_trait[is.na(clinical_trait$days_to_death), 'days_to_last_follow_up']
clinical_trait <- clinical_trait[,-5]
names(clinical_trait) <- c('submitter_id', 'gender', 'overall_survival', 'censoring_status', 'tumor_stage')
clinical_trait <- clinical_trait %>% dplyr::filter(!(is.na(overall_survival))) %>% dplyr::filter(tumor_stage != 'not reported')

7. 通过TCGA_ID将所有样品分组。

思路是：TCGA ID由dash隔开，第四个模块记录了样本的分组。若该模块中编号为01-09则为肿瘤，若编号在10-29则为正常癌旁样本。

#classify the samples by TCGA_IDs (TCGA-XX-XXXX-NN: NN < 10 are cancer samples, otherwise are normal samples).
condition_table <- tidyr::separate(data=data.frame('ids'=colnames(expr_FPKM)[2:length(expr_FPKM)]),
                                  col='ids',
                                  sep='-',
                                  into=c('c1','c2','c3','c4','c5','c6','c7'))
condition_table <- cbind('ids'=colnames(expr_FPKM)[2:length(expr_FPKM)], condition_table) %>%
  cbind('submitter_IDs'=substr(colnames(expr_FPKM)[2:length(expr_FPKM)], 1, 12))
condition_table$c4 <- gsub(condition_table$c4, pattern = '^0..', replacement = 'cancer')
condition_table$c4 <- gsub(condition_table$c4, pattern = '^1..', replacement = 'normal')
condition_table <- condition_table[,c('ids','c4','submitter_IDs')]
names(condition_table) <- c('TCGA_IDs', 'sample_type', 'submitter_id')
condition_table$submitter_id <- as.character(condition_table$submitter_id)
condition_table$TCGA_IDs <- as.character(condition_table$TCGA_IDs)

8. 去除表达矩阵中肿瘤或正常组内可能存在的重复样本，并将cancer样本与临床数据匹配。

由于一个样本可能测序多次，所以可以通过submitter_ID (也就是TCGA_ID的前12位)将组内重复去除。我这里采取的办法是保留第一次在表达矩阵中出现的样本，去除随后出现的重复样本。此外由于此前有部分病人存在信息缺失，已被剔除，现也只保留含有全部信息的病人cancer样本。（normal样本不需匹配临床信息，因为仅用于差异表达分析，不会涉及后续生存分析或WCGNA分析）

注：由于normal组的癌旁组织也来源于病人，所以需要组内分开去除重复，否则来源于统一病人的normal和cancer样本会被去除其一。

#remove the duplicated samples within groups and only maintain the cancer samples with clinical data.
condition_table_cancer <- condition_table[condition_table$sample_type=='cancer',]
condition_table_cancer <- condition_table_cancer[!duplicated(condition_table_cancer$submitter_id),] %>%
  dplyr::inner_join(clinical_trait, by='submitter_id')
condition_table_normal <- condition_table[condition_table$sample_type=='normal',]
condition_table_normal <- condition_table_normal[!duplicated(condition_table_normal$submitter_id),]
condition_table <- rbind(condition_table_cancer[1:3], condition_table_normal)

9. 获得最终可以进行下游分析的规范FPKM表达矩阵。

该矩阵第一列为GENE SYMBOL，第一行为样本TCGA_ID，其余为FPKM表达数据。

#get the final total FPKM expr matrix and expr matrix containing cancer samples for downstream cibersort analyses.
mRNA_exprSet_cancer_only <- mRNA_exprSet[,c('gene_name', condition_table_cancer$TCGA_IDs)]
mRNA_exprSet <- mRNA_exprSet[,c('gene_name', condition_table$TCGA_IDs)]
write.csv(mRNA_exprSet, '/Users/zx/Desktop/TCGA/TCGA-LUSC-white/FPKM/tutorial/FPKM_matrix_LUSC_.csv', quote = F, row.names = F)

到此为止TCGA表达矩阵的数据下载和清洗整合就全部结束了！！！

开不开心！！！（不）

最后输出的表达矩阵就是长这样的：

FPKM矩阵在excel中的部分展示

说在最后：关于counts数据的下载...

几乎一模一样，镜像操作即可！唯一的区别是在初始表达矩阵整合完成后，需要去除dataframe的最后5行，该5行是样本的总结性数据，在后续分析中不会用到，直接去除即可。

某样本中Counts表达矩阵需要删除的部分