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

require(SummarizedExperiment)
require(SingleCellExperiment)
require(singleCellTK)
require(dplyr)
require(Biobase)

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

preprocessMeta <- function(meta, ignore=c("geneSetCollection", "sample")){
  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)
  })
}

sce<-params$object
subTitle <- params$subTitle
studyDesign <- params$studyDesign
samples <- unique(colData(sce)$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)

Empty Drops Detection

description_runEmptyDrops <- descriptionEmptyDrops()

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


emptyData <- c("dropletUtils_emptyDrops_logprob", "dropletUtils_emptyDrops_total")
skipRunEmpty <- any(!emptyData %in% names(colData(sce)))

if (skipRunEmpty) {
  cat("runEmptyDrops did not run successfully. The section of EmptyDrops is skipped")
  plotsEmptyDrops <- NULL
} else {
    plotsEmptyDrops <- tryCatch(
    { plotEmptyDropsResults(sce,sample = colData(sce)$sample,
                                             combinePlot = "none") },
      error = function(e) {
        cat("runEmptyDrops did not run successfully. The section of EmptyDrops is skipped")
        skipRunEmpty <- TRUE
        return(NA)
      }
  )
}

if(!skipRunEmpty) {
  cat(paste0('## ', i, ' \n'))
  cat("\n")
  if(length(samples) == 1) {
    plotsEmptyDrops$scatterEmptyDrops <- list(Sample = list(plotsEmptyDrops$scatterEmptyDrops))
    names(plotsEmptyDrops$scatterEmptyDrops$Sample) <- samples
  }

  cat(description_runEmptyDrops$introduction)
  cat(description_runEmptyDrops$runEmptyDrops)
  cat(description_runEmptyDrops$plotEmptyDropsResults)
  cat(description_runEmptyDrops$plot2)
}

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

  cat("#### Total UMI counts vs Log-Probability \n")
  print(plotsEmptyDrops$scatterEmptyDrops$Sample[[i]])
  cat("\n\n")

  cat("#### Parameters \n\n")
  runEmptyMeta <- sce@metadata$runEmptyDrops
  if(length(samples) == 1) {runEmptyMeta <- list(runEmptyMeta)}

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

description_runBarcodeRank <- descriptionBarcodeRank()

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


barcodeData <- c("dropletUtils_barcodeRank_rank", "dropletUtils_barcodeRank_total", "dropletUtils_barcodeRank_knee", "dropletUtils_barcodeRank_inflection",  "dropletUtils_barcodeRank_fitted") #"sample",

skipBarcode <- (!"runBarcodeRanksMetaOutput" %in% names(S4Vectors::metadata(sce))) || 
              (any(!barcodeData %in% names(S4Vectors::metadata(sce)$runBarcodeRanksMetaOutput[[1]])))

if (skipBarcode) {
  cat("runBarcodeRankDrops did not run successfully. The section of BarcodeRanks is skipped")
  plotsBarcodeRank <- NA
} else {
  plotsBarcodeRank <- tryCatch(
    {plotBarcodeRankScatter(sce,
                            sample = colData(sce)$sample, 
                            title="BarcodeRanks Rank Plot",
                            combinePlot = "none")},
      error = function(e) {
        cat("runBarcodeRankDrops did not run successfully. The section of BarcodeRanks is skipped")
        skipBarcode <- TRUE
        return(NA)
      }
  )
}

if(!skipBarcode) {
  cat(paste0('## ', i, ' \n'))
  cat("\n")
  if(length(samples) == 1) {
    plotsBarcodeRank <- list(Sample = list(plotsBarcodeRank))
    names(plotsBarcodeRank$Sample) <- samples
  }

  #names(plotsBarcodeRank) <- samples
  cat(description_runBarcodeRank$introduction)
  cat(description_runBarcodeRank$runBarcodeRankDrops)
  cat(description_runBarcodeRank$plotBarcodeRankDropsResults)
  cat(description_runBarcodeRank$plot)
}

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

  cat("#### Total UMI counts vs Rank {.tabset .tabset-fade} \n\n")
  print(plotsBarcodeRank$Sample[[samples[i]]])
  cat("\n\n")

  cat("#### Parameters {.tabset .tabset-fade} \n\n")
  runBarcodeMeta <- sce@metadata$runBarcodeRankDrops
  if (length(samples) == 1) { runBarcodeMeta <- list(runBarcodeMeta) }
  x <- preprocessMeta(runBarcodeMeta[[i]])
  cat(t(as.data.frame(x)) %>%
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))
  cat(description_runBarcodeRank$parameter)
  cat("\n\n")  
}
cat("\n")

Summary Statistics {.tabset .tabset-fade}

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

cat("\n")

for (i in seq_along(samples)) {
  cat(paste0("## ", samples[i], " {.tabset .tabset-fade} \n"))
  cat("\n\n")
  sceColData <- colData(sce) %>% as.data.frame()
  sceColData <- sceColData[sceColData$sample == samples[i], ]

  idx <- list(
          'Knee'= sceColData$dropletUtils_BarcodeRank_Knee == 1,
          'Inflection'= sceColData$dropletUtils_BarcodeRank_Inflection == 1, 
          'EmptyDrops' = !is.na(sceColData$dropletUtils_emptyDrops_fdr) &
           (sceColData$dropletUtils_emptyDrops_fdr < 0.01))

  SummaryTable <- lapply(idx, function(x){
  sceColData[x,] %>% 
    dplyr::summarise('Number of cells' = n(), 'Median UMI' = median(sum), 
                     'Median Genes' = median(detected))
  })


  cat(do.call(base::rbind, SummaryTable) %>% 
    knitr::kable(format = "html") %>% kableExtra::kable_styling() %>%
    kableExtra::scroll_box(width = "80%"))

  summary <- "The summary statistics table summarizes QC metrics of droplets that passed different cell detection cutoff. Here it summarizes number of cells passing the cutoff, median of UMI count and median of genes detected per cell."
  cat(summary)
  cat("\n\n") 
}

r cat("\n")

SessionInfo

sessionInfo()

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