# Customizing oncoplots

#### Anand Mayakonda

#### 2022-02-24

    library(maftools)

    #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)

## 0.1 Including Transition/Transversions into oncoplot

    #By default the function plots top20 mutated genes
    oncoplot(maf = laml, draw_titv = TRUE)

## 0.2 Changing colors for variant classifications

    #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)

## 0.3 Including copy number data into oncoplots.

There are two ways one include CN status into MAF. 1. GISTIC results 2. Custom copy number table

### 0.3.1 GISTIC results

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
    )

    oncoplot(maf = laml.plus.gistic, top = 10)

This plot shows frequent deletions in TP53 gene which is located on one of the significantly deleted locus 17p13.2.

### 0.3.2 Custom copy-number table

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('ShallowAmp', 'DeepDel', 'Del', 'Amp'),
      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 ShallowAmp
    #> 2 DNMT3A TCGA-AB-2879        Del
    #> 3 DNMT3A TCGA-AB-2920        Amp
    #> 4 DNMT3A TCGA-AB-2866        Del
    #> 5 DNMT3A TCGA-AB-2892        Del
    #> 6 DNMT3A TCGA-AB-2863 ShallowAmp

    laml.plus.cn = read.maf(maf = laml.maf,
                            cnTable = custom.cn.data,
                            verbose = FALSE)

    oncoplot(maf = laml.plus.cn, top = 5)

## 0.4 Bar plots

`leftBarData`, `rightBarData` and `topBarData` arguments can be used to display additional values as barplots. Below example demonstrates adding gene expression values and mutsig q-values as left and right side bars respectivelly.

    #Selected AML driver genes
    aml_genes = c("TP53", "WT1", "PHF6", "DNMT3A", "DNMT3B", "TET1", "TET2", "IDH1", "IDH2", "FLT3", "KIT", "KRAS", "NRAS", "RUNX1", "CEBPA", "ASXL1", "EZH2", "KDM6A")

    #Variant allele frequcnies (Right bar plot)
    aml_genes_vaf = subsetMaf(maf = laml, genes = aml_genes, fields = "i_TumorVAF_WU", mafObj = FALSE)[,mean(i_TumorVAF_WU, na.rm = TRUE), Hugo_Symbol]
    colnames(aml_genes_vaf)[2] = "VAF"
    head(aml_genes_vaf)
    #>    Hugo_Symbol      VAF
    #> 1:       ASXL1 37.11250
    #> 2:       CEBPA 22.00235
    #> 3:      DNMT3A 43.51556
    #> 4:      DNMT3B 37.14000
    #> 5:        EZH2 68.88500
    #> 6:        FLT3 34.60294

    #MutSig results (Right bar plot)
    laml.mutsig = system.file("extdata", "LAML_sig_genes.txt.gz", package = "maftools")
    laml.mutsig = data.table::fread(input = laml.mutsig)[,.(gene, q)]
    laml.mutsig[,q := -log10(q)] #transoform to log10
    head(laml.mutsig)
    #>      gene        q
    #> 1:   FLT3 12.64176
    #> 2: DNMT3A 12.64176
    #> 3:   NPM1 12.64176
    #> 4:   IDH2 12.64176
    #> 5:   IDH1 12.64176
    #> 6:   TET2 12.64176

    oncoplot(
      maf = laml,
      genes = aml_genes,
      leftBarData = aml_genes_vaf,
      leftBarLims = c(0, 100),
      rightBarData = laml.mutsig,
      rightBarLims = c(0, 20)
    )

## 0.5 Including annotations

Annotations are 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.

    oncoplot(maf = laml, genes = aml_genes, clinicalFeatures = 'FAB_classification')

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, genes = aml_genes,
      clinicalFeatures = 'FAB_classification',
      sortByAnnotation = TRUE,
      annotationColor = fabcolors
    )

## 0.6 Highlighting samples

If you prefer to highlight mutations by a specific attribute, you can use `additionalFeature` argument.

Example: Highlight all mutations where alt allele is C.

    oncoplot(maf = laml, genes = aml_genes,
             additionalFeature = c("Tumor_Seq_Allele2", "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"

## 0.7 Group by Pathways

Genes can be auto grouped based on their Biological processess by setting `pathways = 'auto'` or by providing custom pathway list in the form of a two column tsv file or a data.frame containing gene names and their corresponding pathway.

### 0.7.1 Auto

setting `pathways = 'auto'` draws top 3 most affected pathways

    oncoplot(maf = laml, pathways = "auto", gene_mar = 8, fontSize = 0.6)
    #> Drawing upto top 3 mutated pathways
    #>                                Pathway  N n_affected_genes fraction_affected
    #>  1:                      RTK_signaling 16                6        0.37500000
    #>  2:          Epigenetics_DNA_modifiers  1                1        1.00000000
    #>  3:                    Chromatin_other 14                7        0.50000000
    #>  4:                         Metabolism  2                2        1.00000000
    #>  5:               Transcription_factor 39               11        0.28205128
    #>  6:                   Genome_integrity 14                5        0.35714286
    #>  7:                     MAPK_signaling  9                3        0.33333333
    #>  8:                    Other_signaling 28                5        0.17857143
    #>  9:                           Splicing  6                3        0.50000000
    #> 10:               Histone_modification  3                2        0.66666667
    #> 11:                              Other 22                3        0.13636364
    #> 12: Protein_homeostasis/ubiquitination 15               NA                NA
    #> 13:                   Immune_signaling 10                3        0.30000000
    #> 14:          Chromatin_SWI/SNF_complex  8               NA                NA
    #> 15:                    NOTCH_signaling  1                1        1.00000000
    #> 16:                     PI3K_signaling  9                1        0.11111111
    #> 17:                      RNA_abundance 15                1        0.06666667
    #> 18:            Wnt/B-catenin_signaling  8               NA                NA
    #>     Mutated_samples Fraction_mutated_samples
    #>  1:              63              0.326424870
    #>  2:              48              0.248704663
    #>  3:              45              0.233160622
    #>  4:              37              0.191709845
    #>  5:              35              0.181347150
    #>  6:              28              0.145077720
    #>  7:              25              0.129533679
    #>  8:              13              0.067357513
    #>  9:              11              0.056994819
    #> 10:               4              0.020725389
    #> 11:               4              0.020725389
    #> 12:               4              0.020725389
    #> 13:               3              0.015544041
    #> 14:               2              0.010362694
    #> 15:               1              0.005181347
    #> 16:               1              0.005181347
    #> 17:               1              0.005181347
    #> 18:               0              0.000000000

### 0.7.2 Custom pathways

    pathways = data.frame(
      Genes = c(
        "TP53",
        "WT1",
        "PHF6",
        "DNMT3A",
        "DNMT3B",
        "TET1",
        "TET2",
        "IDH1",
        "IDH2",
        "FLT3",
        "KIT",
        "KRAS",
        "NRAS",
        "RUNX1",
        "CEBPA",
        "ASXL1",
        "EZH2",
        "KDM6A"
      ),
      Pathway = rep(c(
        "TSG", "DNAm", "Signalling", "TFs", "ChromMod"
      ), c(3, 6, 4, 2, 3)),
      stringsAsFactors = FALSE
    )

    head(pathways)
    #>    Genes Pathway
    #> 1   TP53     TSG
    #> 2    WT1     TSG
    #> 3   PHF6     TSG
    #> 4 DNMT3A    DNAm
    #> 5 DNMT3B    DNAm
    #> 6   TET1    DNAm

    oncoplot(maf = laml, pathways = pathways, gene_mar = 8, fontSize = 0.6)

## 0.8 Combining everything

    oncoplot(
      maf = laml.plus.gistic,
      draw_titv = TRUE,
      pathways = pathways,
      clinicalFeatures = c('FAB_classification', 'Overall_Survival_Status'),
      sortByAnnotation = TRUE,
      additionalFeature = c("Tumor_Seq_Allele2", "C"),
      leftBarData = aml_genes_vaf,
      leftBarLims = c(0, 100),
      rightBarData = laml.mutsig[,.(gene, q)],
    )

## 0.9 SessionInfo

    sessionInfo()
    #> R version 4.1.2 (2021-11-01)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server 2012 R2 x64 (build 9600)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.1252
    #> [3] LC_MONETARY=English_United States.1252
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.1252
    #>
    #> attached base packages:
    #> [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    #> [8] methods   base
    #>
    #> other attached packages:
    #>  [1] pheatmap_1.0.12                   doParallel_1.0.17
    #>  [3] iterators_1.0.14                  foreach_1.5.2
    #>  [5] NMF_0.23.0                        Biobase_2.54.0
    #>  [7] cluster_2.1.2                     rngtools_1.5.2
    #>  [9] pkgmaker_0.32.2                   registry_0.5-1
    #> [11] BSgenome.Hsapiens.UCSC.hg19_1.4.3 BSgenome_1.62.0
    #> [13] rtracklayer_1.54.0                Biostrings_2.62.0
    #> [15] XVector_0.34.0                    GenomicRanges_1.46.1
    #> [17] GenomeInfoDb_1.30.1               IRanges_2.28.0
    #> [19] S4Vectors_0.32.3                  BiocGenerics_0.40.0
    #> [21] mclust_5.4.9                      maftools_2.10.05
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] bitops_1.0-7                matrixStats_0.61.0
    #>  [3] RColorBrewer_1.1-2          tools_4.1.2
    #>  [5] bslib_0.3.1                 utf8_1.2.2
    #>  [7] R6_2.5.1                    DBI_1.1.2
    #>  [9] colorspace_2.0-3            withr_2.4.3
    #> [11] tidyselect_1.1.2            curl_4.3.2
    #> [13] compiler_4.1.2              cli_3.2.0
    #> [15] DelayedArray_0.20.0         sass_0.4.0
    #> [17] scales_1.1.1                stringr_1.4.0
    #> [19] digest_0.6.29               Rsamtools_2.10.0
    #> [21] rmarkdown_2.11              R.utils_2.11.0
    #> [23] pkgconfig_2.0.3             htmltools_0.5.2
    #> [25] MatrixGenerics_1.6.0        fastmap_1.1.0
    #> [27] highr_0.9                   rlang_1.0.1
    #> [29] jquerylib_0.1.4             BiocIO_1.4.0
    #> [31] generics_0.1.2              jsonlite_1.8.0
    #> [33] BiocParallel_1.28.3         dplyr_1.0.8
    #> [35] berryFunctions_1.21.0       R.oo_1.24.0
    #> [37] RCurl_1.98-1.6              magrittr_2.0.2
    #> [39] GenomeInfoDbData_1.2.7      Matrix_1.4-0
    #> [41] Rcpp_1.0.8                  munsell_0.5.0
    #> [43] fansi_1.0.2                 abind_1.4-5
    #> [45] lifecycle_1.0.1             R.methodsS3_1.8.1
    #> [47] stringi_1.7.6               yaml_2.3.5
    #> [49] RaggedExperiment_1.18.0     SummarizedExperiment_1.24.0
    #> [51] zlibbioc_1.40.0             plyr_1.8.6
    #> [53] grid_4.1.2                  crayon_1.5.0
    #> [55] lattice_0.20-45             splines_4.1.2
    #> [57] knitr_1.37                  pillar_1.7.0
    #> [59] rjson_0.2.21                reshape2_1.4.4
    #> [61] codetools_0.2-18            XML_3.99-0.8
    #> [63] glue_1.6.2                  evaluate_0.15
    #> [65] data.table_1.14.2           MultiAssayExperiment_1.20.0
    #> [67] vctrs_0.3.8                 gtable_0.3.0
    #> [69] purrr_0.3.4                 assertthat_0.2.1
    #> [71] ggplot2_3.3.5               gridBase_0.4-7
    #> [73] xfun_0.29                   xtable_1.8-4
    #> [75] restfulr_0.0.13             survival_3.2-13
    #> [77] tibble_3.1.6                GenomicAlignments_1.30.0
    #> [79] ellipsis_0.3.2