R数据可视化13:瀑布图/突变图谱
2019-12-06 本文已影响0人
jlyq617
突变相关分析的时候我们经常会选择瀑布图进行展示，瀑布图看起来十分的复杂高端，但是实际上只需要用一个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: 火山图
