nuc_13.2_liana_plots.R

# Copyright (c) [2023] [Ines Rivero Garcia]
# ines.rivero@uni-heidelberg.de

#' In this script we visualise the liana candidate ligands and receptors
library(tidyverse)
library(ggplot2)
library(edgeR)
library(scater)

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

# Load significant CCC results
DBvORAB <- read.csv("results/liana/liana_condition-simplified_DBvORAB_significant_results.csv", 
                    header = TRUE)
DBvControl <- read.csv("results/liana/liana_condition-simplified_DBvControl_significant_results.csv", 
                       header = TRUE)
ORABvControl <- read.csv("results/liana/liana_condition-simplified_ORABvControl_significant_results.csv", 
                         header = TRUE)


total_df <- bind_rows(
    DBvORAB %>% mutate(contrast = "DBvORAB"),
    DBvControl %>% mutate(contrast = "DBvControl"),
    ORABvControl %>% mutate(contrast = "ORABvControl")   
)

summary_df <- total_df %>%
  group_by(contrast) %>%
  summarise(
    total_interactions = n(),
    secreted_ligand = sum(ligand_secreted == "Yes", na.rm = TRUE),
    non_secreted_ligand  = sum(ligand_secreted == "No",  na.rm = TRUE),
    .groups = "drop"
  )


# Barplot of total number of deregulated CCC interactions
pdf("results/liana/liana_condition-simplified_total_interactions_bplot.pdf", 
    height = 4.5, width = 4)
ggplot(summary_df, aes(x = reorder(contrast, -total_interactions), y = total_interactions)) +
  geom_col(color = "black", fill = "grey75") +
  labs(x = "Contrast", y = "Total number of deregulated CCC interactions") +
  theme_bw(base_size = 12) +
  theme(legend.position = "none")
dev.off()

# Barplot of nr deregulated CCC interactions mediated by secreted and non-secreted ligands
pdf("results/liana/liana_condition-simplified_total_interactions_color_secreted_bplot.pdf", 
    height = 4.5, width = 5)
summary_df %>%
    pivot_longer(
        cols = c(secreted_ligand, non_secreted_ligand),
        names_to = "variable",
        values_to = "value"
    ) %>%
    mutate(variable = factor(variable, levels = c("secreted_ligand", "non_secreted_ligand"))) %>%
    ggplot(aes(x = reorder(contrast, -value), y = value, fill = variable)) +
        geom_bar(stat = "identity", color = "black") +
        labs(x = "Contrast", y = "Number of deregulated CCC interactions", fill = "Type of ligand") +
        scale_fill_hue(labels = c("secreted ligand", "non-secreted ligand"))  +
        theme_bw(base_size = 12) +
        theme(legend.position = "bottom")
dev.off()

# Barplot of nr interactions by sender cell
pdf("results/liana/liana_condition-simplified_total_interactions_sender_bplot.pdf", 
    height = 4.5, width = 5)
total_df %>%
    group_by(contrast, source) %>%
    summarise(n = n(), .groups = "drop") %>%
    ggplot(aes(x = reorder(contrast, -n), y = n, fill = source)) +
        geom_bar(stat = "identity", position = "stack") +
        scale_fill_manual(values = c("#e15759", 
                                     "#4e79a7", 
                                     "#59a14f", 
                                     "#b07aa1", 
                                     "#f28e2b", 
                                     "#edc948", 
                                     "#ff9da7", 
                                     "#9c755f")) + 
        labs(x = "Contrast", y = "Number of deregulated CCC interactions", fill = "Sender cell type") +
        theme_bw(base_size = 12)
dev.off()

pdf("results/liana/liana_condition-simplified_total_interactions_sender_pct_bplot.pdf", 
    height = 4.5, width = 5)
total_df %>%
    group_by(contrast, source) %>%
    summarise(n = n(), .groups = "drop") %>%
    ggplot(aes(x = reorder(contrast, -n), y = n, fill = source)) +
        geom_bar(stat = "identity", position = "fill") +
        scale_fill_manual(values = c("#e15759", 
                                     "#4e79a7", 
                                     "#59a14f", 
                                     "#b07aa1", 
                                     "#f28e2b", 
                                     "#edc948", 
                                     "#ff9da7", 
                                     "#9c755f")) + 
        labs(x = "Contrast", y = "Percentage of deregulated CCC interactions", fill = "Sender cell type") +
        theme_bw(base_size = 12)
dev.off()

# Save cell type nr interactions for cytoscape
ORABvControl %>%
    group_by(source, target) %>%
    summarise(n = n(),
              upregulated_interactions = sum(interaction_logFC > 0, na.rm = TRUE),
              downregulated_interactions  = sum(interaction_logFC < 0,  na.rm = TRUE),
              .groups = "drop") %>%
    write.csv("results/liana/liana_condition-simplified_ORABvControl_interactions_network.csv")

DBvControl %>%
    group_by(source, target) %>%
    summarise(n = n(),
              upregulated_interactions = sum(interaction_logFC > 0, na.rm = TRUE),
              downregulated_interactions  = sum(interaction_logFC < 0,  na.rm = TRUE),
              .groups = "drop") %>%
    write.csv("results/liana/liana_condition-simplified_DBvControl_interactions_network.csv")

DBvORAB %>%
    group_by(source, target) %>%
    summarise(n = n(),
              upregulated_interactions = sum(interaction_logFC > 0, na.rm = TRUE),
              downregulated_interactions  = sum(interaction_logFC < 0,  na.rm = TRUE),
              .groups = "drop") %>%
    write.csv("results/liana/liana_condition-simplified_DBvORAB_interactions_network.csv")

# Plot most frequent Fib ligands in ORABvControl
pdf("results/liana/liana_condition-simplified_ORABvControl_fibroblast_ligands.pdf", 
    height = 7, width = 6)
ORABvControl %>%
    filter(source == "Fib") %>%
    group_by(source, ligand_complex, target, interaction_logFC) %>%
    summarise(n = n(), .groups = "drop") %>%
    ggplot(aes(x = target, y = ligand_complex, size = n, color = interaction_logFC)) +
        geom_point() + 
        theme_bw(base_size = 12) +
        scale_color_gradient2(low = "blue", mid = "white", high = "red") +
        labs(y = "ligand", x = "receiver cell type", color = "interaction logFC", size = "nr interactions")
dev.off()

# Plot most frequent Endo ligands in OvSO
pdf("results/liana/liana_condition-simplified_ORABvControl_endo_ligands.pdf", 
    height = 6, width = 6)
ORABvControl %>%
    filter(source == "Endo") %>%
    group_by(source, ligand_complex, target, interaction_logFC) %>%
    summarise(n = n(), .groups = "drop") %>%
    ggplot(aes(x = target, y = ligand_complex, size = n, color = interaction_logFC)) +
        geom_point() + 
        theme_bw(base_size = 12) +
        scale_color_gradient2(low = "blue", mid = "white", high = "red") +
        labs(y = "ligand", x = "receiver cell type", color = "interaction logFC", size = "nr interactions")
dev.off()

# Explore candidate ligands for driving recovery
## Step 1: filter to keep ligands and receptors with logFC > 0 in DB v O.
DBvORAB <- DBvORAB %>%
  filter(interaction_logFC > 0) %>%
  filter(ligand_logFC > 0) %>%
  filter(receptor_logFC > 0)

## Step 2: filter to keep significantly differentially expressed receptors and ligands
DBvORAB <- DBvORAB %>%
  filter(interaction_FDR < 0.05) %>%
  filter(ligand_FDR < 0.05) %>%
  filter(receptor_FDR < 0.05)

## Step 2: filter to keep secreted ligands
DBvORAB <- DBvORAB %>%
  filter(ligand_secreted == "Yes")

## Write candidate table
write.csv(DBvORAB, "results/liana/liana_condition-simplified_DBvORAB_significant_results_candidates.csv")

## Dotplot of candidate interactions
pdf("results/liana/liana_condition-simplified_DBvORAB_significant_results_candidates_dotplot.pdf", 
    height = 4, width = 8)
ggplot(DBvORAB, aes(x = target, y = reorder(interaction, ligand_logCPM))) + 
  geom_point(aes(color = interaction_logCPM, size = -log10(interaction_FDR))) + 
  facet_wrap(~source, ncol = length(unique(DBvORAB$source))) + 
  theme_bw(base_size=12) +
  labs(y = "interaction", 
       x = "target cell type", 
       color = "interaction logCPM", 
       size = "-log10(interaction FDR)")+
  scale_color_gradient(low = "aliceblue", high = "#0571B0") +
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1))
dev.off()

## Explore ligand expression in bulk data
candidate_ligands <- unique(DBvORAB$ligand_complex)

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

coldat <- read_csv("./data/CRC_pseudobulk_coldata.csv") %>%
  dplyr::select(colname, psbulk_n_cells, cell_type, ConditionSimplified, sample_id, Tech_Batch) %>%
  column_to_rownames("colname") %>%
  dplyr::mutate(Condition = as.factor(ConditionSimplified))

pb_data <- pb_data[,rownames(coldat)]
coldat$pseudobulk <- rownames(coldat)

y <- DGEList(counts = pb_data, samples = DataFrame(coldat))
keep <- filterByExpr(y, group = coldat$Condition)
y <- y[keep,]
y <- calcNormFactors(y)
logCPM <- edgeR::cpm(y, log = TRUE)

for(candidate_gene in candidate_ligands){
    mat <- as.data.frame(logCPM[candidate_gene,])
    colnames(mat) <- "candidate_gene"
    mat$sample_id <- rownames(mat)
    mat <- merge(mat, coldat, by.x = "sample_id", by.y = "pseudobulk")
    mat$ConditionSimplified <- factor(mat$ConditionSimplified, levels = c("Control", "ORAB", "DB"))
    
    sender <- DBvORAB %>% filter(ligand_complex == candidate_gene) %>% pull(source) %>% unique()

    for(sendercell in sender){
        pdf(paste0("results/liana/liana_condition-simplified_DBvORAB_candidates_", 
               candidate_gene, 
               "_", 
               sender, 
               "_expression.pdf"),
            height = 2.5, width = 2.5)
        
        print(mat %>% 
            dplyr::filter(cell_type == sendercell) %>%
            ggplot(aes(x = ConditionSimplified, y = candidate_gene, fill = ConditionSimplified)) + 
                geom_boxplot() + 
                geom_point() +
                scale_fill_manual(values = c("Control" = "grey50",
                                             "ORAB" = "#AD1332",
                                             "DB" = "#0571b0")) +
                theme_bw() +
                labs(title = paste0(candidate_gene, " (", sendercell, " expression)"),
                     x = "Condition",
                     y = "Log-norm counts") +
                scale_y_continuous(limits = c(0, round(max(mat$candidate_gene)+1)))+
                theme(plot.title = element_text(hjust = 0.5), 
                      legend.position = "none")
              )
        dev.off()
    }
}


## Explore receptor expression in bulk data
candidate_receptors <- unique(DBvORAB$receptor_complex)
candidate_receptors <- unlist(sapply(candidate_receptor, function(x) unlist(strsplit(x, "_"))))
candidate_receptors <- unique(candidate_receptor)

for(candidate_gene in candidate_receptors){
  mat <- as.data.frame(logCPM[candidate_gene,])
  colnames(mat) <- "candidate_gene"
  mat$sample_id <- rownames(mat)
  mat <- merge(mat, coldat, by.x = "sample_id", by.y = "pseudobulk")
  mat$ConditionSimplified <- factor(mat$ConditionSimplified, levels = c("Control", "ORAB", "DB"))
  
  receiver <- DBvORAB %>% filter(grepl(candidate_gene, receptor_complex)) %>% pull(target) %>% unique()
  
  for(receivercell in receiver){
    pdf(paste0("results/liana/liana_condition-simplified_DBvORAB_candidate-receptor_", 
               candidate_gene, 
               "_", 
               receivercell, 
               "_expression.pdf"),
        height = 2.5, width = 2.5)
    
    print(mat %>% 
            dplyr::filter(cell_type == receivercell) %>%
            ggplot(aes(x = ConditionSimplified, y = candidate_gene, fill = ConditionSimplified)) + 
            geom_boxplot() + 
            geom_point() +
            scale_fill_manual(values = c("Control" = "grey50",
                                         "ORAB" = "#AD1332",
                                         "DB" = "#0571b0")) +
            theme_bw() +
            labs(title = paste0(candidate_gene, " (", receivercell, " expression)"),
                 x = "Condition",
                 y = "Log-norm counts") +
            scale_y_continuous(limits = c(0, round(max(mat$candidate_gene)+1)))+
            theme(plot.title = element_text(hjust = 0.5), 
                  legend.position = "none")
    )
    dev.off()
  }
}


# Plot logFC in sender cell type for each candidate ligand
degs <- read.csv("results/degs/edgeR_res_condition-simplified.csv", header = TRUE, sep = ",")

for(candidate_gene in candidate_ligands){
  sendercell <- DBvORAB %>% filter(ligand_complex == candidate_gene) %>% pull(source) %>% unique()
  degs_filt <- degs %>% filter(gene == candidate_gene) %>% filter(cell_type == sendercell)
  degs_filt <- degs_filt %>% mutate(signif = case_when(FDR < 0.001 ~ "***",
                                                       FDR < 0.01  ~ "**",
                                                       FDR < 0.05  ~ "*",
                                                       TRUE        ~ ""))
  
  pdf(paste0("results/liana/liana_condition-simplified_DBvORAB_candidates_", 
             candidate_gene, 
             "_", 
             sendercell, 
             "_logFC.pdf"),
      height = 2.5, width = 4)
  print(ggplot(degs_filt, aes(x = contrast, y = logFC, fill = logFC)) +
    geom_bar(stat = "identity", width = 0.5) +
    geom_hline(yintercept = 0, color = "black", linewidth = 0.15) +
    geom_text(aes(label = signif, y = ifelse(logFC < 0, logFC + abs(logFC)/2, logFC - logFC/2)), color = "white") +
    coord_flip() +
    scale_fill_gradient2(low = "blue", high = "red", mid = "white") +
    theme_classic() +
    labs(title = paste0(candidate_gene, " secreted by ", sendercell),
         x = "Contrast",
         y = "LogFC") +
    theme(plot.title = element_text(hjust = 0.5, face = "bold")),)
  dev.off()
}

# Plot logFC in receiver cell type for each candidate receptor
for(candidate_gene in candidate_receptors){
  receivercell <- DBvORAB %>% filter(grepl(candidate_gene, receptor_complex)) %>% pull(target) %>% unique()
  degs_filt <- degs %>% filter(gene == candidate_gene) %>% filter(cell_type == receivercell)
  degs_filt <- degs_filt %>% mutate(signif = case_when(FDR < 0.001 ~ "***",
                                                       FDR < 0.01  ~ "**",
                                                       FDR < 0.05  ~ "*",
                                                       TRUE        ~ ""))
  
  pdf(paste0("results/liana/liana_condition-simplified_DBvORAB_candidate-receptors_", 
             candidate_gene, 
             "_", 
             receivercell, 
             "_logFC.pdf"),
      height = 2.5, width = 4)
  print(ggplot(degs_filt, aes(x = contrast, y = logFC, fill = logFC)) +
          geom_bar(stat = "identity", width = 0.5) +
          geom_hline(yintercept = 0, color = "black", linewidth = 0.15) +
          geom_text(aes(label = signif, y = ifelse(logFC < 0, logFC + abs(logFC)/2, logFC - logFC/2)), color = "white") +
          coord_flip() +
          scale_fill_gradient2(low = "blue", high = "red", mid = "white") +
          theme_classic() +
          labs(title = paste0(candidate_gene, " expressed by ", receivercell),
               x = "Contrast",
               y = "LogFC") +
          theme(plot.title = element_text(hjust = 0.5, face = "bold")),)
  dev.off()
}