如何使用TCGAbiolinks进行数据预处理？

引言：在前面我们了解了如何使用TCGAbiolinks检索并获取TCGA数据库的公开数据。今天小编就用前面涉及到的代码，下载今天数据准备需要用到的TCGA样本数据。

一、数据下载阶段

第一步：GDCquery（）筛选我们需要的数据，TCGAbiolinks包下载TCGA数据进行表达差异分析－肝癌案例

library（＂TCGAbiolinks＂）

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

data．category ＝ ＂Transcriptome Profiling＂，

data．type ＝ ＂Gene Expression Quantification＂，

workflow．type ＝ ＂HTSeq － Counts＂）

上图为通过TCGA GDC链接中根据筛选条件查看的符合要求结果。下图为通过GDCquery（）函数中传入对应的参数得到的结果。两者对比，我们可以发现，两者是一模一样的。说明代码执行正确。前面一期中，我们有详细谈及 GDCquery，可做参考。

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

＃getResults（query， rows， cols）根据指定行名或列名从query中获取结果，此处用来获得样本的barcode

＃ 此处共检索出424个barcodes

getResults（）中用到的参数：

参数用法query

来自GDCquery的结果rows用于指定特定的行cols用于指定特定的列

＃ 从samplesDown中筛选出TP（实体肿瘤）样本的barcodes

＃ TCGAquery＿SampleTypes（barcode， typesample）

＃ TP代表PRIMARY SOLID TUMOR；NT－代表Solid Tissue Normal（其他组织样本可参考学习文档）

＃＃此处共检索出371个TP样本barcodes

dataSmTP ＜－ TCGAquery＿SampleTypes（barcode ＝ samplesDown，

typesample ＝ ＂TP＂）

＃ 从samplesDown中筛选出NT（正常组织）样本的barcode

＃此处共检索出50个NT样本barcodes

dataSmNT ＜－ TCGAquery＿SampleTypes（barcode ＝ samplesDown，

typesample ＝ ＂NT＂）

TCGAquery＿SampleTypes中的参数详解：

参数用法barcodeTCGA中的barcodes列表typesample用于指定筛选哪种类型的组织样本，如肿瘤组织“TP”，正常组织“NT”

补充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

＃＃＃设置barcodes参数，筛选符合要求的371个肿瘤样本数据和50正常组织数据

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

data．category ＝ ＂Transcriptome Profiling＂，

data．type ＝ ＂Gene Expression Quantification＂，

workflow．type ＝ ＂HTSeq － Counts＂，

barcode ＝ c（dataSmTP， dataSmNT））

＃barcode参数：根据传入barcodes进行数据过滤

上图为 queryDown＜－GDCquery（）的结果，仅选择了选择371个正常组织和50个肿瘤组织样本。

第二步：GDCdownload（）下载GDCquery（）得到的结果

＃ 下载数据，默认存放位置为当前工作目录下的GDCdata文件夹中。

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

files．per．chunk ＝ 10）

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

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

＃files．per．chunk ＝ NULL：使用API下载大文件的时候，可以把文件分成几个小文件来下载，可以解决下载容易中断的问题。

GDCdownload（query ＝ queryDown）

说明：由于小编前面已经下载过该TCGA数据，所以这里显示的是421个文件已存在。如果还没有下载的话，可能需要根据自己的网速等待一些时间。

显示这样的结果，就算下载成功啦！文件默认保存在 Rstudio默认路径下的GDCdata中。前面就是我们利用第一期知识进行数据下载环节，权当温习功课吧——接下来我们就开始此期的数据处理～～

二、数据处理

第三步：GDCprepare（）将前面GDCquery（）的结果准备成R语言可处理的SE（SummarizedExperiment）文件。

＃读取下载的数据并将其准备到R对象中，在工作目录生成（save＝TRUE）LIHC＿case．rda文件

＃ GDCprepare（）：Prepare GDC data，准备GDC数据，使其可用于R语言中进行分析

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

＂LIHC＿case．rda＂）

GDCprepare（）中的参数：

参数用法query来自GDCquery的结果save是否将结果保存为RData object，默认为TRUEsave．filename文件名，如果没有设置，系统将默认设置directory文件数据的文件夹，默认为“GDCdata”summarizedExperiment是否生成summarizedExperiment对象，默认TRUE

第四步：TCGAanalyze＿Preprocessing（）对数据进行预处理：使用spearman相关系数去除数据中的异常值

＃ 去除dataPrep1中的异常值，dataPrep1数据中含有肿瘤组织和正常组织的数据

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

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

＃ 函数功能描述：Array Array Intensity correlation （AAIC） and correlation boxplot to define outlier

dataPrep2 ＜－ TCGAanalyze＿Preprocessing（object ＝ dataPrep1，

cor．cut ＝ 0．6，

datatype ＝ ＂HTSeq － Counts＂）

＃将预处理后的数据dataPrep2，写入新文件“LIHC＿dataPrep．csv”

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

这里将生成一个array－array intensity correlation（AAIC）相关性热图，如下：

TCGAanalyze＿Preprocessing（）中的参数：

参数用法object来自TCGAprepare的结果cor．cut设置阈值，根据样本中各个样本之间的spearman相关系数进行过滤。默认为0filename设置生成图片文件的名称，默认为PreprocessingOutput．pngwidth生成图片的宽度？？ height生成图片的高度datatype描述RangedSummarizedExperiment 数据类型的字符串

第五步：TCGAtumor＿purity（）筛选肿瘤纯度大于60％的肿瘤barcodes

＃ TCGAtumor＿purity（barcodes， estimate， absolute， lump， ihc， cpe），使用来自5种方法的5个估计值作为阈值对TCGA样本进行过滤，这5个值是estimate， absolute， lump， ihc， cpe，这里设置cpe＝0．6（cpe是派生的共识度量，是将所有方法的标准含量归一化后的均值纯度水平，以使它们具有相等的均值和标准差）

＃筛选肿瘤纯度大于等于60％的样本数据

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

＃ filtered 为被过滤的数据， pure＿barcodes是我们要的肿瘤数据

Purity．LIHC＜－purityDATA＄pure＿barcodes

normal．LIHC＜－purityDATA＄filtered

filtered 为被过滤的数据（为正常组织的数据barcodes）， pure＿barcodes是我们要的肿瘤样本barcodes。

第六步：将肿瘤表达矩阵与正常组织表达矩阵合并，进行基因注释

＃获取肿瘤纯度大于60％的340个肿瘤组织样本＋50个正常组织样本，共计390个样本

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

第七步：进行表达矩阵基因注释

＃基因注释，需要加载“SummarizedExperiment”包，“SummarizedExperiment container”每个由数字或其他模式的类似矩阵的对象表示。行通常表示感兴趣的基因组范围和列代表样品。

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

install．packages（＂BiocManager＂）

＃BiocManager：：install（＂SummarizedExperiment＂） ＃没有的需要执行下载代码

library（＂SummarizedExperiment＂）

rowData（dataPrep1） ＃传入数据dataPrep1必须为SummarizedExperiment对象

＃ 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

＃将结果写入文件“puried．LIHC．cancer．csv”

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

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

第八步：进行表达矩阵标准化和过滤，得到用于差异分析的表达矩阵

｀TCGAanalyze＿Normalization（）｀使用EDASeq软件包标准化mRNA转录本和miRNA。

＃TCGAanalyze＿Normalization（）执行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中的参数：

参数用法tabDFRNAseq表达矩阵，行代表基因，列代表样本geneInfo关于geneLength和gcContent的20531个基因的矩阵，“geneInfoHT”和“geneInfo”可选。method选择标准化的方法，基于’gcContent’ 或 ’geneLength’的标准化方法可选

＃将标准化后的数据再过滤，去除掉表达量较低（count较低）的基因，得到最终的数据

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（）中的参数：

参数用法tabDF数据框或者矩阵，行代表基因，列代表来自TCGA的样本method用于过滤较低count数的基因的方法，有’quantile’， ’varFilter’， ’filter1’， ’filter2’qnt．cut选择均值作为过滤的阈值

最后将过滤后的数据写入文件“TCGA＿LIHC＿final．csv”，就得到我们用于后续差异分析的表达文件：

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

＃保留的是390个样本（前340肿瘤，后50正常组织）

今天的数据预处理就讲到这里，接下来我们将分享：数据分析（差异表达分析、富集分析和聚类分析等）。如果你喜欢的话，就加入我们一起挖数据吧～～