Add meta comparisons

config/config.yaml

@@ -51,15 +51,20 @@ diffexp:
   # model for sleuth differential expression analysis
   models:
     model_X:
-      full: ~condition + batch_effect
-      reduced: ~batch_effect
+      full: ~condition_X + batch_effect_X
+      reduced: ~batch_effect_X
       # Covariate / sample sheet column that shall be used for fold
       # change/effect size based downstream analyses.
-      primary_variable: condition
+      primary_variable: condition_X
       # base level of the primary variable (this should be one of the entries
       # in the primary_variable sample sheet column and will be considered as
       # denominator in the fold change/effect size estimation).
-      base_level: untreated
+      base_level: untreated_X
+    model_Y:
+      full: ~condition_Y + batch_effect_Y
+      reduced: ~batch_effect_Y
+      primary_variable: condition_Y
+      base_level: untreated_Y
   # significance level to use for volcano, ma- and qq-plots
   sig-level:
     volcano-plot: 0.05
@@ -105,6 +110,21 @@ enrichment:
     # the species specified by resources -> ref -> species above
     pathway_database: "reactome"
 
+meta_comparisons:
+  # comparison is only run if set to `true`
+  activate: false
+  # Define any name for comparison
+  model_X_vs_model_Y:
+    items:
+      # Define the two underlying models for the comparison. The models must be defined in the diffexp/models in the config
+      - label: X # Arbitrary label
+        model: model_X # Must match a diffexp model from config
+      - label: Y # Arbitrary label
+        model: model_Y # Must match a diffexp model from config
+    desc: |
+      Comparison between model_X and model_Y effect.
+    label: model_X vs. model_Y effect
+
 bootstrap_plots:
   # desired false discovery rate for bootstrap plots, i.e. a lower FDR will result in fewer boxplots generated
   FDR: 0.01

workflow/envs/pystats.yaml

@@ -0,0 +1,16 @@
+channels:
+  - conda-forge
+  - nodefaults
+dependencies:
+  - polars =0.20.28
+  - pyreadr =0.5
+  - altair =5.2
+  - pyarrow =16.1
+  - vegafusion =1.6
+  - vegafusion-python-embed =1.6
+  - vl-convert-python =1.2
+  - jupyter_core =5.7
+  - ipykernel =6.29
+  - nbconvert =7.14
+  - notebook =7.0
+  - jupyterlab_code_formatter =1.4

workflow/report/meta_comparison_diffexp.rst

@@ -0,0 +1,5 @@
+{{ snakemake.params.desc }}
+Each point represents a gene, x axis shows the {{ snakemake.params.labels[0] }} effect, y-axis the {{ snakemake.params.labels[1] }} effect (both as log2 fold change).
+The color encodes the corresponding q-value.
+By clicking on points, their label can be displayed.
+Holding the Shift key allows to select or deselect labels for multiple genes.

workflow/rules/common.smk

@@ -396,4 +396,17 @@ def all_input(wildcards):
                 ind_transcripts=config["experiment"]["3-prime-rna-seq"]["plot-qc"],
             )
         )
+
+    # meta comparisons
+    if config["meta_comparisons"]["activate"]:
+        wanted_input.extend(
+            directory(
+                expand(
+                    "results/datavzrd-reports/{report_type}_meta_comparison_{meta_comp}",
+                    report_type=["go_terms", "diffexp", "pathways"],
+                    meta_comp=lookup(dpath="meta_comparisons", within=config),
+                )
+            ),
+        )
+
     return wanted_input

workflow/rules/meta_comparisons.smk

@@ -0,0 +1,70 @@
+rule meta_compare_diffex:
+    input:
+        expand(
+            "results/sleuth/diffexp/{model}.genes-representative.diffexp.rds",
+            model=lookup(
+                dpath="meta_comparisons/{meta_comp}/items/*/model", within=config
+            ),
+        ),
+    output:
+        "results/tables/diffexp/meta_compare_{meta_comp}.tsv",
+        "results/meta_comparison/diffexp/{meta_comp}.json",
+    log:
+        notebook="logs/meta_compare_diffexp/{meta_comp}.ipynb",
+    params:
+        desc=lookup(dpath="meta_comparisons/{meta_comp}/desc", within=config),
+        labels=lookup(dpath="meta_comparisons/{meta_comp}/items/*/label", within=config),
+    conda:
+        "../envs/pystats.yaml"
+    notebook:
+        "../scripts/compare_diffexp.py.ipynb"
+
+
+rule meta_compare_enrichment:
+    input:
+        expand(
+            "results/tables/go_terms/{model}.go_term_enrichment.gene_fdr_{gene_fdr}.go_term_sig_study_fdr_{go_term_fdr}.tsv",
+            model=lookup(
+                dpath="meta_comparisons/{meta_comp}/items/*/model", within=config
+            ),
+            gene_fdr=str(config["enrichment"]["goatools"]["fdr_genes"]).replace(
+            ".", "-"
+            ),
+            go_term_fdr=str(config["enrichment"]["goatools"]["fdr_go_terms"]).replace(
+                ".", "-"
+            ),
+        ),
+    output:
+        "results/tables/go_terms/meta_compare_{meta_comp}.tsv",
+        "results/meta_comparison/go_terms/{meta_comp}.json",
+    log:
+        notebook="logs/meta_compare_enrichment/{meta_comp}.ipynb",
+    params:
+        desc=lookup(dpath="meta_comparisons/{meta_comp}/desc", within=config),
+        labels=lookup(dpath="meta_comparisons/{meta_comp}/items/*/label", within=config),
+    conda:
+        "../envs/pystats.yaml"
+    notebook:
+        "../scripts/compare_enrichment.py.ipynb"
+
+
+rule meta_compare_pathways:
+    input:
+        expand(
+            "results/tables/pathways/{model}.pathways.tsv",
+            model=lookup(
+                dpath="meta_comparisons/{meta_comp}/items/*/model", within=config
+            ),
+        ),
+    output:
+        "results/tables/pathways/meta_compare_{meta_comp}.tsv",
+        "results/meta_comparison/pathways/{meta_comp}.json",
+    log:
+        notebook="logs/meta_compare_pathways/{meta_comp}.ipynb",
+    params:
+        desc=lookup(dpath="meta_comparisons/{meta_comp}/desc", within=config),
+        labels=lookup(dpath="meta_comparisons/{meta_comp}/items/*/label", within=config),
+    conda:
+        "../envs/pystats.yaml"
+    notebook:
+        "../scripts/compare_pathways.py.ipynb"