nuc_06_backgroundmarkers.R

# Copyright (c) [2024] [Ricardo O. Ramirez Flores]
# 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)
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()

# Pipeline for differential expression
de_res <- map(cts, function(ct) {
  print(ct)
  ct_meta_data <- coldat %>%
    mutate(test_column = ifelse(cell_type == ct, ct, "rest"))

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

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

  dat <- dat[keep,]

  dat <- calcNormFactors(dat)

  design <- model.matrix(~factor(test_column,
                                 levels = c("rest",ct)), dat$samples)

  colnames(design) <- c("int", ct)

  dat <- estimateDisp(dat, design)

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

  res <- glmQLFTest(fit, coef=ncol(design))

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

  return(de_res)

})

de_res <- de_res %>%
  enframe() %>%
  unnest()

de_res %>%
  dplyr::filter(logFC > 0) %>%
  arrange(name, FDR, - logFC) %>%
  write_csv(file = "data/CRC_cellmarkers.csv")