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

if(!exists("titleMarkerPlots")) titleMarkerPlots <- "Plot Features"
if(!exists("groupTitle")) groupTitle <- "Seurat_louvain_Resolution0.3"
if(!exists("selectedOption")) selectedOption <- "Seurat_louvain_Resolution0.3"
if(!exists("headingMS")) headingMS <- "#"

library(Seurat)
library(dplyr)
library(cowplot)
library(RColorBrewer)
library(ggplot2)
library(knitr)
library(kableExtra)
library(SingleCellExperiment)
library(scater)
library(gridExtra)
library(grid)
library(ggpubr)
library(patchwork)
library(singleCellTK)

cat(headingMS, " Differential Expression (", groupTitle,")\n\n")
cat("Gene markers that are differentially expressed between the previously computed ", selectedOption, " are identified using the Wilcoxon Rank Sum test. Criteria for the identification of these marker genes required them to be detected at a minimum percentage of 25% in either of the two groups of cells (each cluster against all others), and additionally they show 25% log-fold expression difference between two groups.")

data <- runSeuratFindMarkers(data, allGroup = selectedOption, minPCT = 0.25, threshUse = 0.25, verbose = FALSE)
data.markers <- metadata(data)$seuratMarkers

write.csv(data.markers, file = paste0(outputPath, "DEGTable", selectedOption, "-", gsub(" ", "_", Sys.Date()), ".csv"), row.names = FALSE)

topFeaturesTable <- data.frame(data.markers %>% group_by(cluster1) %>% top_n(numTopFeatures, avg_log2FC))
numClusters <- length(unique(data.markers$cluster1))

topFeaturesTable <- data.frame(gene.id = selected.markers, cluster1 = 1)
data.markers <- data.frame(gene.id = selected.markers, cluster1 = 1)

cat(paste0(headingMS, "# ", titleMarkerPlots," {.tabset .tabset-fade}\n\n"))
cat("Gene plots for the top marker genes are visualized below through a feature plot, ridge plot, violin plot, dot plot and a heatmap plot.")

cat(paste0(headingMS, "## Feature Plot {.tabset .tabset-fade .tabset-pills -}\n\n"))
cat("A feature plot visualizes the expression level of a particular gene marker in all cells of the data on a UMAP plot. The feature plots of the top ", numTopFeatures, " marker genes across all clusters are visualized below:")

figHeight <- numTopFeatures * 1.8

template <- paste0(headingMS, "### %s {.tabset .tabset-fade -}
")
template_inside <- paste0(headingMS, "#### %s {-}
")
space <- "

"

for (currentCluster in 0:(numClusters-1)) {
  if(numClusters > 1){
    cat(sprintf(template, paste0("Cluster ", unique(topFeaturesTable$cluster1)[currentCluster+1])))
  }
  selectedFeatures <- topFeaturesTable[topFeaturesTable$cluster1 == unique(data.markers$cluster1)[currentCluster+1], ]$gene.id
  plots <- plotSeuratGenes(data, plotType = "feature", features = selectedFeatures, groupVariable = selectedOption, ncol = 1, cols = c("grey", "blue"))
  p <- singleCellTK:::.ggSCTKTheme(patchwork::wrap_plots(plots, ncol = 2))
  print(p)
  cat(space)
}

cat(paste0(headingMS, "## Ridge Plot {.tabset .tabset-fade .tabset-pills -}\n\n"))
cat("A ridge plot visualizes the expression level of a particular gene marker in all cells of the data separated by clusters in the form of ridges. The ridge plots of the top ", numTopFeatures, " marker genes across all clusters are visualized below:")

for (currentCluster in 0:(numClusters-1)) {
  if(numClusters > 1){
    cat(sprintf(template, paste0("Cluster ", unique(topFeaturesTable$cluster1)[currentCluster+1])))
  }
  selectedFeatures <- topFeaturesTable[topFeaturesTable$cluster1 == unique(data.markers$cluster1)[currentCluster+1], ]$gene.id
  plots <- plotSeuratGenes(data, plotType = "ridge", features = selectedFeatures, groupVariable = selectedOption, ncol = 1, cols = c("grey", "blue"))
  p <- singleCellTK:::.ggSCTKTheme(patchwork::wrap_plots(plots, ncol = 2))
  print(p)
  cat(space)
}

cat(paste0(headingMS, "## Violin Plot {.tabset .tabset-fade .tabset-pills -}\n\n"))
cat("A violin plot visualizes the expression level of a particular gene marker in all cells of the data separated by clusters. The violin plots of the top ", numTopFeatures, " marker genes separated across all clusters are visualized below:")

for (currentCluster in 0:(numClusters-1)) {
  if(numClusters > 1){
    cat(sprintf(template, paste0("Cluster ", unique(topFeaturesTable$cluster1)[currentCluster+1])))
  }
  selectedFeatures <- topFeaturesTable[topFeaturesTable$cluster1 == unique(data.markers$cluster1)[currentCluster+1], ]$gene.id
  plots <- plotSeuratGenes(data, plotType = "violin", features = selectedFeatures, groupVariable = selectedOption, ncol = 1, cols = c("grey", "blue"))
  p <- singleCellTK:::.ggSCTKTheme(patchwork::wrap_plots(plots, ncol = 2))
  print(p)
  cat(space)
}

cat(paste0(headingMS, "## Dot Plot {-}\n\n"))
cat("A dot plot visualizes the change in the expression levels of a particular gene marker in all cells of the data separated by clusters in the form of dots where the color of the dots indicate the expression level of that gene in a cluster and the size of the dot indicates the total percentage of the cells in that cluster.")

plotSeuratGenes(data, plotType = "dot", features = topFeaturesTable$gene.id, groupVariable = selectedOption, ncol = 1, cols = c("grey", "blue")) + theme(axis.text.x = element_text(angle = 90, size = 7), legend.text = element_text(size = 8), legend.key.size = unit(0.4, "cm"), legend.position = "top")

cat(paste0(headingMS, "## Heatmap Plot {-}\n\n"))

plotSeuratGenes(data, plotType = "heatmap", features = topFeaturesTable$gene.id, groupVariable = selectedOption, ncol = 1, cols = c("grey", "blue")) + theme(axis.text.y = element_text(size = 8), legend.text = element_text(size = 8), legend.key.size = unit(0.5, "cm"), legend.position = "top")

cat(paste0(headingMS, "# Table of Top Markers by Cluster\n\n"))
cat("The list of top ", numTopFeatures, " genes for each clusters are listed below:")

topFeaturesTable <- data.frame(data.markers %>% group_by(cluster1) %>% top_n(numTopFeatures, avg_log2FC))

colNames<-colnames(topFeaturesTable)

kable(topFeaturesTable, style = 'html', row.names = FALSE) %>%
  kable_styling(bootstrap_options = "striped") %>%
  scroll_box(width = "100%", height = "500px")

cat("# Session Information\n\n")

sessionInfo()

compbiomed/singleCellTK documentation built on Oct. 27, 2024, 3:26 a.m.