|
上期学习了如何从官方下载最新版的TCGA突变数据,以及突变数据的合并及可视化,详情可戳:
《一文拿捏住突变数据下载、合并、可视化!》
今天我们继续学习TCGA突变数据分析(结合TCGA临床数据),同时对上期批量读取时的报错问题做一个补充。
一. 临床信息处理
- ######临床信息关联:
- #加载往期整合的临床信息矩阵(来自TCGA-KIRC):
- load("KIRC_cl.Rdata")
- View(cl)
复制代码
- #为了关联临床数据和突变数据:将列名"submitter_id"修改为"Tumor_Sample_Barcode";
- colnames(cl)[1]
复制代码
- names(cl)[names(cl) == "submitter_id"] <- "Tumor_Sample_Barcode"
- colnames(cl)[1]
复制代码
二. 突变数据处理
- ######突变数据的合并和读取:
- #批量读入所有maf文件名:
- mafFilePath <- dir(path = "KIRC-maf/",
- pattern = "masked.maf.gz$",
- full.names = T,
- recursive=T)
- head(mafFilePath)
复制代码
- #read.maf函数读取每一个maf文件(来自上期推文):
- mafdata <- lapply(
- mafFilePath,
- function(x){
- read.maf(x,
- clinicalData = cl,#放入临床数据
- isTCGA=TRUE)#如果逻辑值为TRUE,会截取前"Tumor_Sample_Barcode"前12个字符串和Mmaf文件进行匹配;
- }
- )
复制代码
这里可能会出现error,导致读取中断,上期我们没有讲如何处理报错问题,本期提供一种可以跳过error的方法供大家参考:
- #使用try()函数跳过error:
- ##在循环或批量读取中,如果遇到error会立刻终止,且不会返回结果,而try函数允许在error后继续执行代码;
- ##将maf的读取函数read.maf()放在try()中,可以跳过报错,参数silent=TRUE可以选择隐藏/不打印error信息(是否隐藏都不影响结果);
- mafdata <- lapply(
- mafFilePath,
- function(x){
- try(read.maf(x,
- clinicalData = cl,
- isTCGA=TRUE),
- silent=TRUE)}
- )
- View(mafdata)
复制代码
观察lappy后maf文件组成的list,提示“try-error”的子集即为读取时出现报错的maf文件。下面我们需要去掉list中出现“try-error”的子集,并将剩下的maf进行合并。
- #剔除error:
- ##由于在R中没有识别error的class,先自定义一个可以识别错误类型“try-error”的函数;
- is.error <- function(x) inherits(x,"try-error")
- right <- !sapply(mafdata,is.error)
- right[1:70]#简单看一下,逻辑值FALSE为mafdata2中出现error的位置;
复制代码
- Mafdata2 <- mafdata[right] #逻辑值为TRUE的位置为我们所需的;
- #合并maf:
- maf <- merge_mafs(mafdata2)
- maf
复制代码
- save(maf,file = c("KIRC_maf(cl).Rdata"))
复制代码
关于批量读取时报错问题的解决介绍到这里,下面接着上期内容,学习maftools这个包还能完成哪些突变数据的分析。
三. 突变数据分析
1.突变互斥性(exclusive)和共线性(co-occurrence)分析
- ##使用Fisher检验检测重要基因对
- somaticInteractions(maf = maf,
- top = 25,#展示队列中突变前25个基因间的互作
- pvalue = c(0.05, 0.1))
复制代码
- #也可以自定义互作test的基因list:
- gene_list <- c("VHL","PBRM1","TTN","SETD2")
- somaticInteractions(maf = maf, genes = gene_list, pvalue = c(0.05, 0.1))
复制代码
2.瀑布图(oncoplot)补充
- #瀑布图(oncoplot):
- oncoplot(maf = maf, top = 10) #上期瀑布图简单介绍到这里;
复制代码
- #自定义感兴趣的基因:
- genes <- c("VHL","TTN","BAP1","MTOR","ANK3","PCLO")
- oncoplot(maf = maf, genes = genes)
复制代码
- #更改配色:
- vc_cols <- RColorBrewer::brewer.pal(n = 8, name = 'Paired')
- names(vc_cols) <- c(
- 'In_Frame_Ins',
- 'Translation_Start_Site',
- 'Nonstop_Mutation',
- 'Frame_Shift_Ins',
- 'Frame_Shift_Del',
- 'Splice_Site',
- 'In_Frame_Del',
- 'Missense_Mutation'
- )
- oncoplot(maf = maf, colors = vc_cols, top = 5)
复制代码
由于本期我们在读取maf时合并了临床信息,因此有需求的话可以在瀑布图中添加样本注释。
- #添加样本注释:
- #先查看临床信息:
- getClinicalData(x = maf)
复制代码
- #选取注释信息(以癌症分期为例):
- oncoplot(maf = maf, top = 5, clinicalFeatures = 'tumor_stage')
复制代码
- #可以包含多种样本注释(同时添加癌症分期、性别、生死):
- oncoplot(maf = maf, top = 5, clinicalFeatures = c('tumor_stage','gender','status'))
复制代码
3.癌症驱动基因检测
- #基于oncodriveCLUST算法;
- #原理:致癌基因中的大多数变体都在少数特定位点(热点)富集;
- sig <- oncodrive(maf = maf,
- AACol = "HGVSp_Short",
- minMut = 5,
- pvalMethod = 'zscore')
- head(sig)
复制代码
- #散点图展示:
- plotOncodrive(res = sig,
- fdrCutOff = 0.1,
- useFraction = TRUE,
- labelSize = 0.5)
复制代码
散点大小和基因内突变的cluster成正比,x轴显示为在这些cluster中观察到的突变数量(或突变比例),由于预测的两个癌症驱动基因cluster scores相近,在图中挤在一起了。
4.pfam注释/统计
- pfam <- pfamDomains(maf = maf,
- AACol = 'HGVSp_Short',
- top = 10)
- #将pfam 结构域信息添加到氨基酸变化中;
- #圆点代表不同pfam结构域,X轴为出现的突变数量,Y轴以及圆点的大小为影响的基因数;
复制代码
- #Protein summary
- pfam$proteinSummary[,1:7, with = FALSE]
复制代码
- #Domain Summary
- pfam$domainSummary[,1:3, with = FALSE]
复制代码
5.生存分析
通常我们会使用临床特征或基因表达量数据进行生存分析,当然,有突变数据在手也可以将基因突变信息分组;导入的临床数据要确保至少要包含Tumor_Sample_Barcode、二分类事件(1/0) 和时间事件这三列。
- #探究任意基因突变(以BAP1为例):
- mafSurvival(maf = maf,
- genes = 'BAP1', #目标基因
- time = 'OS_time',#生存时间
- Status = 'event', #二分类生死事件(0/1)
- isTCGA = TRUE)
复制代码
可以看到发生BAP1基因突变的病人生存率显著降低。
- #预测生存相关基因集:
- prog_geneset <- survGroup(maf = maf,
- top = 20, #使用前20突变基因
- geneSetSize = 2, #一组的基因数为2
- time = "OS_time",
- Status = "event",
- verbose = FALSE)
- head(prog_geneset)#显示了低存活率相关的基因组合
复制代码
- mafSurvGroup(maf = maf,
- geneSet = c("VHL", "BAP1"),
- time = "OS_time",
- Status = "event")
复制代码
6.比较两个队列(两个maf)
有些时候我们想要比较两个队列的突变数据,比如原发性和复发性 APL的基因突变程度是有差异的, PML和RARA 基因突变发生在复发性APL中,在原发性APL中没有。
我们以KIRC的癌症分期为例,将癌症期和3-41-2期分为两个队列检测差异突变基因。
- #分别提取两个队列Tumor_Sample_Barcode:
- table(cl$tumor_stage)
- cl_12 <- subset(cl, tumor_stage%in%c("i","ii"))$Tumor_Sample_Barcode
- cl_34 <- subset(cl, tumor_stage%in%c("iii","iv"))$Tumor_Sample_Barcode
- #构建两个队列的maf对象:
- maf_cl12 <- subsetMaf(maf = maf, tsb = cl_12, isTCGA = TRUE)
- maf_cl34 <- subsetMaf(maf = maf, tsb = cl_34, isTCGA = TRUE)
- #对两个队列间的所有基因进行Fisher检验(检测差异突变基因):
- comp <- mafCompare(m1 = maf_cl12,
- m2 = maf_cl34,
- m1Name = "stagei&ii",
- m2Name = "stageiii&iv",
- minMut = 5)
- head(comp$results)
复制代码
可以看到两个队列基因发生突变的数量、显著性P值等信息。
- #森林图可视化队列差异突变基因:
- forestPlot(mafCompareRes = comp, pVal = 0.05)
复制代码
- #瀑布图可视化队列:
- genes <- c("BAP1", "ANKS1B", "ZFPM2", "TENM1", "STAG2","RYR3")
- coOncoplot(m1 = maf_cl12,
- m2 = maf_cl34,
- m1Name = 'stagei&ii',
- m2Name = 'stageiii&iv',
- genes = genes,
- removeNonMutated = TRUE) #瀑布图中隐藏无突变样本
复制代码
- #条形图:
- coBarplot(m1 = maf_cl12,
- m2 = maf_cl34,
- m1Name = 'stagei&ii',
- m2Name = 'stageiii&iv',
- genes = genes)
复制代码
7.临床富集分析(基于不同临床特征)
- fab.ce <- clinicalEnrichment(maf = maf,
- clinicalFeature = 'tumor_stage')
- head(fab.ce$groupwise_comparision)
复制代码
- plotEnrichmentResults(enrich_res = fab.ce,
- pVal = 0.03,
- geneFontSize = 0.5,
- annoFontSize = 0.6)
复制代码
8.药物/基因互作分析
- #通过Drug Gene Interaction database分析基因成药性:
- dgi <- drugInteractions(maf = maf, fontSize = 0.75)
- #显示了潜在的可成药基因类别以及其涉及的前5个基因;
复制代码
- #提取已报道/已知的药物-基因互作信息(输入感兴趣基因):
- dnmt3a.dgi <- drugInteractions(genes = "DNMT3A", drugs = TRUE)
- dnmt3a.dgi[,.(Gene, interaction_types, drug_name, drug_claim_name)]
复制代码
9.致癌信号通路分析(针对TCGA中已知的10个致癌信号通路)
- OncogenicPathways(maf = maf)
复制代码
- #瀑布图展示单个pathway:
- PlotOncogenicPathways(maf = maf,
- pathways = "Hippo",
- tsgCol = "red", #红色标注为抑癌基因
- ogCol = "royalblue") #蓝色标注为癌基因
复制代码
好啦,今天的分享就到这里!
参考资料:
https://www.bioconductor.org/pac ... html#1_Introduction
本文作者:基迪奥-喵酱
|
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有帐号?立即注册
x
|