nuc_09.1_estimate_degs.R

# Copyright (c) [2023] [Ricardo O. Ramirez Flores and Ines Rivero-Garcia]
# roramirezf@uni-heidelberg.de

#' In this script we calculate markers of cells
#' using edgeR and pseudobulk profiles of all samples

library(SingleCellExperiment)
library(scater)
library(edgeR)
library(tidyverse)
library(ComplexHeatmap)

setwd("/mnt/sds-hd/sd22b002/projects/ines/heart_revremod")

# Importing pb data
pb_data <- read_csv("data/CRC_pseudobulk.csv") %>%
  column_to_rownames("...1") %>%
  as.matrix() %>%
  t()

# Col-data
coldat <- read_csv("data/CRC_pseudobulk_coldata.csv") %>%
  dplyr::select(colname, psbulk_n_cells, cell_type, ConditionSimplified, sample_id) %>%
  column_to_rownames("colname")

pb_data <- pb_data[,rownames(coldat)]

# Defining cts
cts <- coldat$cell_type %>%
  unique() %>%
  set_names()

mrkr_genes <- read_csv("data/CRC_cellmarkers.csv")
mrkr_genes <- mrkr_genes %>% #dplyr::filter(!name %in% exclude_ct) %>%
  dplyr::filter(FDR < 0.01, logFC > 1) %>%
  dplyr::select(name, gene) %>%
  dplyr::rename("cell_type" = name)

mrkr_genes %>% #dplyr::filter(!name %in% exclude_ct) %>%
  #dplyr::filter(FDR < 0.01, logFC > 1) %>%
  dplyr::filter(gene == "Nppa")


# Pipeline for differential expression
all_de_res <- map(cts, function(ct) {
  print(ct)

  # Banned genes
  b_genes <- mrkr_genes %>%
    dplyr::filter(cell_type != ct) %>%
    pull(gene)

  # Defining the levels and meta-data
  ct_meta_data <- coldat %>%
    dplyr::filter(cell_type == ct) %>%
    dplyr::mutate(ConditionSimplified = as.factor(ConditionSimplified))

  # Processing the input data
  dat <- pb_data[,rownames(ct_meta_data)]

  # Filtering genes

  g_ix <- rownames(dat) %in% b_genes

  dat <- dat[!g_ix,]

  dat <- DGEList(dat, samples = DataFrame(ct_meta_data))

  keep <- filterByExpr(dat, group = ct_meta_data$ConditionSimplified)

  dat <- dat[keep,]

  dat <- calcNormFactors(dat)

  # Defining the contrasts

  design <- model.matrix(~ 0 + ConditionSimplified, ct_meta_data)
  colnames(design) <- ct_meta_data$ConditionSimplified %>% levels()

  dat <- estimateDisp(dat, design)
  fit <- glmQLFit(dat, design, robust=TRUE)

  my_contrasts <- makeContrasts(
    DBvControl = DB - Control,
    ORABvControl = ORAB - Control,
    DBvORAB = DB - ORAB,
    levels = design
  )

  de_res <- map(set_names(colnames(my_contrasts)), function(ctrst) {

    degs <- glmQLFTest(fit, contrast=my_contrasts[,ctrst])

    degs <- degs %>%
      topTags(., n = Inf) %>%
      as.data.frame() %>%
      rownames_to_column("gene")

  }) %>%
    enframe("contrast") %>%
    unnest()

  return(de_res)

})

all_de_res <- all_de_res %>%
  enframe("cell_type") %>%
  unnest()

all_de_res %>%
  dplyr::filter(FDR < 0.05) %>%
  dplyr::filter(gene == "Nppa")

all_de_res %>%
  write_csv(file = "results/degs/edgeR_res_condition-simplified_nobannedgenes.csv")

pdf("results/degs/edgeR_pvals_histogram_condition-simplified.pdf", height = 6, width = 9)
all_de_res %>%
  ggplot(aes(PValue)) +
  geom_histogram() +
  facet_wrap(. ~ contrast)
dev.off()

# Visualise nr DEGs
all_de_res_sign <- all_de_res %>% filter(FDR < 0.05)

summary_counts <- all_de_res_sign %>%
  group_by(cell_type, contrast) %>%
  summarise(
    upregulated = sum(logFC > 0, na.rm = TRUE),
    downregulated = sum(logFC < 0, na.rm = TRUE),
    .groups = "drop"
  )

# Add minus to dowregulated genes for nice plotting
summary_counts_long <- summary_counts
summary_counts_long$downregulated <- -summary_counts_long$downregulated

# Make long
summary_counts_long <- reshape2::melt(summary_counts_long)

# Plot
pdf("results/degs/edgeR_nr_degs_barplot_condition-simplified.pdf", height = 4, width = 8)
summary_counts_long %>%
  ggplot(aes(x = cell_type, y = value, fill = variable)) +
    geom_bar(stat = "identity") +
    theme_bw(base_size = 12) +
    facet_wrap(~ contrast) +
    theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
dev.off()

# Volcano plots
all_de_res$result <- NA

for(i in 1:nrow(all_de_res)){
    if(all_de_res[i, "FDR"] < 0.05){
        if(all_de_res[i, "logFC"] > 0){
            all_de_res[i, "result"] <- "upregulated"
        }else{
            all_de_res[i, "result"] <- "downregulated"
        }
    }else{
        all_de_res[i, "result"] <- "not significant"
    }
}


pdf("results/degs/edgeR_volcanos_condition-simplified.pdf", height = 7, width = 14)
all_de_res %>%
  ggplot(aes(x = logFC, y = -log10(FDR), color = result)) +
    geom_point() +
    geom_vline(xintercept = 0, linewidth = 0.5, linetype = "dashed", color = "grey25") +
    geom_hline(yintercept = -log10(0.05), linewidth = 0.5, linetype = "dashed", color = "grey25") +
    theme_bw(base_size = 14) +
    scale_color_manual(values = c("downregulated" = "#00bfc4", 
                                  "not significant" = "grey", 
                                  "upregulated" = "#f8766d")) +
    facet_grid(contrast ~ cell_type)
dev.off()

pdf("results/degs/edgeR_nr_degs_heatmap_condition-simplified.pdf", height = 5, width = 3)
ggplot(summary_counts_long, aes(x = contrast, y = cell_type, fill = value, label = abs(value))) +
    geom_tile() +
    geom_text() +
    facet_wrap(~variable) +
    scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
    theme_bw()
dev.off()