根据公众号珠江肿瘤的推文 R语言|TCGAbiolinks 包系列（2）——肝癌案例之数据预处理 来学习

肝癌数据下载

第一步：筛选及下载我们想要的数据

#数据筛选
library("TCGAbiolinks")
query <- GDCquery(project = "TCGA-LIHC",
                  data.category = "Transcriptome Profiling",
                  data.type = "Gene Expression Quantification", 
                  workflow.type = "HTSeq - Counts")

samplesDown <- getResults(query,cols=c("cases"))  
#getResults(query, rows, cols)根据指定行名或列名从query中获取结果,此处用来获得样本的barcode
# 此处共检索出424个barcodes
# 从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")
###设置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进行数据过滤

# 下载数据，默认存放位置为当前工作目录下的GDCdata文件夹中。
GDCdownload(queryDown,method = "api", directory = "GDCdata",
            files.per.chunk = 10)
#method ；"API"或者"client"。"API"速度更快，但是容易下载中断。
#directory：下载文件的保存地址GDCdata。
#files.per.chunk = NULL:使用API下载大文件的时候，可以把文件分成几个小文件来下载，可以解决下载容易中断的问题。
GDCdownload(query = queryDown)

数据处理

第二步：读取文件

#读取下载的数据并将其准备到R对象中，在工作目录生成（save=TRUE）LIHC_case.rda文件
# GDCprepare():Prepare GDC data,准备GDC数据，使其可用于R语言中进行分析
dataPrep1 <- GDCprepare(query = queryDown, save = TRUE, save.filename =
                        "LIHC_case.rda")

image.png

第三步：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)相关性热图

image.png

上述参数的含义：

image.png

第四步：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

purityDATA

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

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

#获取肿瘤纯度大于60%的340个肿瘤组织样本+50个正常组织样本,共计390个样本
puried_data <-dataPrep2[,c(Purity.LIHC,normal.LIHC)]

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

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

puried_data

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

`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")

dataNorm

#将标准化后的数据再过滤，去除掉表达量较低（count较低）的基因，得到最终的数据
dataFilt <- TCGAanalyze_Filtering(tabDF = dataNorm,
                                                         method = "quantile", 
                                                          qnt.cut =  0.25)
str(dataFilt)
 #num [1:12985, 1:390] 4600 8157 10713 105 1621 ...
# - attr(*, "dimnames")=List of 2
 # ..$ : chr [1:12985] "A1BG" "A1CF" "A2M" "A4GALT" ...
#  ..$ : chr [1:390] "TCGA-QA-A7B7-01A-11R-A32O-07" "TCGA-DD-A39W-01A-11R-A213-07"     "TCGA-CC-5259-01A-31R-A213-07" "TCGA-CC-A7IF-01A-11R-A33J-07" ...

#保存文件
write.csv(dataFilt,file = "TCGA_LIHC_final.csv",quote = FALSE)  
#保留的是390个样本（前340肿瘤，后50正常组织）

右手的参数

本来想看一下dataFlit长什么样，结果电脑卡住了，暗暗下决心，有钱了要给电脑更新配置。加上家里的网速也真的是慢，跑了大半天才出结果，我太难了。但是学习之路是艰苦奋斗的，学到知识也是真的开心~晚上做生存曲线，嗯。