In compbiomed/singleCellTK: Comprehensive and Interactive Analysis of Single Cell RNA-Seq Data

require(singleCellTK)
require(ggplot2)
require(dplyr)
require(Biobase)
require(cowplot)

docs.base <- paste0("https://www.camplab.net/sctk/v",
                    package.version("singleCellTK"), "/")
docs.qc.Path <- paste0(docs.base, "articles/articles/cmd_qc.html")

metaChecking <- function(metainfo, samples){
  if (any(!samples %in% names(metainfo))) {

    ### if length of metainfo is 1, it indicates that the name is the sample name
    ### This format is compatible with the QC report.
    ### In this case, all samples are ran with single run. Return any sublist is fine. 
    if (length(names(metainfo)) == 1) {
      names(metainfo) <- sample
      return(metainfo)
    }

    ### if QC is ran in R console instead of QC pipeline, we will add place holder.
    ### In this case, different samples also share parameters. Return any sublist is fine 
    newMeta <- list() 
    for (sample in samples) {
      newMeta[[sample]] <- metainfo
    }
    metainfo <- newMeta
  }

  return(metainfo)
}

preprocessMeta <- function(meta, ignore=c("geneSetCollection", "sample", "batch", "batchBackground", "bgBatch")){
  meta <- meta[!names(meta) %in% ignore]
  lapply(meta, function(y){
    while(is.list(y)) {
      if (length(y) != 0){
       y <- y[[1]]
      } else {
       y <- "NULL"
      }
    }

    if(is.vector(y)) {
      y <- paste(y,collapse=' ')
    }

    if(is.null(y)) {
      y <- "NULL"
    }

    return(y)
  })
}

formatPlotList <- function(plotlist, samples) {
  if (length(samples) > 1) {
    return(plotlist)
  }

  plotlist <- list(Sample = list(plotlist))
  names(plotlist$Sample) <- samples
  return(plotlist)
}

sce.qc <- params$object
subTitle <- params$subTitle
studyDesign <- params$studyDesign
reducedDimName <- params$reducedDimName
sceSample <- colData(sce.qc)$sample
samples <- unique(sceSample)

## modify assayType data structure for sce generated by SCTK-QC script pipelinne
if (typeof(S4Vectors::metadata(sce.qc)$assayType) == 'list') {
  assayType <- BiocGenerics::Reduce(dplyr::union, S4Vectors::metadata(sce.qc)$assayType)
  S4Vectors::metadata(sce.qc)$assayType <- assayType
}

sce.qc <- runQuickUMAP(sce.qc, reducedDimName = "QC_UMAP", sample = sceSample)

if (is.null(reducedDimName)) {
  allSampleReducedDim <- "All_UMAP"
  sce.qc <- runQuickUMAP(sce.qc, reducedDimName = allSampleReducedDim, sample = NULL)
} else if (!reducedDimName %in% reducedDimNames(sce.qc)) {
  warning("'reducedDimName' not found in the reducedDimNames of the sce object. ",
  "Generate new reduced dimension reduction for all samples. ")
  allSampleReducedDim <- "All_UMAP"
  sce.qc <- runQuickUMAP(sce.qc, reducedDimName = allSampleReducedDim, sample = NULL)
} else {
  allSampleReducedDim <- reducedDimName
}


### when 'sample' doesnn't exists in colData
if (is.null(samples)) { samples <- "sample" }

---

subtitle: "r subTitle"

r studyDesign

Introduction

cat("\n")

intro <- paste0(
  "Comprehensive quality control (QC) of single-cell RNA-seq data was performed with the [**singleCellTK**](https://github.com/compbiomed/singleCellTK) package. This report contains information about each QC tool and visualization of the QC metrics for each sample. For more information on running this pipeline and performing quality control, see the [**documentation**](",
  docs.qc.Path,
  "). If you use the singleCellTK package for quality control, please include a [**reference**](https://pubmed.ncbi.nlm.nih.gov/35354805/) in your publication."
)

cat(intro)

description_runPerCellQC <- singleCellTK:::descriptionRunPerCellQC()

Summary Statistics {.tabset .tabset-fade}

cat("\n")
cat("# {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
cat("\n")

cat(paste0('## SCTK-QC {.tabset .tabset-fade} \n\n'))
sce.qc <- sampleSummaryStats(inSCE = sce.qc, 
                                   sample = sceSample,
                                   simple = FALSE)

summaryTable <- getSampleSummaryStatsTable(sce.qc, statsName = "qc_table")
cat(apply(summaryTable, 1:2, as.character) %>% 
    knitr::kable(format = "html") %>%
      kableExtra::kable_styling() %>%
      kableExtra::scroll_box(width = "80%"))

cat("\n")
summary <- "The summary statistics table summarizes QC metrics of the cell matrix. This table summarizes the mean and median of UMI counts and median of genes detected per cell, as well as the number and percentages of doublets and estimated ambient RNA scores per dataset."
cat(summary)
cat("\n\n") 

# cellranger
if ("cellranger" %in% listSampleSummaryStatsTables(sce.qc)) {
  cat(paste0('## cellranger {.tabset .tabset-fade} \n\n'))
  crSummaryTable <- getSampleSummaryStatsTable(sce.qc, 'cellranger')
  cat(apply(crSummaryTable, 1:2, as.character) %>% 
    knitr::kable(format = "html") %>%
      kableExtra::kable_styling() %>%
      kableExtra::scroll_box(width = "80%"))
  cat("\n")
  summary <- "The summary statistics table summarizes QC metrics parsed from 'metrics_summary.csv' file if that's available within the cellranger count output folder. This table summarizes the number of cells detected, mean reads / genes per cell, total number of reads, sequencing saturation and other QC metrics of reads alignment"
  cat(summary)
  cat("\n\n") 
}

# starsolo
if ("starsolo" %in% listSampleSummaryStatsTables(sce.qc)) {
  cat(paste0('## starsolo {.tabset .tabset-fade} \n\n'))
  strSummaryTable <- getSampleSummaryStatsTable(sce.qc, 'starsolo')
  cat(apply(strSummaryTable, 1:2, as.character) %>% 
    knitr::kable(format = "html") %>%
      kableExtra::kable_styling() %>%
      kableExtra::scroll_box(width = "80%"))
  cat("\n")
  summary <- "The summary statistics table summarizes QC metrics parsed from 'Summary.csv' file if that's available within the STARsolo output folder. This table summarizes the number of cells detected, mean reads / genes per cell, total number of reads, reads sequencinng quality, sequencing saturation and other QC metrics of reads alignment"
  cat(summary)
  cat("\n\n") 
}
cat("\n")

i="General quality control metrics"

perCellData <- c("sum", "detected", "percent.top_50")
skipCellData <- any(!perCellData %in% names(colData(sce.qc)))

if (skipCellData) {
  #cat("runPerCellQC did not complete successfully. The section of General quality control metrics is skipped")
  plotsQCMetrics <- NULL
} else {
  cat(paste0('# ', i, ' \n'))
  cat("\n")
  plotsQCMetrics<- tryCatch(
    {plotRunPerCellQCResults(sce.qc, 
                             sample = colData(sce.qc)$sample,
                             combinePlot="none",
                             axisSize = NULL,
                             axisLabelSize = NULL,
                             titleSize = NULL)},
    error = function(e) {
      cat("runPerCellQC did not complete successfully. The section of General quality control metrics is skipped")
      skipCellData <<- TRUE
      return(NA)
    }
  )
}

if (!skipCellData) {

  if (length(samples) == 1) {
    plotsQCMetrics <- list(Violin=plotsQCMetrics)
  }
  # names(plotsQCMetrics$Violin)[names(plotsQCMetrics$Violin) == "mito_percent"] <- "MitoGene_Top50_Percent"
  # names(plotsQCMetrics$Violin)[names(plotsQCMetrics$Violin) == "mito_sum"] <- "MitoGene_Sum"
  # names(plotsQCMetrics$Violin)[names(plotsQCMetrics$Violin) == "mito_detected"] <- "MitoGene_Features"
  geneSubsetData <- grep("subsets_.*", names(plotsQCMetrics$Violin), value = TRUE)
  geneSubsetName <- unique(gsub("_sum|_detected|_percent","",
                                gsub("subsets_","",geneSubsetData)))

  cat(description_runPerCellQC$introduction)
  cat(description_runPerCellQC$runPerCellQC)
  cat(description_runPerCellQC$plotRunPerCellQCResults)
  cat(description_runPerCellQC$output)
  cat(description_runPerCellQC$sum)
  cat(description_runPerCellQC$detected)
  cat(description_runPerCellQC$percentTop)
  cat(description_runPerCellQC$subsets)
}

{.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}

cat("\n")
cat(paste0('## Total Counts {.tabset .tabset-fade} \n\n'))
plotsQCMetrics$Violin$Sum
cat("\n\n")
cat(paste0('## Total Features {.tabset .tabset-fade} \n\n'))
plotsQCMetrics$Violin$Detected
cat("\n\n")
cat(paste0('## Percentage of Library Size Occupied by Top 50 Expressed Features {.tabset .tabset-fade} \n\n'))
plotsQCMetrics$Violin$TopPercent
if(!is.null(plotsQCMetrics$Violin$mito_sum)){
  cat("\n\n")
  cat(paste0('## Total Mitochondrial Counts {.tabset .tabset-fade} \n\n'))
  plotsQCMetrics$Violin$mito_sum
}
if(!is.null(plotsQCMetrics$Violin$mito_detected)){
  cat("\n\n")
  cat(paste0('## Total Mitochondrial Features {.tabset .tabset-fade} \n\n'))
  plotsQCMetrics$Violin$mito_detected
}
if(!is.null(plotsQCMetrics$Violin$mito_percent)){
  cat("\n\n")
  cat(paste0('## Percentage of Mitochondrial Counts {.tabset .tabset-fade} \n\n'))
  plotsQCMetrics$Violin$mito_percent
}


if(length(geneSubsetName) > 0){
  for(ix in seq_along(geneSubsetName)){
    cat("\n\n")
    cat(paste0('## Total Counts of Feature Subset "', geneSubsetName[ix], '"{.tabset .tabset-fade} \n\n'))
    print(plotsQCMetrics$Violin[[paste0("subsets_",geneSubsetName[ix],"_sum")]])
    cat("\n\n")
    cat(paste0('## Total Features within Feature Subset "', geneSubsetName[ix], '"{.tabset .tabset-fade} \n\n'))
    print(plotsQCMetrics$Violin[[paste0("subsets_",geneSubsetName[ix],"_detected")]])
    cat("\n\n")
    cat(paste0('## Percentage of Feature Subset "', geneSubsetName[ix], '" Expression{.tabset .tabset-fade} \n\n'))
    print(plotsQCMetrics$Violin[[paste0("subsets_",geneSubsetName[ix],"_percent")]])
  }
}

cat("\n\n")
cat(paste0('## Parameters {.tabset .tabset-fade} \n\n'))

runPerCellMeta <- sce.qc@metadata$scater_addPerCellQC
runPerCellMeta <- metaChecking(runPerCellMeta, samples[[1]])
x <- preprocessMeta(runPerCellMeta[[ samples[[1]] ]])
t(as.data.frame(x)) %>%
  knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
  kableExtra::scroll_box(width = "80%")

cat("\n")
cat(description_runPerCellQC$parameter)

r cat("\n\n")

Doublet Detection

fontSize <- 15
sampleUMAP =  plotSCEDimReduceColData(
    inSCE=sce.qc,
    sample=NULL,
    colorBy="sample",
    conditionClass="factor",
    reducedDimName=allSampleReducedDim,
    xlab=c("QC_UMAP_1"),
    ylab=c("QC_UMAP_2"),
    axisSize = fontSize,
    axisLabelSize = fontSize,
    legendSize = fontSize,
    titleSize = fontSize,
    labelClusters=FALSE,
    title="Sample",
    plotLabels = NULL,
    combinePlot="none"
  )

dbAlgo <- setNames(
  c("scrublet_call", "scDblFinder_doublet_call", "scds_cxds_call", "scds_bcds_call", "scds_hybrid_call"),
  c("Scrublet", "scDblFinder", "Scds_Cxds", "Scds_Bcds", "Scds_Hybrid")
  )

dF <- grep("doubletFinder_doublet_label", names(colData(sce.qc)), value=TRUE)
if (length(dF) != 0) {
  names(dF) <- paste(unlist(strsplit(dF, "_"))[-c(2,3,4)], collapse = "_")
  dbAlgo <- c(dbAlgo, dF)
}

if (length(dbAlgo) != 0) {
  i="Doublet Detection Summary"
  cat(paste0('## ', i, ' \n'))
  cat("## {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
  cat("\n")

  for (algo in dbAlgo) {
    if (!algo %in% names(colData(sce.qc))) {
      next
    }
    algoName <- names(dbAlgo)[dbAlgo == algo]
    cat(paste0("### ", algoName, " {.tabset .tabset-fade} \n"))
    cat("\n\n")

    plotDoublet <- plotSCEDimReduceColData(
        inSCE=sce.qc,
        sample=NULL,
        colorBy=algo,
        conditionClass="factor",
        colorScale = c("lightgray","red"),
        # shape=shape,
        reducedDimName=allSampleReducedDim,
        axisSize = fontSize,
        axisLabelSize = fontSize,
        legendSize = fontSize,
        titleSize = fontSize,
        xlab=c("QC_UMAP_1"),
        ylab=c("QC_UMAP_2"),
        labelClusters=FALSE,
        title=algoName,
        legendTitle = "Doublet Call",
        plotLabels = NULL,
        combinePlot="none"
      )

    print(plot_grid(sampleUMAP, plotDoublet, ncol = 2))
    cat("\n\n")
  }
}

description_Scrublet <- singleCellTK:::descriptionScrublet()

i="Scrublet"

scrubletData <- c("scrublet_score", "scrublet_call")
skipScrublet <- any(!scrubletData %in% names(colData(sce.qc)))

if (skipScrublet) {
  #cat("runScrublet did not complete successfully. The section of Scrublet is skipped")
  plotsScrublet <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
    plotsScrublet<- tryCatch(
      {plotScrubletResults(sce.qc, 
                           reducedDimName="QC_UMAP",
                           sample = colData(sce.qc)$sample,
                           combinePlot = "none",
                           axisSize = NULL,
                           axisLabelSize = NULL,
                           titleSize = NULL)},
    error = function(e) {
      cat("runScrublet did not complete successfully. The section of Scrublet is skipped")
      skipScrublet <<- TRUE
      return(NA)
    }
  )
}

if (!skipScrublet) {
  plotsScrublet <- formatPlotList(plotsScrublet, samples)
  cat(description_Scrublet$introduction)
  cat(description_Scrublet$runScrublet)
  cat(description_Scrublet$plotScrubletResults)
  cat(description_Scrublet$output)
  cat("\n")
}

cat("## {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
cat("\n")
for (sample in samples) {
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n"))
  cat("\n\n")

  cat("#### Scrublet Doublet Assignment {.tabset .tabset-fade} \n\n")
  print(plotsScrublet$Sample[[sample]][["scatter_doubletCall"]])
  cat("\n\n")  

  cat("#### Scrublet Doublet Score {.tabset .tabset-fade} \n")
  print(plotsScrublet$Sample[[sample]][["scatter_doubletScore"]])
  cat("\n\n")

  cat("#### Density Score {.tabset .tabset-fade} \n")
  print(plotsScrublet$Sample[[sample]][["density_doubletScore"]])
  cat("\n\n")

  cat("#### Violin Score {.tabset .tabset-fade} \n ")
  print(plotsScrublet$Sample[[sample]][["violin_doubletScore"]])
  cat("\n\n")

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  runScrubletMeta <- sce.qc@metadata$runScrublet
  runScrubletMeta <- metaChecking(runScrubletMeta, sample)
  # if(length(samples) == 1) {
  #   runScrubletMeta <- list(runScrubletMeta)
  #   names(runScrubletMeta) <- sample
  # }

  x <- preprocessMeta(runScrubletMeta[[sample]])
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))
  cat(description_Scrublet$additionalParam)
  cat("\n\n")
}
cat("\n")

description_DoubletFinder <- singleCellTK:::descriptionDoubletFinder()

i="DoubletFinder"

doubletFinderData <- grep("doubletFinder_doublet", names(colData(sce.qc)), value=TRUE)
skipFinder <- length(doubletFinderData)==0  

if (skipFinder) {
  #cat("runDoubletFinder did not complete successfully. The section of DoubletFinder is skipped")
  plotDoubletFinder <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotDoubletFinder<- tryCatch(
    {plotDoubletFinderResults(inSCE = sce.qc, 
                              reducedDimName = "QC_UMAP",
                              sample = colData(sce.qc)[['sample']],
                              combinePlot = "none",
                              axisSize = NULL,
                              axisLabelSize = NULL,
                              titleSize = NULL)},
    error = function(e) {
      cat("runDoubletFinder did not complete successfully. The section of DoubletFinder is skipped")
      skipFinder <<- TRUE
      return(NA)
    }
  )
}

if (!skipFinder) {
  plotDoubletFinder <- formatPlotList(plotDoubletFinder, samples)
  resolutions <- grep("doubletFinder_doublet_score_resolution_", colnames(colData(sce.qc)), value = TRUE)
  resolutions <- gsub("doubletFinder_doublet_score_resolution_", "", resolutions)

  cat(description_DoubletFinder$introduction)
  cat(description_DoubletFinder$runDoubletFinder)
  cat(description_DoubletFinder$plotDoubletFinderResults)
  cat(description_DoubletFinder$output)
  cat("\n")
}

cat("## {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
cat("\n")
for (sample in samples) {
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n"))
  cat("\n")
  for (resolution in resolutions) {
    cat(paste0("#### Resolution: ", resolution, " {.tabset .tabset-fade} \n"))
    cat("\n")

    cat("##### DoubletFinder Doublet Assignment {.tabset .tabset-fade} \n\n")
    print(plotDoubletFinder$Sample[[sample]][[paste0("Scatter_Call_resolution_", resolution)]])
    cat("\n\n")

    cat("##### DoubletFinder Doublet Score {.tabset .tabset-fade} \n\n")
    print(plotDoubletFinder$Sample[[sample]][[paste0("Scatter_Score_resolution_", resolution)]])
    cat("\n\n")

    cat("##### Density of Doublet Score {.tabset .tabset-fade} \n\n")
    print(plotDoubletFinder$Sample[[sample]][[paste0("Density_resolution_", resolution)]])
    cat("\n\n")

    cat("##### Violin of Doublet Score {.tabset .tabset-fade} \n\n")
    print(plotDoubletFinder$Sample[[sample]][[paste0("violin_resolution_", resolution)]])
    cat("\n\n")

    cat("##### Parameters {.tabset .tabset-fade} \n\n")
    DoubletFinderMeta <- sce.qc@metadata$runDoubletFinder
    DoubletFinderMeta <- metaChecking(DoubletFinderMeta, sample)

    # if(length(samples) == 1) {
    #   DoubletFinderMeta <- list(DoubletFinderMeta)
    #   names(DoubletFinderMeta) <- sample
    # }

    x <- preprocessMeta(DoubletFinderMeta[[sample]])
    cat(t(as.data.frame(x)) %>%
      knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
      kableExtra::scroll_box(width = "80%"))

    cat(description_DoubletFinder$seuratRes)
    cat(description_DoubletFinder$seuratNfeatures)
    cat(description_DoubletFinder$seuratPCs)
    cat(description_DoubletFinder$seuratFormationRate)
    cat("\n\n")
  }
}
cat("\n")

description_ScDblFinder<- singleCellTK:::descriptionScDblFinder()

i="ScDblFinder"
#cat(paste0('## ', i, ' \n'))

scDblFinderData <- c("scDblFinder_doublet_score")
skipScDblFinder <- any(!scDblFinderData %in% names(colData(sce.qc)))

if (skipScDblFinder) {
  #cat("runScDblFinder did not complete successfully. The section of ScDblFinder is skipped")
  plotScDblFinder <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotScDblFinder<- tryCatch(
    {plotScDblFinderResults(inSCE = sce.qc, 
                            sample = colData(sce.qc)[['sample']],
                            combinePlot = "none",
                            reducedDimName = "QC_UMAP",
                            axisSize = NULL,
                            axisLabelSize = NULL,
                            titleSize = NULL)},
    error = function(e) {
      cat("runScDblFinder did not complete successfully. The section of ScDblFinder is skipped")
      skipScDblFinder <<- TRUE
      return(NA)
    }
  )
}

if(!skipScDblFinder) {
  plotScDblFinder <- formatPlotList(plotScDblFinder, samples)
  cat(description_ScDblFinder$introduction)
  cat(description_ScDblFinder$runScDblFinder)
  cat(description_ScDblFinder$plotScDblFinderResults)
  cat(description_ScDblFinder$output)
  cat("\n")
}

cat("## {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
cat("\n")
for (sample in samples){
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n"))
  cat("\n\n")

  cat("#### ScDblFinder Doublet Assignment {.tabset .tabset-fade} \n\n ")
  print(plotScDblFinder$Sample[[sample]][["scatter_doubletCall"]])
  cat("\n\n")

  cat("#### ScDblFinder Doublet Score {.tabset .tabset-fade} \n\n ")
  print(plotScDblFinder$Sample[[sample]][["scatter_doubletScore"]])
  cat("\n\n")

  cat("#### Density Score {.tabset .tabset-fade} \n\n")
  print(plotScDblFinder$Sample[[sample]][["density_doubletScore"]])
  cat("\n\n")

  cat("#### Violin Score {.tabset .tabset-fade} \n\n")
  print(plotScDblFinder$Sample[[sample]][["violin_doubletScore"]])
  cat("\n\n")

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  ScDblFinderMeta <- sce.qc@metadata$runScDblFinder
  ScDblFinderMeta <- metaChecking(ScDblFinderMeta, sample)

  # if(length(samples) == 1) {
  #   DoubletCellMeta <- list(DoubletCellMeta)
  #   names(DoubletCellMeta) <- sample
  # }

  x <- preprocessMeta(ScDblFinderMeta[[sample]])
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))
  cat(description_ScDblFinder$parameter)
  cat("\n\n")
}
cat("\n")

descriprion_CxdsResults<- singleCellTK:::descriptionCXDS()

i="Cxds"
# cat(paste0('## ', i, ' \n'))

cvdsData <- c("scds_cxds_score", "scds_cxds_call")
skipCxds <- any(!cvdsData %in% names(colData(sce.qc)))

if (skipCxds) {
  #cat("runCxds did not complete successfully. The section of Cxds is skipped")
  plotCxds <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotCxds<- tryCatch(
    {plotCxdsResults(inSCE = sce.qc, 
                     sample = colData(sce.qc)[['sample']],
                     combinePlot = "none",
                     reducedDimName = "QC_UMAP",
                     axisSize = NULL,
                     axisLabelSize = NULL,
                     titleSize = NULL)},
    error = function(e) {
      cat("runCxds did not complete successfully. The section of Cxds is skipped")
      skipCxds <<- TRUE
      return(NA)
    }
  )
}

if (!skipCxds) {
  plotCxds <- formatPlotList(plotCxds, samples)
  cat(descriprion_CxdsResults$introduction)
  cat(descriprion_CxdsResults$runCxds)
  cat(descriprion_CxdsResults$plotCxdsResults)
  cat(descriprion_CxdsResults$output)
  cat("\n")
}

cat("## {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
cat("\n")
for (sample in samples){
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n"))
  cat("\n\n")

  cat("#### Cxds Doublet Assignment {.tabset .tabset-fade} \n\n")
  plot(plotCxds$Sample[[sample]][["scatter_doubletCall"]])
  cat("\n\n")

  cat("#### Cxds Doublet Score {.tabset .tabset-fade} \n\n")
  plot(plotCxds$Sample[[sample]][["scatter_doubletScore"]])
  cat("\n\n")

  cat("#### Density Score {.tabset .tabset-fade} \n\n")
  plot(plotCxds$Sample[[sample]][["density_doubletScore"]])
  cat("\n\n")

  cat("#### Violin Score {.tabset .tabset-fade} \n\n")
  plot(plotCxds$Sample[[sample]][["violin_doubletScore"]])
  cat("\n\n")

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  CxdsMeta <- sce.qc@metadata$runCxds
  CxdsMeta <- metaChecking(CxdsMeta, sample)

  # if(length(samples) == 1) {
  #   CxdsMeta <- list(CxdsMeta)
  #   names(CxdsMeta) <- sample
  # }

  x <- preprocessMeta(CxdsMeta[[sample]])
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))

  cat(descriprion_CxdsResults$nTop)
  cat(descriprion_CxdsResults$binThresh)
  cat(descriprion_CxdsResults$verb)
  cat(descriprion_CxdsResults$retRes)
  cat(descriprion_CxdsResults$estNdbl)
  cat("\n\n")

}
cat("\n")

descriprion_BcdsResults<- singleCellTK:::descriptionBCDS()

i="Bcds"
# cat(paste0('## ', i, ' \n'))

bcdsData <- c("scds_bcds_score", "scds_bcds_call")
skipBcds <- any(!bcdsData %in% names(colData(sce.qc)))

if (skipBcds) {
  #cat("runBcds did not complete successfully. The section of Bcds is skipped")
  plotBcds <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotBcds<- tryCatch(
    {plotBcdsResults(inSCE = sce.qc, 
                     sample = colData(sce.qc)[['sample']],
                     combinePlot = "none",
                     reducedDimName = "QC_UMAP",
                     axisSize = NULL,
                     axisLabelSize = NULL,
                     titleSize = NULL)},
    error = function(e) {
      cat("runBcds did not complete successfully. The section of Bcds is skipped")
      skipBcds <<- TRUE
      return(NA)
    }
  )
}

if (!skipBcds) {
  plotBcds <- formatPlotList(plotBcds, samples)
  cat(descriprion_BcdsResults$introduction)
  cat(descriprion_BcdsResults$runBcds)
  cat(descriprion_BcdsResults$plotBcdsResults)
  cat(descriprion_BcdsResults$output)
  cat("\n")
}

cat("## {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
cat("\n")
for (sample in samples){
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n"))
  cat("\n\n")

  cat("#### Bcds Doublet Assignment {.tabset .tabset-fade} \n\n")
  plot(plotBcds$Sample[[sample]][["scatter_doubletCall"]])
  cat("\n\n")

  cat("#### Bcds Doublet Score {.tabset .tabset-fade} \n\n")
  plot(plotBcds$Sample[[sample]][["scatter_doubletScore"]])
  cat("\n\n")

  cat("#### Density Score {.tabset .tabset-fade} \n\n")
  plot(plotBcds$Sample[[sample]][["density_doubletScore"]])
  cat("\n\n")

  cat("#### Violin Score {.tabset .tabset-fade} \n\n")
  plot(plotBcds$Sample[[sample]][["violin_doubletScore"]])
  cat("\n\n")

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  BcdsMeta <- sce.qc@metadata$runBcds
  BcdsMeta <- metaChecking(BcdsMeta, sample)

  # if(length(samples) == 1) {
  #   BcdsMeta <- list(BcdsMeta)
  #   names(BcdsMeta) <- sample
  # }

  x <- preprocessMeta(BcdsMeta[[sample]])
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))

  cat(descriprion_BcdsResults$nTop)
  cat(descriprion_BcdsResults$srat)
  cat(descriprion_BcdsResults$nmax)
  cat(descriprion_BcdsResults$varImp)
  cat("\n\n")

}
cat("\n")

description_ScdsHybrid<- singleCellTK:::descriptionScdsHybrid()

i="ScdsHybrid"
# cat(paste0('## ', i, ' \n'))

hybridData <- c("scds_hybrid_score", "scds_hybrid_call")
skipScdsHybrid <- any(!hybridData %in% names(colData(sce.qc)))

if (skipScdsHybrid) {
  #cat("runCxdsBcdsHybrid did not complete successfully. The section of ScdsHybrid is skipped")
  plotScdsHybrid <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotScdsHybrid<- tryCatch(
    {plotScdsHybridResults(inSCE = sce.qc, 
                           sample = colData(sce.qc)[['sample']],
                           combinePlot = "none",
                           reducedDimName = "QC_UMAP",
                           axisSize = NULL,
                           axisLabelSize = NULL,
                           titleSize = NULL)},
    error = function(e) {
      cat("runCxdsBcdsHybrid did not complete successfully. The section of ScdsHybrid is skipped")
      skipScdsHybrid <<- TRUE
      return(NA)
    }
  )
}

if (!skipScdsHybrid) {
  plotScdsHybrid <- formatPlotList(plotScdsHybrid, samples)
  cat(description_ScdsHybrid$introduction)
  cat(description_ScdsHybrid$runCxdsBcdsHybrid)
  cat(description_ScdsHybrid$plotScdsHybridResults)
  cat(description_ScdsHybrid$output)
  cat("\n")
}

cat("## {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
cat("\n")
for (sample in samples){
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n"))
  cat("\n\n")

  cat("#### ScdsHybrid Doublet Assignment {.tabset .tabset-fade} \n\n")
  plot(plotScdsHybrid$Sample[[sample]][["scatter_doubletCall"]])
  cat("\n\n")

  cat("#### ScdsHybrid Doublet Score {.tabset .tabset-fade} \n\n")
  plot(plotScdsHybrid$Sample[[sample]][["scatter_doubletScore"]])
  cat("\n\n")

  cat("#### Density Score {.tabset .tabset-fade} \n\n")
  plot(plotScdsHybrid$Sample[[sample]][["density_doubletScore"]])
  cat("\n\n")

  cat("#### Violin Score {.tabset .tabset-fade} \n\n")
  plot(plotScdsHybrid$Sample[[sample]][["violin_doubletScore"]])
  cat("\n\n")

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  HybridMeta <- sce.qc@metadata$runCxdsBcdsHybrid
  HybridMeta <- metaChecking(HybridMeta, sample)

  # if(length(samples) == 1) {
  #   HybridMeta <- list(HybridMeta)
  #   names(HybridMeta) <- sample
  # }

  x <- preprocessMeta(HybridMeta[[sample]])
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))

  cat(description_ScdsHybrid$parameters)
  cat("\n\n")
}
cat("\n")

r cat("\n")

Ambient RNA Detection

### get all doublet detection algorithms
fontSize <- 15
amBAlgo <- setNames(
  c("decontX_contamination", "decontX_contamination_bg"),
  c("decontX", "decontX with background")
  )

if (length(amBAlgo) != 0) {
  i="Ambient RNA Detection Summary"
  cat(paste0('## ', i, ' \n'))
  cat("## {.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}")
  cat("\n")

  for (algo in amBAlgo) {
    if (!algo %in% names(colData(sce.qc))) {
      next
    }
    algoName <- names(amBAlgo)[amBAlgo == algo]
    cat(paste0("### ", algoName, " {.tabset .tabset-fade} \n"))
    cat("\n\n")

    plotAmbient <- plotSCEDimReduceColData(
        inSCE=sce.qc,
        sample=NULL,
        colorBy=algo,
        conditionClass="numeric",
        colorScale = c("lightgray","red"),
        reducedDimName=allSampleReducedDim,
        axisSize = fontSize,
        axisLabelSize = fontSize,
        legendSize = fontSize,
        titleSize = fontSize,
        xlab=c("QC_UMAP_1"),
        ylab=c("QC_UMAP_2"),
        labelClusters=FALSE,
        title=algoName,
        legendTitle = "Ambient \nContamination \nScore",
        plotLabels = NULL,
        combinePlot="none"
      )

    print(plot_grid(sampleUMAP, plotAmbient, ncol = 2))
    cat("\n\n")
  }
}

description_DecontX<- singleCellTK:::descriptionDecontX()

i="DecontX"
# cat(paste0('## ', i, ' \n'))

decontXData <- c("decontX_contamination", "decontX_clusters")
skipDecontX <- any(!decontXData %in% names(colData(sce.qc)))

if (skipDecontX) {
  #cat("runDecontX did not complete successfully. The section of DecontX is skipped")
  plotDecontX <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotDecontX<- tryCatch(
    {plotDecontXResults(inSCE = sce.qc, 
                        sample = colData(sce.qc)[['sample']],
                        combinePlot = "none",
                        reducedDimName = "QC_UMAP",
                        axisSize = NULL,
                        axisLabelSize = NULL,
                        titleSize = NULL)},
    error = function(e) {
      cat("runDecontX did not complete successfully. The section of DecontX is skipped")
      skipDecontX <<- TRUE
      return(NA)
    }
  )
}

if (!skipDecontX) {
  plotDecontX <- formatPlotList(plotDecontX, samples)
  cat(description_DecontX$introduction)
  cat(description_DecontX$runDecontX)
  cat(description_DecontX$plotDecontXResults)
  cat(description_DecontX$output)
  cat(description_DecontX$contamination)
  cat(description_DecontX$clustering)
}

{.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}

cat("\n")
for (sample in samples){
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n"))
  cat("\n\n")

  cat("#### DecontX Contamination Score {.tabset .tabset-fade} \n\n")
  plot(plotDecontX$Sample[[sample]][["scatter_decontXContamination"]])
  cat("\n\n")

  cat("#### DecontX Clusters {.tabset .tabset-fade} \n\n")
  plot(plotDecontX$Sample[[sample]][["scatter_decontXClusters"]])
  cat("\n\n")

  cat("#### Density Score {.tabset .tabset-fade} \n\n")
  plot(plotDecontX$Sample[[sample]][["density_decontXContamination"]])
  cat("\n\n")

  cat("#### Violin Score {.tabset .tabset-fade} \n\n")
  plot(plotDecontX$Sample[[sample]][["violin_decontXContamination"]])
  cat("\n\n")

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  DecontXMeta <- sce.qc@metadata$sctk$runDecontX
  DecontXMeta <- metaChecking(DecontXMeta, sample)

  # if(length(samples) == 1) {
  #   DecontXMeta <- list(DecontXMeta)
  #   names(DecontXMeta) <- sample
  # }

  x <- preprocessMeta(DecontXMeta[[sample]])
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))
  cat("\n\n")
}

r cat("\n")

i="DecontX with raw/droplet matrix"
# cat(paste0('## ', i, ' \n'))

decontXData_bg <- c("decontX_contamination_bg", "decontX_clusters_bg")
skipDecontX_bg <- any(!decontXData_bg %in% names(colData(sce.qc)))

if (skipDecontX_bg) {
  #cat("runDecontX did not complete successfully. The section of DecontX is skipped")
  plotDecontX_bg <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotDecontX_bg<- tryCatch(
    {plotDecontXResults(inSCE = sce.qc, 
                        sample = colData(sce.qc)[['sample']],
                        bgResult = TRUE,
                        combinePlot = "none",
                        reducedDimName = "QC_UMAP",
                        axisSize = NULL,
                        axisLabelSize = NULL,
                        titleSize = NULL)},
    error = function(e) {
      cat("runDecontX did not complete successfully. The section of DecontX_bg is skipped")
      skipDecontX_bg <<- TRUE
      return(NA)
    }
  )
}

if (!skipDecontX_bg) {
  plotDecontX_bg <- formatPlotList(plotDecontX_bg, samples)
  cat(description_DecontX$introduction)
  cat(description_DecontX$runDecontX)
  cat(description_DecontX$plotDecontXResults)
  cat(description_DecontX$output)
  cat(description_DecontX$contamination)
  cat(description_DecontX$clustering)
}

{.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}

cat("\n")
for (sample in samples){
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n"))
  cat("\n\n")

  cat("#### DecontX Contamination Score {.tabset .tabset-fade} \n\n")
  plot(plotDecontX_bg$Sample[[sample]][["scatter_decontXContamination"]])
  cat("\n\n")

  cat("#### DecontX Clusters {.tabset .tabset-fade} \n\n")
  plot(plotDecontX_bg$Sample[[sample]][["scatter_decontXClusters"]])
  cat("\n\n")


  cat("#### Density Score {.tabset .tabset-fade} \n\n")
  plot(plotDecontX_bg$Sample[[sample]][["density_decontXContamination"]])
  cat("\n\n")

  cat("#### Violin Score {.tabset .tabset-fade} \n\n")
  plot(plotDecontX_bg$Sample[[sample]][["violin_decontXContamination"]])
  cat("\n\n")

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  DecontXMeta_bg <- sce.qc@metadata$sctk$runDecontX_bg
  DecontXMeta_bg <- metaChecking(DecontXMeta_bg, sample)

  # if(length(samples) == 1) {
  #   DecontXMeta <- list(DecontXMeta)
  #   names(DecontXMeta) <- sample
  # }

  x <- preprocessMeta(DecontXMeta_bg[[sample]])
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))
  cat("\n\n")
}

r cat("\n")

description_SoupX <- singleCellTK:::descriptionSoupX()

i="SoupX"
# cat(paste0('## ', i, ' \n'))

soupXData <- c("soupX_nUMIs", "soupX_contamination")
# soupX cluster label might be supplied by users, thus not checking it.
skipSoupX <- any(!soupXData %in% names(colData(sce.qc)))

if (skipSoupX) {
  #cat("runDecontX did not complete successfully. The section of DecontX is skipped")
  plotSoupX <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotSoupX <- tryCatch(
    {plotSoupXResults(inSCE = sce.qc, 
                      sceSample, 
                      combinePlot = "none")},
    error = function(e) {
      cat("runSoupX did not complete successfully. The section of SoupX is skipped\n\n")
      skipSoupX <<- TRUE
      return(NA)
    }
  )
}

if (!skipSoupX) {
  plotSoupX <- formatPlotList(plotSoupX, samples)
  cat(description_SoupX$introduction)
  cat(description_SoupX$runSoupX)
  cat(description_SoupX$output)
  cat(description_SoupX$contamination)
  cat(description_SoupX$clustering)
  cat(description_SoupX$plotSoupXResults)
}

{.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}

cat("\n\n")
for (sample in samples){
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n\n"))

  cat("#### SoupX Clustering {.tabset .tabset-fade} \n\n")
  plot(plotSoupX$Sample[[sample]][["scatter_soupXClusters"]])
  cat("\n\n")

  cat("#### Soup Fractions {.tabset .tabset-fade} \n\n")

  for (marker in names(plotSoupX$Sample[[sample]])[-1]) {
    cat(paste0("##### ", marker, " {.tabset .tabset-fade} \n\n"))
    plot(plotSoupX$Sample[[sample]][[marker]])
    cat("\n\n")
  }

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  if (length(samples) == 1) {
    SoupXMeta <- getSoupX(sce.qc, soupXSample, background = FALSE)[[1]]$param
    SoupXMeta <- metaChecking(SoupXMeta, soupXSample)
    x <- preprocessMeta(SoupXMeta[[soupXSample]])
  } else {
    SoupXMeta <- getSoupX(sce.qc, sample, background = FALSE)[[1]]$param
    SoupXMeta <- metaChecking(SoupXMeta, sample)
    x <- preprocessMeta(SoupXMeta[[sample]])
  }
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))
  cat("\n\n")
}

r cat("\n")

i="SoupX with raw/droplet matrix"
# cat(paste0('## ', i, ' \n'))

soupXData_bg <- c("soupX_bg_nUMIs", "soupX_bg_contamination")
skipSoupX_bg <- any(!soupXData_bg %in% names(colData(sce.qc)))

if (length(samples) == 1) {
  soupXSample <- "all_cells"
} else {
  soupXSample <- samples
}

if (skipSoupX_bg) {
  plotSoupX <- NULL
} else {
  cat(paste0('## ', i, ' \n'))
  plotSoupX_bg <- tryCatch(
    {plotSoupXResults(inSCE = sce.qc, 
                     soupXSample, 
                     combinePlot = "none",
                     background = TRUE)},
    error = function(e) {
      cat("runSoupX with background did not complete successfully. The section of SoupX with background is skipped")
      skipSoupX_bg <<- TRUE
      return(NA)
    }
  )
}

if (!skipSoupX_bg) {
  plotSoupX_bg <- formatPlotList(plotSoupX_bg, samples)
  cat(description_SoupX$introduction)
  cat(description_SoupX$runSoupX)
  cat(description_SoupX$output)
  cat(description_SoupX$contamination)
  cat(description_SoupX$clustering)
  cat(description_SoupX$plotSoupXResults)
}

{.unlisted .unnumbered .toc-ignore .tabset .tabset-fade}

cat("\n\n")
for (sample in samples){
  cat(paste0("### ", sample, " {.tabset .tabset-fade} \n\n"))

  cat("#### SoupX Clustering {.tabset .tabset-fade} \n\n")
  plot(plotSoupX_bg$Sample[[sample]][["scatter_soupXClusters"]])
  cat("\n\n")

  cat("#### Soup Fractions {.tabset .tabset-fade} \n\n")

  for (marker in names(plotSoupX$Sample[[sample]])[-1]) {
    cat(paste0("##### ", marker, " {.tabset .tabset-fade} \n\n"))
    plot(plotSoupX_bg$Sample[[sample]][[marker]])
    cat("\n\n")
  }

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  if (length(samples) == 1) {
    SoupXMeta <- getSoupX(sce.qc, soupXSample, background = TRUE)[[1]]$param
    SoupXMeta <- metaChecking(SoupXMeta, soupXSample)
    x <- preprocessMeta(SoupXMeta[[soupXSample]])
  } else {
    SoupXMeta <- getSoupX(sce.qc, sample, background = TRUE)[[1]]$param
    SoupXMeta <- metaChecking(SoupXMeta, sample)
    x <- preprocessMeta(SoupXMeta[[sample]])
  }
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))
  cat("\n\n")
}

r cat("\n")

Session Information

Session Information

sessionInfo()

compbiomed/singleCellTK documentation built on May 8, 2024, 6:58 p.m.