查看: 244|回复: 6

[技术前沿] 一文拿捏住突变数据处理、可视化及分析!

[复制链接]

管理员

Rank: 15Rank: 15Rank: 15Rank: 15Rank: 15

主题
666
注册时间
2020.6.16
在线时间
386 小时

发表于 4 天前 | 显示全部楼层 |阅读模式
上期学习了如何从官方下载最新版的TCGA突变数据,以及突变数据的合并及可视化,详情可戳:
《一文拿捏住突变数据下载、合并、可视化!》

今天我们继续学习TCGA突变数据分析(结合TCGA临床数据),同时对上期批量读取时的报错问题做一个补充。

  1. #R包载入:
  2. library(maftools)
复制代码

一. 临床信息处理

  1. ######临床信息关联:
  2. #加载往期整合的临床信息矩阵(来自TCGA-KIRC):
  3. load("KIRC_cl.Rdata")
  4. View(cl)
复制代码


  1. #为了关联临床数据和突变数据:将列名"submitter_id"修改为"Tumor_Sample_Barcode";
  2. colnames(cl)[1]
复制代码


  1. names(cl)[names(cl) == "submitter_id"] <- "Tumor_Sample_Barcode"
  2. colnames(cl)[1]
复制代码


二. 突变数据处理

  1. ######突变数据的合并和读取:
  2. #批量读入所有maf文件名:
  3. mafFilePath <- dir(path = "KIRC-maf/",
  4. pattern = "masked.maf.gz$",
  5. full.names = T,
  6. recursive=T)
  7. head(mafFilePath)
复制代码


  1. #read.maf函数读取每一个maf文件(来自上期推文):
  2. mafdata <- lapply(
  3. mafFilePath,
  4. function(x){
  5. read.maf(x,
  6. clinicalData = cl,#放入临床数据
  7. isTCGA=TRUE)#如果逻辑值为TRUE,会截取前"Tumor_Sample_Barcode"前12个字符串和Mmaf文件进行匹配;
  8. }
  9. )
复制代码


这里可能会出现error,导致读取中断,上期我们没有讲如何处理报错问题,本期提供一种可以跳过error的方法供大家参考:

  1. #使用try()函数跳过error:
  2. ##在循环或批量读取中,如果遇到error会立刻终止,且不会返回结果,而try函数允许在error后继续执行代码;
  3. ##将maf的读取函数read.maf()放在try()中,可以跳过报错,参数silent=TRUE可以选择隐藏/不打印error信息(是否隐藏都不影响结果);
  4. mafdata <- lapply(
  5. mafFilePath,
  6. function(x){
  7. try(read.maf(x,
  8. clinicalData = cl,
  9. isTCGA=TRUE),
  10. silent=TRUE)}
  11. )
  12. View(mafdata)
复制代码


观察lappy后maf文件组成的list,提示“try-error”的子集即为读取时出现报错的maf文件。下面我们需要去掉list中出现“try-error”的子集,并将剩下的maf进行合并。

  1. #剔除error:
  2. ##由于在R中没有识别error的class,先自定义一个可以识别错误类型“try-error”的函数;
  3. is.error <- function(x) inherits(x,"try-error")
  4. right <- !sapply(mafdata,is.error)
  5. right[1:70]#简单看一下,逻辑值FALSE为mafdata2中出现error的位置;
复制代码


  1. Mafdata2 <- mafdata[right] #逻辑值为TRUE的位置为我们所需的;

  2. #合并maf:
  3. maf <- merge_mafs(mafdata2)
  4. maf
复制代码


  1. save(maf,file = c("KIRC_maf(cl).Rdata"))
复制代码

关于批量读取时报错问题的解决介绍到这里,下面接着上期内容,学习maftools这个包还能完成哪些突变数据的分析。

三. 突变数据分析

1.突变互斥性(exclusive)和共线性(co-occurrence)分析

  1. ##使用Fisher检验检测重要基因对
  2. somaticInteractions(maf = maf,
  3. top = 25,#展示队列中突变前25个基因间的互作
  4. pvalue = c(0.05, 0.1))
复制代码



  1. #也可以自定义互作test的基因list:
  2. gene_list <- c("VHL","PBRM1","TTN","SETD2")
  3. somaticInteractions(maf = maf, genes = gene_list, pvalue = c(0.05, 0.1))
复制代码


2.瀑布图(oncoplot)补充

  1. #瀑布图(oncoplot):
  2. oncoplot(maf = maf, top = 10) #上期瀑布图简单介绍到这里;
复制代码


  1. #自定义感兴趣的基因:
  2. genes <- c("VHL","TTN","BAP1","MTOR","ANK3","PCLO")
  3. oncoplot(maf = maf, genes = genes)
复制代码


  1. #更改配色:
  2. vc_cols <- RColorBrewer::brewer.pal(n = 8, name = 'Paired')
  3. names(vc_cols) <- c(
  4. 'In_Frame_Ins',
  5. 'Translation_Start_Site',
  6. 'Nonstop_Mutation',
  7. 'Frame_Shift_Ins',
  8. 'Frame_Shift_Del',
  9. 'Splice_Site',
  10. 'In_Frame_Del',
  11. 'Missense_Mutation'
  12. )
  13. oncoplot(maf = maf, colors = vc_cols, top = 5)
复制代码


由于本期我们在读取maf时合并了临床信息,因此有需求的话可以在瀑布图中添加样本注释。

  1. #添加样本注释:
  2. #先查看临床信息:
  3. getClinicalData(x = maf)
复制代码


  1. #选取注释信息(以癌症分期为例):
  2. oncoplot(maf = maf, top = 5, clinicalFeatures = 'tumor_stage')
复制代码


  1. #可以包含多种样本注释(同时添加癌症分期、性别、生死):
  2. oncoplot(maf = maf, top = 5, clinicalFeatures = c('tumor_stage','gender','status'))
复制代码


3.癌症驱动基因检测

  1. #基于oncodriveCLUST算法;
  2. #原理:致癌基因中的大多数变体都在少数特定位点(热点)富集;
  3. sig <- oncodrive(maf = maf,
  4. AACol = "HGVSp_Short",
  5. minMut = 5,
  6. pvalMethod = 'zscore')
  7. head(sig)
复制代码


  1. #散点图展示:
  2. plotOncodrive(res = sig,
  3. fdrCutOff = 0.1,
  4. useFraction = TRUE,
  5. labelSize = 0.5)
复制代码


散点大小和基因内突变的cluster成正比,x轴显示为在这些cluster中观察到的突变数量(或突变比例),由于预测的两个癌症驱动基因cluster scores相近,在图中挤在一起了。

4.pfam注释/统计

  1. pfam <- pfamDomains(maf = maf,
  2. AACol = 'HGVSp_Short',
  3. top = 10)
  4. #将pfam 结构域信息添加到氨基酸变化中;
  5. #圆点代表不同pfam结构域,X轴为出现的突变数量,Y轴以及圆点的大小为影响的基因数;
复制代码


  1. #Protein summary
  2. pfam$proteinSummary[,1:7, with = FALSE]
复制代码


  1. #Domain Summary
  2. pfam$domainSummary[,1:3, with = FALSE]
复制代码


5.生存分析

通常我们会使用临床特征或基因表达量数据进行生存分析,当然,有突变数据在手也可以将基因突变信息分组;导入的临床数据要确保至少要包含Tumor_Sample_Barcode、二分类事件(1/0) 和时间事件这三列。

  1. #探究任意基因突变(以BAP1为例):
  2. mafSurvival(maf = maf,
  3. genes = 'BAP1', #目标基因
  4. time = 'OS_time',#生存时间
  5. Status = 'event', #二分类生死事件(0/1)
  6. isTCGA = TRUE)
复制代码


可以看到发生BAP1基因突变的病人生存率显著降低。

  1. #预测生存相关基因集:
  2. prog_geneset <- survGroup(maf = maf,
  3. top = 20, #使用前20突变基因
  4. geneSetSize = 2, #一组的基因数为2
  5. time = "OS_time",
  6. Status = "event",
  7. verbose = FALSE)
  8. head(prog_geneset)#显示了低存活率相关的基因组合
复制代码


  1. mafSurvGroup(maf = maf,
  2. geneSet = c("VHL", "BAP1"),
  3. time = "OS_time",
  4. Status = "event")
复制代码


6.比较两个队列(两个maf)

有些时候我们想要比较两个队列的突变数据,比如原发性和复发性 APL的基因突变程度是有差异的, PML和RARA 基因突变发生在复发性APL中,在原发性APL中没有。


我们以KIRC的癌症分期为例,将癌症期和3-41-2期分为两个队列检测差异突变基因。

  1. #分别提取两个队列Tumor_Sample_Barcode:
  2. table(cl$tumor_stage)
  3. cl_12 <- subset(cl, tumor_stage%in%c("i","ii"))$Tumor_Sample_Barcode
  4. cl_34 <- subset(cl, tumor_stage%in%c("iii","iv"))$Tumor_Sample_Barcode

  5. #构建两个队列的maf对象:
  6. maf_cl12 <- subsetMaf(maf = maf, tsb = cl_12, isTCGA = TRUE)
  7. maf_cl34 <- subsetMaf(maf = maf, tsb = cl_34, isTCGA = TRUE)

  8. #对两个队列间的所有基因进行Fisher检验(检测差异突变基因):
  9. comp <- mafCompare(m1 = maf_cl12,
  10. m2 = maf_cl34,
  11. m1Name = "stagei&ii",
  12. m2Name = "stageiii&iv",
  13. minMut = 5)
  14. head(comp$results)
复制代码


可以看到两个队列基因发生突变的数量、显著性P值等信息。

  1. #森林图可视化队列差异突变基因:
  2. forestPlot(mafCompareRes = comp, pVal = 0.05)
复制代码


  1. #瀑布图可视化队列:
  2. genes <- c("BAP1", "ANKS1B", "ZFPM2", "TENM1", "STAG2","RYR3")
  3. coOncoplot(m1 = maf_cl12,
  4. m2 = maf_cl34,
  5. m1Name = 'stagei&ii',
  6. m2Name = 'stageiii&iv',
  7. genes = genes,
  8. removeNonMutated = TRUE) #瀑布图中隐藏无突变样本
复制代码


  1. #条形图:
  2. coBarplot(m1 = maf_cl12,
  3. m2 = maf_cl34,
  4. m1Name = 'stagei&ii',
  5. m2Name = 'stageiii&iv',
  6. genes = genes)
复制代码


7.临床富集分析(基于不同临床特征)

  1. fab.ce <- clinicalEnrichment(maf = maf,
  2. clinicalFeature = 'tumor_stage')
  3. head(fab.ce$groupwise_comparision)
复制代码


  1. plotEnrichmentResults(enrich_res = fab.ce,
  2. pVal = 0.03,
  3. geneFontSize = 0.5,
  4. annoFontSize = 0.6)
复制代码


8.药物/基因互作分析

  1. #通过Drug Gene Interaction database分析基因成药性:
  2. dgi <- drugInteractions(maf = maf, fontSize = 0.75)
  3. #显示了潜在的可成药基因类别以及其涉及的前5个基因;
复制代码


  1. #提取已报道/已知的药物-基因互作信息(输入感兴趣基因):
  2. dnmt3a.dgi <- drugInteractions(genes = "DNMT3A", drugs = TRUE)
  3. dnmt3a.dgi[,.(Gene, interaction_types, drug_name, drug_claim_name)]
复制代码


9.致癌信号通路分析(针对TCGA中已知的10个致癌信号通路)

  1. OncogenicPathways(maf = maf)
复制代码


  1. #瀑布图展示单个pathway:
  2. PlotOncogenicPathways(maf = maf,
  3. pathways = "Hippo",
  4. tsgCol = "red", #红色标注为抑癌基因
  5. ogCol = "royalblue") #蓝色标注为癌基因
复制代码


好啦,今天的分享就到这里!

参考资料:
https://www.bioconductor.org/pac ... html#1_Introduction


本文作者:基迪奥-喵酱

本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
新的一天加油!
回复

使用道具 举报

迅猛龙

Rank: 8Rank: 8

主题
3
注册时间
2021.6.22
在线时间
54 小时

发表于 3 天前 | 显示全部楼层
值得学习
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
101 小时

发表于 3 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

钵水母

Rank: 3Rank: 3

主题
0
注册时间
2021.4.14
在线时间
14 小时

发表于 3 天前 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
101 小时

发表于 前天 21:19 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

功夫熊猫

Rank: 10Rank: 10Rank: 10

主题
0
注册时间
2020.11.21
在线时间
101 小时

发表于 昨天 12:04 | 显示全部楼层
新的一天加油!
回复

使用道具 举报

帝王蝶

Rank: 4

主题
0
注册时间
2016.10.5
在线时间
43 小时

发表于 昨天 14:29 | 显示全部楼层
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

快速回复 返回顶部 返回列表