如何使用TCGAbiolinks進(jìn)行數(shù)據(jù)預(yù)處理？

引言：在前面我們了解了如何使用TCGAbiolinks檢索并獲取TCGA數(shù)據(jù)庫(kù)的公開數(shù)據(jù)。今天小編就用前面涉及到的代碼，下載今天數(shù)據(jù)準(zhǔn)備需要用到的TCGA樣本數(shù)據(jù)。

一、數(shù)據(jù)下載階段

第一步：GDCquery（）篩選我們需要的數(shù)據(jù)，TCGAbiolinks包下載TCGA數(shù)據(jù)進(jìn)行表達(dá)差異分析－肝癌案例

library（＂TCGAbiolinks＂）

query ＜－ GDCquery（project ＝ ＂TCGA－LIHC＂，

data．category ＝ ＂Transcriptome Profiling＂，

data．type ＝ ＂Gene Expression Quantification＂，

workflow．type ＝ ＂HTSeq － Counts＂）

上圖為通過TCGA GDC鏈接中根據(jù)篩選條件查看的符合要求結(jié)果。下圖為通過GDCquery（）函數(shù)中傳入對(duì)應(yīng)的參數(shù)得到的結(jié)果。兩者對(duì)比，我們可以發(fā)現(xiàn)，兩者是一模一樣的。說明代碼執(zhí)行正確。前面一期中，我們有詳細(xì)談及 GDCquery，可做參考。

samplesDown ＜－ getResults（query，cols＝c（＂cases＂））

＃getResults（query， rows， cols）根據(jù)指定行名或列名從query中獲取結(jié)果，此處用來獲得樣本的barcode

＃ 此處共檢索出424個(gè)barcodes

getResults（）中用到的參數(shù)：

參數(shù)用法query

來自GDCquery的結(jié)果rows用于指定特定的行cols用于指定特定的列

＃ 從samplesDown中篩選出TP（實(shí)體腫瘤）樣本的barcodes

＃ TCGAquery＿SampleTypes（barcode， typesample）

＃ TP代表PRIMARY SOLID TUMOR；NT－代表Solid Tissue Normal（其他組織樣本可參考學(xué)習(xí)文檔）

＃＃此處共檢索出371個(gè)TP樣本barcodes

dataSmTP ＜－ TCGAquery＿SampleTypes（barcode ＝ samplesDown，

typesample ＝ ＂TP＂）

＃ 從samplesDown中篩選出NT（正常組織）樣本的barcode

＃此處共檢索出50個(gè)NT樣本barcodes

dataSmNT ＜－ TCGAquery＿SampleTypes（barcode ＝ samplesDown，

typesample ＝ ＂NT＂）

TCGAquery＿SampleTypes中的參數(shù)詳解：

參數(shù)用法barcodeTCGA中的barcodes列表typesample用于指定篩選哪種類型的組織樣本，如腫瘤組織“TP”，正常組織“NT”

補(bǔ)充TCGA中的組織樣本類型：

TPPRIMARY SOLID TUMORTMMetastaticTRRECURRENT SOLID TUMORTAMAdditional MetastaticTBPrimary Blood Derived Cancer－Peripheral BloodTHOCHuman Tumor Original CellsTRBMRecurrent Blood Derived Cancer－Bone MarrowTBM Primary Blood Derived Cancer－Bone MarrowTAPAdditional－New PrimaryNB Blood Derived Normal NTSolid Tissue NormalNBCBuccal Cell Normal？？？NEBVEBV Immortalized NormalNBMBone Marrow Normal

＃＃＃設(shè)置barcodes參數(shù)，篩選符合要求的371個(gè)腫瘤樣本數(shù)據(jù)和50正常組織數(shù)據(jù)

queryDown ＜－ GDCquery（project ＝ ＂TCGA－LIHC＂，

data．category ＝ ＂Transcriptome Profiling＂，

data．type ＝ ＂Gene Expression Quantification＂，

workflow．type ＝ ＂HTSeq － Counts＂，

barcode ＝ c（dataSmTP， dataSmNT））

＃barcode參數(shù)：根據(jù)傳入barcodes進(jìn)行數(shù)據(jù)過濾

上圖為 queryDown＜－GDCquery（）的結(jié)果，僅選擇了選擇371個(gè)正常組織和50個(gè)腫瘤組織樣本。

第二步：GDCdownload（）下載GDCquery（）得到的結(jié)果

＃ 下載數(shù)據(jù)，默認(rèn)存放位置為當(dāng)前工作目錄下的GDCdata文件夾中。

GDCdownload（queryDown，method ＝ ＂api＂， directory ＝ ＂GDCdata＂，

files．per．chunk ＝ 10）

＃method ；＂API＂或者＂client＂。＂API＂速度更快，但是容易下載中斷。

＃directory：下載文件的保存地址。Default： GDCdata。

＃files．per．chunk ＝ NULL：使用API下載大文件的時(shí)候，可以把文件分成幾個(gè)小文件來下載，可以解決下載容易中斷的問題。

GDCdownload（query ＝ queryDown）

說明：由于小編前面已經(jīng)下載過該TCGA數(shù)據(jù)，所以這里顯示的是421個(gè)文件已存在。如果還沒有下載的話，可能需要根據(jù)自己的網(wǎng)速等待一些時(shí)間。

顯示這樣的結(jié)果，就算下載成功啦！文件默認(rèn)保存在 Rstudio默認(rèn)路徑下的GDCdata中。前面就是我們利用第一期知識(shí)進(jìn)行數(shù)據(jù)下載環(huán)節(jié)，權(quán)當(dāng)溫習(xí)功課吧——接下來我們就開始此期的數(shù)據(jù)處理～～

二、數(shù)據(jù)處理

第三步：GDCprepare（）將前面GDCquery（）的結(jié)果準(zhǔn)備成R語言可處理的SE（SummarizedExperiment）文件。

＃讀取下載的數(shù)據(jù)并將其準(zhǔn)備到R對(duì)象中，在工作目錄生成（save＝TRUE）LIHC＿case．rda文件

＃ GDCprepare（）：Prepare GDC data，準(zhǔn)備GDC數(shù)據(jù)，使其可用于R語言中進(jìn)行分析

dataPrep1 ＜－ GDCprepare（query ＝ queryDown， save ＝ TRUE， save．filename ＝

＂LIHC＿case．rda＂）

GDCprepare（）中的參數(shù)：

參數(shù)用法query來自GDCquery的結(jié)果save是否將結(jié)果保存為RData object，默認(rèn)為TRUEsave．filename文件名，如果沒有設(shè)置，系統(tǒng)將默認(rèn)設(shè)置directory文件數(shù)據(jù)的文件夾，默認(rèn)為“GDCdata”summarizedExperiment是否生成summarizedExperiment對(duì)象，默認(rèn)TRUE

第四步：TCGAanalyze＿Preprocessing（）對(duì)數(shù)據(jù)進(jìn)行預(yù)處理：使用spearman相關(guān)系數(shù)去除數(shù)據(jù)中的異常值

＃ 去除dataPrep1中的異常值，dataPrep1數(shù)據(jù)中含有腫瘤組織和正常組織的數(shù)據(jù)

＃ TCGAanalyze＿Preprocessing（object， cor．cut ＝ 0， filename ＝ NULL，

width ＝ 1000， height ＝ 1000， datatype ＝ names（assays（object））［1］）

＃ 函數(shù)功能描述：Array Array Intensity correlation （AAIC） and correlation boxplot to define outlier

dataPrep2 ＜－ TCGAanalyze＿Preprocessing（object ＝ dataPrep1，

cor．cut ＝ 0．6，

datatype ＝ ＂HTSeq － Counts＂）

＃將預(yù)處理后的數(shù)據(jù)dataPrep2，寫入新文件“LIHC＿dataPrep．csv”

write．csv（dataPrep2，file ＝ ＂LIHC＿dataPrep．csv＂，quote ＝ FALSE）

這里將生成一個(gè)array－array intensity correlation（AAIC）相關(guān)性熱圖，如下：

TCGAanalyze＿Preprocessing（）中的參數(shù)：

參數(shù)用法object來自TCGAprepare的結(jié)果cor．cut設(shè)置閾值，根據(jù)樣本中各個(gè)樣本之間的spearman相關(guān)系數(shù)進(jìn)行過濾。默認(rèn)為0filename設(shè)置生成圖片文件的名稱，默認(rèn)為PreprocessingOutput．pngwidth生成圖片的寬度？？ height生成圖片的高度datatype描述RangedSummarizedExperiment 數(shù)據(jù)類型的字符串

第五步：TCGAtumor＿purity（）篩選腫瘤純度大于60％的腫瘤barcodes

＃ TCGAtumor＿purity（barcodes， estimate， absolute， lump， ihc， cpe），使用來自5種方法的5個(gè)估計(jì)值作為閾值對(duì)TCGA樣本進(jìn)行過濾，這5個(gè)值是estimate， absolute， lump， ihc， cpe，這里設(shè)置cpe＝0．6（cpe是派生的共識(shí)度量，是將所有方法的標(biāo)準(zhǔn)含量歸一化后的均值純度水平，以使它們具有相等的均值和標(biāo)準(zhǔn)差）

＃篩選腫瘤純度大于等于60％的樣本數(shù)據(jù)

purityDATA ＜－ TCGAtumor＿purity（colnames（dataPrep1）， 0， 0， 0， 0， 0．6）

＃ filtered 為被過濾的數(shù)據(jù)， pure＿barcodes是我們要的腫瘤數(shù)據(jù)

Purity．LIHC＜－purityDATA＄pure＿barcodes

normal．LIHC＜－purityDATA＄filtered

filtered 為被過濾的數(shù)據(jù)（為正常組織的數(shù)據(jù)barcodes）， pure＿barcodes是我們要的腫瘤樣本barcodes。

第六步：將腫瘤表達(dá)矩陣與正常組織表達(dá)矩陣合并，進(jìn)行基因注釋

＃獲取腫瘤純度大于60％的340個(gè)腫瘤組織樣本＋50個(gè)正常組織樣本，共計(jì)390個(gè)樣本

puried＿data ＜－dataPrep2［，c（Purity．LIHC，normal．LIHC）］

第七步：進(jìn)行表達(dá)矩陣基因注釋

?；蜃⑨?，需要加載“SummarizedExperiment”包，“SummarizedExperiment container”每個(gè)由數(shù)字或其他模式的類似矩陣的對(duì)象表示。行通常表示感興趣的基因組范圍和列代表樣品。

＃if （！requireNamespace（＂BiocManager＂， quietly ＝ TRUE））

install．packages（＂BiocManager＂）

＃BiocManager：：install（＂SummarizedExperiment＂） ＃沒有的需要執(zhí)行下載代碼

library（＂SummarizedExperiment＂）

rowData（dataPrep1） ＃傳入數(shù)據(jù)dataPrep1必須為SummarizedExperiment對(duì)象

＃ DataFrame with 56512 rows and 3 columns

＃ ensembl＿gene＿id external＿gene＿name original＿ensembl＿gene＿id

＃ ＜character＞ ＜character＞ ＜character＞

＃ ENSG00000000003 ENSG00000000003 TSPAN6 ENSG00000000003．13

＃ ENSG00000000005 ENSG00000000005 TNMD ENSG00000000005．5

＃ ENSG00000000419 ENSG00000000419 DPM1 ENSG00000000419．11

＃ ENSG00000000457 ENSG00000000457 SCYL3 ENSG00000000457．12

＃將結(jié)果寫入文件“puried．LIHC．cancer．csv”

rownames（puried＿data）＜－rowData（dataPrep1）＄external＿gene＿name

write．csv（puried＿data，file ＝ ＂puried．LIHC．csv＂，quote ＝ FALSE）

第八步：進(jìn)行表達(dá)矩陣標(biāo)準(zhǔn)化和過濾，得到用于差異分析的表達(dá)矩陣

｀TCGAanalyze＿Normalization（）｀使用EDASeq軟件包標(biāo)準(zhǔn)化mRNA轉(zhuǎn)錄本和miRNA。

＃TCGAanalyze＿Normalization（）執(zhí)行EDASeq包中的如下功能：

1． EDASeq：：newSeqExpressionSet

2． EDASeq：：withinLaneNormalization

3． EDASeq：：betweenLaneNormalization

4． EDASeq：：counts

dataNorm ＜－ TCGAanalyze＿Normalization（tabDF ＝ puried＿data，

geneInfo ＝ geneInfo，

method ＝ ＂gcContent＂）

TCGAanalyze＿Normalization中的參數(shù)：

參數(shù)用法tabDFRNAseq表達(dá)矩陣，行代表基因，列代表樣本geneInfo關(guān)于geneLength和gcContent的20531個(gè)基因的矩陣，“geneInfoHT”和“geneInfo”可選。method選擇標(biāo)準(zhǔn)化的方法，基于’gcContent’ 或 ’geneLength’的標(biāo)準(zhǔn)化方法可選

＃將標(biāo)準(zhǔn)化后的數(shù)據(jù)再過濾，去除掉表達(dá)量較低（count較低）的基因，得到最終的數(shù)據(jù)

dataFilt ＜－ TCGAanalyze＿Filtering（tabDF ＝ dataNorm，

method ＝ ＂quantile＂，

qnt．cut ＝ 0．25）

str（dataFilt）

＃num ［1：13083， 1：340］ 274 2432 60347 1012 1947 ．．．

＃－ attr（＊， ＂dimnames＂）＝List of 2

＃ ．．＄ ： chr ［1：13083］ ＂A1BG＂ ＂A1CF＂ ＂A2M＂ ＂A4GALT＂ ．．．

＃ ．．＄ ： chr ［1：390］ ＂TCGA－DD－AAD5－01A－11R－A41C－07＂ ＂TCGA－DD－A4NO－01A－11R－A28V－07＂ ＂TCGA－EP－A2KA－01A－11R－A180－07＂ ＂TCGA－DD－AACP－01A－11R－A41C－07＂ ．．．

TCGAanalyze＿Filtering（）中的參數(shù)：

參數(shù)用法tabDF數(shù)據(jù)框或者矩陣，行代表基因，列代表來自TCGA的樣本method用于過濾較低count數(shù)的基因的方法，有’quantile’， ’varFilter’， ’filter1’， ’filter2’qnt．cut選擇均值作為過濾的閾值

最后將過濾后的數(shù)據(jù)寫入文件“TCGA＿LIHC＿final．csv”，就得到我們用于后續(xù)差異分析的表達(dá)文件：

write．csv（dataFilt，file ＝ ＂TCGA＿LIHC＿final．csv＂，quote ＝ FALSE）

＃保留的是390個(gè)樣本（前340腫瘤，后50正常組織）

今天的數(shù)據(jù)預(yù)處理就講到這里，接下來我們將分享：數(shù)據(jù)分析（差異表達(dá)分析、富集分析和聚類分析等）。如果你喜歡的話，就加入我們一起挖數(shù)據(jù)吧～～