GenVisR绘制瀑布图/突变图谱

一 R包的下载

下载地址：http://bioconductor.org/packages/release/bioc/html/GenVisR.html

学习过程中参考的两篇文章（感谢）：

https://www.jianshu.com/p/726310b02a56

http://www.jintiankansha.me/t/PyrlMUw0di

我下载的是R3.5.3-64，只要R3.3.0以上即可。然后打开R，运行下面两步：

都是在windows上面操作的，linux画图还得弄交互器，不然展示不方便。

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

biocLite("GenVisR") ####需要很长时间，耐心等待。

然后他可能会提示：Update all/some/none? [a/s/n]: 输入a即可（全部更新）

然后调用一下，安装成功！

library(GenVisR)

二瀑布图的绘制

输入的文件是一个表格，包含了总共23个样品的突变信息，MGI格式，如下：

其中必须用到的只有三列，即，样品名称，基因名，aa突变位点类别。

读取文件（指定1,2,8,12,13这几列）：

第一列，样品名；第二列，染色体名；第八列，基因名；12列，突变区域；13列，aa突变情况。

mutation_data=read.table(file="E:/liufei/R/total_sample_plot", header=TRUE, sep="\t")

mutation=mutation_data[,c(1,2,8,12,13)]

设定列名：

colnames(mutation)[c(1,2,3,4,5)]=c("sample","chr","gene_name","trv_type","amino.acid.change")

画图：

这里必须要加fileType="MGI"，否则R会认为输入是一个MAF格式（默认），导致无法绘图。

铛当当当~画出来啦~

waterfall(mutation, fileType="MGI")

我们再来把图形充实一下：

列出X轴（mainXlabel=TRUE，默认是FALSE）：

waterfall(mutation, fileType="MGI", mainXlabel=TRUE)

方格中显示氨基酸替换，并设置字体大小为1.8。

waterfall(mutation, fileType="MGI", mainXlabel=TRUE, mainLabelCol="amino.acid.change", mainLabelSize=1.8)

还可以设置方格中为染色体名称：

waterfall(mutation, fileType="MGI", mainXlabel=TRUE, mainLabelCol="chr", mainLabelSize=1.8)

指定基因和样品作图：

waterfall(mutation, fileType="MGI", mainXlabel=TRUE, mainLabelCol="amino.acid.change", mainLabelSize=1.8,plotSamples=c("1801167_N_T","1801169_N_T","1801942_N_T","1801943_N_T","1801944_N_T"),plotGenes=c("MACF1","NBPF10","LINC00842","CLCNKA","SRGAP2B","OR2T10","C2orf71","HEATR5B","ANKRD36C","PCDP1"))

指定基因和样品并按照设定的基因顺序和样品顺序作图：

waterfall(mutation, fileType="MGI", mainXlabel=TRUE, mainLabelCol="amino.acid.change", mainLabelSize=1.8,plotSamples=c("1801167_N_T","1801169_N_T","1801942_N_T","1801943_N_T","1801944_N_T"),plotGenes=c("MACF1","NBPF10","LINC00842","CLCNKA","SRGAP2B","OR2T10","C2orf71","HEATR5B","ANKRD36C","PCDP1"),sampOrder=c("1801167_N_T","1801169_N_T","1801942_N_T","1801943_N_T","1801944_N_T"),geneOrder=c("MACF1","NBPF10","LINC00842","CLCNKA","SRGAP2B","OR2T10","C2orf71","HEATR5B","ANKRD36C","PCDP1"))

设置该突变在样品中出现频率的阈值，低于此频率的将被过滤掉。

waterfall(mutation, fileType="MGI", mainXlabel=TRUE, mainLabelCol="amino.acid.change", mainLabelSize=1.8,mainRecurCutoff = 0.5)

基本上主要能用到的示例都列举了，如果还需要其他情况示例，查看帮助文档。

下面介绍一下GenVisR主要常用的参数。

三 GenVisR参数解析

1.fileType

fileType的格式可以有三种：MAF（默认）；MGI； Custom

◾若读入格式为MAF，则必须包含以下几列信息："Tumor_Sample_Barcode", "Hugo_Symbol", "Variant_Classification"；

◾若读入格式为MGI，则必须包含以下几列信息："sample","gene_name","trv_type"

◾若读入格式为自定义格式（Custom），则必须包含以下几列信息："sample", "gene", "variant_class"

详细请看：下载后R包中自带的介绍（GenVisR_intro.html）。

2.mainXlabel

mainXlabel：布尔型，是否显示X轴的sample name；默认为FALSE

mainXlabel=TRUE 表示列出X轴的样品名。

3.mainLabelCol

mainLabelCol设置方格中显示内容，可以是氨基酸变化情况（ mainLabelCol="amino.acid.change"），也可以是染色体情况（ mainLabelCol="chr"）。需要注意的是，这里mainLabelCol=后面跟的是列名。

4.mainLabelSize

mainLabelSize设置方格中字体的大小，如果样品数目较多，可以设置小一些。

5.plotSamples

plotSamples是指定绘制的样品名（中将图形），如果不指定，会将全部样品都绘制。如果指定，图形中只出现指定样品名的突变情况。

6.sampOrder

sampOrder指定样品名（X轴）排列顺序，从左到右。

7.plotGenes

plotGenes是指定绘制的基因名称（左侧图形）

8.geneOrder

geneOrder是指定基因的排序，从上到下。

9.clinData

clinData包含临床信息的data frame，默认为NULL；若进行自定义，列名必须为 "sample", "variable","value";此外，clinData 必须以”long” format的形式存在。

其他的参数可以参见一篇博客，介绍的很详细：

http://www.jintiankansha.me/t/PyrlMUw0di

下面主要介绍一下这个图形各个部分的含义.

四图形解析

图形主要分为上/左/中三个部分：

上：

每个样品对应的总体基因的突变密度。

红色是同义突变。

蓝色是非同义突变。

该图最上面的panel表示mutations in sample/coverage space∗1,000,000，也被称为Mutation Burden，这个coverage space可以自定义（默认是44100000），这里规定每个样本的coverage space是相同的，如果需要的话可以自定义每个样本的coverage space进行计算，不想显示这部分的话设置plotMutBurden = FALSE。

左：

特定基因（这里是mainRecurCutoff大于0.5的基因，自己可以设定mainRecurCutoff）在所有样本中的突变频率。

左侧panel表示每个突变位点在样本中出现的频率，该值为0.5表示有一般样本携带此突变。

中：

特定基因（这里是mainRecurCutoff大于0.5的基因，自己可以设定mainRecurCutoff）在对应特定样本的突变类型。

中间panel表示每个样本每个基因携带的突变信息，颜色表示该样本指定基因的突变类型（总共18种）。

五增加临床信息

在收集样本的时候，往往也会记录样本的情况，比如样本的年龄和得的癌症的类型等等。

这里我们来生成一个临床数据，我们以乳腺癌的4种分子分型LumA，LumB，HER2和Basal加上normal（正常，对照）作为我们样本的疾病情况，并把病人年龄分为20-30、31-50、51-60、和61+四个阶段。

读入临床数据：

clinical=read.table(file="E:/liufei/R/clin", header=TRUE, sep="\t")

这一步是官方给定的，是为了生成随机的clinical（样品名和突变文件mutation一致）。

但是我刚刚给定了clinical，所以这几步可以忽略。

subtype <- c("lumA", "lumB", "her2", "basal", "normal")
subtype <- sample(subtype, 23, replace = TRUE)
age <- c("20-30", "31-50", "51-60", "61+")
age <- sample(age, 23, replace = TRUE)
sample <- as.character(unique(mutation$sample))
clinical <- as.data.frame(cbind(sample, age,subtype))

clinical 是wide格式

wide格式转化成long格式：

用clinical的melt即可。

library(reshape2)

clinical <- melt(clinical, id.vars = c("sample"))

这时候，clinical变成了long格式。

画图：

waterfall(mutation, fileType="MGI",clinDat = clinical, clinVarCol = c(lumA = "blue4", lumB = "deepskyblue", her2 = "hotpink2", basal = "firebrick2", normal = "green4", `20-30` = "#ddd1e7", `31-50` = "#bba3d0", `51-60` = "#9975b9", `61+` = "#7647a2"), mainRecurCutoff = 0.05, maxGenes = 25, clinLegCol = 2, clinVarOrder = c("lumA", "lumB", "her2", "basal", "normal", "20-30", "31-50", "51-60", "61+"),mainXlabel=TRUE,mainLabelCol="amino.acid.change", mainLabelSize=2)

这样，图形下方就展示出临床信息啦。

再来看一下这里面几个参数：

clinLegCol：临床信息legend分为多少列。只有当clinData存在时才有效；

clinVarOrder：按指定顺序plot临床信息；

clinVarCol：根据clinData中的”variable”指定颜色；

如果参数改成：

waterfall(mutation, fileType="MGI",clinDat = clinical,clinVarOrder = c("lumA", "lumB", "her2", "basal", "normal", "20-30", "31-50", "51-60", "61+"),clinLegCol = 2)