从高通量基因组技术直观地可视化和解读数据仍然具有挑战性。R中的基因组可视化(GenVisR)试图通过提供高度可定制的出版物质量图形来减轻这一负担，支持多个物种，并主要关注cohort level(即多个样本/患者)。GenVisR试图保持高度的灵活性，同时利用ggplot2和bioconductor的能力来实现这一目标。

部分功能跟maftools的功能很相似。

安装：从Bioconductor安装：

if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("GenVisR")
library(GenVisR)

功能：

GenVisR主要有以下13个功能，我将用两个数据集去测试（一个从TCGA下载，一个是自己的内部数据）。

1. Waterfall (mutation overview graphic)
2. lolliplot (mutation hotspot graphic)
3. genCov (sequence coverage graphic)
4. TvTi (transition/transversion graphic)
5. cnSpec(copy altered cohort graphic)
6. cnView(copy altered single sample graphic)
7. covBars(sequencing coverage cohort)
8. cnFreq(proportional copy number alterations)
9. ideoView(ideogram graphic)
10. lohSpec(Loss of Heterozygosity Spectrum)
11. lohView(Loss of Heterozygosity View)
12. compldent(snp identity graphic)
13. geneViz(Transcript Represenation)

1. Waterfall plot

首先需要一个MAF文件或MGI文件作为输入， 文件中的突变类型或者叫“Variant_Classification”包含下列字段：

image.png

可以下载一个TCGA的某一个肿瘤类型的MAF文件看看，再用自己跑出来的单个样本的MAF文件试试。

使用TCGAbiolinks下载：随便选择一个肿瘤，就用COAD。
BiocManager::install("TCGAbiolinks")
library(TCGAbiolinks)

# 下载TCGA的MAF
maf_file = GDCquery_Maf("COAD",save.csv = T, directory = "GDCdata", pipelines = "varscan")
waterfall(maf_file, fileType = "MAF")

有报错：

image.png

报错内容是：该MAF中Variant_Classification中有不认识的字段。

看看到底是哪个不认识：

image.png

肉眼比较一下，发现TCGA结果中多出来了一个Splice_Region是GenVisR不认识的，看看属于Splice_Region的突变多不多呢？

image.png

还挺多，有1929个，删掉试试：

maf_file = subset(maf_file, Variant_Classification != "Splice_Region")
waterfall(maf_file, fileType = "MAF")

gene太多了导致左侧都堆在一起了。

image.png

下载的TCGA COAD肿瘤的MAF中有399个样本，再来看看我自己一个样本的MAF文件结果：

my_maf = read.csv("mysample.variants.funcotated.without.header.MAF.xls",header = T, sep = "\t")
waterfall(maf_file, fileType = "MAF")

image.png

虽然可以设置plotGenes来指定特定的基因，但好像不能想maftools那样指定top基因

# Plot only the specified genes 
waterfall(brcaMAF, plotGenes = c("PIK3CA", "TP53", "USH2A", "MLL3", "BRCA1"))

如果有临床信息（clinical data）可以通过一些代码附加上，并画出如下图：

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

# Melt the clinical data into 'long' format.
library(reshape2)
clinical <- melt(clinical, id.vars = c("sample"))

# Run waterfall
waterfall(brcaMAF, 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"), plotGenes = c("PIK3CA",
    "TP53", "USH2A", "MLL3", "BRCA1"), clinLegCol = 2, clinVarOrder = c("lumA", "lumB",
    "her2", "basal", "normal", "20-30", "31-50", "51-60", "61+"))

image.png

但是这里都没有Clinical Data,就不演示了。

GenVisR 基因组数据可视化实战(二)
GenVisR 基因组数据可视化实战(三)
GenVisR 基因组数据可视化实战 (四)