nuc_12.6_cell_states_funcomics.R

# Copyright (c) [2023] [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(tidyverse)
library(decoupleR)
library(biomaRt)
library(ComplexHeatmap)

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

#' Basic function to convert human to mouse gene names
convertMouseGeneList <- function(x){
  require("biomaRt")
  human = useMart("ensembl", dataset = "hsapiens_gene_ensembl",
                  host = "https://dec2021.archive.ensembl.org/")
  mouse = useMart("ensembl", dataset = "mmusculus_gene_ensembl",
                  host = "https://dec2021.archive.ensembl.org/")

  genesV2 = getLDS(attributes = c("mgi_symbol"),
                   filters = "mgi_symbol",
                   values = x , mart = mouse,
                   attributesL = c("hgnc_symbol"),
                   martL = human, uniqueRows=T)

  humanx <- unique(genesV2[, 2])

  return(genesV2)
}

# Now we take the differential stats
cts <- c("CM", "Endo", "Fib")

for(celltype in cts){
    all_de_res <- read_csv(file = paste0("results/cell_states/", 
                                         celltype, 
                                         "_condition-simplified_cellstates_clustermarkers.csv"))
    all_de_res$cell_type <- celltype

    all_de_res <- all_de_res %>%
      dplyr::mutate(sign_stat = sign(logFC) * F) %>%
      dplyr::select(cell_type, name, gene, sign_stat) %>%
      pivot_wider(names_from = name,
                  values_from = sign_stat,
                  values_fill = 0) %>%
      group_by(cell_type) %>%
      nest() %>%
      mutate(data = map(data, ~ .x %>%
                          unique() %>%
                          column_to_rownames("gene") %>%
                          as.matrix())) %>%
      deframe()

    all_de_res <- map(all_de_res, function(pb_data) {

      mice_genes <- convertMouseGeneList(x = rownames(pb_data))

      # Duplicated mice and human genes, are deleted
      repeated_m <- mice_genes$MGI.symbol[which(duplicated(mice_genes$MGI.symbol))] %>% unique()
      repeated_h <- mice_genes$MGI.symbol[which(duplicated(mice_genes$HGNC.symbol))] %>% unique()
      out_genes <- unique(repeated_m, repeated_h)

      # useful genes
      useful_genes <- rownames(pb_data)[!(rownames(pb_data) %in% out_genes)]

      mice_genes <- mice_genes %>%
        dplyr::filter(MGI.symbol %in% useful_genes) %>%
        dplyr::filter(MGI.symbol %in% rownames(pb_data)) %>%
        column_to_rownames("MGI.symbol")

      mice_genes <- na.omit(mice_genes[useful_genes, , drop = F])

      mice_genes <- mice_genes[rownames(mice_genes) %in% rownames(pb_data), , drop =F]
    
      # Finally filter and rename
      pb_data <- pb_data[rownames(mice_genes),]
      rownames(pb_data) <- mice_genes$HGNC.symbol

      return(pb_data)

    })

    all_de_res <- map(all_de_res, function(x)

      x[!duplicated(rownames(x)),]

      )

    # 2 Get functions for enrichment -----------------------------

    msigdb <- decoupleR::get_resource("MSigDB")

    msigdb_hallmarks <- msigdb %>%
      dplyr::filter(collection == "hallmark") %>%
      dplyr::select(genesymbol, geneset) %>%
      dplyr::mutate(weight = 1) %>%
      dplyr::rename("source" = geneset,
                "target" = genesymbol) %>%
      unique()

    gsets <- list("CM" = read_csv("data/gsets/pathophysiol_processed/cm.csv"),
                  "Fib" = read_csv("data/gsets/pathophysiol_processed/mural.csv"),
                  "Endo" = read_csv("data/gsets/pathophysiol_processed/vasc.csv"),
                  "progeny" = decoupleR::get_progeny(top = 500),
                  "collectri" = decoupleR::get_collectri(organism = "human"),
                  "HMARKS" = msigdb_hallmarks)
    # Rename collectri MOR to weight for compatibility
    colnames(gsets$collectri) <- c("source", "target", "weight")

    # I will define a function that performs enrichment of needed genesets
    get_enriched_terms <- function(ct, gset) {

      ct_load_list <- all_de_res[ct]
      net <- gsets[gset] %>%
        bind_rows() %>%
        unique()

      enrich_ct <- map(ct_load_list, function(ct_load) {

        enrich_res <- decoupleR::run_ulm(mat = ct_load,
                                         network = net,
                                         .source = "source",
                                         .target = "target",
                                         .mor = "weight",
                                         minsize = 4)


        enrich_res <- enrich_res %>%
          dplyr::select(-statistic) %>%
          group_by(condition) %>%
          dplyr::mutate(adj_p = p.adjust(p_value)) %>%
          dplyr::arrange(adj_p, desc(abs(score)))


      })

      return(enrich_ct)

    }

    # Then we define a plotting function that takes this list
    plot_sign <- function(enrich_ct, pval_thrs = 0.05, top_n = 10){

      sign_targets <- map(enrich_ct, function(ct_load) {

        ct_load %>%
          dplyr::filter(adj_p <= pval_thrs) %>%
          dplyr::filter(score > 0) %>%
          dplyr::arrange(adj_p, desc(abs(score))) %>%
          group_by(condition) %>%
          dplyr::slice(1:top_n) %>%
          pull(source) %>%
          unique()

      }) %>% unlist() %>%
        unique()
    
      plt_mat <- map(enrich_ct, ~ .x %>% dplyr::filter(source %in% sign_targets) %>%
                       dplyr::select(-p_value) %>%
                       dplyr::select(-adj_p) %>%
                       pivot_wider(names_from = condition, values_from = score))

      plt_mat <- plt_mat[[1]] %>%
        column_to_rownames("source") %>%
        as.matrix()

      max_score <- max(plt_mat, na.rm = T)

      col_split <- rep(names(enrich_ct), 1, each = 2)

      col_fun_fact <- circlize::colorRamp2(seq((max_score + 0.5) *
                                                 -1, max_score + 0.5, length = 50),
                                           grDevices::hcl.colors(50,"Blue-Red 2",
                                                             rev = T))

      hmap_plot <- Heatmap(plt_mat, name = "estimate",
                           cluster_columns = F,
                           row_names_gp = gpar(fontsize = 7),
                           column_names_gp = gpar(fontsize = 7),
                           column_title_gp = gpar(fontsize = 7),
                           show_row_dend = F,
                           col = col_fun_fact)

      return(hmap_plot)

    }

    # Make final function

    enrich_factors <- function(ct, gset, pval_thrs = 0.05, top_n = 15) {

      test_enrich_ct <- get_enriched_terms(ct, gset)

      plt_h <- plot_sign(enrich_ct = test_enrich_ct,
                         pval_thrs = pval_thrs,
                         top_n = top_n)

      return(list(data = test_enrich_ct, plt = plt_h))

    }

    # Generate results
    result <- enrich_factors(ct = celltype,
                           gset = c(celltype, "progeny", "HMARKS"),
                           pval_thrs = 0.05,
                           top_n = 10)
    write_csv(result$data[[1]], paste0("./results/cell_states/", 
                                       celltype, 
                                       "_condition-simplified_cellstates_funcomics.csv"))
    
    # Generate collectri results
    result_collectri <- enrich_factors(ct = celltype,
                                       gset = "collectri",
                                       pval_thrs = 0.05,
                                       top_n = 10)
    write_csv(result_collectri$data[[1]], paste0("./results/cell_states/",
                                       celltype,
                                       "_condition-simplified_cellstates_funcomics_collectri.csv"))
}