TCGA数据下载和整理的三种方法

2021-06-17 本文已影响0人 Hayley笔记

TCGA数据下载方法：gdc-client，Xena和gdcRNAtools

TCGA癌症类型

需要下载的数据

组学信息（样本）： $\color{#ea4335}{counts文件}$ （存储单个病人表达数据，需整理为表达矩阵） + $\color{#ea4335}{json文件}$ （存储样本文件的详细信息，可以为RNA-seq，miRNA-seq，exon/CNV等等）

临床信息（患者）： $\color{#ea4335}{xml文件}$ （存储单个病人的临床信息，需整理为临床信息表格，包含了患者ID，生存/死亡，分期，性别/种族/年龄，死亡时间等信息。）

1. gdc-client（官方下载工具，优先选择⚠️）

1.1 网页选择数据下载manifest清单文件（作为gdc下载的参数之一）

在Repository勾选自己需要的case和file类型。以CHOL为例：
（这里选择了一个project中的所有文件，也可以更自由的勾选想要的文件）

选择count文件，case-Project选择TCGA-CHOL(这里没显示出来)

选择clinical文件，case-Project选择TCGA-CHOL

选择好了分别点击Manifest下载

下载的两个manifest文件：

id是将来下载下来的文件夹的名字；filename是每个个文件的名字；md5是一串编码，可以告诉你下载的数据是否完整；size是文件大小。

将下载的两个文件放在工作目录下

1.2 使用gdc-client下载count文件和xml文件

将gdc-client(mac)或gdc-client.exe(windows)放在工作目录下

library(stringr)
proj="TCGA-CHOL"
if(!dir.exists("clinical"))dir.create("clinical") #新建文件夹存放要下载的clinical数据
if(!dir.exists("expdata"))dir.create("expdata") #新建文件夹存放要下载的count数据
dir() #列出目录下的文件

#下面两行命令在terminal完成
./gdc-client download -m gdc_manifest_cl.2020-03-23.txt -d clinical
./gdc-client download -m gdc_manifest_expdata.2020-03-23.txt -d expdata

length(dir("./clinical/"))
length(dir("./expdata/"))

可以看到，下载的文件是按样本存放的，每个文件单独有一个文件夹，而我们需要得到的是表格，需要将他们批量读入R语言并整理。

1.3 整理临床信息

library(XML)

处理单个文件（示例）

result <- xmlParse("./clinical/142aea0e-7a7b-4ac4-9dbb-0f62e2379599/nationwidechildrens.org_clinical.TCGA-W5-AA2O.xml")
rootnode <- xmlRoot(result) #查看这个xml文件有多少个节点
xmlSize(rootnode)
[1] 2
#print(rootnode[1])  ##查看各段（节点）中的信息，第一段一般都是头文件信息
#print(rootnode[2])  ##我们需要的患者信息在第二段
xmldataframe <- xmlToDataFrame(rootnode[2]) #把xml文件转成data.frame
head(t(xmlToDataFrame(rootnode[2])))
View(xmldataframe)
# 生成只有一行的单个患者的信息数据框

批量处理文件，并生成临床信息表格

xmls = dir("clinical/",pattern = "*.xml$",recursive = T) #列出所有的xml文件
cl = list()
for(i in 1:length(xmls)){
  result <- xmlParse(paste0("clinical/",xmls[[i]]))
  rootnode <- xmlRoot(result)
  cl[[i]] = xmlToDataFrame(rootnode[2])
} #生成包含了所有单个转换好的文件数据框的列表
clinical <- do.call(rbind,cl) #⚠️把cl列表中的每一个元素（数据框）按行合并
View(clinical)

1.4 整理表达矩阵

1.4.1 探索数据

先任选两个counts文件读取，并观察geneid的顺序是否一致。

options(stringsAsFactors = F)
x = read.table("expdata/03aee74e-4e37-4a58-a720-c90e807d2f40/be5bc6a0-9720-47ac-953e-fa8d0c32cd82.htseq.counts.gz",
               row.names = 1,sep = "\t")
x2 = read.table("expdata/10d08172-48d2-49e7-b760-721163492cc1/c1071bcd-5a0c-4e09-a578-fc4b6dbe26ad.htseq.counts.gz",
                row.names = 1,sep = "\t")
identical(rownames(x),rownames(x2))

由此可知，他们的geneid顺序是一致的，可以直接cbind，不会导致顺序错乱。

1.4.2 批量读取所有的counts.gz文件，并按列合并得到count矩阵

count_files = dir("expdata/",pattern = "*.htseq.counts.gz$",recursive = T)

exp = list()
for(i in 1:length(count_files)){
  exp[[i]] <- read.table(paste0("expdata/",count_files[[i]]),row.names = 1,sep = "\t")
}
  
exp <- do.call(cbind,exp)
dim(exp)
# [1] 60488    45
exp[1:4,1:4]
#                      V2 V2.1 V2.2 V2.3
# ENSG00000000003.13 2504  226 4107 9646
# ENSG00000000005.5     0    5    0    1
# ENSG00000000419.11 1272 1146  741 1266
# ENSG00000000457.12  504  602  312 1317

发现问题：这样产生出来的表达矩阵没有列名。
解决办法：找到一个文件名与样本ID一一对应的文件：cart-json文件。

1.4.3 下载并处理cart-json文件

下载文件
选择好需要下载的count文件后点击加入cart，再点击cart，选择metadata即可下载。

寻找对应关系

count矩阵exp的列所对应的样本名file_name与count_files的顺序一致，要将其转换为样本id，也就是在json中寻找file_name和sampleID的对应关系

meta <- jsonlite::fromJSON("metadata.cart.2020-03-23.json")
View(meta)

colnames(meta)
ids <- meta$associated_entities; class(ids) #sampleID储存在meta$associated_entities这一列中。
# ⚠️ meta是个data.frame，而meta$associated_entities是个list，比较神奇
ids[[1]] 
ids[[1]][,1] #提出data.frame的第一列，也就是我们需要的sampleID（单个）

可以看到，meta$associated_entities是个列表，这个列表里包含数据框，数据框的第一列内容就是tcga样本id。（注意，换了数据需要自己探索存放在哪一列。不一定是完全一样的，需要确认清楚。）

生成对应关系矩阵

ID = sapply(ids,function(x){x[,1]})
file2id = data.frame(file_name = meta$file_name,
                     ID = ID)
View(file2id) #得到file_name和sampleID的对应矩阵

1.4.4 修改矩阵列名

文件名与TCGA样本ID的对应关系已经得到，接下来是将其添加到表达矩阵中，成为行名。需要找到读取文件的顺序，一一对应修改。

head(file2id$file_name,2)
# [1] "4696ce44-29bf-41ea-b866-ddb17c376e94.htseq.counts.gz"
# [2] "46fd54e8-5ab7-43f1-bb88-01a163ff121f.htseq.counts.gz"
head(count_files,2)
# [1] "03aee74e-4e37-4a58-a720-c90e807d2f40/be5bc6a0-9720-47ac-953e-fa8d0c32cd82.htseq.counts.gz"
# [2] "0d2c466e-d8b8-4b8f-9909-0be2175fa6a0/39fae157-a126-4635-a212-137065a398f9.htseq.counts.gz"
# 可以看到文件名还是有一些不一样，count_files的文件名包含了文件夹的名字，需要使用str_split分开，后面的才是和file2id中和TCGA样本ID对应的file_name
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]
table(count_files2 %in% file2id$file_name)

count_files2的顺序就是列名的顺序，根据它来调整file2id的顺序。此处需要再次理解一下match函数

file2id = file2id[match(count_files2,file2id$file_name),]
colnames(exp) = file2id$ID
View(exp)

1.5 表达矩阵过滤和添加分组信息

表达矩阵整理完成,需要过滤一下那些在很多样本里表达量都为0的基因。过滤标准不唯一。

dim(exp)
#exp = exp[rowSums(exp)>0,] #这个标准太宽松
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ] #至少在10个样本中有表达（总共36个样本）
dim(exp)
exp[1:4,1:4]
exp = as.matrix(exp)

添加分组信息
样本ID第1-12位是相应的患者ID，第14和15位，<10是tumor，>=10是normal。

table(str_sub(colnames(exp),14,15))
Group = ifelse(as.numeric(str_sub(colnames(exp),14,15)) < 10,'tumor','normal')
Group = factor(Group,levels = c("normal","tumor"))
table(Group)
save(exp,clinical,Group,proj,file = paste0(proj,"_gdc.Rdata"))

2 Xena（网页工具）

网址：https://xenabrowser.net/datapages/

Xena简化了gdc-client的下载方法，打包了TCGA中的每一个project，允许单个project去下载相应的数据。
下载count数据，临床信息和survival信息（这里是代码下载方法，也可以在网页直接下载）
Xena网站把临床信息和survival分开了，因此需要单独下载

if(T){
  download.file(url = "https://gdc.xenahubs.net/download/TCGA-CHOL.htseq_counts.tsv.gz",destfile = "counts.tsv.gz") 
  download.file(url = "https://gdc.xenahubs.net/download/TCGA-CHOL.GDC_phenotype.tsv.gz",destfile = "phenotype.tsv.gz")
  download.file(url = "https://gdc.xenahubs.net/download/TCGA-CHOL.survival.tsv.gz",destfile = "survival.tsv.gz")
}
# 这里下载的是CHOL的数据，需要下载别的数据的时候，只需要把CHOL换成别的癌症的缩写即可。
# download.file需要2个参数，url设置下载的地址，destfile设置下载下来文件要叫什么名字。

需要注意的是，这里下载的count数据是经过了log转化的，因此在分析前需要先进行逆转log操作。

dat = read.table("counts.tsv.gz",check.names = F,row.names = 1,header = T) #压缩格式可以直接读取
#逆转log
dat = as.matrix(2^dat - 1)
dat[1:4,1:4]
as.character(dat[1:100,1:10]) #逆转之后矩阵中会出现一些小数，有小数在进行Deseq2运算的时候会报错

# 用apply转换为整数矩阵
exp = apply(dat, 2, as.integer)
exp[1:4,1:4] #行名消失
rownames(exp) = rownames(dat)

clinical = read.table("phenotype.tsv.gz",fill = T,header = T,sep = "\t")
surv = read.table("survival.tsv.gz",header = T)
clinical[1:4,1:4]
surv[1:4,1:4]

3 gdcRNAtools（R包，为GDC数据定制）

一站式分析完整教程：http://bioconductor.org/packages/devel/bioc/vignettes/GDCRNATools/inst/doc/GDCRNATools.html

代码来自2021生信技能树数据挖掘课