2024-06-17

不可多得的GISTIC2.0教学之你不能错过的R包maftools结果可视化

大海哥生信果 2023-09-11 19:00:19

转自公众号：生信滩
http://mp.weixin.qq.com/s?__biz=MzkzMzUwODcwMg==&mid=2247485109&idx=1&sn=e48e66fe78e00023cafc7f9e9bd08aab

点击蓝字关注我们

小伙伴们大家好，我还是那个大海哥，今天大海哥要分享的还是接着前几天分享的GISTIC2.0教学，不过今天是这个教学的最后一篇！终于要完结了，希望你们可以把这个技术牢记于心哦！不说题外话了~今天大海哥要用到一个R包，名字叫做maftools，那么分析之前，我们先了解一下什么是maftools？

“maftools” 包提供了一系列用于分析和展示MAF文件中突变数据的功能，使研究人员能够更好地理解和解释肿瘤中的遗传变异。

以下是一些”maftools” 包提供的主要功能：

数据加载与处理：

“maftools” 允许用户加载MAF文件，并对数据进行必要的预处理和过滤，以便后续分析。

突变频率分析：用户可以使用”maftools” 包绘制突变的频率分布直方图、柱状图等，帮助他们了解每个基因的突变频率和分布情况。

突变类型分析：

该包能够绘制不同突变类型（例如错义突变、无义突变、插入、缺失等）的分布图，以及这些突变类型在不同样本中的分布情况。

基因级别分析：

用户可以通过”maftools” 包在基因级别上进行分析，包括绘制基因突变频率、突变谱、热图等，有助于了解哪些基因在样本中最常见的突变。

生存分析：

该包允许用户基于MAF数据进行生存分析，探究突变与患者生存率之间的关系。

复杂变异分析：

“maftools” 可以帮助用户识别并绘制出复杂变异模式，比如多基因突变、合并突变等。

突变注释可视化：

用户可以通过”maftools” 将突变注释信息（如功能影响、COSMIC注释等）与突变数据结合在一起进行可视化，从而更好地理解突变的功能和意义。

互作网络分析：

该包还支持构建基于突变数据的蛋白质互作网络，以揭示突变在细胞信号通路中的影响。

总之，”maftools” 是一个强大的工具，那么说了这么多，让我们开始吧！

#上次我们说到，通过在线分析，得到了下图中的一堆文件

#现在让我们来处理吧!library(tidyverse)# 建议下载github最新版本的maftools包library(maftools)#读入GISTIC输出的文件laml.gistic = readGistic(gisticAllLesionsFile = './Gistic 2.0 results/all_lesions.conf_99.txt',                   gisticAmpGenesFile = './Gistic 2.0 results/amp_genes.conf_99.txt',                    gisticDelGenesFile = './Gistic 2.0 results/del_genes.conf_99.txt',                    gisticScoresFile = './Gistic 2.0 results/scores.gistic', isTCGA = T)

#然后稍等片刻，这一步加载需要消耗一丢丢时间~#然后开始可视化重要的结果##绘图#genome plot 部分突出了明显的扩增和删失区域gisticChromPlot(gistic = laml.gistic, ref.build = 'hg38') ## 参考基因组选hg38

#Bubble plot 绘制显著改变的cytobands,以及改变的样本数和它所包含的基因数的函数。每一个气泡的大小是根据-log10转化q值。gisticBubblePlot(gistic = laml.gistic)

#oncoplot  gisticOncoPlot(gistic = laml.gistic,                sortByAnnotation = TRUE, top =20)

#我们还可以挑选出q值小于0.25的基因(sig_cytoband <- laml.gistic@cytoband.summary %>% filter(qvalues<0.25) %>% .$Unique_Name)#看看扩展基因有多少table(grepl(pattern = "AP",sig_cytoband))## ## FALSE  TRUE ##    34     3## 挑选q值小于0.25的Amplification cytobandsig_AP_cytoband <- laml.gistic@cytoband.summary %>%   filter(qvalues<0.25) %>%   .$Unique_Name %>%   .[grepl(pattern = "AP",.)]## 将Amplification 的基因挑选出来sig_AP_gene <- laml.gistic@data %>% filter(Cytoband %in% sig_AP_cytoband) %>% .$Hugo_Symbol unique(sig_AP_gene)#同样的挑选出丢失基因## 将Deletion cytoband 的基因挑选出来sig_DP_gene <- laml.gistic@data %>% filter(Cytoband %in% sig_DP_cytoband) %>% .$Hugo_Symbol unique(sig_DP_gene)#找到了基因之后就可以“为所欲为了….”拷贝数缺失区域的基因那么多，来个GO、kegg、GSEA这些分析是不是都可以。#简单富集分析一下library(clusterProfiler)library(org.Hs.eg.db)library(topGO)sig_DP_entrezId = mapIds(x = org.Hs.eg.db,                       keys = unique(sig_DP_gene),                       keytype = "SYMBOL",                       column = "ENTREZID")sig_DP_entrezId <- na.omit(sig_DP_entrezId)
go_ALL <- enrichGO(gene = sig_DP_entrezId,                       OrgDb = org.Hs.eg.db,                       keyType = "ENTREZID",                       ont = "ALL",                       pvalueCutoff = 0.5,                       qvalueCutoff = 0.5)

##分析完成后，作图### 数目太多，我只要各个条目的前6个BP_df<- go_ALL@result %>% arrange(desc(Count)) %>% filter(ONTOLOGY=="BP") %>% head()CC_df<- go_ALL@result %>% arrange(desc(Count)) %>% filter(ONTOLOGY=="CC") %>% head()MF_df<- go_ALL@result %>% arrange(desc(Count)) %>% filter(ONTOLOGY=="MF") %>% head()go_df <- rbind(BP_df,CC_df,MF_df)go_df$go_term_order=factor(x = c(1:nrow(go_df)),labels = go_df$Description)library(ggsci)ggplot(data=go_df, aes(x=go_term_order,y=Count, fill=ONTOLOGY)) +  geom_bar(stat="identity", width=0.8)  +   scale_fill_jama() +   theme_classic() +  xlab("GO term") + ylab("Number of Genes") + labs(title = "The Most Enriched GO Terms")+   theme(axis.text.x=element_text(face = "bold", color="gray50",angle = 60,vjust = 1, hjust = 1 ))