feat!: allow for multiple gene_lists

.test/config/config.yaml

@@ -44,7 +44,9 @@ diffexp:
     qq-plot: 0.05
   genes_of_interest:
     activate: true
-    genelist: "resources/gene_list.tsv"
+    gene_lists:
+      gene_list_1: "resources/gene_list.tsv"
+      gene_list_2: "resources/gene_list.tsv"
 
 diffsplice:
   activate: false

.test/resources/gene_list.tsv

@@ -1,18 +1 @@
-STAT1
-IRF1
-HLA-A
-HLA-DRB1
-TYR
-PMEL
-DCT
-MLANA
-MITF
-CDK2
-SOX10
-ERBB3
-LEF1
-CTNNB1
-CDH1
-FN1
-NGFR
-AXL
+STAT1	IRF1	HLA-A	HLA-DRB1	TYR	PMEL	DCT	MLANA	MITF	CDK2	SOX10

.test/resources/gene_list_2.tsv

@@ -0,0 +1 @@
+MLANA	MITF	CDK2	SOX10	ERBB3	LEF1	CTNNB1	CDH1	FN1	NGFR	AXL	

.test/three_prime/config/config.yaml

@@ -44,7 +44,8 @@ diffexp:
     qq-plot: 0.05
   genes_of_interest:
     activate: false
-    genelist: "resources/gene_list.tsv"
+    gene_lists:
+      gene_list_1: "resources/gene_list.tsv"
 
 diffsplice:
   activate: false

config/config.yaml

@@ -136,13 +136,19 @@ diffexp:
     ma-plot: 0.05
     qq-plot: 0.05
   # heatmap and bootstrap plots for given set of genes:
-  # If you want a heatmap and bootstrap plots for a particular set of genes,
-  # you can set `activate: true` and provide a genelist file. In this file,
-  # list all your HGNC gene symbols of interest in one line, separated by
-  # tabulators (tabs).
+  # If you want a heatmap and bootstrap plots for particular sets of genes,
+  # you can set `activate: true` and provide one or more gene list files. In
+  # those files, list all your HGNC gene symbols of interest in one line,
+  # separated by tabulators (tabs). The workflow will generate one heatmap
+  # plot per file, and a bootstrap plot for each of the genes contained in
+  # any of the provided files.
   genes_of_interest:
     activate: false
-    genelist: "config/gene_list.tsv"
+    gene_lists:
+      # use a descriptive gene list name, as this will be used as an infix in
+      # the filename of the resulting heatmap plot
+      gene_list_1: "config/gene_list.tsv"
+
 
 diffsplice:
   # isoformSwitchAnalyzer

workflow/rules/common.smk

@@ -242,15 +242,6 @@ def kallisto_params(wildcards, input):
     return extra
 
 
-def input_genelist(predef_genelist):
-    if config["diffexp"]["genes_of_interest"]["activate"] == True:
-        predef_genelist = config["diffexp"]["genes_of_interest"]["genelist"]
-    else:
-        predef_genelist = []
-
-    return predef_genelist
-
-
 def all_input(wildcards):
     """
     Function defining all requested inputs for the rule all (below).
@@ -339,9 +330,10 @@ def all_input(wildcards):
         wanted_input.extend(
             expand(
                 [
-                    "results/plots/diffexp-heatmap/{model}.diffexp-heatmap.{mode}.pdf",
+                    "results/plots/diffexp-heatmap/{model}.{gene_list}.diffexp-heatmap.{mode}.pdf",
                 ],
                 model=config["diffexp"]["models"],
+                gene_list=lookup(within=config, dpath("diffexp/gene_of_interest/gene_lists")),
                 mode=["predefined"],
             )
         )

workflow/rules/diffexp.smk

@@ -146,7 +146,7 @@ rule plot_bootstrap:
         color_by=config["bootstrap_plots"]["color_by"],
         fdr=config["bootstrap_plots"]["FDR"],
         top_n=config["bootstrap_plots"]["top_n"],
-        genes=config["diffexp"]["genes_of_interest"],
+        gene_lists=lookup(within=config, dpath="diffexp/genes_of_interest/gene_lists"),
     log:
         "logs/plots/bootstrap/{model}/{model}.plot_bootstrap.log",
     script:
@@ -258,18 +258,21 @@ rule tpm_matrix:
 rule plot_diffexp_heatmap:
     input:
         logcountmatrix_file="results/tables/logcount-matrix/{model}.logcount-matrix.tsv",
-        predef_genelist=input_genelist,
+        predef_gene_list=lookup(within=config, dpath="diffexp/genes_of_interest/gene_lists/{gene_list}"),
     output:
         diffexp_heatmap=report(
-            "results/plots/diffexp-heatmap/{model}.diffexp-heatmap.{mode}.pdf",
+            "results/plots/diffexp-heatmap/{model}.{gene_list}.diffexp-heatmap.{mode}.pdf",
             caption="../report/plot-heatmap.rst",
             category="Heatmaps",
-            labels={"model": "{model}-{mode}"},
+            labels={
+                "model": "{model}-{mode}",
+                "gene list": "{gene_list}",
+            },
         ),
     params:
         model=get_model,
     log:
-        "logs/plots/diffexp-heatmap/{model}.diffexp-heatmap.{mode}.log",
+        "logs/plots/diffexp-heatmap/{model}.{gene_list}.diffexp-heatmap.{mode}.log",
     conda:
         "../envs/heatmap.yaml"
     script:

workflow/schemas/config.schema.yaml

@@ -100,10 +100,13 @@ properties:
         properties:
           activate: 
             type: boolean
-          genelist:
-            type: string
+          gene_lists:
+            type: object
+            patternProperties:
+              "^.+$":
+                type: string
         required:
-          - genelist
+          - gene_lists
     required:
       - models
       - genes_of_interest

workflow/scripts/plot-bootstrap.R

@@ -20,11 +20,19 @@ top_genes <- results %>%
 
 
 if (snakemake@params[["genes"]]$activate == TRUE) {
-    gene_table <- read.csv(snakemake@params[["genes"]]$genelist,
-        sep = "\t")
-    names(gene_table) <- c("ext_gene")
-    genes_of_interest <- tibble(gene_table) %>%
-        distinct(ext_gene)
+  genes = c()
+  for (gene_list in unname(snakemake@params[["gene_lists"]])) {
+    new_genes <- read_tsv(
+      gene_list,
+      col_names = FALSE
+    )
+    genes <- union(genes, new_genes |> unlist() |> unname())
+  }
+  genes_of_interest <- enframe(
+    genes,
+    name = NULL,
+    value = "ext_gene"
+  )
 } else {
     # "genes" is null, if the list provided in config.yaml is empty
     genes_of_interest <- tibble(ext_gene = character())

workflow/scripts/plot_diffexp_heatmap.R

@@ -30,7 +30,7 @@ if (snakemake@wildcards[["mode"]] == "topn") {
 } else if (snakemake@wildcards[["mode"]] == "predefined") {
     # Adding gene list to the variable
     predefine_genelist <-
-        read.table(snakemake@input[["predef_genelist"]],
+        read.table(snakemake@input[["predef_gene_list"]],
             sep = "\t"
         )
     selectedgenes <-