• 发文
  • 评论
  • 微博
  • 空间
  • 微信

生信代码:数据预处理(TCGAbiolinks包)

科研菌 2021-01-07 09:32 发文

引言:在前面我们了解了如何使用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正常组织)

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

声明:本文为OFweek维科号作者发布,不代表OFweek维科号立场。如有侵权或其他问题,请及时联系我们举报。
2
评论

评论

    相关阅读

    暂无数据

    科研菌

    每天为您提供高质量的原创科研论文...

    举报文章问题

    ×
    • 营销广告
    • 重复、旧闻
    • 格式问题
    • 低俗
    • 标题夸张
    • 与事实不符
    • 疑似抄袭
    • 我有话要说
    确定 取消

    举报评论问题

    ×
    • 淫秽色情
    • 营销广告
    • 恶意攻击谩骂
    • 我要吐槽
    确定 取消

    用户登录×

    请输入用户名/手机/邮箱

    请输入密码