突变信息瀑布图

生信喵

背景

用R代码画出文章PMID: 28622513里的oncoprint图，又叫瀑布图。

使用场景

展示多个样本中多个基因的突变情况，包括但不限于SNP、indel、CNV。

尤其是几十上百的大样本量的情况下，一目了然，看出哪个基因在哪个人群里突变多，以及多种突变类型的分布。

场景一：作为全基因组测序或全外显子组测序文章的第一个图。

场景二：展示感兴趣的癌症类型里，某一个通路的基因突变情况。

下载TCGA数据

需要下载突变信息和临床信息，作为输入数据。

#source("https://bioconductor.org/biocLite.R")
#biocLite("TCGAbiolinks")
library(TCGAbiolinks)

#查看癌症项目名缩写
cancerName<-data.frame(id = TCGAbiolinks:::getGDCprojects()$project_id, name = TCGAbiolinks:::getGDCprojects()$name)

cancerName前几行如下：

下载临床信息

参数 project = 后面写你要看的癌症名称缩写

query <- GDCquery(
    project = "TCGA-LIHC",
    data.category = "Clinical",
    data.format = "bcr xml"
)
GDCdownload(query)


cliquery <- GDCprepare_clinic(query, clinical.info = "patient")
colnames(cliquery)[1] <- "Tumor_Sample_Barcode"
str(cliquery)

cliquery前几行如下：

理论上，$符号后面的单词都可以作为分类信息画出来。

history_hepato_carcinoma_risk_factors里有50种，需要提取其中的Hepatitis B、Hepatitis C和cirrhosis。

cliquery$HBV <- 0
cliquery$HBV[grep(pattern = "Hepatitis B",cliquery$history_hepato_carcinoma_risk_factors)] <- 1

cliquery$HCV <- 0
cliquery$HCV[grep(pattern = "Hepatitis C",cliquery$history_hepato_carcinoma_risk_factors)] <- 1

cliquery$cirrhosis <- 0
cliquery$cirrhosis[grep(pattern = "cirrhosis",cliquery$history_hepato_carcinoma_risk_factors)] <- 1
table(cliquery$HBV)
#   0   1 
# 270 107 
table(cliquery$HCV)
# 0   1 
# 321  56 
table(cliquery$cirrhosis)
#   0   1 
# 373   4

后面我们将以 neoplasm_histologic_grade，gender，race_list，HCV，HBV，cirrhosis作为分类画图。

下载突变数据

之前的TCGA的MAF文件是可以下载的，每个癌症包含4种软件得到的突变文件：4种找mutation的pipeline任选其一：muse, varscan2, somaticsniper, mutect2。

后来就改版了，不让你随便下载了。但其实还是可以下载的，只不过没有那么多选择了！

现在的情况是每个样本都是一个单独的maf文件，需要下载后自己整理，就像整理表达矩阵那样。

但是现在我们有TCGAbiolinks，根本不需要自己动手，直接三步走即可得到我们需要的MAF文件。

snp <- GDCquery(
    project = "TCGA-LIHC", 
    data.category = "Simple Nucleotide Variation",
    data.type = "Masked Somatic Mutation",
    access = "open"
)
  
GDCdownload(snp)
GDCprepare(snp, save = T,save.filename = "TCGA-LIHC_SNP.Rda")

这样得到的这个Rda文件其实是一个数据框，不过由于内容和之前的MAF文件一模一样，所以也是可以直接用maftools读取使用的。

maftools是一个非常强大的突变数据可视化和分析的R包，这个包在bioconductor上，需要的自行安装。

开始画图

#source("https://bioconductor.org/biocLite.R")
#biocLite("maftools")
library(maftools)

load(file = "TCGA-LIHC_SNP.Rda")
maf.lihc <- data # 数据框
maf <- read.maf(maf = maf.lihc, clinicalData = cliquery, isTCGA = T)
# -Validating
# -Silent variants: 11113 
# -Summarizing
# --Possible FLAGS among top ten genes:
#   TTN
#   MUC16
# -Processing clinical data
# -Finished in 2.100s elapsed (2.020s cpu) 

data前几行如下：

用默认配色方案展示突变频率排名前20的基因

pdf("oncoplotTop20.pdf",width = 12,height = 12)
oncoplot(maf = maf,
         top = 20,
         fontSize = 0.7,legendFontSize = 1,annotationFontSize = 0.7,
         clinicalFeatures = c("race_list","neoplasm_histologic_grade","gender","HCV","HBV","cirrhosis")
         )
dev.off()

用默认配色方案展示一个列表里的基因

把你感兴趣的基因名存到 easy_input.txt文件里，每个基因为一行。就可以只画这几个基因，其中没发生突变的基因会被忽略。

geneName<-read.table("easy_input.txt",header = F,as.is = T)
genelist<-geneName$V1

pdf("oncoplotLIST.pdf",width = 12,height = 12)
oncoplot(maf = maf,
         genes = genelist,
         fontSize = 0.7,legendFontSize = 1,annotationFontSize = 0.7,
         clinicalFeatures = c("race_list","neoplasm_histologic_grade","gender","HCV","HBV","cirrhosis") 
         )
dev.off()

重新配色展示突变频率排名前20的基因

默认的颜色不够好看，还可以改配色

#突变
col = RColorBrewer::brewer.pal(n = 10, name = 'Paired')
names(col) = c('Frame_Shift_Del','Missense_Mutation', 'Nonsense_Mutation', 'Frame_Shift_Ins','In_Frame_Ins', 'Splice_Site', 'In_Frame_Del','Nonstop_Mutation','Translation_Start_Site','Multi_Hit')

#人种
racecolors = RColorBrewer::brewer.pal(n = 4,name = 'Spectral')
names(racecolors) = c("ASIAN", "WHITE", "BLACK_OR_AFRICAN_AMERICAN",  "AMERICAN_INDIAN_OR_ALASKA_NATIVE")

annocolors <- list(
  race_list = c(ASIAN = "#D7191C", 
                WHITE = "#FDAE61", 
                `BLACK OR AFRICAN AMERICAN` = "#ABDDA4", 
                `AMERICAN INDIAN OR ALASKA NATIVE` = "#2B83BA"),
  HCV = c(`0` = "grey", `1` = "black"),
  HBV = c(`0` = "grey", `1` = "black"),
  cirrhosis = c(`0` = "grey", `1` = "black"),
  # 添加缺失的临床特征颜色
  gender = c(Male = "blue", Female = "pink"),  # 名称必须与MAF中的性别标签完全一致
  neoplasm_histologic_grade = c(G1 = "green", G2 = "yellow", G3 = "red")  # 示例颜色
)

#画图
pdf("oncoplotTop20_col.pdf",width = 12,height = 12)
oncoplot(maf = maf,
         colors = col,#给突变配色
         annotationColor = annocolors,#给临床信息配色
         top = 20,
         fontSize = 0.7,legendFontSize = 1,annotationFontSize = 0.7,
         clinicalFeatures = c("race_list","neoplasm_histologic_grade","gender","HCV","HBV","cirrhosis"),
         writeMatrix =T)
dev.off()

生信喵

文中所用数据可以关注公众号“生信喵实验柴”

发送关键词“20250330”获取