#path to TCGA LAML MAF file
laml.maf = system.file('extdata', 'tcga_laml.maf.gz', package = 'maftools')
#clinical information containing survival information and histology. This is optional
laml.clin = system.file('extdata', 'tcga_laml_annot.tsv', package = 'maftools')
laml = read.maf(maf = laml.maf,
clinicalData = laml.clin,
verbose = FALSE)
#One can use any colors, here in this example color palette from RColorBrewer package is used
vc_cols = RColorBrewer::brewer.pal(n = 8, name = 'Paired')
names(vc_cols) = c(
'Frame_Shift_Del',
'Missense_Mutation',
'Nonsense_Mutation',
'Multi_Hit',
'Frame_Shift_Ins',
'In_Frame_Ins',
'Splice_Site',
'In_Frame_Del'
)
print(vc_cols)
#> Frame_Shift_Del Missense_Mutation Nonsense_Mutation Multi_Hit
#> "#A6CEE3" "#1F78B4" "#B2DF8A" "#33A02C"
#> Frame_Shift_Ins In_Frame_Ins Splice_Site In_Frame_Del
#> "#FB9A99" "#E31A1C" "#FDBF6F" "#FF7F00"
oncoplot(maf = laml, colors = vc_cols, top = 10)
There are two ways one include CN status into MAF. 1. GISTIC results 2. Custom copy number table
Most widely used tool for copy number analysis from large scale studies is GISTIC and we can simultaneously read gistic results along with MAF. GISTIC generates numerous files but we need mainly four files all_lesions.conf_XX.txt
, amp_genes.conf_XX.txt
, del_genes.conf_XX.txt
, scores.gistic
where XX is confidence level. These files contain significantly altered genomic regions along with amplified and deleted genes respectively.
#GISTIC results LAML
all.lesions =
system.file("extdata", "all_lesions.conf_99.txt", package = "maftools")
amp.genes =
system.file("extdata", "amp_genes.conf_99.txt", package = "maftools")
del.genes =
system.file("extdata", "del_genes.conf_99.txt", package = "maftools")
scores.gis =
system.file("extdata", "scores.gistic", package = "maftools")
#Read GISTIC results along with MAF
laml.plus.gistic = read.maf(
maf = laml.maf,
gisticAllLesionsFile = all.lesions,
gisticAmpGenesFile = amp.genes,
gisticDelGenesFile = del.genes,
gisticScoresFile = scores.gis,
isTCGA = TRUE,
verbose = FALSE,
clinicalData = laml.clin
)
This plot shows frequent deletions in TP53 gene which is located on one of the significantly deleted locus 17p13.2.
In case there is no GISTIC results available, one can generate a table containing CN status for known genes in known samples. This can be easily created and read along with MAF file.
For example lets create a dummy CN alterations for DNMT3A
in random 20 samples.
set.seed(seed = 1024)
barcodes = as.character(getSampleSummary(x = laml)[,Tumor_Sample_Barcode])
#Random 20 samples
dummy.samples = sample(x = barcodes,
size = 20,
replace = FALSE)
#Genarate random CN status for above samples
cn.status = sample(
x = c('Amp', 'Del'),
size = length(dummy.samples),
replace = TRUE
)
custom.cn.data = data.frame(
Gene = "DNMT3A",
Sample_name = dummy.samples,
CN = cn.status,
stringsAsFactors = FALSE
)
head(custom.cn.data)
#> Gene Sample_name CN
#> 1 DNMT3A TCGA-AB-2898 Amp
#> 2 DNMT3A TCGA-AB-2879 Amp
#> 3 DNMT3A TCGA-AB-2920 Del
#> 4 DNMT3A TCGA-AB-2866 Amp
#> 5 DNMT3A TCGA-AB-2892 Amp
#> 6 DNMT3A TCGA-AB-2863 Amp
laml.plus.cn = read.maf(maf = laml.maf,
cnTable = custom.cn.data,
verbose = FALSE)
oncoplot(maf = laml.plus.cn, top = 5)
This data should be a data.frame or a tsv file with two required columns titled gene
and q
.
For example, including mutsig q
values into oncoplot.
Similar to significance values included as right bar plot, it’s also possible to include expression (or any sort of continuous) data to left side of the plot
#Dummy expression values for top 20 genes
set.seed(seed = 1024)
exprs_tbl = data.frame(genes = getGeneSummary(x = laml)[1:20, Hugo_Symbol],
exprn = rnorm(n = 10, mean = 12, sd = 5))
head(exprs_tbl)
#> genes exprn
#> 1 FLT3 8.106686
#> 2 DNMT3A 10.052618
#> 3 NPM1 1.831008
#> 4 IDH2 7.088134
#> 5 IDH1 13.239450
#> 6 TET2 1.480677
oncoplot(maf = laml, exprsTbl = exprs_tbl)
Annotations are usually stored in clinical.data
slot of MAF.
getClinicalData(x = laml)
#> Tumor_Sample_Barcode FAB_classification days_to_last_followup
#> 1: TCGA-AB-2802 M4 365
#> 2: TCGA-AB-2803 M3 792
#> 3: TCGA-AB-2804 M3 2557
#> 4: TCGA-AB-2805 M0 577
#> 5: TCGA-AB-2806 M1 945
#> ---
#> 189: TCGA-AB-3007 M3 1581
#> 190: TCGA-AB-3008 M1 822
#> 191: TCGA-AB-3009 M4 577
#> 192: TCGA-AB-3011 M1 1885
#> 193: TCGA-AB-3012 M3 1887
#> Overall_Survival_Status
#> 1: 1
#> 2: 1
#> 3: 0
#> 4: 1
#> 5: 1
#> ---
#> 189: 0
#> 190: 1
#> 191: 1
#> 192: 0
#> 193: 0
Include FAB_classification
from clinical data as one of the sample annotations.
More than one annotations can be included by passing them to the argument clinicalFeatures
. Above plot can be further enhanced by sorting according to annotations. Custom colors can be specified as a list of named vectors for each levels.
#Color coding for FAB classification
fabcolors = RColorBrewer::brewer.pal(n = 8,name = 'Spectral')
names(fabcolors) = c("M0", "M1", "M2", "M3", "M4", "M5", "M6", "M7")
fabcolors = list(FAB_classification = fabcolors)
print(fabcolors)
#> $FAB_classification
#> M0 M1 M2 M3 M4 M5 M6 M7
#> "#D53E4F" "#F46D43" "#FDAE61" "#FEE08B" "#E6F598" "#ABDDA4" "#66C2A5" "#3288BD"
oncoplot(
maf = laml,
clinicalFeatures = 'FAB_classification',
sortByAnnotation = TRUE,
annotationColor = fabcolors
)
If you prefer to highlight mutations by a specific attribute, you can use additionalFeature
argument.
Example: Highlight all mutations where alt allele is C.
Note that first argument (Tumor_Seq_Allele2) must a be column in MAF file, and second argument (C) is a value in that column. If you want to know what columns are present in the MAF file, use getFields
.
getFields(x = laml)
#> [1] "Hugo_Symbol" "Entrez_Gene_Id" "Center"
#> [4] "NCBI_Build" "Chromosome" "Start_Position"
#> [7] "End_Position" "Strand" "Variant_Classification"
#> [10] "Variant_Type" "Reference_Allele" "Tumor_Seq_Allele1"
#> [13] "Tumor_Seq_Allele2" "Tumor_Sample_Barcode" "Protein_Change"
#> [16] "i_TumorVAF_WU" "i_transcript_name"
sessionInfo()
#> R version 3.6.1 (2019-07-05)
#> Platform: x86_64-pc-linux-gnu (64-bit)
#> Running under: Ubuntu 18.04.3 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.10-bioc/R/lib/libRblas.so
#> LAPACK: /home/biocbuild/bbs-3.10-bioc/R/lib/libRlapack.so
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> attached base packages:
#> [1] stats4 parallel stats graphics grDevices utils datasets
#> [8] methods base
#>
#> other attached packages:
#> [1] pheatmap_1.0.12 doParallel_1.0.15
#> [3] iterators_1.0.12 foreach_1.4.7
#> [5] NMF_0.21.0 bigmemory_4.5.33
#> [7] Biobase_2.46.0 cluster_2.1.0
#> [9] rngtools_1.4 pkgmaker_0.27
#> [11] registry_0.5-1 BSgenome.Hsapiens.UCSC.hg19_1.4.0
#> [13] BSgenome_1.54.0 rtracklayer_1.46.0
#> [15] Biostrings_2.54.0 XVector_0.26.0
#> [17] GenomicRanges_1.38.0 GenomeInfoDb_1.22.0
#> [19] IRanges_2.20.1 S4Vectors_0.24.1
#> [21] BiocGenerics_0.32.0 maftools_2.2.10
#>
#> loaded via a namespace (and not attached):
#> [1] splines_3.6.1 R.utils_2.9.2
#> [3] assertthat_0.2.1 GenomeInfoDbData_1.2.2
#> [5] Rsamtools_2.2.1 yaml_2.2.0
#> [7] pillar_1.4.2 lattice_0.20-38
#> [9] glue_1.3.1 digest_0.6.23
#> [11] RColorBrewer_1.1-2 colorspace_1.4-1
#> [13] R.oo_1.23.0 plyr_1.8.5
#> [15] htmltools_0.4.0 Matrix_1.2-18
#> [17] XML_3.98-1.20 pkgconfig_2.0.3
#> [19] bibtex_0.4.2 zlibbioc_1.32.0
#> [21] purrr_0.3.3 xtable_1.8-4
#> [23] scales_1.1.0 BiocParallel_1.20.0
#> [25] tibble_2.1.3 ggplot2_3.2.1
#> [27] withr_2.1.2 SummarizedExperiment_1.16.0
#> [29] lazyeval_0.2.2 survival_3.1-8
#> [31] magrittr_1.5 crayon_1.3.4
#> [33] mclust_5.4.5 evaluate_0.14
#> [35] R.methodsS3_1.7.1 tools_3.6.1
#> [37] data.table_1.12.8 lifecycle_0.1.0
#> [39] matrixStats_0.55.0 gridBase_0.4-7
#> [41] stringr_1.4.0 munsell_0.5.0
#> [43] DelayedArray_0.12.0 compiler_3.6.1
#> [45] rlang_0.4.2 grid_3.6.1
#> [47] RCurl_1.95-4.12 bigmemory.sri_0.1.3
#> [49] bitops_1.0-6 rmarkdown_2.0
#> [51] gtable_0.3.0 codetools_0.2-16
#> [53] reshape2_1.4.3 R6_2.4.1
#> [55] GenomicAlignments_1.22.1 knitr_1.26
#> [57] dplyr_0.8.3 stringi_1.4.3
#> [59] Rcpp_1.0.3 wordcloud_2.6
#> [61] tidyselect_0.2.5 xfun_0.11