Manuscript figures v1 (#63)

* Figure 1 UMAPs. Color scheme. Other setup #49 #55

* Enrichment plotting improvements, keeping filtering to keep top terms #49 #52 #55

* Split into smk rulefiles per dataset (rather than per figure) #49 #52 #55

* Enrichment plots: Select only terms from bone marrow or PBMC. Sort by developmental hierarchy. Two terms per cell type to allow for more overlap #49 #52 #55

* Add modality in title of each plot #49 #52 #55

* Enrichment terms abbreviations #49 #52 #55

* Don't remove levels, since this make duplicate terms #49 #52 #55

* First draft logFC heatmap #49 #52 #55

* Enrichment heatmap labelling improvements & removal of tiny dots #49 #52 #55

* Differentially accessible genes (max region per gene) #49 $52 #55

* Differentially accessible regions #49 $52 #55

* Enrichment analysis axis label improvement #49 #52 #55

* Add downregulated differential features #49 #52 #55

* Labelling improvements #49 #52 #55

* Final heatmap manually ordered by cell type #49 #52 #55

* Data driven reordering with dendsort #49 #52 #55

* Improve cell type colors. Add clustering methods to ggheatmap from ggalign. #49 #52 #55

* Patchworked differential heatmap for RNA #49 #52 #55

* Handcrafted marker labelling function for differential heatmap #49 #52 #55

* Adapt differential heatmap for ATAC (mapping regions to genenames, keeping duplicates, subsetting markers) #49 #52 #55

* Test marker annotation passed #49 #52 #55

* differential heatmap: colors on top, text labels, marker annot for RNA good spacing, labels not cut off #49 #52 #55

* If too many markers, choose unique top markers per cell type #49 #52 #55

* Cuttree into clusters of genes, then sort by x axis order of cell types, to get data driven diagonal heatmap #49 #52 #55

* Add title to differential heatmap #49 #52 #55

* Better layout for heatmap #49 #52 #55

* Already nice crossprediction plot with automatic layout #49 #52 #55

* Already nice crossprediction plot with hierarchical layout #49 #52 #55

* Add lines for missing edges #49 #52 #55

* Color additional edges compared to Corces grey. Add legend #49 #52 #55

* Modualrize crossprediction plot and apply to ATAC data #49 #52 #55

* Code layout #49 #52 #55

* EP vs TA scatterplot v1 #49 #52 #55

* Rasterize epi potential scatterplots #56

* Umap for integrated and unintegrated #56

* CFA plots #56

* Figure 2 fix text and image size #52

* Figure 3 fix text and image size #56

* First draft INT GSEA heatmap #56

* INT GSEA heatmap improvements #56

* INT GSEA heatmap for Reactome in addition to GOBP #56

* Reload RNA figures script that was deleted by accident

* TF plot Mono v1 #56

* TF plot HSC and Mono v2 #56

* Prioritise papalexi TFs #56

* Color code TF EP TA overlap between cell types. Boldness code overlap with Papalexi dataset. #56

* Papalexi UMAPs and first draft cross prediction #61

* Papalexi crossprediction plot #61

* Lollipop plot for TA KO enrichment. Better titles. #61

* Prepare for file syncing to local

* Cell cycle umaps #61

* Bubble plot for all TF enrichments in TA #61

* Differential heatmap KO signatures #61

* GSEA bubble plot Papalexi #61

* Snakemakify figures.smk #62

* Update CorcesINT unsupervised annot to match main #62

* Bug fixes: Wrong paths etc #62

---------

Co-authored-by: Stephan Reichl <53785552+sreichl@users.noreply.github.com>

Author: bednarsky

.gitignore

@@ -156,6 +156,14 @@ metadata/
 */.snakemake/*
 .snakemake/
 
+paper
+
 # macOS specific files
 *DS_Store*
-._*
+._*
+
+# vscode
+settings.json
+
+# just a helper, don't want to sync this
+rsync_input_files_to_local.sh

config/config.yaml

@@ -8,8 +8,8 @@ project_name: "MrBiomics"
 
 #### Datasets and Workflows to include ###
 workflows:
-    quickstart:
-        enrichment_analysis: "config/quickstart/quickstart_enrichment_analysis_config.yaml"
+    # quickstart:
+    #     enrichment_analysis: "config/quickstart/quickstart_enrichment_analysis_config.yaml"
     digits:
         unsupervised_analysis: "config/digits/digits_unsupervised_analysis_config.yaml"
     CorcesRNA:

workflow/Snakefile

@@ -38,7 +38,7 @@ include: os.path.join("rules", "CorcesATAC.smk")
 include: os.path.join("rules", "CorcesINT.smk")
 include: os.path.join("rules", "Papalexi2021scCRISPR.smk")
 include: os.path.join("rules", "docs.smk")
-
+include: os.path.join("rules", "figures.smk")
 
 #### Target Rule ####
 rule all:
@@ -81,10 +81,15 @@ rule all:
         rules.Papalexi2021scCRISPR_enrichment_analysis_all.input,
         rules.Papalexi2021scCRISPR_functional_similarity.output, # custom rule
         #### Custom rules used only for docs/wiki (keep commented) ####
-        rules.CorcesRNA_plots.output,
-        rules.CorcesATAC_plots.output,
-        # rules.Papalexi2021scCRISPR_plots.output,
+        # rules.CorcesRNA_plots.output,
+        # rules.CorcesATAC_plots.output,
         # rules.CorcesINT_plots.output,
+        # rules.Papalexi2021scCRISPR_plots.output,
+        #### Figures ####
+        rules.CorcesRNA_figures.output,
+        rules.CorcesATAC_figures.output,
+        rules.CorcesINT_figures.output,
+        rules.Papalexi_figures.output,
     default_target: True # important for module usage
     threads: 1
     resources:

workflow/envs/ggplot.yaml

@@ -1,30 +1,34 @@
+name: ggplot
 channels:
   - conda-forge
   - bioconda
   - nodefaults
 dependencies:
-  - r-ggplot2
-  - r-tidyr
-  - r-tibble
-  - r-patchwork
-  - r-scales
-  - r-dplyr
-  - r-reshape2
-  - r-pheatmap
-  - r-ggrepel
-  - r-stringr
-  - r-gridGraphics
-  - r-ggplotify
-  - r-directlabels
-  - r-data.table
-  - r-svglite
-  - r-sna
-  - r-RColorBrewer
-  - r-Hmisc
-  - r-cairo
-  - cairo
-  - pango
-  - libpng
-  - libjpeg-turbo
-  - libtiff
-  - gtk2
+  - r-ggplot2=3.5
+  - r-ggdendro=0.2
+  - r-ggrepel=0.9
+  - r-ggrastr=1.0
+  - r-dendsort=0.3
+  - r-tidyr=1.3
+  - r-tibble=3.3
+  - r-patchwork=1.3
+  - r-scales=1.4
+  - r-dplyr=1.1
+  - r-reshape2=1.4
+  - r-stringr=1.5
+  - r-gridgraphics=0.5
+  - r-ggplotify=0.1
+  - r-directlabels=2025.6
+  - r-data.table=1.17
+  - r-svglite=2.1
+  - r-sna=2.8
+  - r-rcolorbrewer=1.1
+  - r-hmisc=5.2
+  - r-cairo=1.6
+  - cairo=1.18
+  - pango=1.56
+  - libpng=1.6
+  - libjpeg-turbo=3.1
+  - libtiff=4.7
+  - gtk2=2.24
+  - r-ggtext=0.1

workflow/rules/figures.smk

@@ -1,18 +1,26 @@
 # Snakemake file for creating the figure for the manuscript
-
-# TODO: think of sensible place to put plots? maybe MrBiomics/docs/figure1 or MrBiomics/results/figures?
-
-#### CorcesRNA - Figures (custom rule) ####
+#### CorcesRNA - Figure panels (custom rule) ####
 rule CorcesRNA_figures:
     input:
+        umap_coords = os.path.join("results/CorcesRNA/unsupervised_analysis/normCQN_integrated/UMAP/UMAP_correlation_15_0.1_2_data.csv"),
+        dea_ova = os.path.join("results/CorcesRNA/dea_limma/normCQN_OvA_cell_type/results.csv"),
         enrichment_results = os.path.join("results/CorcesRNA/enrichment_analysis/cell_types/preranked_GSEApy/Azimuth_2023/cell_types_Azimuth_2023_all.csv"),
         crossprediction_adj_mtx = os.path.join("results/CorcesRNA/special_analyses/crossprediction/adjacency_matrix.csv"),
     output:
-        enrichment_plot = os.path.join("docs/CorcesRNA/enrichment_analysis.pdf"),
-        crossprediction_plot = os.path.join("docs/CorcesRNA/crossprediction_plot.pdf"),
+        umap_plot = os.path.join("paper/CorcesRNA/umap.pdf"),
+        dea_heatmap_plot = os.path.join("paper/CorcesRNA/differential_heatmap.pdf"),
+        enrichment_plot = os.path.join("paper/CorcesRNA/enrichment.pdf"),
+        crossprediction_plot = os.path.join("paper/CorcesRNA/crossprediction.pdf"),
+        crossprediction_coordinates = os.path.join("paper/CorcesRNA/crossprediction_coordinates.csv"),
     params:
         figure_theme_path = workflow.source_path("../scripts/figure_theme.R"),
-        adj_p = 0.05, # unused
+        figure_utils_path = workflow.source_path("../scripts/figure_utils.R"),
+        # enrichment analysis
+        fdr_threshold = 0.05,
+        log2FC_threshold = 3,
+        # lineage reconstructions
+        lineage_tree_cut_off = 0.07,
+        hierarchy_coordinates = True
     resources:
         mem_mb=config.get("mem", "16000"),
     threads: config.get("threads", 1)
@@ -22,3 +30,113 @@ rule CorcesRNA_figures:
         os.path.join("logs","rules","CorcesRNA_figures.log"),
     script:
         "../scripts/CorcesRNA_figures.R"
+
+#### CorcesATAC - Figure panels (custom rule) ####
+rule CorcesATAC_figures:
+    input:
+        umap_coords = os.path.join("results/CorcesATAC/unsupervised_analysis/normCQN_integrated/UMAP/UMAP_correlation_15_0.1_2_data.csv"),
+        dea_ova = os.path.join("results/CorcesATAC/dea_limma/normCQN_OvA_cell_type/results.csv"),
+        enrichment_results = os.path.join("results/CorcesATAC/enrichment_analysis/cell_types_up/GREAT/Azimuth_2023/cell_types_up_Azimuth_2023_all.csv"),
+        crossprediction_adj_mtx = os.path.join("results/CorcesATAC/special_analyses/crossprediction/adjacency_matrix.csv"),
+    output:
+        umap_plot = os.path.join("paper/CorcesATAC/umap.pdf"),
+        dea_heatmap_plot = os.path.join("paper/CorcesATAC/differential_heatmap.pdf"),
+        enrichment_plot = os.path.join("paper/CorcesATAC/enrichment.pdf"),
+        crossprediction_plot = os.path.join("paper/CorcesATAC/crossprediction.pdf"),
+        crossprediction_coordinates = os.path.join("paper/CorcesATAC/crossprediction_coordinates.csv"),
+    params:
+        figure_theme_path = workflow.source_path("../scripts/figure_theme.R"),
+        figure_utils_path = workflow.source_path("../scripts/figure_utils.R"),
+        fdr_threshold = 0.05,
+        log2FC_threshold = 3,
+        lineage_tree_cut_off = 0.2,
+        hierarchy_coordinates = True
+    resources:
+        mem_mb=config.get("mem", "16000"),
+    threads: config.get("threads", 1)
+    conda:
+        "../envs/ggplot.yaml"
+    log:
+        os.path.join("logs","rules","CorcesATAC_figures.log"),
+    script:
+        "../scripts/CorcesATAC_figures.R"
+
+# #### CorcesINT - Figure panels (custom rule) ####
+rule CorcesINT_figures:
+    input:
+        unintegrated_umap_coords = os.path.join("results/CorcesINT/unsupervised_analysis/normupperquartile/UMAP/UMAP_correlation_15_0.1_2_data.csv"),
+        integrated_umap_coords = os.path.join("results/CorcesINT/unsupervised_analysis/normupperquartile_integrated/UMAP/UMAP_correlation_15_0.1_2_data.csv"),
+        unintegrated_cfa = os.path.join("results/CorcesINT/spilterlize_integrate/all/normupperquartile_CFA.csv"),
+        integrated_cfa = os.path.join("results/CorcesINT/spilterlize_integrate/all/normupperquartile_integrated_CFA.csv"),
+        norm_counts = os.path.join("results/CorcesINT/spilterlize_integrate/all/normupperquartile_integrated.csv"),
+        metadata = os.path.join("results/CorcesINT/spilterlize_integrate/all/annotation.csv"),
+        dea_results = os.path.join("results/CorcesINT/dea_limma/normupperquartile_integrated/results.csv"),
+        gene_annotation = os.path.join("results/CorcesRNA/rnaseq_pipeline/counts/gene_annotation.csv"),
+        GO_enrichment_results = os.path.join("results/CorcesINT/enrichment_analysis/cell_types/preranked_GSEApy/GO_Biological_Process_2025/cell_types_GO_Biological_Process_2025_all.csv"),
+        Reactome_enrichment_results = os.path.join("results/CorcesINT/enrichment_analysis/cell_types/preranked_GSEApy/ReactomePathways/cell_types_ReactomePathways_all.csv"),
+        Mono_TF_EP = os.path.join("results/CorcesINT/enrichment_analysis/Mono_EP/RcisTarget/hg38_500bp_up_100bp_down_v10clust/Mono_EP_hg38_500bp_up_100bp_down_v10clust.csv"),
+        Mono_TF_TA = os.path.join("results/CorcesINT/enrichment_analysis/Mono_TA/RcisTarget/hg38_500bp_up_100bp_down_v10clust/Mono_TA_hg38_500bp_up_100bp_down_v10clust.csv"),
+        HSC_TF_EP = os.path.join("results/CorcesINT/enrichment_analysis/HSC_EP/RcisTarget/hg38_500bp_up_100bp_down_v10clust/HSC_EP_hg38_500bp_up_100bp_down_v10clust.csv"),
+        HSC_TF_TA = os.path.join("results/CorcesINT/enrichment_analysis/HSC_TA/RcisTarget/hg38_500bp_up_100bp_down_v10clust/HSC_TA_hg38_500bp_up_100bp_down_v10clust.csv"),
+    output:
+        unintegrated_cfa_plot = os.path.join("paper/CorcesINT/unintegrated_cfa.pdf"),
+        integrated_cfa_plot = os.path.join("paper/CorcesINT/integrated_cfa.pdf"),
+        integrated_umap_plot = os.path.join("paper/CorcesINT/integrated_umap.pdf"),
+        unintegrated_umap_plot = os.path.join("paper/CorcesINT/unintegrated_umap.pdf"),
+        epigenetic_scatter_dir = directory(os.path.join("paper/CorcesINT/correlation_plots")),
+        GO_enrichment_plot = os.path.join("paper/CorcesINT/GO_enrichment.pdf"),
+        Reactome_enrichment_plot = os.path.join("paper/CorcesINT/Reactome_enrichment.pdf"),
+        TF_plot = os.path.join("paper/CorcesINT/TF_lollipop.pdf"),
+    params:
+        figure_theme_path = workflow.source_path("../scripts/figure_theme.R"),
+        figure_utils_path = workflow.source_path("../scripts/figure_utils.R"),
+        # thresholds for EP/TA categorization
+        adjp_th = 0.05,
+        fdr_threshold = 0.05,
+        lfc_th = 1,
+        ave_expr_th = 0,
+        max_genes_tf_plot = 25
+    resources:
+        mem_mb=config.get("mem", "16000"),
+    threads: config.get("threads", 1)
+    conda:
+        "../envs/ggplot.yaml"
+    log:
+        os.path.join("logs","rules","CorcesINT_figures.log"),
+    script:
+        "../scripts/CorcesINT_figures.R"
+
+#### Papalexi - Figure panels (custom rule) ####
+rule Papalexi_figures:
+    input:
+        CORRECTED_umap_coords = os.path.join("results/Papalexi2021scCRISPR/unsupervised_analysis/merged_CORRECTED/UMAP/UMAP_correlation_10_0.1_2_data.csv"),
+        CORRECTED_metadata = os.path.join("results/Papalexi2021scCRISPR/scrnaseq_processing_seurat/merged/CORRECTED/metadata.csv"),
+        MIXSCAPE_umap_coords = os.path.join("results/Papalexi2021scCRISPR/unsupervised_analysis/merged_MIXSCAPE_LDA/UMAP/UMAP_correlation_10_0.1_2_data.csv"),
+        MIXSCAPE_metadata = os.path.join("results/Papalexi2021scCRISPR/mixscape_seurat/merged/FILTERED_metadata.csv"),
+        KO_crossprediction_adj_mtx = os.path.join("results/Papalexi2021scCRISPR/special_analyses/crossprediction/adjacency_matrix.csv"),
+        SPI1_TA_results = os.path.join("results/Papalexi2021scCRISPR/enrichment_analysis/SPI1/preranked_GSEApy/Corces_TA_signatures/SPI1_Corces_TA_signatures.csv"),
+        KO_DEA_results = os.path.join("results/Papalexi2021scCRISPR/dea_seurat/KO_mixscape/results.csv"),
+        KO_enrichment_results_GOBP = os.path.join("results/Papalexi2021scCRISPR/enrichment_analysis/KO/preranked_GSEApy/GO_Biological_Process_2025/KO_GO_Biological_Process_2025_all.csv"),
+        KO_enrichment_results_Reactome = os.path.join("results/Papalexi2021scCRISPR/enrichment_analysis/KO/preranked_GSEApy/ReactomePathways/KO_ReactomePathways_all.csv"),
+    output:
+        umap_corrected_fig = os.path.join("paper/Papalexi/umap_CORRECTED.pdf"),
+        umap_lda_fig = os.path.join("paper/Papalexi/umap_LDA.pdf"),
+        crossprediction_fig = os.path.join("paper/Papalexi/crossprediction.pdf"),
+        spi1_ta_lollipop_fig = os.path.join("paper/Papalexi/SPI1_TA_lollipop.pdf"),
+        ko_DEA_heatmap_fig = os.path.join("paper/Papalexi/KO_DEA_heatmap.pdf"),
+        ko_enrichment_GOBP_fig = os.path.join("paper/Papalexi/KO_enrichment_GOBP_bubbleplot.pdf"),
+    params:
+        figure_theme_path = workflow.source_path("../scripts/figure_theme.R"),
+        figure_utils_path = workflow.source_path("../scripts/figure_utils.R"),
+        ko_column = "gene",
+        phase_column = "Phase",
+        fdr_threshold = 0.05
+    resources:
+        mem_mb=config.get("mem", "16000"),
+    threads: config.get("threads", 1)
+    conda:
+        "../envs/ggplot.yaml"
+    log:
+        os.path.join("logs","rules","Papalexi_figures.log"),
+    script:
+        "../scripts/Papalexi_figures.R"

workflow/scripts/CorcesATAC_figures.R

@@ -0,0 +1,101 @@
+# set correct working directory -> project/repo folder (only for dev, not prod)
+# getwd()
+# setwd('..')
+# getwd()
+
+#### libraries
+# source this for libraries, MrBiomics theme and utility functions
+source(snakemake@params[["figure_theme_path"]])
+source(snakemake@params[["figure_utils_path"]])
+
+# FIXME snakemakeify
+# ## libraries are loaded in the source files
+# source(snakemake@params[["figure_theme_path"]])
+# ## input
+# enrichment_results_path <- snakemake@input[["enrichment_results"]]
+# ## output
+# enrichment_plot_path <- snakemake@output[["enrichment_plot"]]
+# ## params
+# fdr_threshold <- snakemake@params[["fdr_threshold"]]
+
+# inputs
+CorcesATAC_umap_coords_path <- snakemake@input[["umap_coords"]]
+CorcesATAC_dea_OvA_path <- snakemake@input[["dea_ova"]]
+CorcesATAC_enrichment_results_path <- snakemake@input[["enrichment_results"]]
+CorcesATAC_crossprediction_adj_mtx_path <- snakemake@input[["crossprediction_adj_mtx"]]
+
+# outputs
+atac_umap_path <- snakemake@output[["umap_plot"]]
+atac_dea_heatmap_path <- snakemake@output[["dea_heatmap_plot"]]
+atac_enrichment_path <- snakemake@output[["enrichment_plot"]]
+atac_crossprediction_path <- snakemake@output[["crossprediction_plot"]]
+atac_crossprediction_coordinates_path <- snakemake@output[["crossprediction_coordinates"]]
+
+# params
+fdr_threshold <- snakemake@params[["fdr_threshold"]]
+log2FC_threshold <- snakemake@params[["log2FC_threshold"]]
+lineage_tree_cut_off <- snakemake@params[["lineage_tree_cut_off"]]
+hierarchy_coordinates <- snakemake@params[["hierarchy_coordinates"]]
+######### UMAPs (from unsupervised analysis) ############
+# Create UMAP plots
+atac_umap_plot <- umap_plot(CorcesATAC_umap_coords_path, atac_umap_path, title = "ATAC")
+
+
+######### DEA HEATMAP ############
+atac_dea_heatmap_plot <- plot_differential_features_heatmap(
+    dea_results_path = CorcesATAC_dea_OvA_path,
+    fig_path = atac_dea_heatmap_path,
+    fdr_threshold = fdr_threshold,
+    log2FC_threshold = log2FC_threshold,
+    title = "ATAC",
+    feature = 'Regions',
+    ct_clst_dist = "euclidean",
+    ct_clst_method = "ward.D2",
+    feature_clst_dist = "maximum",
+    feature_clst_method = "ward.D2",
+    q_mask = 0.025,
+    label_box_size_factor = 1
+)
+
+
+######### ENRICHMENT HEATMAP ############
+create_atac_enrichment_df <- function(enrichment_results_path, fdr_threshold = 0.05) {
+    # Load enrichment analysis result
+    df <- data.frame(fread(file.path(enrichment_results_path), header=TRUE))
+    
+    # Adapt for ATAC data
+    df_formatted <- df %>%
+        rename(Term = description, statistic = p_adjust, score = fold_enrichment) %>%
+        mutate(name = sub("_up$", "", name),
+                score = ifelse(is.infinite(log2(score)), NaN, log2(score))) %>%
+        mutate(name = recode(name, !!!DATA_TO_CELL_TYPE_COLORS_MAPPING))
+    
+    return(df_formatted)
+}
+
+atac_df_formatted <- create_atac_enrichment_df(CorcesATAC_enrichment_results_path, fdr_threshold)
+atac_heatmap_df <- prepare_for_heatmap(df_formatted = atac_df_formatted,
+                                      fdr_threshold = fdr_threshold,
+                                      tissues_to_keep = c("PBMC", "Bone Marrow"))
+atac_enrichment_plot <- plot_enrichment_heatmap(
+    heatmap_df = atac_heatmap_df, 
+    fig_path = atac_enrichment_path,
+    fill_lab = "log2(FE)",
+    size_lab = "-log10(p-adj.)",
+    title = "ATAC",
+    ylabel = "Enrichment term\n(GREAT, Azimuth 2023)"
+)
+
+
+######### Lineage reconstruction using crossprediction ############
+# plot adapted from: https://gist.github.com/dsparks/4331058
+
+
+atac_crosspred_p <- plot_crossprediction_from_adjacency(
+    adjacency_matrix_path = CorcesATAC_crossprediction_adj_mtx_path,
+    fig_path = atac_crossprediction_path,
+    coordinates_out_path = atac_crossprediction_coordinates_path,
+    lineage_tree_cut_off = lineage_tree_cut_off,
+    hierarchy_coordinates = hierarchy_coordinates,
+    modality_label = "ATAC",
+)