TCGA數(shù)據(jù)下載方法：gdc-client纯命，Xena和gdcRNAtools
需要下載的數(shù)據(jù)
- 組學(xué)信息（樣本）：
（存儲(chǔ)單個(gè)病人表達(dá)數(shù)據(jù)，需整理為表達(dá)矩陣） +
（存儲(chǔ)樣本文件的詳細(xì)信息，可以為RNA-seq，miRNA-seq，exon/CNV等等）
- 臨床信息（患者）：
（存儲(chǔ)單個(gè)病人的臨床信息小泉，需整理為臨床信息表格，包含了患者ID冕杠，生存/死亡微姊，分期，性別/種族/年齡分预，死亡時(shí)間等信息兢交。）
1. gdc-client（官方下載工具，優(yōu)先選擇??）
1.1 網(wǎng)頁(yè)選擇數(shù)據(jù)下載manifest清單文件（作為gdc下載的參數(shù)之一）
在Repository勾選自己需要的case和file類型笼痹。以CHOL為例：
（這里選擇了一個(gè)project中的所有文件配喳，也可以更自由的勾選想要的文件）
選擇count文件酪穿，case-Project選擇TCGA-CHOL(這里沒(méi)顯示出來(lái))
選擇clinical文件，case-Project選擇TCGA-CHOL
選擇好了分別點(diǎn)擊Manifest下載
下載的兩個(gè)manifest文件：
將下載的兩個(gè)文件放在工作目錄下
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數(shù)據(jù)
if(!dir.exists("expdata"))dir.create("expdata") #新建文件夾存放要下載的count數(shù)據(jù)
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/"))
可以看到，下載的文件是按樣本存放的泌绣，每個(gè)文件單獨(dú)有一個(gè)文件夾钮追，而我們需要得到的是表格，需要將他們批量讀入R語(yǔ)言并整理赞别。
1.3 整理臨床信息
library(XML)
處理單個(gè)文件（示例）
result <- xmlParse("./clinical/142aea0e-7a7b-4ac4-9dbb-0f62e2379599/nationwidechildrens.org_clinical.TCGA-W5-AA2O.xml")
rootnode <- xmlRoot(result) #查看這個(gè)xml文件有多少個(gè)節(jié)點(diǎn)
xmlSize(rootnode)
[1] 2
#print(rootnode[1]) ##查看各段（節(jié)點(diǎn)）中的信息，第一段一般都是頭文件信息
#print(rootnode[2]) ##我們需要的患者信息在第二段
xmldataframe <- xmlToDataFrame(rootnode[2]) #把xml文件轉(zhuǎn)成data.frame
head(t(xmlToDataFrame(rootnode[2])))
View(xmldataframe)
# 生成只有一行的單個(gè)患者的信息數(shù)據(jù)框
批量處理文件配乓，并生成臨床信息表格
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])
} #生成包含了所有單個(gè)轉(zhuǎn)換好的文件數(shù)據(jù)框的列表
clinical <- do.call(rbind,cl) #??把cl列表中的每一個(gè)元素（數(shù)據(jù)框）按行合并
View(clinical)
1.4 整理表達(dá)矩陣
1.4.1 探索數(shù)據(jù)
先任選兩個(gè)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，不會(huì)導(dǎo)致順序錯(cuò)亂腰埂。
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
發(fā)現(xiàn)問(wèn)題：這樣產(chǎn)生出來(lái)的表達(dá)矩陣沒(méi)有列名。
解決辦法：找到一個(gè)文件名與樣本ID一一對(duì)應(yīng)的文件：cart-json文件屿笼。
1.4.3 下載并處理cart-json文件
- 下載文件
選擇好需要下載的count文件后點(diǎn)擊加入cart牺荠，再點(diǎn)擊cart，選擇metadata即可下載驴一。
- 尋找對(duì)應(yīng)關(guān)系
count矩陣exp的列所對(duì)應(yīng)的樣本名file_name與count_files的順序一致休雌，要將其轉(zhuǎn)換為樣本id，也就是在json中尋找file_name和sampleID的對(duì)應(yīng)關(guān)系
meta <- jsonlite::fromJSON("metadata.cart.2020-03-23.json")
View(meta)
colnames(meta)
ids <- meta$associated_entities; class(ids) #sampleID儲(chǔ)存在meta$associated_entities這一列中肝断。
# ?? meta是個(gè)data.frame杈曲，而meta$associated_entities是個(gè)list，比較神奇
ids[[1]]
ids[[1]][,1] #提出data.frame的第一列胸懈，也就是我們需要的sampleID（單個(gè)）
可以看到担扑，meta$associated_entities是個(gè)列表，這個(gè)列表里包含數(shù)據(jù)框趣钱，數(shù)據(jù)框的第一列內(nèi)容就是tcga樣本id涌献。（注意，換了數(shù)據(jù)需要自己探索存放在哪一列首有。不一定是完全一樣的洁奈，需要確認(rèn)清楚间唉。）
- 生成對(duì)應(yīng)關(guān)系矩陣
ID = sapply(ids,function(x){x[,1]})
file2id = data.frame(file_name = meta$file_name,
ID = ID)
View(file2id) #得到file_name和sampleID的對(duì)應(yīng)矩陣
1.4.4 修改矩陣列名
文件名與TCGA樣本ID的對(duì)應(yīng)關(guān)系已經(jīng)得到，接下來(lái)是將其添加到表達(dá)矩陣中利术，成為行名呈野。需要找到讀取文件的順序，一一對(duì)應(yīng)修改印叁。
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分開(kāi)轮蜕，后面的才是和file2id中和TCGA樣本ID對(duì)應(yīng)的file_name
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]
table(count_files2 %in% file2id$file_name)
count_files2的順序就是列名的順序昨悼，根據(jù)它來(lái)調(diào)整file2id的順序。此處需要再次理解一下match函數(shù)
file2id = file2id[match(count_files2,file2id$file_name),]
colnames(exp) = file2id$ID
View(exp)
1.5 表達(dá)矩陣過(guò)濾和添加分組信息
表達(dá)矩陣整理完成,需要過(guò)濾一下那些在很多樣本里表達(dá)量都為0的基因跃洛。過(guò)濾標(biāo)準(zhǔn)不唯一率触。
dim(exp)
#exp = exp[rowSums(exp)>0,] #這個(gè)標(biāo)準(zhǔn)太寬松
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ] #至少在10個(gè)樣本中有表達(dá)（總共36個(gè)樣本）
dim(exp)
exp[1:4,1:4]
exp = as.matrix(exp)
添加分組信息
樣本ID第1-12位是相應(yīng)的患者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（網(wǎng)頁(yè)工具）
網(wǎng)址：https://xenabrowser.net/datapages/
Xena簡(jiǎn)化了gdc-client的下載方法细燎，打包了TCGA中的每一個(gè)project两曼，允許單個(gè)project去下載相應(yīng)的數(shù)據(jù)。
下載count數(shù)據(jù)玻驻，臨床信息和survival信息（這里是代碼下載方法悼凑，也可以在網(wǎng)頁(yè)直接下載）
Xena網(wǎng)站把臨床信息和survival分開(kāi)了，因此需要單獨(dú)下載
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的數(shù)據(jù)璧瞬，需要下載別的數(shù)據(jù)的時(shí)候户辫，只需要把CHOL換成別的癌癥的縮寫(xiě)即可。
# download.file需要2個(gè)參數(shù)嗤锉，url設(shè)置下載的地址寸莫，destfile設(shè)置下載下來(lái)文件要叫什么名字。
需要注意的是档冬，這里下載的count數(shù)據(jù)是經(jīng)過(guò)了log轉(zhuǎn)化的膘茎，因此在分析前需要先進(jìn)行逆轉(zhuǎn)log操作。
dat = read.table("counts.tsv.gz",check.names = F,row.names = 1,header = T) #壓縮格式可以直接讀取
#逆轉(zhuǎn)log
dat = as.matrix(2^dat - 1)
dat[1:4,1:4]
as.character(dat[1:100,1:10]) #逆轉(zhuǎn)之后矩陣中會(huì)出現(xiàn)一些小數(shù)酷誓，有小數(shù)在進(jìn)行Deseq2運(yùn)算的時(shí)候會(huì)報(bào)錯(cuò)
# 用apply轉(zhuǎn)換為整數(shù)矩陣
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數(shù)據(jù)定制）
一站式分析完整教程：http://bioconductor.org/packages/devel/bioc/vignettes/GDCRNATools/inst/doc/GDCRNATools.html
代碼來(lái)自2021生信技能樹(shù)數(shù)據(jù)挖掘課
