see ../../../validation_schulteschrepping/differential_expression/

library(fgsea)

## Loading required package: Rcpp

library(tidyverse)

## ── Attaching packages ────────────────────────────────── tidyverse 1.3.0 ──

## ✔ ggplot2 3.3.0     ✔ purrr   0.3.3
## ✔ tibble  2.1.3     ✔ dplyr   0.8.5
## ✔ tidyr   1.0.2     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0

## ── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

Healthy vs COVID

DE_RES_IN_PATH <- "data/externalData/SchulteSchrepping/differential_expression_cohort1_no_cutoffs/sample_groups/t0_plus_healthy/results/healthy_vs_covid/toptables/COVID-Healthy/id.celltype/8_pDCs--model@healthy_vs_covid--coef@COVID-Healthy--toptab.tsv"

LE_IN_PATH <- "input/genesets/pdc_apoptosis_le/pDC_PC1_covid_LE.rds"

#FIG_OUT_PATH <- "plots/validations/SchulteSchrepping/pdc_apoptosis_no_cell_number_cutoff/pdc_apoptosis_cohort1_using_our_LE_no_cutoff_covid-healthy.pdf"
#dir.create(dirname(FIG_OUT_PATH), recursive = TRUE)

toptab <- read_tsv(DE_RES_IN_PATH)

## Parsed with column specification:
## cols(
##   gene = col_character(),
##   logFC = col_double(),
##   AveExpr = col_double(),
##   t = col_double(),
##   P.Value = col_double(),
##   adj.P.Val = col_double(),
##   B = col_double()
## )

genes <- readRDS(LE_IN_PATH)

stat_vec <- toptab %>% select(gene, t) %>% deframe()

enrich_res <- fgsea(pathways = list(pdc_apoptosis_le = genes),
                    stats = stat_vec, nperm = 10000)

pval <- round(enrich_res$pval, 5)

p <- plotEnrichment(pathway = genes, stats = stat_vec) +
             annotate("text",  x=Inf, y = Inf, label = paste("p =", pval), vjust=1, hjust=1) +
             ggtitle("cohort1 pDC apoptosis using LE from our data - no cell # cutoff")
#ggsave(plot = p, filename = FIG_OUT_PATH, width = 4, height = 4)
print(p)

Severity

DE_RES_IN_PATH <- "data/externalData/SchulteSchrepping/differential_expression_cohort1_no_cutoffs/sample_groups/t0_covid_only/results/severity/toptables/severe-mild/id.celltype/8_pDCs--model@severity--coef@severe-mild--toptab.tsv"

#FIG_OUT_PATH <- "plots/validations/SchulteSchrepping/pdc_apoptosis_no_cell_number_cutoff/pdc_apoptosis_cohort1_using_our_LE_no_cutoff_severe-mild.pdf"
#dir.create(dirname(FIG_OUT_PATH), recursive = TRUE)

toptab <- read_tsv(DE_RES_IN_PATH)

## Parsed with column specification:
## cols(
##   gene = col_character(),
##   logFC = col_double(),
##   AveExpr = col_double(),
##   t = col_double(),
##   P.Value = col_double(),
##   adj.P.Val = col_double(),
##   B = col_double()
## )

stat_vec <- toptab %>% select(gene, t) %>% deframe()

enrich_res <- fgsea(pathways = list(pdc_apoptosis_le = genes),
                    stats = stat_vec, nperm = 10000)

## Warning in fgsea(pathways = list(pdc_apoptosis_le = genes), stats = stat_vec, : There are ties in the preranked stats (0.02% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.

pval <- round(enrich_res$pval, 5)

p <- plotEnrichment(pathway = genes, stats = stat_vec) +
             annotate("text",  x=Inf, y = Inf, label = paste("p =", pval), vjust=1, hjust=1) +
             ggtitle("cohort1 pDC apoptosis severe vs mild enrichment using LE from our data \nno cell# cutoff")
#ggsave(plot = p, filename = FIG_OUT_PATH, height = 4, width = 4)
print(p)

GSEA of PDC apoptosis leading edge genes in schulte schrepping et al. 2020 data

Healthy vs COVID

Severity