R数据可视化13:瀑布图/突变图谱

突变相关分析的时候我们经常会选择瀑布图进行展示，瀑布图看起来十分的复杂高端，但是实际上只需要用一个GenVisR的R包即可解决。今天就让我们来看一下如何做瀑布图吧。另外，因为本人不是做该方向的，如果有描述不准确的地方，请大家指正。

什么是瀑布图 Waterfall Plot

Wiki上介绍的瀑布图分为两种，一类是2D形式，另一类是3D形式。我们简单介绍一下2D形式的瀑布图。该类瀑布图用于描述一系列中间正值或负值如何影响初始值。通常，初始值和最终值（端点）由整列表示，而中间值则显示为基于上一列的值开始的浮动列。这些列可以用不同的颜色标记，以区分正值和负值。

传统瀑布图的例子

可以看到该例子展示了获利能力的分析。但是用于展示突变的瀑布图和传统的瀑布图并不完全一样，不过他们的展现形式很像。

SNP的瀑布图

在SNP的瀑布图中，横轴表示的是不同的样本，纵轴是基因，填充则代表该基因发生突变，不同的颜色代表不同的突变情况。上面的柱状图是对于每个样本突变情况的统计。

所以从该图我们既能够获得每个样本的具体信息，也能够从全局了解这一组样本的整体情况，很好地展示了突变的情况。

怎么做瀑布图

本次作图我们使用一个叫做GenVisR的R包,该包除了提供瀑布图还提供了多种其他形式较为复杂的、用于展现多个样本突变情况的数据图（见下图），具体的作图方法大家可以参考GenVisR使用手册。

GenVisR提供的可视化方式

今天我们要使用该包提供的一个叫做brcaMAF的数据表，通过名字也可以看出这是乳腺癌的数据，该数据包含50个样本，来源于TCGA，格式为MAF格式。

MAF格式是Mutation Annotation Format的缩写，是一个以制表符分隔的文本文件，其中聚合了来自VCF文件的突变信息，该文件格式标准由TCGA制定，包含了一些关于突变的常见信息，进一步的具体信息详见MAF格式介绍

1）需要什么格式的数据
我们首先来看一下brcaMAF数据的情况，可以看到该数据包括了55列信息，如Hugo_Symbol、Chromosome等等,一共观察到了2773个突变。

colnames(brcaMAF)
 [1] "Hugo_Symbol"                   "Entrez_Gene_Id"                "Center"                        "NCBI_Build"                    "Chromosome"                    "Start_Position"               
 [7] "End_Position"                  "Strand"                        "Variant_Classification"        "Variant_Type"                  "Reference_Allele"              "Tumor_Seq_Allele1"            
[13] "Tumor_Seq_Allele2"             "dbSNP_RS"                      "dbSNP_Val_Status"              "Tumor_Sample_Barcode"          "Matched_Norm_Sample_Barcode"   "Match_Norm_Seq_Allele1"       
[19] "Match_Norm_Seq_Allele2"        "Tumor_Validation_Allele1"      "Tumor_Validation_Allele2"      "Match_Norm_Validation_Allele1" "Match_Norm_Validation_Allele2" "Verification_Status"          
[25] "Validation_Status"             "Mutation_Status"               "Sequencing_Phase"              "Sequence_Source"               "Validation_Method"             "Score"                        
[31] "BAM_File"                      "Sequencer"                     "Tumor_Sample_UUID"             "Matched_Norm_Sample_UUID"      "chromosome_name_WU"            "start_WU"                     
[37] "stop_WU"                       "reference_WU"                  "variant_WU"                    "type_WU"                       "gene_name_WU"                  "transcript_name_WU"           
[43] "transcript_species_WU"         "transcript_source_WU"          "transcript_version_WU"         "strand_WU"                     "transcript_status_WU"          "trv_type_WU"                  
[49] "c_position_WU"                 "amino_acid_change_WU"          "ucsc_cons_WU"                  "domain_WU"                     "all_domains_WU"                "deletion_substructures_WU"    
[55] "transcript_error"      
head(brcaMAF)
  Hugo_Symbol Entrez_Gene_Id           Center NCBI_Build Chromosome Start_Position End_Position Strand Variant_Classification Variant_Type Reference_Allele Tumor_Seq_Allele1 Tumor_Seq_Allele2 dbSNP_RS dbSNP_Val_Status
1       A2ML1         144568 genome.wustl.edu         37         12        8994108      8994108      +      Missense_Mutation          SNP                G                 G                 C    novel                 
2       AADAC             13 genome.wustl.edu         37          3      151545656    151545656      +      Missense_Mutation          SNP                A                 A                 G    novel                 
3       AADAT          51166 genome.wustl.edu         37          4      170991750    170991750      +                 Silent          SNP                G                 G                 A    novel                 
4        AASS          10157 genome.wustl.edu         37          7      121756793    121756793      +      Missense_Mutation          SNP                G                 G                 A    novel                 
5        ABAT              0 genome.wustl.edu         37         16        8857982      8857982      +                 Silent          SNP                G                 G                 A    novel                 
6       ABCA3             21 genome.wustl.edu         37         16        2335631      2335631      +      Missense_Mutation          SNP                C                 T                 T    novel                 
          Tumor_Sample_Barcode  Matched_Norm_Sample_Barcode Match_Norm_Seq_Allele1 Match_Norm_Seq_Allele2 Tumor_Validation_Allele1 Tumor_Validation_Allele2 Match_Norm_Validation_Allele1 Match_Norm_Validation_Allele2
1 TCGA-A1-A0SO-01A-22D-A099-09 TCGA-A1-A0SO-10A-03D-A099-09                      G                      G                        G                        C                             G                             G
2 TCGA-A2-A0EU-01A-22W-A071-09 TCGA-A2-A0EU-10A-01W-A071-09                      A                      A                                                                                                              
3 TCGA-A2-A0ER-01A-21W-A050-09 TCGA-A2-A0ER-10A-01W-A055-09                      G                      G                                                                                                              
4 TCGA-A2-A0EN-01A-13D-A099-09 TCGA-A2-A0EN-10A-01D-A099-09                      G                      G                        G                        A                             G                             G
5 TCGA-A1-A0SI-01A-11D-A142-09 TCGA-A1-A0SI-10B-01D-A142-09                      G                      G                                                                                                              
6 TCGA-A2-A0D0-01A-11W-A019-09 TCGA-A2-A0D0-10A-01W-A021-09                      C                      C                                                                                                              
  Verification_Status Validation_Status Mutation_Status Sequencing_Phase Sequence_Source Validation_Method Score BAM_File      Sequencer                    Tumor_Sample_UUID             Matched_Norm_Sample_UUID
1             Unknown             Valid         Somatic         Phase_IV             WXS Illumina_WXS_gDNA     1    dbGAP Illumina GAIIx b3568259-c63c-4eb1-bbc7-af711ddd33db 17ba8cdb-e35b-4496-a787-d1a7ee7d4a1e
2             Unknown          Untested         Somatic         Phase_IV             WXS              none     1    dbGAP Illumina GAIIx de30da8f-903f-428e-a63d-59625fc858a9 1583a7c5-c835-44fa-918a-1448abf6533d
3             Unknown          Untested         Somatic         Phase_IV             WXS              none     1    dbGAP Illumina GAIIx 31ed187e-9bfe-4ca3-8cbb-10c1e0184331 2bc2fdaf-fb2f-4bfd-9e20-e20edff6633a
4             Unknown             Valid         Somatic         Phase_IV             WXS Illumina_WXS_gDNA     1    dbGAP Illumina GAIIx 12362ad7-6866-4e7a-9ec6-8a0a68df8896 ad478c68-a18b-4529-ad7a-86039e6da6b1
5             Unknown          Untested         Somatic         Phase_IV             WXS              none     1    dbGAP Illumina GAIIx e218c272-a7e1-4bc9-b8c5-d2d1c903550f fbcab9dc-4a6b-4928-9459-699c9932e3e1
6             Unknown          Untested         Somatic         Phase_IV             WXS              none     1    dbGAP Illumina GAIIx 3f20d0fe-aaa1-40f1-b2c1-7f070f93aef5 bbf1c43d-d7b3-4574-a074-d22ad537829c
  chromosome_name_WU  start_WU   stop_WU reference_WU variant_WU type_WU gene_name_WU transcript_name_WU transcript_species_WU transcript_source_WU transcript_version_WU strand_WU transcript_status_WU trv_type_WU
1                 12   8994108   8994108            G          C     SNP        A2ML1        NM_144670.3                 human              genbank                58_37c         1            validated    missense
2                  3 151545656 151545656            A          G     SNP        AADAC        NM_001086.2                 human              genbank                58_37c         1             reviewed    missense
3                  4 170991750 170991750            G          A     SNP        AADAT        NM_016228.3                 human              genbank                58_37c        -1             reviewed      silent
4                  7 121756793 121756793            G          A     SNP         AASS        NM_005763.3                 human              genbank                58_37c        -1             reviewed    missense
5                 16   8857982   8857982            G          A     SNP         ABAT        NM_000663.4                 human              genbank                58_37c         1             reviewed      silent
6                 16   2335631   2335631            C          T     SNP        ABCA3        NM_001089.2                 human              genbank                58_37c        -1             reviewed    missense
  c_position_WU amino_acid_change_WU ucsc_cons_WU                                                    domain_WU
1        c.1224              p.W408C        0.995                                                         NULL
2         c.896              p.N299S        0.000      HMMPfam_Abhydrolase_3,superfamily_alpha/beta-Hydrolases
3         c.708               p.L236        1.000 HMMPfam_Aminotran_1_2,superfamily_PLP-dependent transferases
4         c.788              p.T263M        1.000                                          HMMPfam_AlaDh_PNT_C
5         c.423               p.E141        0.987     HMMPfam_Aminotran_3,superfamily_PyrdxlP-dep_Trfase_major
6        c.3295             p.D1099N        0.980                                                         NULL
                                                                                                                                                                                                                                               all_domains_WU
1             HMMPfam_A2M,HMMPfam_A2M_N,superfamily_Terpenoid cyclases/Protein prenyltransferases,HMMPfam_A2M_recep,superfamily_Alpha-macroglobulin receptor domain,HMMPfam_A2M_N_2,HMMPfam_A2M_comp,HMMPfam_Thiol-ester_cl,PatternScan_ALPHA_2_MACROGLOBULIN
2                                                                                                                                                                         PatternScan_LIPASE_GDXG_SER,HMMPfam_Abhydrolase_3,superfamily_alpha/beta-Hydrolases
3                                                                                                                                                                                                HMMPfam_Aminotran_1_2,superfamily_PLP-dependent transferases
4 HMMPfam_Saccharop_dh,HMMPfam_AlaDh_PNT_C,HMMPfam_AlaDh_PNT_N,superfamily_NAD(P)-binding Rossmann-fold domains,superfamily_Formate/glycerate dehydrogenase catalytic domain-like,superfamily_Glyceraldehyde-3-phosphate dehydrogenase-like C-terminal domain
5                                                                                                                                                                    HMMPfam_Aminotran_3,PatternScan_AA_TRANSFER_CLASS_3,superfamily_PyrdxlP-dep_Trfase_major
6                                                                                                                            HMMPfam_ABC_tran,HMMSmart_SM00382,PatternScan_ABC_TRANSPORTER_1,superfamily_P-loop containing nucleoside triphosphate hydrolases
  deletion_substructures_WU transcript_error
1                         -        no_errors
2                         -        no_errors
3                         -        no_errors
4                         -        no_errors
5                         -        no_errors
6                         -        no_errors

那么我们的MAF文件也需要那么多信息吗？并非如此，和很多其他作图所需数据一样，其中有一些信息是必须提供的，另外一些是非必须的。
具体地分为三类情况，瀑布图地功能提供了三种数据格式的选择：
1.MAF
必须包括列名为"Tumor_Sample_Barcode","Hugo_Symbol","Variant_Classification"的信息
2.MGI
必须包括列名为"sample","gene_name","trv_type"的信息
3.Custom
必须包括列名为"sample", "gene", "variant_class"的信息

MGI也是一种以制表符分割的文本文件，具体的可以见链接MGI格式介绍

2）如何作图
waterfall函数有很多参数，可以根据需求展示突变信息，那么下面就来一步作图，我们展示几种常用的参数用途，其他更多具体参数的意义可以查看帮助?waterfall（本文的代码来源GenVisR官方手册）

# 最基本的作图
waterfall(brcaMAF, fileType="MAF")
# 展示至少在6%的样本中存在的突变
waterfall(brcaMAF, mainRecurCutoff = 0.06)
#特定基因的突变图谱
waterfall(brcaMAF, plotGenes = c("PIK3CA", "TP53", "USH2A", "MLL3", "BRCA1"))

waterfall

我们常常需要结合样本的临床信息分析突变情况，那么要怎样同时展示样本的临床信息呢？

#建立临床信息
#分组情况为5组
subtype <- c("lumA", "lumB", "her2", "basal", "normal")
#随机放回的分配组别
subtype <- sample(subtype, 50, replace = TRUE)
#年龄分为4组
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))
head(clinical)
                        sample subtype   age
1 TCGA-A1-A0SO-01A-22D-A099-09   basal 51-60
2 TCGA-A2-A0EU-01A-22W-A071-09   basal 31-50
3 TCGA-A2-A0ER-01A-21W-A050-09  normal   61+
4 TCGA-A2-A0EN-01A-13D-A099-09  normal 31-50
5 TCGA-A1-A0SI-01A-11D-A142-09    lumA 31-50
6 TCGA-A2-A0D0-01A-11W-A019-09  normal   61+

# Melt the clinical data into 'long' format.
library(reshape2)
#整理数据表
clinical <- melt(clinical, id.vars = c("sample"))
head(clinical)
                        sample variable  value
1 TCGA-A1-A0SO-01A-22D-A099-09  subtype  basal
2 TCGA-A2-A0EU-01A-22W-A071-09  subtype  basal
3 TCGA-A2-A0ER-01A-21W-A050-09  subtype normal
4 TCGA-A2-A0EN-01A-13D-A099-09  subtype normal
5 TCGA-A1-A0SI-01A-11D-A142-09  subtype   lumA
6 TCGA-A2-A0D0-01A-11W-A019-09  subtype normal

# 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+"))#设定顺序

添加临床信息的waterfall

自己的数据可以整理成MAF格式，也可以选择Custom格式，要注意的是MAF格式和MGI格式对mutation type的类别名字有固定要求，如果你的mutation命名方式或者有不在下列类型中的突变类型，请选择Custom类别，该作图方式对mutation type的类别名字没有限制。

MAF	MGI
Nonsense_Mutation	nonsense
Frame_Shift_Ins	frame_shift_del
Frame_Shift_Del	frame_shift_ins
Translation_Start_Site	splice_site_del
Splice_Site	splice_site_ins
Nonstop_Mutation	splice_site
In_Frame_Ins	nonstop
In_Frame_Del	in_frame_del
Missense_Mutation	in_frame_ins
5’Flank	missense
3’Flank	splice_region_del
5’UTR	splice_region_ins
3’UTR	splice_region
RNA	5_prime_flanking_region
Intron	3_prime_flanking_region
IGR	3_prime_untranslated_region
Silent	5_prime_untranslated_region
Targeted_Region	rna
intronic
silent

那么关于瀑布图的分享就到这里啦～

往期R数据可视化分享
R数据可视化12: 曼哈顿图
 R数据可视化11: 相关性图
 R数据可视化10: 蜜蜂图 Beeswarm
R数据可视化9: 棒棒糖图 Lollipop Chart
R数据可视化8: 金字塔图和偏差图
 R数据可视化7: 气泡图 Bubble Plot
R数据可视化6: 面积图 Area Chart
R数据可视化5: 热图 Heatmap
R数据可视化4: PCA和PCoA图
 R数据可视化3: 直方/条形图
 R数据可视化2: 箱形图 Boxplot
R数据可视化1: 火山图