TCGA學(xué)習(xí)01：數(shù)據(jù)下載與整理 - 簡書
TCGA學(xué)習(xí)02：差異分析 - 簡書
TCGA學(xué)習(xí)03：生存分析 - 簡書
TCGA學(xué)習(xí)04：建模預(yù)測-cox回歸 - 簡書
TCGA學(xué)習(xí)04：建模預(yù)測-lasso回歸 - 簡書
TCGA學(xué)習(xí)04：建模預(yù)測-隨機(jī)森林&向量機(jī) - 簡書

前言交代

1、學(xué)習(xí)參考

之前參加了生信技能樹花花老師的TCGA數(shù)據(jù)挖掘試講課，收獲很多，最近整理一下上課筆記，同時參考了老師的簡書相關(guān)教程。生信入門的朋友也可微信加入生信星球公眾號，個人覺得很好的一個學(xué)習(xí)平臺。一起來學(xué)習(xí)吧~

2、筆記內(nèi)容

（1）在TCGA公共數(shù)據(jù)庫，下載某一類癌癥的某一種RNA的count數(shù)據(jù)，以及相應(yīng)的病人臨床信息；
（2）根據(jù)count數(shù)據(jù)進(jìn)行差異分析，結(jié)合病人臨床信息進(jìn)行生存分析以及建模預(yù)測

課程還介紹有幾點(diǎn)更深入的研究，就不記錄了，可以參看老師的簡書教程。

3、學(xué)習(xí)目標(biāo)

TCGA的大致流程；R操作.....

第一步：數(shù)據(jù)下載與整理

1、下載

1.1 官網(wǎng)下載“腳本”文件

• （1）進(jìn)入GDC官網(wǎng)的repository條目https://portal.gdc.cancer.gov/repository

• （2）先在case選項(xiàng)中分別選擇Primary Site（癌癥位置）、Programme（選TCGA），還可以再選Project，我這里沒選
  出于小樣本量考慮，選了testis的site，后來查了下原來是.......

  
  
  case

• （3）基于選擇的case，在file選項(xiàng)中選擇miRNAcount數(shù)據(jù)（因?yàn)樾。?，注意改下名字，避免沖突，下同。

  
  
  miRNA count數(shù)據(jù)

• （4）基于選擇的case，在file選項(xiàng)中選擇病人臨床信息數(shù)據(jù)，注意要是bcr xml格式

  
  
  image.png

這里注意一下count矩陣信息有157個樣本；臨床信息則只有135份，對差異分析應(yīng)該沒影響；不知道會對后面生存分析等會不會有影響，存疑。
還有會注意到上面下載的兩個文件只有幾十KB大小，因?yàn)閮H僅是類似下載“腳本”的文件，還需要利用GDC提供的程序下載到本地。

• （5）下載程序有不同系統(tǒng)版本，在https://gdc.cancer.gov/access-data/gdc-data-transfer-tool選擇適合自己系統(tǒng)的，一般window的R中提供的terminal終端使用window版本即可。下載代碼如下：

mkdir mirna clinical
head gdc_manifest.2020-05-02.clinical.txt
./gdc-client download -m gdc_manifest.2020-05-02.mirna.txt -d mirna/
# -d表示儲存路徑
./gdc-client download -m gdc_manifest.2020-05-02.clinical.txt -d clinical/
# 切換到R
length(dir("./clinical/"))  #統(tǒng)計(jì)子目錄數(shù)
length(dir("./mirna/"))

• 以上兩種數(shù)據(jù)文件就下載好了，但需要整理成我們想要的表達(dá)矩陣與臨床信息的文件類型。

老師教程還提供了另外兩種下載方式，基于R包實(shí)現(xiàn)的。但GDC下載還是比較常用的。

2、整理

2.1 整理臨床信息

library(XML)  #clinical臨床信息為xml格式，我們目的轉(zhuǎn)換成df
result <- xmlParse("./clinical/0f133012-23ef-4237-acfd-47b132b99775/nationwidechildrens.org_clinical.TCGA-W5-AA2Q.xml")
rootnode <- xmlRoot(result)
rootsize <- xmlSize(rootnode)  #結(jié)果顯示有兩個節(jié)點(diǎn)
print(rootnode[2])  #病人信息在第二部分
xmldataframe <- xmlToDataFrame(rootnode[2])  #轉(zhuǎn)換成df，發(fā)現(xiàn)是長長的一行
head(t(xmlToDataFrame(rootnode[2])))

臨床信息很多

xmls = dir("clinical/",pattern = "*.xml$",recursive = T)
#目的是取出所有的xml文件，pattern = "*.xml$"用正則表達(dá)式匹配；recursive = T表示遞歸模式。
head(xmls) #是帶路徑的

head(xmls)

td = function(x){
  result <- xmlParse(file.path("clinical/",x))  #補(bǔ)成全路徑
  rootnode <- xmlRoot(result)
  xmldataframe <- xmlToDataFrame(rootnode[2])[c(
    'bcr_patient_barcode',
    'vital_status',
    'days_to_death',
    'days_to_last_followup',
    'race_list',
    'days_to_birth',
    'gender' ,
    'stage_event'
  )]
  return(t(xmldataframe))  #返回每個病人的信息df
}  #編寫小函數(shù)，批量處理
cl = lapply(xmls,td)  # xmls即為td函數(shù)的x，lapply逐個處理，結(jié)果儲存在cl列表中
do.call(cbind,cl)

cl_df <- t(do.call(cbind,cl))  #將所有病人信息合并
cl_df[1:3,1:3]
class(cl_df) #為矩陣
clinical = data.frame(cl_df)
rownames(clinical) <- stringr::str_to_upper(clinical$bcr_patient_barcode)
clinical <- clinical[,-1]
dim(clinical)
clinical[1:4,1:4]

最終的clinical數(shù)據(jù)

2.2 整理表達(dá)矩陣

length(dir("./mirna/"))   #157
x = read.table("mirna/01b6ce76-61a8-4bc4-b0d7-f69d4ce187d6/3fa67c62-6856-42fe-a9f1-fa6831f42b53.mirbase21.mirnas.quantification.txt"
               ,header = TRUE)
head(x)  

原始每個樣本表達(dá)矩陣信息

mis = dir("mirna/",pattern = "*tification.txt$",recursive = T)
mis
length(mis)
#同樣先取名
ex = function(x){
  result <- read.table(file.path("mirna/",x),sep = "\t",header = TRUE)[,1:2]
  return(result)
}

#批量讀取函數(shù)
mi = lapply(mis,ex)   #結(jié)果是第一列為基因名，第二列為count數(shù)的兩列數(shù)據(jù)
mi_df <- t(do.call(cbind,mi))  #合并并轉(zhuǎn)置，因?yàn)?5個病人所以90行
dim(mi_df)
tmp = mi_df[1:4,1:4]
tmp

tmp

colnames(mi_df) <- mi_df[1,]   #添加列名
mi_df[1:4,1:5]
#奇數(shù)列（基因名）是重復(fù)多余的，只保留偶數(shù)列
mi_df <- mi_df[seq(2,nrow(mi_df),2),]
dim(mi_df)
mi_df[1:4,1:4]
#轉(zhuǎn)為數(shù)值型
mi_df <- apply(mi_df, 2, as.numeric)
mi_df[1:4,1:4]

mi_df

• 表達(dá)矩陣基本成形，但會發(fā)現(xiàn)致命的缺點(diǎn)--沒有行名，即病人信息，目前只知道原始矩陣的文件名（見上head(xmls)圖），但不是我們需要的病人ID類型，因此還要回GDC網(wǎng)站下載相應(yīng)的json文件，里面包含有名稱對應(yīng)關(guān)系。

• 下載方法（如下圖）：選擇count或臨床信息時，全選加入購物車，點(diǎn)擊購物車，然后點(diǎn)擊Metadata即可。這里因?yàn)槭墙o表達(dá)矩陣加行名，所以選擇count數(shù)據(jù)的157個樣本到購物車。

  
  
  下載json文件

#接下來利用json文件給表達(dá)矩陣加上行名----
dim(mi_df)
meta <- jsonlite::fromJSON("metadata.cart.2020-05-02.157.json")
#觀察到associated_entities列中每一格都是一個df，里面存放著病人的各種相關(guān)ID；
#其中entity_submitter_id即我們想要的,meta有157行，就對上了
meta$associated_entities[[1]]  

meta里的id

entity <- meta$associated_entities  
#取第四列，為列表，包含4個元素，每個元素為一個df
class(entity[157])
entity[[157]][4]
class(entity)
jh = function(x){
  as.character(x[4])
}   #取df中的entity_submitter_id，即我們最終想給矩陣加上的樣本名
jh(entity[[1]])
ID = sapply(entity,jh)   #取得了所有的病人ID
head(ID)

image.png

options(stringsAsFactors = F)
file2id = data.frame(file_name = meta$file_name,  
                     #meta中file_name列即為我們下的mrna的文件名字
                     ID = ID)   
#file2id表格157行2列，即儲存著我們想要的文件名與ID的對應(yīng)關(guān)系
head(mis)  #mis文件是我們下載mrna的gz文件的遞歸路徑，注意是有順序的

head(mis)

#注意mis的樣本順序與之前表達(dá)矩陣樣本順序相同
mis2 = stringr::str_split(mis,"/",simplify = T)[,2]   #取文件名
mis2[1] %in% file2id$file_name #序列匹配與否
#接下來就要把file2id列表按mis2排序，再把病人ID加到表達(dá)矩陣上

#match(A,B)  
head(match(mis2,file2id$file_name))
#排序可以驗(yàn)證一下，match返回值的第一個是38，意思mis2的第一個元素是file2id$file_name的第38個元素。
mis2[1] == file2id$file_name[38]
row_tcga = file2id[match(mis2,file2id$file_name),]
#調(diào)整file2id行順序以適應(yīng)我們下載的miran數(shù)據(jù)
#如上match(mis2,file2id$file_name)第一個值為4；就把第四行放到第一行
rownames(mi_df) = row_tcga$ID #最后一步 
#給mi_df添加行名
mi_df[1:4,1:4]

mi_df[1:4,1:4]

expr = t(mi_df)
dim(expr)
expr = expr[apply(expr, 1, function(x) {
  sum(x > 1) > 9  #過濾掉低表達(dá)基因
}), ]  
dim(expr) 
group_list=ifelse(as.numeric(substr(colnames(expr),14,15)) < 10,'tumor','normal')
# 01轉(zhuǎn)為tumor，大于10轉(zhuǎn)為normal
table(group_list)  #只有tumor樣本
# group_list
# tumor 
# 157 
group_list <- factor(group_list,levels = c("normal","tumor"))  #分組因子

• 為什么01轉(zhuǎn)為tumor，大于10轉(zhuǎn)為normal呢？參考癌癥類型和樣本代號詳解TCGA一文中，Sample號從01-29的，其中01-09是tumor，也就是癌癥樣本；其中10-29是normal，也就是癌旁。
  
  TCGA編號
• 不過遺憾的是我嘗試找的兩個小樣本癌中都是只有tumor樣本(<10)，就是說不能做差異分析了，就暫時先用花花老師的數(shù)據(jù)吧。

save(expr,clinical,group_list,file = "gdc.Rdata")

• 綜上，我們一共得到三個文件，用于接下來的分析

  
  
  表達(dá)矩陣&樣本信息&分組信息
• 不足之處：（1）表達(dá)矩陣與樣本信息數(shù)量不同；（2）表達(dá)矩陣全是tumor樣本。