rmd
插入代碼塊ctrl+alt+I

級(jí)別標(biāo)題：＃.png

TCGA 癌癥類型.png

1.數(shù)據(jù)下載
gdc-client 官方下載工具
UCSC Xena 瀏覽器 分類打包，直接下載
gdcRNAtools基于gdc-client下載并簡(jiǎn)化整理，用R語(yǔ)言完成

image.png

image.png

image.png

image.png

title: "gdc-client數(shù)據(jù)下載"
author: "Sun Xiaojie"
output: html_document
editor_options:
chunk_output_type: console

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
knitr::opts_chunk$set(fig.width = 6,fig.height = 6,collapse = TRUE)
knitr::opts_chunk$set(message = FALSE)

本文的內(nèi)容是用GDC下載并整理表達(dá)矩陣和臨床信息數(shù)據(jù)。

1.從網(wǎng)頁(yè)選擇數(shù)據(jù)，下載manifest文件

數(shù)據(jù)存放網(wǎng)站：https://portal.gdc.cancer.gov/

在Repository勾選自己需要的case和file類型。以CHOL為例：
case-Project選擇TCGA-CHOL。

2.從網(wǎng)頁(yè)下載數(shù)據(jù)
示例miRNA-Seq，clinical數(shù)據(jù)和joson沒有區(qū)別
gdc tcga網(wǎng)站
case——program??TCGA——project??選擇目標(biāo)癌癥

files——Data Category??clinical——date format??bcr xml
下載manifes清單

files——data category??transcriptome profoling ——date type ??gene Expression qualification——experimental strategy ??RNA-Seq——workflow type??HTseq-counts
manifest進(jìn)行下載，儲(chǔ)存至新建文件夾

file-選擇如圖：

左右分別是expdata 和clinical的樣本選擇截圖。選好后，點(diǎn)擊右側(cè)manifest鍵下載對(duì)應(yīng)的清單文件。

2.使用gdc-client工具下載

注意：

將gdc-client(mac)或gdc-client.exe(windows)放在工作目錄下；

將manifest文件放在工作目錄下。

options(stringsAsFactors = F)
library(stringr)
cancer_type="TCGA-CHOL"  #定義變量
if(!dir.exists("clinical"))dir.create("clinical")   #判斷文件是否存在，否則新建
if(!dir.exists("expdata"))dir.create("expdata")
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語(yǔ)言并整理。

3.整理臨床信息

library(XML)  #加載讀取xml的包
result <- xmlParse("./clinical/142aea0e-7a7b-4ac4-9dbb-0f62e2379599/nationwidechildrens.org_clinical.TCGA-W5-AA2O.xml")  #

rootnode <- xmlRoot(result)  #提取xml文件節(jié)點(diǎn)
xmlSize(rootnode)

print(rootnode[1])  #對(duì)節(jié)點(diǎn)輸出，看看各自節(jié)點(diǎn)信息
#print(rootnode[2])  #患者有用的臨床信息都在第二個(gè)節(jié)點(diǎn)
xmldataframe <- xmlToDataFrame(rootnode[2])  #轉(zhuǎn)成數(shù)據(jù)框
head(t(xmlToDataFrame(rootnode[2]))) #轉(zhuǎn)置數(shù)據(jù)框，按行排列

xmls = dir("clinical/",pattern = "*.xml$",recursive = T)
cl = list()   #讀取"clinical/"文件夾下所有這樣格式的xml，所有子文件夾也要掃描
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)
clinical[1:3,1:3]

4.整理表達(dá)矩陣

探索數(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")  #患者1的基因表達(dá)矩陣
x2 = read.table("expdata/10d08172-48d2-49e7-b760-721163492cc1/c1071bcd-5a0c-4e09-a578-fc4b6dbe26ad.htseq.counts.gz",
                row.names = 1,sep = "\t")  #患者2的基因表達(dá)矩陣
identical(rownames(x),rownames(x2)) #驗(yàn)證患者基因表達(dá)順序是否一致

由此可知，他們的geneid順序是一致的，可以直接cbind，不會(huì)導(dǎo)致順序錯(cuò)亂。

批量讀取所有的counts.gz文件（將不同患者的表達(dá)矩陣所處文件名匯總到一個(gè)列表）

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

前面是文件夾的名字/后面是文件的名字

image.png

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)
exp[1:4,1:4]

image.png

image.png

發(fā)現(xiàn)問題：這樣產(chǎn)生出來的表達(dá)矩陣沒有列名（也就是不同患者id信息無法顯示，之后就與臨床信息無法對(duì)應(yīng)）

解決辦法：找到一個(gè)文件名與樣本ID一一對(duì)應(yīng)的文件。cart-json文件。

meta <- jsonlite::fromJSON("metadata.cart.2020-03-23.json")
colnames(meta)
ids <- meta$associated_entities;class(ids)
ids[[1]]
ids[[1]][,1]

可以看到，meta$associated_entities是個(gè)列表，這個(gè)列表里包含數(shù)據(jù)框，數(shù)據(jù)框的第一列內(nèi)容就是tcga樣本id。

image.png

注意，換了數(shù)據(jù)需要自己探索存放在哪一列。不一定是完全一樣的，需要確認(rèn)清楚。
批量提取associated_entities（ids）列表元素里面的基因名

ID = sapply(ids, function(x){x[,1]})  
#當(dāng)列表元素里只有一個(gè)值時(shí)， sapply將列表轉(zhuǎn)化為帶名字的向量 
#創(chuàng)造函數(shù)function(x){運(yùn)算方式}
#取ids列表每個(gè)元素里面的第一列組成向量
file2id = data.frame(file_name = meta$file_name,
                     ID = ID) #以基因表達(dá)矩陣為模版，根據(jù)文件夾與臨床樣本ID對(duì)應(yīng)關(guān)系，調(diào)整臨床信息數(shù)據(jù)

image.png

image.png

文件名與TCGA樣本ID的對(duì)應(yīng)關(guān)系已經(jīng)得到，接下來是將其添加到表達(dá)矩陣中，成為行名。需要找到讀取文件的順序，一一對(duì)應(yīng)修改。

head(file2id$file_name,2)
head(count_files,2)
count_files2 = stringr::str_split(count_files,"/",simplify = T)[,2]
table(count_files2 %in% file2id$file_name)
#檢查基因表達(dá)矩陣文件名與臨床信息文件名是否完全一致

文件名格式并不完全一樣，截取count_files“/”后面的名字組成新向量

image.png

count_files2（表達(dá)矩陣）的順序就是列名的順序，根據(jù)它來調(diào)整file2id（樣本）的順序。此處需要再次理解一下match函數(shù)。

file2id = file2id[match(count_files2,file2id$file_name),]
colnames(exp) = file2id$ID
exp[1:4,1:4]

表達(dá)矩陣整理完成,需要過濾一下那些在很多樣本里表達(dá)量都為0的基因。過濾標(biāo)準(zhǔn)不唯一。

dim(exp)
#exp = exp[rowSums(exp)>0,]  
 #篩選條件所有樣本A基因表達(dá)量總和>0則保留此行
exp = exp[apply(exp, 1, function(x) sum(x > 1) > 9), ]
#對(duì)exp這個(gè)矩陣1（行）/2（列）執(zhí)行后面的函數(shù)，A基因表達(dá)了>1的樣本總數(shù)>9的行保留，>9這個(gè)條件可根據(jù)情況改變
dim(exp)
exp[1:4,1:4]
exp = as.matrix(exp)

分組信息

image.png

根據(jù)樣本ID的第14-15位，給樣本分組（tumor和normal）

table(str_sub(colnames(exp),14,15))  #根據(jù)前圖判斷對(duì)照和腫瘤，截取tacg基因名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,cancer_type,file = paste0(cancer_type,"_gdc.Rdata"))