【工具】TCGAbiolinks分析TCGA數據（DEA篇）

首先，你需要下載TCGAbiolinks

> source("https://bioconductor.org/biocLite.R")

> biocLite("TCGAbiolinks")

記住要Citation哦

「TCGAbiolinks: an R/Bioconductor package for integrative analysis of TCGA

data.」 Nucleic acids research (2015): gkv1507. (Colaprico,

Antonio and Silva, Tiago C. and Olsen, Catharina and Garofano, Luciano and

Cava, Claudia and Garolini, Davide and Sabedot, Thais S. and Malta, Tathiane M.

and Pagnotta, Stefano M. and Castiglioni, Isabella and Ceccarelli, Michele and

Bontempi, Gianluca and Noushmehr, Houtan 2016)

下載數據

可供下載的數據有兩種，一個是Harmonized 數據；另一個是Legacy數據

兩者數據有什麼區別呢？

Legacy數據是TCGA原始數據，Harmonized是將不同項目的原始數據進行整合、處理過的數據。

這裡我們只考慮Harmonized數據

用到的Function是GDCquery(project=」」,

data.category=」」, data.type=」」,

workflow.type=」」)

下圖是一個在RNA層面上的兩種數據，Harmonized 和Legacy，的樹狀圖總結。

先看一下自己所在位置，好知道過會兒下載到哪裡了:

getwd()

接下來開始進行下載前的query：

> library(TCGAbiolinks)

> query <- GDCquery(project = "TCGA-COAD",

data.category = "Transcriptome Profiling",

data.type = "Gene Expression Quantification",

workflow.type = "HTSeq - FPKM-UQ")

其中，project的選項有如下：

你可以把例子中的project隨意替換成你想要研究的腫瘤。

接下來，你需要下載query

> GDCdownload(query)

之後，將下載好的query轉換成一個SummerizedExperiment的文件，這個以rda為後綴的文件是一個總結性文件，有了它，我們可以不再需要之前下載的raw數據，所以後面的remove.files.prepared可以選擇True，這樣會把之前下載的大量文件刪除，當然也可以留著不刪除（即default）。

> dataCOAD <- GDCprepare(query, save = TRUE,

save.filename = "dataCOAD_summerizedExperiment.rda",

remove.files.prepared = TRUE)

可以看一看rda文件，用到的package是SummarizedExperiment

> library(SummarizedExperiment)

> samples.information=colData(dataCOAD)

數據準備好了，我們接下來開始進行DEA分析。

所謂DEA，也就是Differential Expression Analysis，將Tumor組和對照組進行比較。

首先，將剛才GDCprepare好的數據進行normalization，用normalization()

這裡注意geneInfo=geneInfoHT，default其實是geneInfo，但由於我們前面選擇的是HTseq，所以要選擇geneInfoHT

> dataNorm <- TCGAanalyze_Normalization(tabDF = dataCOAD, geneInfo = geneInfoHT)

之後，常規選擇，用Filtering()

> dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,

method ="quantile",

qnt.cut = 0.25)

接著，定義對照組（這裡的對照組是Solid normal tissue），用到SampleType()

> samplesNT <- TCGAquery_SampleTypes(barcode = colnames(dataFilt),

typesample = c("NT"))

定義腫瘤組，用SampleType()

> samplesTP <- TCGAquery_SampleTypes(barcode = colnames(dataFilt),

typesample = c("TP"))

進行DEA分析，用到DEA()

> dataDEGs <- TCGAanalyze_DEA(mat1 =dataFilt[,samplesNT],

mat2 = dataFilt[,samplesTP],

Cond1type = "Normal",

Cond2type = "Tumor",

fdr.cut = 0.01 ,

logFC.cut = 1,

method = "glmLRT")

最後，將分析好的數據整入進一個表格里，用到LevelTab()

> dataDEGsFiltLevel <- TCGAanalyze_LevelTab(dataDEGs,"Tumor","Normal",

dataFilt[,samplesTP],dataFilt[,samplesNT])

將表格保存到一個csv的文件

write.csv(dataDEGsFiltLevel,file="DEA_COAD.csv")

打開文件，你將得到一個這樣的表格

DEA分析完畢了，是不是很簡單呢？

PS:

如果運行途中總是奇怪地出現錯誤，就輸入rm(list = ls())，清空前面的數據，從頭開始運行。

另外，對於RNA-seq來說，有些TCGA的數據是只有腫瘤的數據，而沒有對照組。

就RNA-seq而言，對照組一般是Solid Tissue Normal，而不是血，原因是由於血和腫瘤的RNA差異太大了。然而Solid Tissue Normal占的是少數，有些時候你會發現這種腫瘤根本就沒有。這種情況下，就沒法進行DEA分析了。

另外，教大家一個辨別Tissue Type的方法：

• TCGA-12-4567-01-blah-blah --> 這是Normal
• TCGA-12-4567-11-blah-blah --> 這是tumor

注意黑體的部分。01-09是tumor；10-19是Normal；20-29是Control

推薦閱讀：