生信代码：数据预处理（TCGAbiolinks包）

科研菌

发布于 2021-01-12 11:54:18

6.9K1

发布于 2021-01-12 11:54:18

文章被收录于专栏：科研菌

引言：在前面我们了解了如何使用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中的参数详解：

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

补充TCGA中的组织样本类型：

TP	PRIMARY SOLID TUMOR	TM	Metastatic
TR	RECURRENT SOLID TUMOR	TAM	Additional Metastatic
TB	Primary Blood Derived Cancer-Peripheral Blood	THOC	Human Tumor Original Cells
TRBM	Recurrent Blood Derived Cancer-Bone Marrow	TBM	Primary Blood Derived Cancer-Bone Marrow
TAP	Additional-New Primary	NB	Blood Derived Normal
NT	Solid Tissue Normal	NBC	Buccal Cell Normal‍‍‍
NEBV	EBV Immortalized Normal	NBM	Bone 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，默认为TRUE
save.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相关系数进行过滤。默认为0
filename	设置生成图片文件的名称，默认为PreprocessingOutput.png
width	生成图片的宽度‍‍
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中的参数：

参数	用法
tabDF	RNAseq表达矩阵，行代表基因，列代表样本
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正常组织）

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

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2021-01-06，如有侵权请联系 cloudcommunity@tencent.com 删除

数据处理

本文分享自科研菌微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

数据处理

登录后参与评论

0 条评论

热度