SignIT: Mutation Signature Analysis in Individual Tumours with MCMC

This vignette demonstrates the use of SignIT to decipher mutation signatures in individual tumours, explore signature bleed interactions, and perform temporal dissection of mutation signatures. SignIT is available for download and installation from https://www.github.com/eyzhao/SignIT. See installation instructions in the README.

SignIT on a Patient’s Mutation Catalog

data(example_catalog)
pander(head(example_catalog))

data(example_variants)
pander(head(example_variants, 10))

example_catalog_2 <- with(
  example_variants,
  mutations_to_catalog(chr, pos, ref, alt)
)

## No BSgenome object provided. Using default: BSgenome.Hsapiens.UCSC.hg19

## 
## Attaching package: 'BiocGenerics'

## The following objects are masked from 'package:parallel':
## 
##     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
##     clusterExport, clusterMap, parApply, parCapply, parLapply,
##     parLapplyLB, parRapply, parSapply, parSapplyLB

## The following objects are masked from 'package:dplyr':
## 
##     combine, intersect, setdiff, union

## The following objects are masked from 'package:stats':
## 
##     IQR, mad, sd, var, xtabs

## The following objects are masked from 'package:base':
## 
##     anyDuplicated, append, as.data.frame, cbind, colMeans,
##     colnames, colSums, do.call, duplicated, eval, evalq, Filter,
##     Find, get, grep, grepl, intersect, is.unsorted, lapply,
##     lengths, Map, mapply, match, mget, order, paste, pmax,
##     pmax.int, pmin, pmin.int, Position, rank, rbind, Reduce,
##     rowMeans, rownames, rowSums, sapply, setdiff, sort, table,
##     tapply, union, unique, unsplit, which, which.max, which.min

## 
## Attaching package: 'S4Vectors'

## The following objects are masked from 'package:dplyr':
## 
##     first, rename

## The following object is masked from 'package:tidyr':
## 
##     expand

## The following object is masked from 'package:base':
## 
##     expand.grid

## 
## Attaching package: 'IRanges'

## The following objects are masked from 'package:dplyr':
## 
##     collapse, desc, slice

## The following object is masked from 'package:purrr':
## 
##     reduce

## 
## Attaching package: 'XVector'

## The following object is masked from 'package:purrr':
## 
##     compact

## 
## Attaching package: 'Biostrings'

## The following object is masked from 'package:base':
## 
##     strsplit

## Chromosome levels match after prepending "chr". Proceeding...

## Warning in get_trinucleotide(chr, pos, ref, genome = genome): "chr" was
## automatically prepended to all chromosome names to match genome seqlevels

pander(head(example_catalog_2))

patient_signature_exposures <- get_exposures(
  mutation_catalog = example_catalog,
  n_cores = 1
)

## Sampling

## 
## SAMPLING FOR MODEL 'signature_model' NOW (CHAIN 1).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 400 [  0%]  (Warmup)
## Iteration:  40 / 400 [ 10%]  (Warmup)
## Iteration:  80 / 400 [ 20%]  (Warmup)
## Iteration: 120 / 400 [ 30%]  (Warmup)
## Iteration: 160 / 400 [ 40%]  (Warmup)
## Iteration: 200 / 400 [ 50%]  (Warmup)
## Iteration: 201 / 400 [ 50%]  (Sampling)
## Iteration: 240 / 400 [ 60%]  (Sampling)
## Iteration: 280 / 400 [ 70%]  (Sampling)
## Iteration: 320 / 400 [ 80%]  (Sampling)
## Iteration: 360 / 400 [ 90%]  (Sampling)
## Iteration: 400 / 400 [100%]  (Sampling)
## 
##  Elapsed Time: 9.77 seconds (Warm-up)
##                4.31 seconds (Sampling)
##                14.08 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'signature_model' NOW (CHAIN 2).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 400 [  0%]  (Warmup)
## Iteration:  40 / 400 [ 10%]  (Warmup)
## Iteration:  80 / 400 [ 20%]  (Warmup)
## Iteration: 120 / 400 [ 30%]  (Warmup)
## Iteration: 160 / 400 [ 40%]  (Warmup)
## Iteration: 200 / 400 [ 50%]  (Warmup)
## Iteration: 201 / 400 [ 50%]  (Sampling)
## Iteration: 240 / 400 [ 60%]  (Sampling)
## Iteration: 280 / 400 [ 70%]  (Sampling)
## Iteration: 320 / 400 [ 80%]  (Sampling)
## Iteration: 360 / 400 [ 90%]  (Sampling)
## Iteration: 400 / 400 [100%]  (Sampling)
## 
##  Elapsed Time: 10.67 seconds (Warm-up)
##                4.71 seconds (Sampling)
##                15.38 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'signature_model' NOW (CHAIN 3).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 400 [  0%]  (Warmup)
## Iteration:  40 / 400 [ 10%]  (Warmup)
## Iteration:  80 / 400 [ 20%]  (Warmup)
## Iteration: 120 / 400 [ 30%]  (Warmup)
## Iteration: 160 / 400 [ 40%]  (Warmup)
## Iteration: 200 / 400 [ 50%]  (Warmup)
## Iteration: 201 / 400 [ 50%]  (Sampling)
## Iteration: 240 / 400 [ 60%]  (Sampling)
## Iteration: 280 / 400 [ 70%]  (Sampling)
## Iteration: 320 / 400 [ 80%]  (Sampling)
## Iteration: 360 / 400 [ 90%]  (Sampling)
## Iteration: 400 / 400 [100%]  (Sampling)
## 
##  Elapsed Time: 10.13 seconds (Warm-up)
##                3.44 seconds (Sampling)
##                13.57 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'signature_model' NOW (CHAIN 4).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 400 [  0%]  (Warmup)
## Iteration:  40 / 400 [ 10%]  (Warmup)
## Iteration:  80 / 400 [ 20%]  (Warmup)
## Iteration: 120 / 400 [ 30%]  (Warmup)
## Iteration: 160 / 400 [ 40%]  (Warmup)
## Iteration: 200 / 400 [ 50%]  (Warmup)
## Iteration: 201 / 400 [ 50%]  (Sampling)
## Iteration: 240 / 400 [ 60%]  (Sampling)
## Iteration: 280 / 400 [ 70%]  (Sampling)
## Iteration: 320 / 400 [ 80%]  (Sampling)
## Iteration: 360 / 400 [ 90%]  (Sampling)
## Iteration: 400 / 400 [100%]  (Sampling)
## 
##  Elapsed Time: 9.54 seconds (Warm-up)
##                3.63 seconds (Sampling)
##                13.17 seconds (Total)

patient_exposure_posterior_plot <- plot_exposure_posteriors(patient_signature_exposures)
print(patient_exposure_posterior_plot)

patient_exposure_nnls_plot <- plot_nnls_solution(patient_signature_exposures)
print(patient_exposure_nnls_plot)

plot_two_signature_hexplot(patient_signature_exposures, 'Signature 3', 'Signature 8')

plot_signature_pairwise_bleed(patient_signature_exposures)

plot_exposure_posteriors_bleed(patient_signature_exposures, min_bleed = 0.25)

summary_table <- get_exposure_summary_table(patient_signature_exposures)
pander(summary_table)

Structural Variant Signatures

Mutation types are not fixed. As long as you provide a reference mutation signature matrix which matches the mutation types of your mutation catalog, you can decipher mutation signatures. We will briefly demonstrate an example using structural variant signatures.

The structural variant reference signatures from Nik-Zainal et al. (2016) can be accessed as follows

sv_reference <- get_reference_signatures(signature_set = 'nikzainal_sv')
pander(head(sv_reference))

We then use a mutation catalog with matching mutation types

data(example_sv_catalog)
pander(example_sv_catalog)

Next, we perform SignIT analysis as before.

sv_exposures <- get_exposures(example_sv_catalog, reference_signatures = sv_reference, n_cores = 1)

## Sampling

## 
## SAMPLING FOR MODEL 'signature_model' NOW (CHAIN 1).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 400 [  0%]  (Warmup)
## Iteration:  40 / 400 [ 10%]  (Warmup)
## Iteration:  80 / 400 [ 20%]  (Warmup)
## Iteration: 120 / 400 [ 30%]  (Warmup)
## Iteration: 160 / 400 [ 40%]  (Warmup)
## Iteration: 200 / 400 [ 50%]  (Warmup)
## Iteration: 201 / 400 [ 50%]  (Sampling)
## Iteration: 240 / 400 [ 60%]  (Sampling)
## Iteration: 280 / 400 [ 70%]  (Sampling)
## Iteration: 320 / 400 [ 80%]  (Sampling)
## Iteration: 360 / 400 [ 90%]  (Sampling)
## Iteration: 400 / 400 [100%]  (Sampling)
## 
##  Elapsed Time: 0.44 seconds (Warm-up)
##                0.33 seconds (Sampling)
##                0.77 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'signature_model' NOW (CHAIN 2).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 400 [  0%]  (Warmup)
## Iteration:  40 / 400 [ 10%]  (Warmup)
## Iteration:  80 / 400 [ 20%]  (Warmup)
## Iteration: 120 / 400 [ 30%]  (Warmup)
## Iteration: 160 / 400 [ 40%]  (Warmup)
## Iteration: 200 / 400 [ 50%]  (Warmup)
## Iteration: 201 / 400 [ 50%]  (Sampling)
## Iteration: 240 / 400 [ 60%]  (Sampling)
## Iteration: 280 / 400 [ 70%]  (Sampling)
## Iteration: 320 / 400 [ 80%]  (Sampling)
## Iteration: 360 / 400 [ 90%]  (Sampling)
## Iteration: 400 / 400 [100%]  (Sampling)
## 
##  Elapsed Time: 0.54 seconds (Warm-up)
##                0.58 seconds (Sampling)
##                1.12 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'signature_model' NOW (CHAIN 3).
## 
## Gradient evaluation took 0 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 400 [  0%]  (Warmup)
## Iteration:  40 / 400 [ 10%]  (Warmup)
## Iteration:  80 / 400 [ 20%]  (Warmup)
## Iteration: 120 / 400 [ 30%]  (Warmup)
## Iteration: 160 / 400 [ 40%]  (Warmup)
## Iteration: 200 / 400 [ 50%]  (Warmup)
## Iteration: 201 / 400 [ 50%]  (Sampling)
## Iteration: 240 / 400 [ 60%]  (Sampling)
## Iteration: 280 / 400 [ 70%]  (Sampling)
## Iteration: 320 / 400 [ 80%]  (Sampling)
## Iteration: 360 / 400 [ 90%]  (Sampling)
## Iteration: 400 / 400 [100%]  (Sampling)
## 
##  Elapsed Time: 0.41 seconds (Warm-up)
##                0.33 seconds (Sampling)
##                0.74 seconds (Total)
## 
## 
## SAMPLING FOR MODEL 'signature_model' NOW (CHAIN 4).
## 
## Gradient evaluation took 0.01 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 100 seconds.
## Adjust your expectations accordingly!
## 
## 
## Iteration:   1 / 400 [  0%]  (Warmup)
## Iteration:  40 / 400 [ 10%]  (Warmup)
## Iteration:  80 / 400 [ 20%]  (Warmup)
## Iteration: 120 / 400 [ 30%]  (Warmup)
## Iteration: 160 / 400 [ 40%]  (Warmup)
## Iteration: 200 / 400 [ 50%]  (Warmup)
## Iteration: 201 / 400 [ 50%]  (Sampling)
## Iteration: 240 / 400 [ 60%]  (Sampling)
## Iteration: 280 / 400 [ 70%]  (Sampling)
## Iteration: 320 / 400 [ 80%]  (Sampling)
## Iteration: 360 / 400 [ 90%]  (Sampling)
## Iteration: 400 / 400 [100%]  (Sampling)
## 
##  Elapsed Time: 0.42 seconds (Warm-up)
##                0.32 seconds (Sampling)
##                0.74 seconds (Total)

And we can use the same plotting functions,

plot_exposure_posteriors_bleed(sv_exposures, signature_trim = 'Rearrangement Signature')

And the same summarization function

pander(get_exposure_summary_table(sv_exposures))

Temporally Dissecting Mutation Signatures

We next explore how to perform temporal dissection of mutation signatures. The input data here is a data.frame or tibble of mutations, with some additional metadata. The following glimpse shows the necessary formatting of the input. Feel free to load data(example_variants) to see the complete input data for yourself.

glimpse(example_variants)

## Observations: 9,506
## Variables: 15
## $ chr            <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ pos            <dbl> 137744, 1479032, 2038218, 2615086, 2615087, 262...
## $ ref            <chr> "G", "T", "C", "C", "T", "C", "G", "A", "C", "A...
## $ alt            <chr> "C", "A", "T", "T", "C", "A", "A", "G", "T", "G...
## $ A              <int> 0, 8, 0, 0, 0, 24, 6, 87, 0, 81, 0, 4, 0, 25, 0...
## $ C              <int> 9, 3, 57, 41, 5, 191, 0, 0, 84, 0, 38, 15, 48, ...
## $ G              <int> 77, 0, 1, 0, 0, 0, 33, 7, 0, 10, 0, 0, 34, 0, 2...
## $ T              <int> 0, 104, 7, 5, 41, 0, 0, 0, 7, 0, 30, 0, 0, 27, ...
## $ total_depth    <int> 86, 115, 65, 46, 46, 215, 39, 94, 91, 91, 68, 1...
## $ alt_depth      <int> 9, 8, 7, 5, 5, 24, 6, 7, 7, 10, 38, 4, 48, 25, ...
## $ cnv_start      <int> 10022, 582430, 582430, 582430, 582430, 582430, ...
## $ cnv_end        <int> 582429, 8975823, 8975823, 8975823, 8975823, 897...
## $ tumour_copy    <int> 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,...
## $ tumour_content <dbl> 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63, 0.63,...
## $ normal_copy    <int> 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2,...

Signature timing is performed using get_population_signatures. Note that this process can take quite some time depending on your machine. We are using method = ‘vb’, which employs the variational Bayes solution and is much more computationally economical.

Note that if the parameter n_populations is not provided, an automatic process of model selection is performed to attempt to estimate the optimal number of signatures.

A preferred approach (if computational resources are available) is to perform complete population signature inference for five models: 1-5 populations. Then, compute the WAIC of each model with compute_population_signatures_waic(population_signatures) and choose the model which minimizes it. WAIC is a preferred model selection metric for Bayesian inference.

population_signatures <- get_population_signatures(
  example_variants,
  n_populations = 2,
  method = 'vb'
)

## Chromosome levels match after prepending "chr". Proceeding...

## Warning in get_trinucleotide(chr, pos, ref, genome = genome): "chr" was
## automatically prepended to all chromosome names to match genome seqlevels

## No BSgenome object provided. Using default: BSgenome.Hsapiens.UCSC.hg19
## Chromosome levels match after prepending "chr". Proceeding...

## Warning in get_trinucleotide(chr, pos, ref, genome = genome): "chr" was
## automatically prepended to all chromosome names to match genome seqlevels

## Extracting parameters from population model with 2 populations

## Subsetting reference signatures on request

## Sampling

## Running model with provided reference signatures: 13 signatures, 9506 mutations, and 2 clusters

## Sampling

## ------------------------------------------------------------
## EXPERIMENTAL ALGORITHM:
##   This procedure has not been thoroughly tested and may be unstable
##   or buggy. The interface is subject to change.
## ------------------------------------------------------------
## 
## 
## 
## Gradient evaluation took 2.35 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 23500 seconds.
## Adjust your expectations accordingly!
## 
## 
## Begin eta adaptation.
## Iteration:   1 / 250 [  0%]  (Adaptation)
## Iteration:  50 / 250 [ 20%]  (Adaptation)
## Iteration: 100 / 250 [ 40%]  (Adaptation)
## Iteration: 150 / 250 [ 60%]  (Adaptation)
## Iteration: 200 / 250 [ 80%]  (Adaptation)
## Success! Found best value [eta = 1] earlier than expected.
## 
## Begin stochastic gradient ascent.
##   iter       ELBO   delta_ELBO_mean   delta_ELBO_med   notes 
##    100     -8e+04             1.000            1.000
##    200     -8e+04             0.500            1.000
##    300     -8e+04             0.334            0.001   MEDIAN ELBO CONVERGED
## 
## Drawing a sample of size 1000 from the approximate posterior... 
## COMPLETED.

Now, let’s visualize the results.

plot_population_signatures(population_signatures)

## Warning in function_list[[k]](value): NAs introduced by coercion

## Warning in function_list[[k]](value): NAs introduced by coercion

This analysis is revealing two tumour subpopulations. The red higher-prevalence subpopulation is inferred to be “earlier”. The Proportion refers to the fraction of mutations contributed from each subpopulation.

We see here that Signature 3 exhibits early bias while signatures 16, 17, and 30 have late bias.

For quantitative analysis workflows, use the following to provide a summary of these results:

pander(summarise_population_signatures(population_signatures))

## Warning in function_list[[k]](value): NAs introduced by coercion

## Warning in function_list[[k]](value): NAs introduced by coercion

Session Info

The below output specifies the packages used in this analysis, and their version numbers.

print(sessionInfo())

## R version 3.4.1 (2017-06-30)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS release 6.7 (Final)
## 
## Matrix products: default
## BLAS: /projects/ezhao_prj/dependencies/miniconda3/envs/dependencies/lib/R/lib/libRblas.so
## LAPACK: /projects/ezhao_prj/dependencies/miniconda3/envs/dependencies/lib/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
##  [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] BSgenome.Hsapiens.UCSC.hg19_1.4.0 BSgenome_1.44.2                  
##  [3] rtracklayer_1.36.6                Biostrings_2.46.0                
##  [5] XVector_0.18.0                    GenomicRanges_1.30.0             
##  [7] GenomeInfoDb_1.14.0               IRanges_2.12.0                   
##  [9] S4Vectors_0.16.0                  BiocGenerics_0.24.0              
## [11] bindrcpp_0.2                      pander_0.6.1                     
## [13] ggraph_1.0.1                      hexbin_1.27.1                    
## [15] cowplot_0.9.2                     rstan_2.17.2                     
## [17] StanHeaders_2.17.2                dplyr_0.7.4                      
## [19] purrr_0.2.4                       readr_1.1.1                      
## [21] tidyr_0.8.0                       tibble_1.3.4                     
## [23] ggplot2_2.2.1                     tidyverse_1.1.1                  
## [25] nnls_1.4                          signit_1.0.1                     
## [27] Rcpp_0.12.13                     
## 
## loaded via a namespace (and not attached):
##  [1] nlme_3.1-131               bitops_1.0-6              
##  [3] matrixStats_0.52.2         lubridate_1.6.0           
##  [5] RColorBrewer_1.1-2         doParallel_1.0.11         
##  [7] httr_1.3.1                 rprojroot_1.3-1           
##  [9] tools_3.4.1                backports_1.1.2           
## [11] R6_2.2.2                   lazyeval_0.2.1            
## [13] colorspace_1.3-2           tidyselect_0.2.3          
## [15] gridExtra_2.3              mnormt_1.5-5              
## [17] compiler_3.4.1             rvest_0.3.2               
## [19] Biobase_2.38.0             xml2_1.1.1                
## [21] DelayedArray_0.4.1         labeling_0.3              
## [23] scales_0.5.0               psych_1.7.3.21            
## [25] stringr_1.2.0              digest_0.6.13             
## [27] Rsamtools_1.30.0           foreign_0.8-67            
## [29] rmarkdown_1.8.5            pkgconfig_2.0.1           
## [31] htmltools_0.3.6            readxl_1.0.0              
## [33] rlang_0.1.4                bindr_0.1                 
## [35] jsonlite_1.5               BiocParallel_1.12.0       
## [37] inline_0.3.14              RCurl_1.95-4.8            
## [39] magrittr_1.5               GenomeInfoDbData_0.99.1   
## [41] Matrix_1.2-7.1             munsell_0.4.3             
## [43] viridis_0.5.0              stringi_1.1.6             
## [45] yaml_2.1.16                MASS_7.3-47               
## [47] SummarizedExperiment_1.8.0 zlibbioc_1.24.0           
## [49] plyr_1.8.4                 grid_3.4.1                
## [51] ggrepel_0.7.0              forcats_0.2.0             
## [53] udunits2_0.13              lattice_0.20-34           
## [55] haven_1.0.0                hms_0.4.1                 
## [57] knitr_1.17                 igraph_1.1.2              
## [59] reshape2_1.4.3             codetools_0.2-15          
## [61] rmutil_1.1.0               XML_3.98-1.6              
## [63] glue_1.2.0                 evaluate_0.10.1           
## [65] modelr_0.1.0               tweenr_0.1.5              
## [67] foreach_1.4.4              cellranger_1.1.0          
## [69] gtable_0.2.0               assertthat_0.2.0          
## [71] ggforce_0.1.1              broom_0.4.2               
## [73] viridisLite_0.3.0          iterators_1.0.9           
## [75] GenomicAlignments_1.12.2   units_0.4-4