inst/doc/MesKit.R

In MesKit: A tool kit for dissecting cancer evolution from multi-region derived tumor biopsies via somatic alterations

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(echo = TRUE)

## ----echo=FALSE, paged.print=FALSE, rownames.print=FALSE----------------------
MAFtable <- read.table(system.file("extdata", "CRC_HZ.maf", package = "MesKit"), header=TRUE)
extractLines <- rbind(MAFtable[1, ], MAFtable[6600, ])
extractLines <- rbind(extractLines, MAFtable[15000, ])
data.frame(extractLines, row.names = NULL)

## ----echo=FALSE, paged.print=FALSE, rownames.print=FALSE----------------------
ClinInfo <- read.table(system.file("extdata", "CRC_HZ.clin.txt", package = "MesKit"), header = TRUE)
ClinInfo[1:5, ]

## ----echo=FALSE, paged.print=FALSE, rownames.print=FALSE----------------------
ccfInfo <- read.table(system.file("extdata", "CRC_HZ.ccf.tsv", package = "MesKit"), header = TRUE)
ccfInfo[1:5, ]

## ----echo=FALSE, paged.print=FALSE, rownames.print=FALSE----------------------
segInfo <- read.table(system.file("extdata", "CRC_HZ.seg.txt", package = "MesKit"), header = TRUE)
segInfo[1:5, ]

## ---- eval=FALSE--------------------------------------------------------------
#  # Installation of MesKit requires Bioconductor version 3.12 or higher
#  if (!requireNamespace("BiocManager", quietly = TRUE)){
#    install.packages("BiocManager")
#  }
#  # The following initializes usage of Bioc 3.12
#  BiocManager::install(version = "3.12")
#  BiocManager::install("MesKit")

## ----eval=FALSE---------------------------------------------------------------
#  if (!requireNamespace("devtools", quietly = TRUE)) {
#    install.packages("devtools")
#  }
#  devtools::install_github("Niinleslie/MesKit")

## ----message=FALSE,warning=FALSE----------------------------------------------
library(MesKit)

## -----------------------------------------------------------------------------
maf.File <- system.file("extdata/", "CRC_HZ.maf", package = "MesKit")
ccf.File <- system.file("extdata/", "CRC_HZ.ccf.tsv", package = "MesKit")
clin.File <- system.file("extdata", "CRC_HZ.clin.txt", package = "MesKit")
# Maf object with CCF information
maf <- readMaf(mafFile = maf.File,
               ccfFile = ccf.File,
               clinicalFile  = clin.File,
               refBuild = "hg19")  

## ----message=FALSE, fig.align='left', fig.width=11, fig.height=6.5------------
# Driver genes of CRC collected from [IntOGen] (https://www.intogen.org/search) (v.2020.2)
driverGene.file <- system.file("extdata/", "IntOGen-DriverGenes_COREAD.tsv", package = "MesKit")
driverGene <- as.character(read.table(driverGene.file)$V1)
mut.class <- classifyMut(maf, class =  "SP", patient.id = 'V402')
head(mut.class)

## ----message=FALSE, fig.align='left', fig.width=11, fig.height=6.5------------
plotMutProfile(maf, class =  "SP", geneList = driverGene, use.tumorSampleLabel = TRUE)

## ----message=FALSE------------------------------------------------------------
# Read segment file
segCN <- system.file("extdata", "CRC_HZ.seg.txt", package = "MesKit")
# Read gistic output files
all.lesions <- system.file("extdata", "COREAD_all_lesions.conf_99.txt", package = "MesKit")
amp.genes <- system.file("extdata", "COREAD_amp_genes.conf_99.txt", package = "MesKit")
del.genes <- system.file("extdata", "COREAD_del_genes.conf_99.txt", package = "MesKit")
seg <- readSegment(segFile = segCN, gisticAllLesionsFile = all.lesions,
                   gisticAmpGenesFile = amp.genes, gisticDelGenesFile = del.genes)
seg$V402[1:5, ]

## ---- fig.width=10, fig.align='left', fig.height=5----------------------------
plotCNA(seg, patient.id = c("V402", "V750", "V824"), use.tumorSampleLabel = TRUE)

## ----warning=FALSE, fig.width=8, fig.height=5---------------------------------
# calculate MATH score of each sample
mathScore(maf, patient.id = 'V402')

## ----fig.height=5, fig.width=6, message=FALSE---------------------------------
AUC <- ccfAUC(maf, patient.id = 'V402', use.tumorSampleLabel = TRUE)
names(AUC)

## ----fig.height=4, fig.width=4.5, message=FALSE, fig.align='left'-------------
# show cumulative density plot of CCF
AUC$CCF.density.plot

## ----message=FALSE, fig.height=4.5, fig.width=12------------------------------
mut.cluster <- mutCluster(maf, patient.id = 'V402', 
                          use.ccf = TRUE, use.tumorSampleLabel = TRUE)
clusterPlots <- mut.cluster$cluster.plot
cowplot::plot_grid(plotlist = clusterPlots[1:6])

## ----fig.align='left', fig.width=5, fig.height=4.5, fig.asp=1.0---------------
# calculate the Fst of brain metastasis from V402
calFst(maf, patient.id = 'V402', plot = TRUE, use.tumorSampleLabel = TRUE, 
       withinTumor  = TRUE, number.cex = 10)[["V402_BM"]]

## ----fig.align='left', fig.width=5, fig.height=4.5, fig.asp=1.0---------------
# calculate the Nei's genetic distance of brain metastasis from V402
calNeiDist(maf, patient.id = 'V402', use.tumorSampleLabel = TRUE, 
           withinTumor  = TRUE, number.cex = 10)[["V402_BM"]]

## ----fig.align='left', fig.width=4, fig.height=4.5, message=FALSE, warning=FALSE----
ccf.list <- compareCCF(maf, pairByTumor = TRUE, min.ccf = 0.02,
                       use.adjVAF = TRUE, use.indel = FALSE)
V402_P_BM <- ccf.list$V402$`P-BM`
# visualize via smoothScatter R package
graphics::smoothScatter(matrix(c(V402_P_BM[, 3], V402_P_BM[, 4]),ncol = 2),
                xlim = c(0, 1), ylim = c(0, 1),
                colramp = colorRampPalette(c("white", RColorBrewer::brewer.pal(9, "BuPu"))),
                xlab = "P", ylab = "BM")
  
## show driver genes
gene.idx <- which(V402_P_BM$Hugo_Symbol %in% driverGene) 
points(V402_P_BM[gene.idx, 3:4], cex = 0.6, col = 2, pch = 2)
text(V402_P_BM[gene.idx, 3:4], cex = 0.7, pos = 1,
       V402_P_BM$Hugo_Symbol[gene.idx])
title("V402 JSI = 0.341", cex.main = 1.5)

## ----fig.align='left', fig.width=7, fig.height=7, fig.asp=1.0-----------------
JSI.res <- calJSI(maf, patient.id = 'V402', pairByTumor = TRUE, min.ccf = 0.02, 
                  use.adjVAF = TRUE, use.indel = FALSE, use.tumorSampleLabel = TRUE)
names(JSI.res)

## ----fig.align='left', fig.width=7, fig.height=7, fig.asp=1.0-----------------
# show the JSI result
JSI.res$JSI.multi
JSI.res$JSI.pair

## ----message=FALSE, warning=FALSE, fig.height=4, fig.width=4------------------
neutralResult <- testNeutral(maf, min.mut.count = 10, patient.id = 'V402', use.tumorSampleLabel = TRUE)
neutralResult$neutrality.metrics
neutralResult$model.fitting.plot$P_1

## -----------------------------------------------------------------------------
phyloTree <- getPhyloTree(maf, patient.id = "V402", method = "NJ", min.vaf = 0.06)

## ---- message=FALSE, warning=FALSE, fig.height=4, fig.width=10----------------
tree.NJ <- getPhyloTree(maf, patient.id = 'V402', method = "NJ")
tree.MP <- getPhyloTree(maf, patient.id = 'V402', method = "MP")
# compare phylogenetic trees constructed by two approaches
compareTree(tree.NJ, tree.MP, plot = TRUE, use.tumorSampleLabel = TRUE)

## ---- message=FALSE, warning=FALSE, fig.width=4.5, fig.height=4.5-------------
library(org.Hs.eg.db)
library(clusterProfiler)
# Pathway enrichment analysis
V402.branches <- getMutBranches(phyloTree)
# pathway enrichment for private mutated genes of the primary tumor in patient V402
V402_Public <- V402.branches[V402.branches$Mutation_Type == "Private_P", ]
geneIDs = suppressMessages(bitr(V402_Public$Hugo_Symbol, fromType="SYMBOL", 
                              toType=c("ENTREZID"), OrgDb="org.Hs.eg.db"))
KEGG_V402_Private_P  = enrichKEGG(
  gene = geneIDs$ENTREZID,
  organism = 'hsa',
  keyType = 'kegg',
  pvalueCutoff = 0.05,
)
dotplot(KEGG_V402_Private_P)

## ----message=FALSE, warning = FALSE-------------------------------------------
# load the genome reference
library(BSgenome.Hsapiens.UCSC.hg19)

## ---- warning = FALSE, message=FALSE------------------------------------------
mutClass <- mutTrunkBranch(phyloTree, CT = TRUE, plot = TRUE)
names(mutClass)

## ---- warning = FALSE, message=FALSE, fig.height=4.5, fig.width=4.5-----------
mutClass$mutTrunkBranch.res
mutClass$mutTrunkBranch.plot

## ---- warning = FALSE, message=FALSE, fig.height=2.5, fig.width=8-------------
trimatrix_V402 <- triMatrix(phyloTree, level = 5)
# Visualize the 96 trinucleodide mutational profile
plotMutSigProfile(trimatrix_V402)[[1]]


## ---- fig.height=5.2, fig.width=8---------------------------------------------
# reconstruct mutational profile of V402 using COSMIC V2 signatures
fit_V402 <- fitSignatures(trimatrix_V402, signaturesRef = "cosmic_v2")
# Compare the reconstructed mutational profile with the original mutational profile
plotMutSigProfile(fit_V402)[[1]]

## ---- warning = FALSE, message=FALSE, fig.height=4,5, fig.width=7.5-----------
# Below plot shows cosine similarities between the mutational profile of each group and COSMIC signatures
library(ComplexHeatmap)
ComplexHeatmap::Heatmap(fit_V402$V402$cosine.similarity, name = "Cosine similarity")

## ---- fig.align='left', fig.width=12, fig.height=6, message=FALSE, warning=FALSE----
# A phylogenetic tree along with binary and CCF heatmap of mutations 
phylotree_V402 <- plotPhyloTree(phyloTree, use.tumorSampleLabel = TRUE)
binary_heatmap_V402 <- mutHeatmap(maf, min.ccf = 0.04, use.ccf = FALSE, patient.id = "V402", use.tumorSampleLabel = TRUE)
CCF_heatmap_V402 <- mutHeatmap(maf, use.ccf = TRUE, patient.id = "V402", 
                                 min.ccf = 0.04, use.tumorSampleLabel = TRUE)
cowplot::plot_grid(phylotree_V402, binary_heatmap_V402,
                   CCF_heatmap_V402, nrow = 1, rel_widths = c(1.5, 1, 1))

## -----------------------------------------------------------------------------
sessionInfo()