4) Differential analysis

This script performs differential expressed genes (DEG) analysis and Gene Set Enrichment Analysis (GSEA) per cell type/cluster. Cluster free analysis is covered in script 5). More info about the package cacoa here: https://github.com/kharchenkolab/cacoa pre-print: https://www.biorxiv.org/content/10.1101/2022.03.15.484475v1.full.pdf Vignette: https://pklab.med.harvard.edu/viktor/cacoa/walkthrough_short.html

Setup

#Load helper functions
source("/people/nrq364/NeuSiC/scRNA_helper_new.R")

devtools::load_all("/people/nrq364/cacoa_dev")
#library(cacoa)
library(magrittr)
library(conos)
library(pagoda2)
library(qs)
library(ggplot2)
library(ggrastr)
library(openxlsx)
library(EnhancedVolcano)
library(cowplot)

1. Read in data

Read in data generated in notebook 3.

cao <- qread("cao.qs", nthreads = 10)

2. Differentially expressed genes (DEGs) analysis

DESeq2 uses samplewise pseudobulk gene expression.

• Larger clusters generally lead to more DEGs due to increased statistical power and higher overall counts in the pseudobulk approach.
• Smaller clusters may yield fewer DEGs because of lower counts, higher variability, and reduced statistical power.

cao$estimateDEPerCellType(n.cores=50, test = "DESeq2.Wald", verbose = T, independent.filtering = T) 

2.1 Save DEG result as table

Extract DEG results from cao object. DEG is a list. Contains all genes.

DEG <- cao$test.results$de %>% lapply("[[", "res")

Save as an excel file, with one sheet per cell type.

write.xlsx(x = DEG, file= "DEG.xlsx")

Save only DE genes with padj < 0.05

DEG_sig <- lapply(DEG, function(df) df[df$padj < 0.05, ])

Save as excel file with one sheet per cell type.

write.xlsx(x = DEG_sig, file="DEG_sig.xlsx")

2.2 Visualization of DE genes as Volcano plot

Use the package EnhancedVolcano to plot volcano plots of cell types of interest. Default cut-offs: p.adj < 0.05, log2FC > |1.5|

  EnhancedVolcano(DEG$CGE_IPC,
    lab = rownames(DEG$CGE_IPC),
    x = 'log2FoldChange',
    y = 'padj',
    title = 'CGE IPCs',
    pCutoff = 0.05,
    FCcutoff = 1.5,
    pointSize = 2.0,
    labSize = 4.0,
    colAlpha = 0.8,
    legendPosition = 'right',
    legendLabSize = 12,
    legendIconSize = 2.0,
    drawConnectors = TRUE,
    widthConnectors = 0.75, ylim = c(0,6), subtitleLabSize = 0)

Plot all cell types.

volcano_list <- lapply(names(DEG), function(cell_type) {
  EnhancedVolcano(DEG[[cell_type]],
                  lab = rownames(DEG[[cell_type]]),
                  x = 'log2FoldChange',
                  y = 'pvalue',
                  title = cell_type,
                  pCutoff = 0.05,
                  FCcutoff = 1.5,
                  pointSize = 2.0,
                  labSize = 4.0,
                  colAlpha = 0.8,
                  drawConnectors = TRUE,
                  widthConnectors = 0.75,
                  titleLabSize = 8,      # title font
                  subtitleLabSize = 0,   # subtitle font
                  captionLabSize = 4,     # caption font
                  axisLabSize = 8, 
                  ylim = c(0,6)) + theme(
    legend.position = "none")
})

library(patchwork)
wrap_plots(volcano_list, ncol = 4)

3. Gene set enrichment analysis (GSEA)

GSEA takes the whole ranked list of DE genes, so n.top.genes parameter is not relevant here. BH method used for multiple comparison correction.

cao$estimateOntology(type="GSEA", org.db=org.Mm.eg.db::org.Mm.eg.db)

3.1 Save GSEA result as table

Use a function to save terms with p.adj < 0.05 as a .tsv file and read it into R. It is one huge table with all cell types and all 3 subtypes of gene ontologies (biological process = BP, molecular function = MF, cellular compartment = CC).

cao$saveOntologyAsTable("GSEA.tsv", name="GSEA")
GSEA <- read.table("GSEA.tsv", sep="\t", header = T)

Transform GSEA table into a list of tables (one table per cell type).

cell_types <- unique(GSEA$CellType)
GSEA_list <- lapply(cell_types, {function(x) GSEA[GSEA$CellType == x, ]})
names(GSEA_list) <- cell_types

write.xlsx(x = GSEA_list, file= "GSEA.xlsx")

Explanation of result table

• Description: Gene set name
• setSize: number of genes with gene-level statistic values. To be more specific, if we input a pathway gene set with 58 genes (e.g. HALLMARK_MYC_TARGETS_V2 ), while our gene list only have 54 of it, so the result will only show setSize = 54
• enrichmentScore: also called ES, same as in Broad GSEA implementation. It reflects the degree to which a gene set is overrepresented at the top or bottom of a ranked list of genes.
• NES: normalized enrichment score is the primary statistic for examining gene set enrichment results. By normalizing the enrichment score, GSEA accounts for differences in gene set size and in correlations between gene sets and the expression dataset. Therefore, NES can be used to compare analysis results across gene sets. A positive normalized enrichment scores (NES) will indicate that genes in set S will be mostly represented at the top of your list (logFC > 0 or up-regulated genes).
• rank: The position in the ranked list at which the maximum enrichment score occurred. If gene sets achieve the maximum enrichment score near the top or bottom of the ranked list, the rank at max is either very small or very large.
• leading_edge: includes three statistics
  • tags: The percentage of gene hits before (for positive ES) or after (for negative ES) the peak in the running enrichment score, which indicates the percentage of genes contributing to the enrichment score.
  • list: The percentage of genes in the ranked gene list before (for positive ES) or after (for negative ES) the peak in the running enrichment score. This gives an indication of where in the list the enrichment score is attained.
  • signal: If the gene set is entirely within the first N positions in the list, then the signal strength is maximal or 100%. If the gene set is spread throughout the list, then the signal strength decreases towards 0%.

3.2 Visualization of GSEA results as heatmap

As default only GO terms of BP (biological process) subtype are plotted. Similar GO terms are clustered together to reduce redundancy. As a default we plot the top 20 up- and downregulated GO terms.

cao$plotOntologyHeatmap(name="GSEA", genes="down", subtype="BP", top.n = 20, cluster = F)

Loading required package: DOSE

Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
DOSE v4.0.1 Learn more at https://yulab-smu.top/contribution-knowledge-mining/

Please cite:

Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology
Semantic and Enrichment analysis. Bioinformatics. 2015, 31(4):608-609

cao$plotOntologyHeatmap(name="GSEA", genes="up", subtype="BP", top.n = 20, cluster = F)