## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
    collapse = TRUE,
    comment = "#>",
    fig.width = 7, # Width of plots in inches
    fig.height = 5 # Height of plots in inches
)

## ----install, eval=F----------------------------------------------------------
# BiocManager::install("scafari")

## ----loading , eval=T, warning=FALSE------------------------------------------
library(scafari)

## ----reading, message=FALSE---------------------------------------------------
library(SingleCellExperiment)

# Load .h5 file
h5_file_path <- system.file("extdata", "demo.h5", package = "scafari")
h5 <- h5ToSce(h5_file_path)
sce <- h5$sce_amp
se.var <- h5$se_var

## ----metadata-----------------------------------------------------------------
# Get metadata
metadata(sce)

## ----loglopplot, message=FALSE------------------------------------------------
logLogPlot(sce)

## ----normalize_rc-------------------------------------------------------------
sce <- normalizeReadCounts(sce = sce)

## ----known_canon, eval=TRUE---------------------------------------------------
library(R.utils)

temp_dir <- tempdir()
known_canon_path <- file.path(
    temp_dir,
    "UCSC_hg19_knownCanonical_goldenPath.txt"
)

if (!file.exists(known_canon_path)) {
    url <- paste0(
        "http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/",
        "knownCanonical.txt.gz"
    )
    destfile <- file.path(
        temp_dir,
        "UCSC_hg19_knownCanonical_goldenPath.txt.gz"
    )

    # Download the file
    download.file(url, destfile)

    # Decompress the file
    R.utils::gunzip(destfile, remove = FALSE)
}

## ----annotate_amp-------------------------------------------------------------
try(annotateAmplicons(sce = sce, known.canon = known_canon_path))

## ----plot_amp_dist, warning=FALSE---------------------------------------------
plotAmpliconDistribution(sce)

## ----plot_norm_rc-------------------------------------------------------------
plotNormalizedReadCounts(sce = sce)

## ----plot_panel_uniformity, warning=FALSE-------------------------------------
plotPanelUniformity(sce, interactive = FALSE)

## ----filter_variants, message = FALSE-----------------------------------------
library(SummarizedExperiment)

filteres <- filterVariants(
    depth.threshold = 10,
    genotype.quality.threshold = 30,
    vaf.ref = 5,
    vaf.het = 35,
    vaf.hom = 95,
    min.cell = 50,
    min.mut.cell = 1,
    se.var = se.var,
    sce = sce,
    shiny = FALSE
)
se.f <-
    SummarizedExperiment(
        assays =
            list(
                VAF = t(filteres$vaf.matrix.filtered),
                Genotype = t(filteres$genotype.matrix.filtered),
                Genoqual = t(filteres$genoqual.matrix.filtered)
            ),
        rowData = filteres$variant.ids.filtered,
        colData = filteres$cells.keep
    )

# Filter out cells in sce object
# Find the indices of the columns to keep
indices_to_keep <- match(filteres$cells.keep,
    SummarizedExperiment::colData(sce)[[1]],
    nomatch = 0
)

# Subset the SCE using these indices
sce_filtered <- sce[, indices_to_keep]
SingleCellExperiment::altExp(sce_filtered, "variants") <- se.f

## ----annotate_variants--------------------------------------------------------
suppressPackageStartupMessages(library(DT))
suppressPackageStartupMessages(library(dplyr))

# Load sce object with filtered variants
# Check metadata
metadata(altExp(sce_filtered))

# Annotate variants
sce_filtered <- annotateVariants(sce = sce_filtered)

# Check metadata. Now a annotation flag is present
metadata(altExp(sce_filtered))

# The annotations are stored in the rowData
rowData(altExp(sce_filtered)) %>%
    as.data.frame() %>%
    head() %>%
    datatable()

## ----plot_variant_heatmap, fig.height=8---------------------------------------
plotVariantHeatmap(sce = sce_filtered)

## ----plot_gt_quality----------------------------------------------------------
plotGenotypequalityPerGenotype(sce = sce_filtered)

## ----variants_of_interest-----------------------------------------------------
variants.of.interest <- c(
    "KIT:chr4:55599436:T/C",
    "TET2:chr4:106158216:G/A",
    "FLT3:chr13:28610183:A/G",
    "TP53:chr17:7577427:G/A"
)

## ----cluster_variants, fig.show='hide'----------------------------------------
plot <- clusterVariantSelection(
    sce = sce_filtered,
    variants.of.interest = variants.of.interest,
    n.clust = 3
)

names(plot)

## ----show_clusters------------------------------------------------------------
plot$clusterplot

## ----plot_clusters_vaf, warning=FALSE-----------------------------------------
plotClusterVAF(
    sce = sce_filtered,
    variants.of.interest = variants.of.interest,
    gg.clust = plot$clusterplot
)

## ----plot_clusters_vaf_map, warning=FALSE-------------------------------------
plotClusterVAFMap(
    sce = sce_filtered,
    variants.of.interest = variants.of.interest,
    gg.clust = plot$clusterplot
)

## ----plot_clusters_gt, warning=FALSE------------------------------------------
plotClusterGenotype(
    sce = sce_filtered,
    variants.of.interest = variants.of.interest,
    gg.clust = plot$clusterplot
)

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