一文拿捏住突变数据处理、可视化及分析！

基迪奥-Jt桃

上期学习了如何从官方下载最新版的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. 
   
3. #合并maf：
   
4. maf <- merge_mafs(mafdata2)
   
5. 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. 
   
6. #构建两个队列的maf对象：
   
7. maf_cl12 <- subsetMaf(maf = maf, tsb = cl_12, isTCGA = TRUE)
   
8. maf_cl34 <- subsetMaf(maf = maf, tsb = cl_34, isTCGA = TRUE)
   
9. 
   
10. #对两个队列间的所有基因进行Fisher检验(检测差异突变基因):
    
11. comp <- mafCompare(m1 = maf_cl12,
    
12. m2 = maf_cl34,
    
13. m1Name = "stagei&ii",
    
14. m2Name = "stageiii&iv",
    
15. minMut = 5)
    
16. 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


本文作者：基迪奥-喵酱

yyn1996

值得学习

2012csc

Juncau

2012csc

2012csc

wang649547