nuc_12.3_cell_states_CODA.R

# Copyright (c) [2024] [Ricardo O. Ramirez Flores]
# roramirezf@uni-heidelberg.de

#' Perform compositional analysis of the data
library(tidyverse)
library(compositions)
library(ggpubr)
library(ggrepel)

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

cts <- c("CM", "Endo", "Fib")

for(celltype in cts){
    print(celltype)
    
    # Load data
    meta_data <- read_csv(paste0("results/cell_states/", 
                                 celltype, 
                                 "_condition-simplified_cellstates_pseudobulk_coldata.csv")) %>%
        dplyr::select(colname, psbulk_n_cells, cell_states, ConditionSimplified, Tech_Batch, sample_id) %>%
        column_to_rownames("colname") %>%
        group_by(sample_id, ConditionSimplified, Tech_Batch) %>%
        dplyr::summarise(total_cells = sum(psbulk_n_cells)) %>%
        unique() %>%
        as_tibble()
    
    # Load col data
    clr_props_mat <- read_csv(paste0("results/cell_states/", 
                                     celltype, 
                                     "_condition-simplified_cellstates_pseudobulk_coldata.csv")) %>%
        dplyr::select(colname, psbulk_n_cells, cell_states, ConditionSimplified, Tech_Batch, sample_id) %>%
        column_to_rownames("colname") %>%
        group_by(sample_id) %>%
        dplyr::mutate(total_cells = sum(psbulk_n_cells)) %>%
        dplyr::select(-Tech_Batch) %>%
        dplyr::mutate(ct_props = psbulk_n_cells/total_cells) %>%
        dplyr::select(sample_id, cell_states, ct_props) %>%
        pivot_wider(names_from = cell_states, values_from = ct_props) %>%
        column_to_rownames("sample_id") %>%
        as.matrix() %>%
        compositions::acomp() %>%
        compositions::clr()
    
    # Test if there are differences in composition based on the condition
    long_clrs <- clr_props_mat %>%
        as.data.frame() %>%
        rownames_to_column("sample_id") %>%
        pivot_longer(-sample_id, names_to = "cell_states") %>%
        left_join(meta_data, by = "sample_id")

    aov_clrs <- long_clrs %>%
        group_by(cell_states) %>%
        nest() %>%
        dplyr::mutate(cond_test = map(data, function(dat){
            aov(value~ConditionSimplified, data = dat) %>%
            broom::tidy()
        })) %>%
        dplyr::select(cell_states, cond_test) %>%
        ungroup() %>%
        unnest() %>%
        dplyr::filter(term == "ConditionSimplified")

    # Adjust pvalue
    aov_clrs$p.adj <- p.adjust(aov_clrs$p.value, method = "BH")

    # Print result
    print(aov_clrs)

    # Save results
    write_csv(aov_clrs, paste0("results/cell_states/", 
                               celltype, 
                               "_condition-simplified_cellstates_aov_clr_res.csv"))
  
    # Visualization
    long_clrs$ConditionSimplified <- factor(long_clrs$ConditionSimplified, levels = c("Control", "ORAB", "DB"))
    ct_plts <- long_clrs %>%
        ggplot(aes(x = ConditionSimplified, y = value, fill = ConditionSimplified)) +
        geom_boxplot() +
        geom_point() +
        theme_bw(base_size = 12) +
        scale_fill_manual(values = c('ORAB'='#ad1332', 'DB'='#0571b0', 'Control' = 'grey50'))+
        ylab("Centered log of cell type proportions") +
        theme(axis.text.x = element_text(angle =90, hjust = 1, vjust = 0.5)) +
        facet_wrap(.~cell_states,scales = "free_y")

    pdf(paste0("results/cell_states/", 
               celltype, 
               "_condition-simplified_cellstates_clr_dists.pdf"),height = 4.5,width = 9)
    plot(ct_plts)
    dev.off()

    # Select cells with significant ANOVA to perform Tukey test on them.
    significant_cells <- aov_clrs %>% filter(p.adj < 0.05) %>% pull(cell_states)

    if(length(significant_cells) > 0){
        # Tukey test for selected contrasts
        tukey_results <- vector(mode = "list", length = length(significant_cells))
        for(i in 1:length(significant_cells)){
            ct <- unique(significant_cells)[i]
            model <- aov(value~ConditionSimplified, data=long_clrs[long_clrs$cell_states == ct,])
            tukey_results[[i]] <- TukeyHSD(model, conf.level=.95)$ConditionSimplified
            names(tukey_results)[i] <- ct
        }
        print(tukey_results)
        
        # Convert the list of matrices into one tidy data frame for plotting
        pvals_df <- map_dfr(names(tukey_results), function(cell_states) {
            df <- as.data.frame(tukey_results[[cell_states]])
            df <- df %>%
                mutate(comparison = rownames(tukey_results[[cell_states]]), cell_states = cell_states) %>%
                separate(comparison, into = c("group1", "group2"), sep = "-") %>%
                rename(diff = diff, lwr = lwr, upr = upr,p.adj = `p adj`)
            return(df)
        })
        rownames(pvals_df) <- 1:nrow(pvals_df)
        pvals_df <- pvals_df %>%
            mutate(signif_label = case_when(p.adj < 0.05 ~ "*", TRUE ~ "ns"))
        sig_pvals_df <- pvals_df %>% filter(p.adj < 0.05)
        sig_pvals_df <- sig_pvals_df %>%
            group_by(cell_states) %>%
            mutate(y_position = max(long_clrs$value[long_clrs$cell_states == unique(cell_states)]) * 1.05 + 0.2 * row_number()) %>%
            ungroup()

        # Visualisation
        ct_plts_signif <- long_clrs %>%
            ggplot(aes(x = ConditionSimplified, y = value, fill = ConditionSimplified)) +
                geom_boxplot(color="black") +
                geom_point(color = "black") +
                theme_bw() +
                scale_fill_manual(values = c("ORAB"="#ad1332", "DB"="#0571b0", "Control"="grey50"))+
                ylab("Centered log of cell type proportions") +
                theme(axis.text.x = element_text(angle =90, hjust = 1, vjust = 0.5)) +
                facet_wrap(.~cell_states, nrow = 2, scales= "free_y") +
                theme(legend.position = "none") +
                stat_pvalue_manual(data = sig_pvals_df,
                               label = "signif_label",
                               xmin = "group1",
                               xmax = "group2",
                               y.position = "y_position",
                               tip.length = 0.01,
                               bracket.size = 0.3,
                               size = 4,
                               inherit.aes = FALSE) +
                scale_y_continuous(expand = expansion(mult = c(0.05, 0.1)))

        pdf(paste0("results/cell_states/", 
                   celltype, 
                   "_condition-simplified_cellstates_clr_dists_signif.pdf"),
            height = 4, width = 8)
        plot(ct_plts_signif)
        dev.off()
    }

    # Visualise actual proportions with Tukey significances
    props_mat <- read_csv(paste0("results/cell_states/", 
                                     celltype, 
                                     "_condition-simplified_cellstates_pseudobulk_coldata.csv")) %>%
        dplyr::select(colname, psbulk_n_cells, cell_states, ConditionSimplified, Tech_Batch, sample_id) %>%
        column_to_rownames("colname") %>%
        group_by(sample_id) %>%
        dplyr::mutate(total_cells = sum(psbulk_n_cells)) %>%
        dplyr::select(-Tech_Batch) %>%
        dplyr::mutate(ct_props = psbulk_n_cells/total_cells) %>%
        dplyr::select(sample_id, cell_states, ct_props) %>%
        pivot_wider(names_from = cell_states, values_from = ct_props) %>%
        column_to_rownames("sample_id") %>%
        as.matrix() #%>%
        #compositions::acomp() %>%
        #compositions::clr()

    long_props <- props_mat %>%
        as.data.frame() %>%
        rownames_to_column("sample_id") %>%
        pivot_longer(-sample_id, names_to = "cell_states") %>%
        left_join(meta_data, by = "sample_id")

    long_props$ConditionSimplified <- factor(long_props$ConditionSimplified, levels = c("Control", "ORAB", "DB"))

    # Not all pseudobulk samples contain nuclei from all cell states, therefore remove NA's
    long_props <- long_props[complete.cases(long_props),]

    # Change y_position of statistical significance marks
    sig_pvals_df <- sig_pvals_df %>%
        group_by(cell_states) %>%
        mutate(y_position = max(long_props$value[long_props$cell_states == unique(cell_states)])+0.05*row_number()) %>%
        ungroup()                                        

    prop_plts <- long_props %>%
        ggplot(aes(x = ConditionSimplified, y = value, color = ConditionSimplified, 
                   fill = ConditionSimplified)) +
            geom_boxplot(alpha = 0.25) +
            geom_point() +
            theme_bw() +
            scale_color_manual(values = c("ORAB"="#AD1332", "DB"="#0571B0", "Control"="grey50"))+
            scale_fill_manual(values = c("ORAB"="#AD1332", "DB"="#0571B0", "Control"="grey50"))+
            ylab("Cell state proportions") +
            theme(axis.text.x = element_text(angle =90, hjust = 1, vjust = 0.5)) +
            facet_wrap(.~cell_states, nrow = 2) +
            #scale_y_continuous(limits = c(0, 1)) +
            theme(legend.position = "none") +
            stat_pvalue_manual(data = sig_pvals_df,
                               label = "signif_label",
                               xmin = "group1",
                               xmax = "group2",
                               y.position = "y_position",
                               tip.length = 0.01,
                               bracket.size = 0.3,
                               size = 4,
                               inherit.aes = FALSE)

    pdf(paste0("results/cell_states/", 
               celltype, 
               "_condition-simplified_cellstates_proportions_plot.pdf"),
              height = 4, width = 8)
    plot(prop_plts)
    dev.off()

    
    # PCA
    pca_comps <- prcomp(clr_props_mat) %>% summary()
    manifold <- pca_comps$x
    expl_var <- pca_comps$importance["Proportion of Variance",] %>%
      enframe("PC", "expl_var")
    loadings <- pca_comps$rotation

    long_manifold <- manifold %>%
        as.data.frame() %>%
        rownames_to_column("sample_id") %>%
        pivot_longer(-sample_id, names_to = "PC") %>%
        left_join(meta_data, by = "sample_id")

    pca_plt <- long_manifold %>%
        dplyr::filter(PC %in% c("PC1", "PC2")) %>%
        pivot_wider(names_from = "PC", values_from = "value") %>%
        ggplot(aes(x = PC1, y = PC2, color = ConditionSimplified)) +
        theme_bw(base_size = 12) +
        geom_point(size = 2) +
        scale_color_manual(values = c('ORAB'='#ad1332', 'DB'='#0571b0', 'Control' = 'grey50'))
    
    pdf(paste0("results/cell_states/", 
               celltype, 
               "_condition-simplified_cellstates_clr_pca.pdf"),height = 3,width = 3)
    plot(pca_plt)
    dev.off()

    # Test that there are no compositional differences between batches
    clr_props_mat <- read_csv(paste0("results/cell_states/", 
                                     celltype, 
                                     "_condition-simplified_cellstates_pseudobulk_coldata.csv")) %>%
      dplyr::select(colname, psbulk_n_cells, cell_states, ConditionSimplified, Tech_Batch, sample_id) %>%
      column_to_rownames("colname") %>%
      group_by(sample_id) %>%
      dplyr::mutate(total_cells = sum(psbulk_n_cells)) %>%
      dplyr::select(-ConditionSimplified) %>%
      dplyr::mutate(ct_props = psbulk_n_cells/total_cells) %>%
      dplyr::select(sample_id, cell_states, ct_props) %>%
      pivot_wider(names_from = cell_states, values_from = ct_props) %>%
      column_to_rownames("sample_id") %>%
      as.matrix() %>%
      compositions::acomp() %>%
      compositions::clr()

    long_clrs <- clr_props_mat %>%
      as.data.frame() %>%
      rownames_to_column("sample_id") %>%
      pivot_longer(-sample_id, names_to = "cell_states") %>%
      left_join(meta_data, by = "sample_id")

    aov_clrs <- long_clrs %>%
      group_by(cell_states) %>%
      nest() %>%
      dplyr::mutate(cond_test = map(data, function(dat){
          aov(value~Tech_Batch, data = dat) %>%
          broom::tidy()
      })) %>%
      dplyr::select(cell_states, cond_test) %>%
      ungroup() %>%
      unnest() %>%
      dplyr::filter(term == "Tech_Batch")

    # Adjust pvalue
    aov_clrs$p.adj <- p.adjust(aov_clrs$p.value, method = "BH")  

    print(aov_clrs)
}