Introduction to SigsPack

Loading the package

library(SigsPack)

Loading a VCF

In most cases you will want to analyse real data, e.g. fit a sample’s mutational profile to known mutational signatures. You can easily load your data and use it with the package.

if (require(BSgenome.Hsapiens.UCSC.hg19)) {
  sample <- vcf2mut_cat(
    system.file("extdata", "example.vcf.gz", package="SigsPack"),
    BSgenome.Hsapiens.UCSC.hg19
  )
}
#> Loading required package: BSgenome.Hsapiens.UCSC.hg19
#> Loading required package: BSgenome
#> Loading required package: BiocGenerics
#> Loading required package: parallel
#> 
#> 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 object is masked from 'package:SigsPack':
#> 
#>     normalize
#> The following objects are masked from 'package:stats':
#> 
#>     IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#> 
#>     Filter, Find, Map, Position, Reduce, anyDuplicated, append,
#>     as.data.frame, basename, cbind, colnames, dirname, do.call,
#>     duplicated, eval, evalq, get, grep, grepl, intersect,
#>     is.unsorted, lapply, mapply, match, mget, order, paste, pmax,
#>     pmax.int, pmin, pmin.int, rank, rbind, rownames, sapply,
#>     setdiff, sort, table, tapply, union, unique, unsplit, which,
#>     which.max, which.min
#> Loading required package: S4Vectors
#> Loading required package: stats4
#> 
#> Attaching package: 'S4Vectors'
#> The following object is masked from 'package:base':
#> 
#>     expand.grid
#> Loading required package: IRanges
#> Loading required package: GenomeInfoDb
#> Loading required package: GenomicRanges
#> Loading required package: Biostrings
#> Loading required package: XVector
#> 
#> Attaching package: 'Biostrings'
#> The following object is masked from 'package:base':
#> 
#>     strsplit
#> Loading required package: rtracklayer

It will be transformed into a ‘mutational profile’ sorting all single nucleotide variants according to their trinucleotide context and mutation. The result is a 96 by 1 matrix that can be used as input to the analysis functions of the package as well as in most other packages of this field.

Simulating data

Instead of using real-life data it is also possible to generate synthetic data that can be used to analyse signatures stability. The following code generates 10 mutational catalogues with exposure to the COSMIC signatures 2, 6, 15 and 27. The catalogues consist of mutation counts sampled from a distribution that is a linear combination of the aforementioned signatures. The weight of each of these signatures can optionally be specified, too. For convenience, the COSMIC reference signature matrix is included in the package and used as default in the functions. However, it is also possible to use all functions with a custom signature matrix.

data("cosmicSigs")

cats <- create_mut_catalogues(10, 500, P=cosmicSigs, sig_set = c(2,6,15,27))
knitr::kable(head(cats[['catalogues']]))

Estimating signature exposures and bootstrapping samples

The function bootstrap_mut_catalogues bootstraps a sample returning a specified amount of re-samples which can be used to gain a better variability estimation of the sample’s signature exposure.

The signature exposure is calculated using quadratic programming. This can be done on one or several samples at once. The COSMIC signatures have been included in the package, and are used by default. However, it is possible to use your own signature matrix, or use a sub-set of COSMIC signatures, instead of the whole matrix.

reps <- bootstrap_mut_catalogues(n = 1000, original = cats[["catalogues"]][,1])

# using only signatures 4, 17, 23 and 30 for signature estimation
sigs <- signature_exposure(reps, sig_set = c(4,17,23,30))

print(sigs$exposures[,1])
#>  Signature 1  Signature 2  Signature 3  Signature 4  Signature 5 
#>  0.000000000  0.000000000  0.000000000  0.429266785  0.000000000 
#>  Signature 6  Signature 7  Signature 8  Signature 9 Signature 10 
#>  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000 
#> Signature 11 Signature 12 Signature 13 Signature 14 Signature 15 
#>  0.000000000  0.000000000  0.000000000  0.000000000  0.000000000 
#> Signature 16 Signature 17 Signature 18 Signature 19 Signature 20 
#>  0.000000000  0.061009737  0.000000000  0.000000000  0.000000000 
#> Signature 21 Signature 22 Signature 23 Signature 24 Signature 25 
#>  0.000000000  0.000000000  0.006687235  0.000000000  0.000000000 
#> Signature 26 Signature 27 Signature 28 Signature 29 Signature 30 
#>  0.000000000  0.000000000  0.000000000  0.000000000  0.503036242

With one command you can bootstrap a mutational catalogue and obtain a table and a plot illustrating the results of the signature estimation for this sample and the bootstrapped re-samples. The plot shows the distribution of estimated signature exposure for all the re-samples, highlighting the one of the original mutational catalogue, thus providing insights on the reliability of the estimates.

report <- summarize_exposures(reps[,1])

knitr::kable(
  head(report)
  )

Tri-nucleotide contexts and normalization

Accurate exposures estimation requires matching tri-nucleotide frequencies between observations and the signature matrix. The COSMIC signature matrix provided in the package is relative to the whole genome (GRCh37) tri-nucleotide frequencies. So if you want to fit those signatures to exome data, the data need to be normalized to match the signatures prior to exposures estimation.

Let’s say, that the vcf we derived the mutational catalogue from contained exome data. In the case of this example, we want to scale our exome data to the genome ‘space’ to match the COSMIC reference signatures. Hence, we need both the tri-nucleotide distribution of the human genome (GRCh37) and the exome that our data were collected on.

Extract the tri-nucleotide context frequencies from a genome (BSgenome object) or exome and use them to normalize the data.

if (require(BSgenome.Hsapiens.UCSC.hg19)){
  genome_contexts <- get_context_freq(BSgenome.Hsapiens.UCSC.hg19)
  exome_contexts <- get_context_freq(
    BSgenome.Hsapiens.UCSC.hg19,
    system.file("extdata", "example.bed.gz", package="SigsPack")
  )
  normalized_mut_cat <- SigsPack::normalize(sample, exome_contexts, hg19context_freq)
}

The normalization returns frequencies, not count numbers. If you want to use them for exposure estimation or for bootstrapping, the catalogue size must be scaled, or input together with the normalized catalogue.

if (require(BSgenome.Hsapiens.UCSC.hg19)) {
  sigs_norm <- signature_exposure(sum(sample) * normalized_mut_cat)
  report_norm <- summarize_exposures(normalized_mut_cat, m=sum(sample))
  reps_norm <- bootstrap_mut_catalogues(
      n=1000,
      original=normalized_mut_cat, 
      m=sum(sample)
  )
}

sessionInfo

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] BSgenome.Hsapiens.UCSC.hg19_1.4.0 BSgenome_1.54.0                  
#>  [3] rtracklayer_1.46.0                Biostrings_2.54.0                
#>  [5] XVector_0.26.0                    GenomicRanges_1.38.0             
#>  [7] GenomeInfoDb_1.22.0               IRanges_2.20.0                   
#>  [9] S4Vectors_0.24.0                  BiocGenerics_0.32.0              
#> [11] SigsPack_1.0.0                   
#> 
#> loaded via a namespace (and not attached):
#>  [1] Rcpp_1.0.2                  lattice_0.20-38            
#>  [3] prettyunits_1.0.2           Rsamtools_2.2.0            
#>  [5] assertthat_0.2.1            zeallot_0.1.0              
#>  [7] digest_0.6.22               BiocFileCache_1.10.0       
#>  [9] R6_2.4.0                    backports_1.1.5            
#> [11] RSQLite_2.1.2               evaluate_0.14              
#> [13] highr_0.8                   httr_1.4.1                 
#> [15] pillar_1.4.2                zlibbioc_1.32.0            
#> [17] rlang_0.4.1                 GenomicFeatures_1.38.0     
#> [19] progress_1.2.2              curl_4.2                   
#> [21] blob_1.2.0                  Matrix_1.2-17              
#> [23] rmarkdown_1.16              BiocParallel_1.20.0        
#> [25] stringr_1.4.0               RCurl_1.95-4.12            
#> [27] bit_1.1-14                  biomaRt_2.42.0             
#> [29] DelayedArray_0.12.0         compiler_3.6.1             
#> [31] xfun_0.10                   pkgconfig_2.0.3            
#> [33] askpass_1.1                 htmltools_0.4.0            
#> [35] tidyselect_0.2.5            openssl_1.4.1              
#> [37] SummarizedExperiment_1.16.0 tibble_2.1.3               
#> [39] GenomeInfoDbData_1.2.2      quadprog_1.5-7             
#> [41] matrixStats_0.55.0          XML_3.98-1.20              
#> [43] crayon_1.3.4                dplyr_0.8.3                
#> [45] dbplyr_1.4.2                rappdirs_0.3.1             
#> [47] GenomicAlignments_1.22.0    bitops_1.0-6               
#> [49] grid_3.6.1                  DBI_1.0.0                  
#> [51] magrittr_1.5                stringi_1.4.3              
#> [53] vctrs_0.2.0                 tools_3.6.1                
#> [55] bit64_0.9-7                 glue_1.3.1                 
#> [57] Biobase_2.46.0              purrr_0.3.3                
#> [59] hms_0.5.1                   yaml_2.2.0                 
#> [61] AnnotationDbi_1.48.0        memoise_1.1.0              
#> [63] knitr_1.25                  VariantAnnotation_1.32.0