Computation of phylogenetic trees and clustering of mutations

Before you begin

library(mitoClone2)

You should have identified true somatic variants in the mitochondrial (and nuclear) genome. The remaining vignettes of this package document how to get there. Here, we start with count matrices of the alternative and the reference alleles, across a number of sites of interest. Such data is available from two patients (P1, P2) as part of this package. Only P1 is analyzed as part of this vignette but the commented code for P2 should be fully functional.

load(system.file("data/M_P1.RData", package = "mitoClone2"))
load(system.file("data/N_P1.RData", package = "mitoClone2"))
P1 <- mutationCallsFromMatrix(as.matrix(M_P1), as.matrix(N_P1))

A first important step is to decide which mutations to include in the clustering. The default is to use all mutations that are covered in at least 20% of the cells, but this assignment can be changed manually.

Compute a phylogenetic tree

The next step is to run SCITE or PhISCS to compute the most likely phylogenetic tree. PhISCS is bundled in this package, but the package needs to be run in an environment where gurobi and the gurobipy python package are available. For example, you could set up a conda environment that contains this package. SCITE does not require any additional licenses but may need to be manually compiled.

## this next step takes approx 4.1 minutes to run
tmpd <- tempdir()
dir.create(paste0(tmpd, "/p1"))
P1 <- varCluster(P1, tempfolder = paste0(tmpd, "/p1"), method = "SCITE")
#> Warning in max(which(x == 1)): no non-missing arguments to max; returning -Inf

This step can take a while to run. It computes a likely phylogenetic tree of all the mutations. In the case of SCITE, multiple equally likely tree can be produced. In it’s current state, this package simply selects the first tree from the list.

Identify clones and assign cells to clones

In many cases, the order of the leaves on these trees is arbitrary, because mutations systematically co-occur. We therefore cluster the mutations into clones. In detail, we take every every branch on the tree and then shuffle the order of mutations in that branch while re-calculating the likelihood. If swapping nodes leads to small changes in the likelihood, these nodes are then merged into a “clone”. The parameter min.lik that controls the merging is set arbitrarily, see below for more information.

P1 <- clusterMetaclones(P1, min.lik = 1)
#> Warning in mutcalls@mut2clone[branches[[i]]] <-
#> as.integer(max(mutcalls@mut2clone) + : number of items to replace is not a
#> multiple of replacement length

This step also assigns each cell to the most likely clone, and provides an estimate of the likelihood. Refer to help(mutationCalls) for more info on how these results are stored.

Finally, the clustering can be plotted.

plotClones(P1)

Parameter choice

The parameter min.lik that controls the merging is set arbitrarily. In practice, the goal of these analyses is to group mutations into clones for subsequent analyses (such as differential expression analyses) and it may make sense to overwrite the result of clusterMetaclones manually; for example, if a subclone defned on a mitochondrial mutation only should be treated as part of a more clearly defined upstream clone for differential expression analysis.

To overwrite the result of clusterMetaclones, first retrieve the assignment of mutations to clones:

m2c <- mitoClone2:::getMut2Clone(P1)
print(m2c)
#>   X1282GA   X2537GA   X3335TC   X3350TC   X5492TC X7527GDEL   X8167TC  X11196GA 
#>         2         3         5         5         1         4         5         2 
#>  X11559GA  X14462GA  X16233AG      EAPP     SRSF2     CEBPA      KLF7      root 
#>         4         3         5         4         4         3         2         5

## To e.g. treat the mt:2537G>A and mt:14462:G>A mutations as a
## subclone distinct from CEBPA, we can assign them a new clonal
## identity
m2c[c("X2537GA", "X14462GA")] <- as.integer(6)

P1.new <- mitoClone2:::overwriteMetaclones(P1, m2c)
plotClones(P1.new)

Session information

#> R version 4.1.0 (2021-05-18)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows Server x64 (build 17763)
#> 
#> 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] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] mitoClone2_0.99.20
#> 
#> loaded via a namespace (and not attached):
#>   [1] bitops_1.0-7                matrixStats_0.60.0         
#>   [3] bit64_4.0.5                 RColorBrewer_1.1-2         
#>   [5] filelock_1.0.2              progress_1.2.2             
#>   [7] httr_1.4.2                  GenomeInfoDb_1.29.3        
#>   [9] tools_4.1.0                 bslib_0.2.5.1              
#>  [11] utf8_1.2.2                  R6_2.5.0                   
#>  [13] DBI_1.1.1                   BiocGenerics_0.39.1        
#>  [15] colorspace_2.0-2            tidyselect_1.1.1           
#>  [17] prettyunits_1.1.1           bit_4.0.4                  
#>  [19] curl_4.3.2                  compiler_4.1.0             
#>  [21] Biobase_2.53.0              formatR_1.11               
#>  [23] xml2_1.3.2                  DelayedArray_0.19.1        
#>  [25] rtracklayer_1.53.0          sass_0.4.0                 
#>  [27] scales_1.1.1                deepSNV_1.39.0             
#>  [29] rappdirs_0.3.3              stringr_1.4.0              
#>  [31] digest_0.6.27               Rsamtools_2.9.1            
#>  [33] rmarkdown_2.9               XVector_0.33.0             
#>  [35] pkgconfig_2.0.3             htmltools_0.5.1.1          
#>  [37] MatrixGenerics_1.5.1        highr_0.9                  
#>  [39] dbplyr_2.1.1                fastmap_1.1.0              
#>  [41] BSgenome_1.61.0             rlang_0.4.11               
#>  [43] RSQLite_2.2.7               VGAM_1.1-5                 
#>  [45] farver_2.1.0                jquerylib_0.1.4            
#>  [47] BiocIO_1.3.0                generics_0.1.0             
#>  [49] jsonlite_1.7.2              BiocParallel_1.27.2        
#>  [51] dplyr_1.0.7                 VariantAnnotation_1.39.0   
#>  [53] RCurl_1.98-1.3              magrittr_2.0.1             
#>  [55] GenomeInfoDbData_1.2.6      Matrix_1.3-4               
#>  [57] Rcpp_1.0.7                  munsell_0.5.0              
#>  [59] S4Vectors_0.31.0            fansi_0.5.0                
#>  [61] lifecycle_1.0.0             stringi_1.7.3              
#>  [63] yaml_2.2.1                  SummarizedExperiment_1.23.1
#>  [65] zlibbioc_1.39.0             BiocFileCache_2.1.1        
#>  [67] grid_4.1.0                  blob_1.2.2                 
#>  [69] parallel_4.1.0              crayon_1.4.1               
#>  [71] lattice_0.20-44             Biostrings_2.61.1          
#>  [73] splines_4.1.0               GenomicFeatures_1.45.0     
#>  [75] hms_1.1.0                   KEGGREST_1.33.0            
#>  [77] knitr_1.33                  pillar_1.6.1               
#>  [79] GenomicRanges_1.45.0        rjson_0.2.20               
#>  [81] biomaRt_2.49.2              stats4_4.1.0               
#>  [83] XML_3.99-0.6                glue_1.4.2                 
#>  [85] evaluate_0.14               png_0.1-7                  
#>  [87] vctrs_0.3.8                 gtable_0.3.0               
#>  [89] purrr_0.3.4                 assertthat_0.2.1           
#>  [91] cachem_1.0.5                ggplot2_3.3.5              
#>  [93] xfun_0.24                   restfulr_0.0.13            
#>  [95] pheatmap_1.0.12             tibble_3.1.3               
#>  [97] GenomicAlignments_1.29.0    AnnotationDbi_1.55.1       
#>  [99] memoise_2.0.0               IRanges_2.27.0             
#> [101] Rhtslib_1.25.1              ellipsis_0.3.2