---

title: "Analysis of scRNA-seq data (compilation)"

date: "Last modified: `r format(Sys.time(), '%B %d, %Y')`"

tags: [scRNA-seq, melanoma, immunotherapy, PBMCs] 

output:

  html_document:

    theme: flatly

    highlight: zenburn

    toc: true

    number_sections: false

    toc_depth: 2

    toc_float:

      collapsed: false

      smooth_scroll: true

    code_folding: hide

    self_contained: yes

---

# Project Description

This project aims to identify therapeutic targets for melanoma patients based on scRNAseq data analysis.  

## Datasets

Disease datasets:  

  * Dataset 1 (Treatment Response Dataset)

    - Contains treatment response data

    - Includes responder/non-responder classifications

    - Focus on immunotherapy outcomes

  * Dataset 2 (Tumor Cell Dataset)

    - Characterizes tumor cell heterogeneity

    - Includes cell type annotations

    - Focus on tumor cell populations

  * Dataset 3 (Tumor Microenvironment Dataset)

    - Maps the tumor microenvironment

    - Contains malignant/non-malignant classifications

    - Includes detailed immune cell typing

Control datasets:  

  * Dataset 4 (Control Melanocyte Dataset)

    - Normal melanocyte control data

    - Reference for non-malignant state

    - Baseline expression profiles

  * Dataset 5 (Control Immune Dataset)

    - Normal immune cell control data

    - Reference for immune cell states

    - Baseline immune profiles

## Processes/Analyses

The following processing or analytical steps were conducted:  

  0. Environment initialization

  1. scRNAseq data consolidation and quality control  

  2. Cell clustering based on UMAP dimensionality reduction  

  3. Biomarker identification based on UMAP cluster groups  

  4. Differential gene expression analysis  

  5. Gene set enrichment analysis  

  6. Cell classification using Random Forests  

  7. Output saving (if prompted)  

# 0. Environment initialization

Install required packages.

```{r install_packages, results="hide", message=F, warning=F, error=F}

# Initiate Renv (for reproducibility)

renv::init()

# Define packages to install

pkg.list = c('renv', 'svDialogs', 'vroom', 'dplyr', 'DT', 'openxlsx', 'progress',

             'ggplot2', 'ggrepel', 'ggvenn', 'ggfortify', 'pheatmap', 

             'patchwork', 'Seurat', 'harmony', 'MAST', 'enrichR', 'UpSetR', 

             'singleCellNet')

# Install MAST (if needed)

if ('MAST' %in% pkg.install) {

  if (!requireNamespace('BiocManager', quietly=T)) {

    install.packages('BiocManager')

  }

  BiocManager::install('MAST')

  pkg.install <- pkg.install[!('MAST' %in% pkg.install)]

}

# Install SingleCellNet (if needed)

if ('singleCellNet' %in% pkg.install) {

  install.packages("devtools")

  devtools::install_github("pcahan1/singleCellNet")

  pkg.install <- pkg.install[!('singleCellNet' %in% pkg.install)]

}

# Define packages not already installed

pkg.install <- pkg.list[!(pkg.list %in% installed.packages()[, 'Package'])]

# Install uninstalled packages

if (length(pkg.install) > 0) {

  install.packages(pkg.install)

}

```

Load installed packages.  

```{r load_packages, results="hide", message=F, warning=F, error=F}

# Load packages

library(renv)           # for version control + reproducibility

library(svDialogs)      # for prompting user-input

library(vroom)          # for quickly reading data

library(dplyr)          # for data processing

library(DT)             # to display datatables

library(rio)            # to load all worksheets in a workbook

library(openxlsx)       # to write data to excel

library(progress)       # to display progress bar

library(ggplot2)        # for data visualization

library(ggrepel)        # to use geom_pont_repel()

library(ggvenn)         # to visualize venn diagrams

library(ggfortify)      # to visualize PCA plots

library(pheatmap)       # to visualize heatmaps

library(patchwork)      # for combining plots

library(Seurat)         # for scRNA-seq analysis

library(harmony)        # to integration of scRNA-seq data

library(MAST)           # for scRNAseq DEG analysis

library(enrichR)        # to perform GSEA

library(UpSetR)         # to visualize multiple set comparisons

library(singleCellNet)  # for RF implementation

```

## Custom functions

`qc_filter()`: ilter data based on QC for scRNA-seq data. 

```{r qc_filter, results="hide", message=F, warning=F, error=F}

qc_filter <- function(obj, feat.t=c(200, 2500), pct.mt.t=5, var.method='vst', 

                      feat.n=2000, qc.plot=T, top.n=10, title='') {

  ############################ FUNCTION DESCRIPTION ############################

  # feat.t = lower and upper limits on unique gene counts

  # pct.mt.t = threshold of level in mitochondrial contamination

  # var.method = method for selecting highly variable genes

  # feat.n = number of variable genes to select

  # qc.plot = boolean whether to generate plots to decide downstream analyses

  # title = string to use for plot title

  ############################## BEGIN FUNCTION ################################

  # determine percentage of mitochondrial contamination

  obj[['pct.mt']] <- PercentageFeatureSet(obj, pattern='^MT-')

  # filter + nomalize + scale data

  obj <- obj %>% 

    subset(., subset=(nFeature_RNA > feat.t[1]) & (nFeature_RNA < feat.t[2]) & 

             (pct.mt < pct.mt.t)) %>% NormalizeData(.) %>% 

    FindVariableFeatures(., selection.method=var.method) %>% 

    ScaleData(.) %>% RunPCA(., features=VariableFeatures(object=.))

  # generate follow-up QC plots (if prompted)

  if (qc.plot) { 

    p1 <- VariableFeaturePlot(obj) %>% 

      LabelPoints(plot=., points=head(VariableFeatures(obj), top.n), repel=T)

    p2 <- ElbowPlot(obj)

    plot(p1 + p2 + plot_annotation(title=title))

  }

  # return output

  return(obj)

}

```

`de_analyze()`: conduct differential expression (DE) analysis based on unpaired Student's t-test. 

```{r de_analyze, results="hide", message=F, warning=F, error=F}

de_analyze <- function(m1, m2, alt='two.sided', paired=F, var.equal=F, 

                       adj.method='bonferroni', t=0.05, de.plot=F, title='') { 

  ############################ FUNCTION DESCRIPTION ############################

  # m1, m2 = expression matrices to compare

  # alt, paired, var.equal = arguments for t.test() function

  # adj.method = method for calculating adjusted p-value

  # t = threshold for significance 

  # de.plot = boolean whether to generate a volcano plot

  ############################## BEGIN FUNCTION ################################

  # make sure two matrices have same number of rows

  if (nrow(m1) != nrow(m2)) { 

    stop('Row length does not match between the provided matrices.')

  }

  # make sure gene names align between matrices

  if (!all(rownames(m1) == rownames(m2))) { 

    stop('Gene names do not align between provided matrices.')

    }

  # instantiate output variable

  results <- data.frame(gene=rownames(m1), 

                        t.stat=vector(mode='numeric', length=nrow(m1)), 

                        p.val=vector(mode='numeric', length=nrow(m1)))

  # conduct unpaired t-test with unequal variance for each gene

  pb <- progress_bar$new(

    format='  analyzing [:bar] :percent time left: :eta', total=nrow(m1))

  for (i in 1:nrow(m1)) { 

    pb$tick()

    x <- m1[i, ]; y <- m2[i, ]

    r <- t.test(x, y, alternative=alt, paired=paired, var.equal=var.equal)

    results$t.stat[i] <- r$statistic

    results$p.val[i] <- r$p.value

  }

  # determine adjusted p-values

  results$q.val <- p.adjust(results$p.val, method=adj.method)

  # add additional fields

  results <- results %>%

    mutate(Significance=case_when(q.val < t & t.stat > 0 ~ 'Up',

                                  q.val < t & t.stat < 0 ~ 'Down',

                                  T ~ 'NS')) %>% arrange(q.val)

  # generate volcano plot (if prompted)

  if (de.plot) { 

    p <- results %>% arrange(t.stat) %>% ggplot(data=., 

           aes(x=t.stat, y=-log10(q.val), col=Significance, label=gene)) +

      geom_point() + geom_text_repel() + theme_minimal() + 

      scale_color_manual(values=c('blue', 'black', 'red')) + ggtitle(title)

    plot(p)

  }

  # return output

  return(results)

}

```

## Additional settings.  

```{r settings}

# Adjust system settings

Sys.setenv("VROOM_CONNECTION_SIZE" = 131072 * 2)

# Save plots? (default: F)

save.plots <- dlgInput('Save all plots? (T/F)', F)$res

# Set seed (for reproducibility)

set.seed(123)

```

# 1. scRNAseq data consolidation and quality control

Load all scRNA-seq files or load already saved workspace.  

```{r load_data, warning=F, message=F}

f <- list.files()

if (any(endsWith(f, '.RData'))) {

  load(f[endsWith(f, '.RData')][1])

} else { 

  # Define variables to store data

  obj <- list(); meta <- list()

  # Dataset 1: Treatment Response Dataset

  meta$dataset1 <- read.delim(file='data/dataset1/metadata.txt')

  x <- vroom('data/dataset1/counts.txt', col_names=F)

  x <- x[, 1:nrow(meta$dataset1)] %>% as.matrix(.)

  y <- read.delim('data/dataset1/genes.txt', header=F)

  rownames(x) <- make.unique(y$V1)

  colnames(x) <- meta$dataset1$ID

  obj$dataset1 <- CreateSeuratObject(x, project='dataset1',

                                      min.cells=3, min.features=200)

  obj$dataset1$Treatment <- as.factor(meta$dataset1$Treatment)

  obj$dataset1$Response <- as.factor(meta$dataset1$Response)

  rm(x, y); gc()

  # Dataset 2: Tumor Cell Dataset

  meta$dataset2 <- read.csv('data/dataset2/metadata.csv')

  x <- vroom('data/dataset2/expression_matrix.csv')

  y <- as.matrix(x[, -1]); rownames(y) <- x$...1

  obj$dataset2 <- CreateSeuratObject(y, project='dataset2',

                                      min.cells=3, min.features=200)

  obj$dataset2$type <- as.factor(meta$dataset2$cell.types)

  # Dataset 3: Tumor Microenvironment Dataset

  x <- vroom('data/dataset3/expression_matrix.txt')

  y <- as.matrix(x[-c(1:3), -1]); rownames(y) <- x$Cell[-c(1:3)]

  obj$dataset3 <- CreateSeuratObject(y, project='dataset3',

                                     min.cells=3, min.features=200)

  meta$dataset3 <- data.frame(Tumor=t(x[1, -1]),

                              Malignant=case_when(x[2, -1] == 1 ~ 'No',

                                                  x[2, -1] == 2 ~ 'Yes',

                                                  x[2, -1] == 0 ~ 'Unresolved'),

                              NMType=case_when(x[3, -1] == 1 ~ 'T',

                                               x[3, -1] == 2 ~ 'B',

                                               x[3, -1] == 3 ~ 'Macro',

                                               x[3, -1] == 4 ~ 'Endo',

                                               x[3, -1] == 5 ~ 'CAF',

                                               x[3, -1] == 6 ~ 'NK',

                                               x[3, -1] == 0 ~ 'NA'))

  obj$dataset3$Malignant <- meta$dataset3$Malignant

  obj$dataset3$NMType <- meta$dataset3$NMType

  # Dataset 4: Control Melanocyte Dataset

  x <- vroom('data/dataset4/expression_matrix.csv')

  y <- as.matrix(x[, -1]); rownames(y) <- x$...1

  obj$dataset4 <- CreateSeuratObject(y, project='dataset4',

                                      min.cells=3, min.features=200)

  # Dataset 5: Control Immune Dataset

  meta$dataset5 <- read.table('data/dataset5/metadata.tsv',

                               sep='\t', header=T)

  x <- vroom('data/dataset5/matrix.tsv')

  y <- as.matrix(x[, -1]); rownames(y) <- x$Gene.Name

  obj$dataset5 <- CreateSeuratObject(y, project='dataset5',

                                      min.cells=3, min.features=200)

  obj$dataset5$sample <- as.factor(meta$dataset5$sample)

  # Clear unnecessary memory

  rm(f, x, y, meta); gc()

}

```

Combine all datasets together into single Seurat object. Also apply `harmony` to remove clustering bias based on dataset source.  

```{r combine_data, warning=F, message=F} 

# Combine all datasets

if (!exists('data.all')) { 

  data.all <- list()

}

if (!('raw' %in% names(data.all))) { 

  data.all$raw <- merge(obj$dataset1, obj$dataset2) %>%

    merge(., obj$dataset3) %>% merge(., obj$dataset4) %>% merge(., obj$dataset5)

}

# Add source information

if (!('Source' %in% names(data.all$raw@meta.data))) { 

  data.all$raw$Source <- c(rep(names(obj)[1], ncol(obj$dataset1)),

                           rep(names(obj)[2], ncol(obj$dataset2)),

                           rep(names(obj)[3], ncol(obj$dataset3)),

                           rep(names(obj)[4], ncol(obj$dataset4)),

                           rep(names(obj)[5], ncol(obj$dataset5)))

}

# Apply QC

if (!('proc' %in% names(data.all))) { 

  data.all$proc <- qc_filter(data.all$raw)

}

# Visualize integration results

p1 <- DimPlot(object=data.all$proc, reduction='pca', pt.size=0.1, 

              group.by='Source')

p2 <- VlnPlot(object=data.all$proc, features='PC_1', pt.size=0.1, 

              group.by='Source')

p1 + p2

if (save.plots) { 

  ggsave('plots/QC_no_harmony.pdf', width=10, height=5, units='in')

}

# Apply harmony (to remove clustering based on dataset source)

if (!('harmony' %in% names(data.all$proc@reductions))) { 

  data.all$proc <- data.all$proc %>% RunHarmony('Source', plot_convergence=T)

}

# Visualize integration results (after harmony)

p1 <- DimPlot(object=data.all$proc, reduction='harmony', pt.size=0.1, 

              group.by='Source')

p2 <- VlnPlot(object=data.all$proc, features='harmony_1', pt.size=0.1, 

              group.by='Source')

p1 + p2

if (save.plots) { 

  ggsave('plots/QC_w_harmony.pdf', width=10, height=5, units='in')

}

```

# 2. Cell clustering based on UMAP dimensionality reduction

Cluster cells based on uniform manifold approximation and projection (UMAP).  

```{r UMAP_clustering, warning=F, message=T}

# Visualize UMAP plots (runtime: ~5 minutes)

if (!('umap' %in% names(data.all$proc@reductions))) { 

  data.all$proc <- data.all$proc %>% 

      RunUMAP(reduction='harmony', dims=1:20) %>% 

      FindNeighbors(reduction='harmony', dims=1:20) %>% 

      FindClusters(resolution=0.5) %>% identity()

}

#   by dataset source

p <- DimPlot(data.all$proc, reduction='umap', group.by='Source', pt.size=0.1, 

        split.by='Source') + ggtitle('UMAP split by dataset source'); plot(p)

#   by cluster (unlabeled)

p <- DimPlot(data.all$proc, reduction='umap', label=T, pt.size=0.1) + 

  ggtitle('UMAP of combined scRNA-seq data (unlabeled)'); plot(p)

if (save.plots) { 

  ggsave('plots/UMAP_unlabeled.pdf', width=10, height=5, units='in')

}

#   by cluster (labeled)

lab <- c('Tumor cell', 'Technical error 1', 'Cytotoxic CD8 T cell 1', 'B cell', 

         'Melanocytes 1', 'Macrophage', 'Cytotoxic CD8 T cell 2', 'Technical error 2', 

         'Memory T cell', 'Melanocytes 2', 'Dysfunctional CD8 T cell', 

         'Cytotoxic CD8 T cell 3', '? 1', 'Melanocytes 3', 'Activated cell', 

         'Exhausted T cell', 'Cytotoxic T cell', '? 2', '? 3', 'Melanocytes 4', '? 4')

names(lab) <- levels(data.all$proc)

plot.data <- data.all$proc %>% RenameIdents(., lab)

p <- DimPlot(plot.data, reduction='umap', label=T, pt.size=0.1) + 

  ggtitle('UMAP of combined scRNA-seq data (labeled)'); plot(p)

if (save.plots) { 

  ggsave('plots/UMAP_labeled.pdf', width=10, height=5, units='in')

}

```

# 3. Biomarker identification based on UMAP cluster groups

Determine biomarkers based on UMAP cluster assignment. 

```{r biomarkers, warning=F, message=T}

# Find all biomarkers based on clustering (runtime: ~30 minutes)

if (!('bm' %in% names(data.all))) { 

  data.all$bm <- FindAllMarkers(data.all$proc, min.pct=0.25, logfc.threshold=0.25)

}

# View table of top 3 biomarkers for each cluster

d <- data.all$bm %>% group_by(cluster) %>% slice_max(n=3, order_by=avg_log2FC)

datatable(d)

# Visualize ridge plots based on biomarkers of interest

ridge.feat <- dlg_list(sort(unique(d$gene)), multiple=T)$res

p <- RidgePlot(data.all$proc, features=ridge.feat, 

          ncol=ceiling(length(ridge.feat) / 2)); plot(p)

if (save.plots) { 

  ggsave('plots/biomarker_ridge_plots.pdf', width=10, height=10, units='in')

}

```

# 4. Differential gene expression analysis

Determine differentially expressed genes (DEGs) for the following group comparisons:  

  * Disease vs. control (tumor)  

  * Disease vs. control (immune cells, bulk)  

  * Disease vs. control (immune cells, cluster-based)  

  * Responder vs. non-responder (specific to GSE120575)  

Of note, DEGs were determined using a custom function for unpaired Student's t-test and [MAST](https://bioconductor.org/packages/release/bioc/html/MAST.html). 

```{r DE_initialization, warning=F, message=T, eval=!('data' %in% de)}

# Instantiate DE variable + relevant fields

de <- list(); de$ttest <- list(); de$mast <- list()

# Define data to work with

de$data <- GetAssayData(object=data.all$proc, slot='data')

```

## Case: disease vs. control (tumor)

Comparison between benign vs. tumor skin cells. 

Estimated runtime: ~15 minutes. 

```{r DE_tumor, warning=F, message=T, eval=!('tumor' %in% de$ttest)}

# Start timer

t.start <- Sys.time()

# Define cell groups to compare

ix <- intersect(rownames(obj$dataset3), rownames(obj$dataset2)) %>%

  intersect(., rownames(obj$dataset4))

jx1 <- data.all$proc$Source == names(obj)[3] | data.all$proc$Malignant %in% 'Yes'

jx2 <- data.all$proc$Source == names(obj)[4]

# DE analysis (t-test)

x <- de$data[ix, jx1]; y <- de$data[ix, jx2]

de$ttest$tumor <- de_analyze(x, y) %>% na.omit

# DE analysis (MAST)

x <- cbind(de$data[ix, jx1], de$data[ix, jx2]) %>% CreateSeuratObject(.)

Idents(object=x, cells=1:sum(jx1)) <- 'Disease'

Idents(object=x, cells=sum(jx1)+1:ncol(x)) <- 'Control'

de$mast$tumor <- FindMarkers(object=x, ident.1='Disease', ident.2='Control', 

                             test.use='MAST')

# End timer + log time elapsed

t.end <- Sys.time()

t.end - t.start

```

## Case: disease vs. control (immune cells, bulk)

Bulk comparison between healthy vs. diseased immune cells. 

Estimated runtime: ~20 minutes. 

```{r DE_immune_bulk,warning=F, message=T, eval=!('immune' %in% de$ttest)}

# Start timer

t.start <- Sys.time()

# Define cell groups to compare

ix <- intersect(rownames(obj$dataset1), rownames(obj$dataset2)) %>% 

  intersect(., rownames(obj$dataset5))

jx1 <- data.all$proc$Source == names(obj)[1] | data.all$proc$Malignant %in% 'No'

jx2 <- data.all$proc$Source == names(obj)[5]

# DE analysis (t-test)

x <- de$data[ix, jx1]; y <- de$data[ix, jx2]

de$ttest$immune <- de_analyze(x, y) %>% na.omit()

# DE analysis (MAST)

x <- cbind(de$data[ix, jx1], de$data[ix, jx2]) %>% CreateSeuratObject(.)

Idents(object=x, cells=1:sum(jx1)) <- 'Tumor'

Idents(object=x, cells=sum(jx1)+1:ncol(x)) <- 'Benign'

de$mast$immune <- FindMarkers(object=x, ident.1='Tumor', ident.2='Benign', 

                              test.use='MAST')

# End timer + log time elapsed

t.end <- Sys.time()

t.end - t.start

```

## Case: disease vs. control (immune cells, cluster-based)

Immune cell-specific comparisons between healthy vs. diseased cells. 

Estimated runtime: ~15 minutes. 

```{r DE_immune_cluster, warning=F, message=T, eval=!('B' %in% de$ttest)}

# Start timer

t.start <- Sys.time()

# Define labeled scRNAseq data

lab <- c('Tumor cell', 'Technical error 1', 'Cytotoxic CD8 T cell 1', 'B cell', 

         'Melanocytes 1', 'Macrophage', 'Cytotoxic CD8 T cell 2', 'Technical error 2', 

         'Memory T cell', 'Melanocytes 2', 'Dysfunctional CD8 T cell', 

         'Cytotoxic CD8 T cell 3', '? 1', 'Melanocytes 3', 'Activated cell', 

         'Exhausted T cell', 'Cytotoxic T cell', '? 2', '? 3', 'Melanocytes 4', '? 4')

names(lab) <- levels(data.all$proc)

plot.data <- data.all$proc %>% RenameIdents(., lab)

# Iterate through each immune cell group (runtime: ~5 minutes)

immune.cells <- c('B', 'CD8', 'Macrophage', 'Memory', 'Activated')

for (i in 1:length(immune.cells)) { 

  # Define row + column indices

  ix <- intersect(rownames(obj$dataset1), rownames(obj$dataset2)) %>% 

    intersect(., rownames(obj$dataset5))

  c <- grepl(immune.cells[i], Idents(plot.data))

  jx1 <- (plot.data$Source == names(obj)[1] | 

            data.all$proc$Malignant %in% 'No') & c

  jx2 <- data.all$proc$Source == names(obj)[5] & c

  # # Print number of genes + cells

  # cat(sprintf('No. of genes: %d\tNo. of G1: %d\tNo. of G2: %d\n', 

  #             length(ix), sum(jx1), sum(jx2)))

  # DE analysis (t-test)

  x <- de$data[ix, jx1]; y <- de$data[ix, jx2]

  de$ttest[[immune.cells[i]]] <- de_analyze(x, y) %>% na.omit()

  # DE analysis (MAST)

  x <- cbind(de$data[ix, jx1], de$data[ix, jx2]) %>% CreateSeuratObject(.)

  Idents(object=x, cells=1:sum(jx1)) <- 'Tumor'

  Idents(object=x, cells=sum(jx1)+1:ncol(x)) <- 'Benign'

  de$mast[[immune.cells[i]]] <- FindMarkers(object=x, ident.1='Tumor',

                                            ident.2='Benign', test.use='MAST')

}

# End timer + log time elapsed

t.end <- Sys.time()

t.end - t.start

```

## Case: responder vs. non-responder (GSE120575)

Comparison between responder vs. non-responder cells to immunotherapy. 

Estimated runtime: ~1 hour.  

```{r DE_response, warning=F, message=T, eval=!('response' %in% de$ttest)}

# Start timer

t.start <- Sys.time()

# Define cell groups to compare

ix <- rownames(obj$dataset1)

jx1 <- data.all$proc$Response %in% 'Responder'

jx2 <- data.all$proc$Response %in% 'Non-responder'

# DE analysis (t-test)

x <- de$data[ix, jx1]; y <- de$data[ix, jx2]

de$ttest$response <- de_analyze(x, y) %>% na.omit

# DE analysis (MAST)

x <- cbind(de$data[ix, jx1], de$data[ix, jx2]) %>% CreateSeuratObject(.)

Idents(object=x, cells=1:sum(jx1)) <- 'Responder'

Idents(object=x, cells=sum(jx1)+1:ncol(x)) <- 'Non-responder'

de$mast$response <- FindMarkers(object=x, ident.1='Responder', 

                                ident.2='Non-responder', test.use='MAST')

# End timer + log time elapsed

t.end <- Sys.time()

t.end - t.start

```

## Visualizations

**Description**: visual inspection of DE analysis results. 

### T-test vs. MAST DEGs

Visual comparison of DEGs determined with t-test vs. MAST (venn diagrams).

```{r DE_venn, warning=F, message=T}

# Compare DEGs b/t the two methods

f <- names(de$mast); p <- list()

for (i in 1:length(f)) { 

  x <- de$ttest[[f[i]]]; y <- de$mast[[f[i]]]

  g <- intersect(x$gene[x$q.val < 0.05], row.names(y)[y$p_val_adj < 0.05])

  # g <- intersect(x$gene[x$q.val < 0.05], y$gene[y$p_val_adj < 0.05]) 

  cat(sprintf('No. of intersecting DEGs for %s: %d\n', f[i], length(g)))

  ggvenn(list('t-test'=x$gene[x$q.val < 0.05], 

              'MAST'=row.names(y)[y$p_val_adj < 0.05])) + ggtitle(f[i])

  ggsave(paste0('plots/DE_venn_', f[i], '.pdf'), width=5, height=5, units='in')

}

```

### Responder vs. non-responder DEGs

Clustergram (or heatmap) of DEGs identified for responder vs. non-responder comparison. 

```{r DE_response_heat, warning=F, message=T}

# Define data to plot

n <- 50

g <- de$mast$response$gene[1:n]

jx1 <- which(data.all$proc$Response %in% 'Responder')[1:(n/2)]

jx2 <- which(data.all$proc$Response %in% 'Non-responder')[1:(n/2)]

x <- de$data[g, c(jx1, jx2)] %>% as.data.frame(.)

# Specify heatmap arguments

#   color palette

col.pal <- colorRampPalette(colors=c('white', 'red'), space='Lab')(100)

#   column annotations (response)

col.annot <- data.frame(Response=rep(c('Responder', 'Non-responder'), each=n/2), 

                        row.names=colnames(x))

# Heatmap

p <- pheatmap(x, cluster_rows=T, cluster_cols=F, scale='none', color=col.pal, 

              annotation_col=col.annot, cellheight=(500/n), cellwidth=(500/n), 

              show_rownames=T, show_colnames=F, gaps_col=n/2, 

              main='DEG heatmap (response)')

# Save plot

tiff('plots/DE_heatmap_response.tiff', units='in', width=11, height=8, res=300)

print({p}); dev.off()

```

PCA of DEGs identified for responder vs. non-responder comparison. 

```{r DE_response_PCA, warning=F, message=T}

# Define data to plot

g <- de$mast$response$gene[de$mast$response$p_val_adj < 0.05]

jx <- data.all$proc$Response %in% c('Responder', 'Non-responder')

x <- de$data[g, jx] %>% as.data.frame(.) %>% t(.)

m <- data.frame(Response=data.all$proc$Response[jx])

# Apply PCA

response_pca <- prcomp(x, scale.=T)

# Visualize PCA plot

autoplot(response_pca, data=m, colour='Response') + ggtitle('DEG PCA (response)')

# Save plot

ggsave('plots/DE_PCA_response.tiff', units='in', width=7, height=5, dpi=300)

```

### Tumor vs. benign DEGs

Clustergram (or heatmap) of DEGs identified for tumor vs. benign comparison. 

```{r DE_tumor_heat, warning=F, message=T}

# Define data to plot

n <- 50

g <- de$mast$tumor$gene[1:n]

jx1 <- which(data.all$proc$Source == names(obj)[3] | 

  data.all$proc$Malignant %in% 'Yes')[1:(n/2)]

jx2 <- which(data.all$proc$Source == names(obj)[4])[1:(n/2)]

x <- de$data[g, c(jx1, jx2)] %>% as.data.frame(.)

# Specify heatmap arguments

#   color palette

col.pal <- colorRampPalette(colors=c('white', 'red'), space='Lab')(100)

#   column annotations (response)

col.annot <- data.frame(Malignant=rep(c('Yes', 'No'), each=n/2), 

                        row.names=colnames(x))

# Heatmap

p <- pheatmap(x, cluster_rows=T, cluster_cols=F, scale='none', color=col.pal, 

              annotation_col=col.annot, cellheight=(500/n), cellwidth=(500/n), 

              show_rownames=T, show_colnames=F, gaps_col=n/2, 

              main='DEG heatmap (tumor)')

# Save plot

tiff('plots/DE_heatmap_tumor.tiff', units='in', width=11, height=8, res=300)

print({p}); dev.off()

```

PCA of DEGs identified for tumor vs. benign comparison. 

```{r DE_tumor_PCA, warning=F, message=T}

# Define data to plot

g <- de$mast$tumor$gene[1:100]

jx1 <- data.all$proc$Source == names(obj)[3] | data.all$proc$Malignant %in% 'Yes'

jx2 <- data.all$proc$Source == names(obj)[4]

x <- cbind(de$data[g, jx1], de$data[g, jx2]) %>% as.data.frame(.) %>% t(.)

m <- data.frame(Malignant=rep(c('Yes', 'No'), c(sum(jx1), sum(jx2))))

# Apply PCA

pca.tumor <- prcomp(x, scale.=T)

# Visualize PCA plot

autoplot(pca.tumor, data=m, colour='Malignant') + ggtitle('DEG PCA (tumor)')

# Save plot

ggsave('plots/DE_PCA_tumor.tiff', units='in', width=7, height=5, dpi=300)

```

# 5. Gene set enrichment analysis

Determine enriched pathways based on the DEG results (from MAST) using [enrichR](https://cran.r-project.org/web/packages/enrichR/vignettes/enrichR.html)..

Total runtime: ~2 minutes.  

```{r GSEA_analysis, warning=F, message=T}

# Start timer

t.start <- Sys.time()

# Instantiate GSEA variable

gsea <- list()

# Ensure comparison to human genes

setEnrichrSite('Enrichr')

# Select databases to inquire

dbs <- dlg_list(sort(listEnrichrDbs()$libraryName), multiple=T)$res

# Loop through all DE fields

f <- names(de$mast)

for (i in 1:length(f)) {

  x <- de$mast[[f[i]]]

  y <- enrichr(x$gene[x$p_val_adj < 0.05], dbs)

  z <- y$GO_Biological_Process_2021

  n <- sum(z$Adjusted.P.value < 0.05)

  gsea[[f[i]]] <- z[z$Adjusted.P.value < 0.05, ]

  cat(sprintf('No. of enriched pathways for %s: %d\n', f[i], n))

  plotEnrich(z, orderBy='Adjusted.P.value') + 

    ggtitle(sprintf('Case: %s (Total = %d)', f[i], n))

  ggsave(paste0('plots/GSEA_', f[i], '.pdf'), width=7, height=5, units='in')

}

# End timer + log time elapsed

t.end <- Sys.time()

t.end - t.start

```

Visualize intersecting pathways between comparisons.  

```{r GSEA_intersect, warning=F, message=T}

# Define plots data

df <- list('tumor'=gsea$tumor$Term, 

           'immune'=gsea$immune$Term, 

           'B'=gsea$B$Term, 

           'CD8'=gsea$CD8$Term, 

           'Macrophage'=gsea$Macrophage$Term, 

           'Memory'=gsea$Memory$Term, 

           'Activated T'=gsea$Activated$Term, 

           'response'=gsea$response$Term)

# Generate and save UpSet plot

pdf(file='plots/GSEA_intersect.pdf', width=11, height=8)

upset(fromList(df), sets=names(df), order.by='freq')

```

# 6. Cell classification using Random Forests

Apply Random Forests to construct a cell type classifier using [SingleCellNet](https://github.com/pcahan1/singleCellNet).

Define training + testing data. 

```{r ML_train_test, warning=F, message=T}

# Pull required data

s <- extractSeurat(data.all$proc, exp_slot_name='counts')

rm(data.all); gc()

st <- s$sampTab %>% mutate(ID=rownames(.))

exp <- s$expDat %>% as(., 'dgCMatrix')

# Define label field

label <- dlg_list(title='Choose label: ', sort(colnames(st)), multiple=F)$res

ID <- dlg_list(title='Choose cell ID: ', sort(colnames(st)), multiple=F)$res

# Re-order by label field

ix <- order(st[[label]])

st <- st[ix, ]

exp <- exp[, ix]

# Train/Test split

stList1 <- splitCommon(sampTab=st, ncells=100, dLevel=label)

stTrain <- stList1[[1]]

expTrain <- exp[, rownames(stTrain)]

stList2 <- splitCommon(sampTab=stList1[[2]], ncells=100, dLevel=label)

stTest <- stList2[[1]]

expTest <- exp[, rownames(stTest)]

```

Apply SCN (runtime: ~5 minutes). 

```{r SCN, warning=F, message=T}

# Start timer

t.start <- Sys.time()

# Train model

x <- scn_train(stTrain=stTrain, expTrain=expTrain, nTopGenes=10, nRand=70, 

               nTrees=1000, nTopGenePairs=25, dLevel=label)

# Test model

y <- scn_predict(cnProc=x[['cnProc']], expDat=expTest, nrand=50)

# End timer + log time elapsed

t.end <- Sys.time()

t.end - t.start

```

Model assessment and visualization of results.  

```{r ML_evaluation, warning=F, message=F}

# Assess model

z <- assess_comm(ct_scores=y, stTrain=stTrain, stQuery=stTest, dLevelSID=ID,

                 classTrain=label, classQuery=label)

plot_PRs(z) + ggtitle('Model performance (test set)')

if (save.plots) {

  ggsave('plots/ML_performance.pdf', width=5, height=5, units='in')

}

# Visualize results

#   classification results

nrand <- 50

sla <- as.vector(stTest[[label]])

names(sla) <- as.vector(rownames(stTest))

slaRand <- rep('rand', nrand)

names(slaRand) <- paste('rand_', 1:nrand, sep='')

sla <- append(sla, slaRand)

p <- sc_hmClass(classMat=y, grps=sla, max=300, isBig=T)

if (save.plots) {

  pdf('plots/ML_classification.pdf', width=7, height=5,

      title='Classification results (test set)')

  p

}

#   attribution plot

plot_attr(classRes=y, sampTab=stTest, nrand=nrand, sid=ID, dLevel=label) + 

  xlab('Predicted group') + ylab('True class ratio') + 

  ggtitle('Attribution plot (test set)')

if (save.plots) {

  ggsave('plots/ML_attribution.pdf', width=7, height=5, units='in')

}

```

# 7. Output saving (if prompted)

```{r save_outputs, warning=F, message=T}

# Save outputs

save.data <- dlgInput('Save all results? (T/F)', F)$res %>% as.logical(.)

if (save.data) { 

  # Biomarker results

  xl.list <- list('biomarkers'=data.all$bm)

  write.xlsx(xl.list, 'results/biomarkers.xlsx', row.names=F)

  # DEG results

  write.xlsx(de$ttest, 'results/DE_ttest.xlsx', row.names=F)

  write.xlsx(de$mast, 'results/DE_MAST.xlsx', row.names=T)

  # GSEA results

  write.xlsx(gsea, 'results/GSEA.xlsx', row.names=F)

}

# Save workspace

save.wksp <- dlgInput('Save R workspace? (T/F)', F)$res %>% as.logical(.)

if (save.wksp) { 

  # Save workspace (runtime: ~2 hours)

  save.image(file='scRNAseq_wksp.RData')

}

# Save project state (for reproducibility)

renv::snapshot()

```

# Visualize data

```{r visualize_data, warning=F, message=F, fig.width=12, fig.height=6}

# Create UMAP plots

p1 <- DimPlot(data.all$proc, reduction='umap', group.by='Source',

              label=T, repel=T) +

  ggtitle('Dataset Source')

p2 <- DimPlot(data.all$proc, reduction='harmony', group.by='Source',

              label=T, repel=T) +

  ggtitle('Dataset Source (Harmony)')

# Combine plots

p1 + p2

# Save plot

ggsave('plots/UMAP_integration.pdf', width=12, height=6)

```

# Visualize tumor vs melanocyte comparison

```{r visualize_tumor, warning=F, message=F, fig.width=12, fig.height=6}

# Create volcano plot

p1 <- EnhancedVolcano(tumor.vs.melanocyte,

                      lab=rownames(tumor.vs.melanocyte),

                      x='avg_log2FC', y='p_val_adj',

                      pCutoff=0.05, FCcutoff=1,

                      title='Tumor vs Melanocyte DEGs')

# Create heatmap

top.genes <- rownames(tumor.vs.melanocyte)[tumor.vs.melanocyte$p_val_adj < 0.05]

top.genes <- top.genes[order(abs(tumor.vs.melanocyte$avg_log2FC[

  tumor.vs.melanocyte$p_val_adj < 0.05]), decreasing=T)][1:50]

plot.data <- subset(data.all$proc,

                    subset=Source %in% c(names(obj)[3:4]))

plot.data <- ScaleData(plot.data, features=top.genes)

p2 <- DoHeatmap(plot.data, features=top.genes,

                group.by='Source') +

  ggtitle('Top 50 Tumor vs Melanocyte DEGs')

# Combine plots

p1 + p2

# Save plot

ggsave('plots/tumor_vs_melanocyte.pdf', width=12, height=6)

```

# Visualize immune vs non-malignant comparison

```{r visualize_immune, warning=F, message=F, fig.width=12, fig.height=6}

# Create volcano plot

p1 <- EnhancedVolcano(immune.vs.nonmalignant,

                      lab=rownames(immune.vs.nonmalignant),

                      x='avg_log2FC', y='p_val_adj',

                      pCutoff=0.05, FCcutoff=1,

                      title='Immune vs Non-malignant DEGs')

# Create heatmap

top.genes <- rownames(immune.vs.nonmalignant)[

  immune.vs.nonmalignant$p_val_adj < 0.05]

top.genes <- top.genes[order(abs(immune.vs.nonmalignant$avg_log2FC[

  immune.vs.nonmalignant$p_val_adj < 0.05]), decreasing=T)][1:50]

plot.data <- subset(data.all$proc,

                    subset=Source %in% c(names(obj)[c(1,2,5)]))

plot.data <- ScaleData(plot.data, features=top.genes)

p2 <- DoHeatmap(plot.data, features=top.genes,

                group.by='Source') +

  ggtitle('Top 50 Immune vs Non-malignant DEGs')

# Combine plots

p1 + p2

# Save plot

ggsave('plots/immune_vs_nonmalignant.pdf', width=12, height=6)

```

# Visualize immune cell type comparisons

```{r visualize_immune_types, warning=F, message=F, fig.width=12, fig.height=6}

# Create plots for each cell type

for (i in 1:length(immune.cells)) {

  # Read results

  immune.cell.type <- read.csv(paste0('results/',

                                     immune.cells[i],

                                     '_cell_DEGs.csv'))

  # Create volcano plot

  p1 <- EnhancedVolcano(immune.cell.type,

                        lab=rownames(immune.cell.type),

                        x='avg_log2FC', y='p_val_adj',

                        pCutoff=0.05, FCcutoff=1,

                        title=paste0(immune.cells[i], ' Cell DEGs'))

  # Create heatmap

  top.genes <- rownames(immune.cell.type)[immune.cell.type$p_val_adj < 0.05]

  top.genes <- top.genes[order(abs(immune.cell.type$avg_log2FC[

    immune.cell.type$p_val_adj < 0.05]), decreasing=T)][1:50]

  plot.data <- subset(data.all$proc,

                      subset=Source %in% c(names(obj)[c(1:3,5)]))

  plot.data <- ScaleData(plot.data, features=top.genes)

  p2 <- DoHeatmap(plot.data, features=top.genes,

                  group.by='Source') +

    ggtitle(paste0('Top 50 ', immune.cells[i], ' Cell DEGs'))

  # Combine plots

  p <- p1 + p2

  # Save plot

  ggsave(paste0('plots/', immune.cells[i], '_cell_DEGs.pdf'),

         plot=p, width=12, height=6)

}

```

# Visualize responder vs non-responder comparison

```{r visualize_response, warning=F, message=F, fig.width=12, fig.height=6}

# Create volcano plot

p1 <- EnhancedVolcano(response.vs.nonresponse,

                      lab=rownames(response.vs.nonresponse),

                      x='avg_log2FC', y='p_val_adj',

                      pCutoff=0.05, FCcutoff=1,

                      title='Response vs Non-response DEGs')

# Create heatmap

top.genes <- rownames(response.vs.nonresponse)[

  response.vs.nonresponse$p_val_adj < 0.05]

top.genes <- top.genes[order(abs(response.vs.nonresponse$avg_log2FC[

  response.vs.nonresponse$p_val_adj < 0.05]), decreasing=T)][1:50]

plot.data <- subset(data.all$proc,

                    subset=Source %in% names(obj)[1])

plot.data <- ScaleData(plot.data, features=top.genes)

p2 <- DoHeatmap(plot.data, features=top.genes,

                group.by='Response') +

  ggtitle('Top 50 Response vs Non-response DEGs')

# Combine plots

p1 + p2

# Save plot

ggsave('plots/response_vs_nonresponse.pdf', width=12, height=6)

```