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

@@ -32,6 +32,7 @@ scatter:
 
 diffexp:
   exclude:
+    - SRR8309099_2
   models:
     model_X:
       full: ~condition
@@ -44,7 +45,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

.test/three_prime/config/samples.tsv

@@ -4,4 +4,5 @@ SRR8309094	Control
 SRR8309095	Treated
 SRR8309097	Treated
 SRR8309098	Control
-SRR8309099	Treated
+SRR8309099	Treated
+SRR8309099_2	Treated

.test/three_prime/config/units.tsv

@@ -4,4 +4,5 @@ SRR8309094	u1	430	43	quant_seq_test_data/SRR8309094.fastq.gz
 SRR8309095	u1	430	43	quant_seq_test_data/SRR8309095.fastq.gz					
 SRR8309097	u1	430	43	quant_seq_test_data/SRR8309097.fastq.gz					
 SRR8309098	u1	430	43	quant_seq_test_data/SRR8309098.fastq.gz					
-SRR8309099	u1	430	43	quant_seq_test_data/SRR8309099.fastq.gz					
+SRR8309099	u1	430	43	quant_seq_test_data/SRR8309099.fastq.gz					
+SRR8309099_2	u1	430	43	quant_seq_test_data/SRR8309099.fastq.gz					

config/config.yaml

@@ -31,11 +31,14 @@ resources:
     species: homo_sapiens
     # ensembl release version:
     # Update this to the latest working version, when you first set up a new
-    # analysis on a dataset. Later, it only makes sense to update (or downgrade)
-    # the release versions if either (a) the version you are using consistently
-    # fails to download (some Ensembl release versions are just broken) or
-    # (b) you know that a newer version will include changes that will fix some
-    # error or adds transcripts that will be relevant to your analysis.
+    # analysis on a dataset. You can usually find the latest release in the
+    # Ensembl blog, by looking at the latest posts of the release category:
+    # https://www.ensembl.info/category/01-release/
+    # Later, it only makes sense to update (or downgrade) the release versions
+    # if either (a) the version you are using consistently fails to download
+    # (some Ensembl release versions are just broken) or (b) you know that a
+    # newer version will include changes that will fix some error or adds
+    # transcripts that will be relevant to your analysis.
     release: "114"
     # genome build:
     # Usually, this should just be the main build listed in:
@@ -91,8 +94,11 @@ scatter:
 
 diffexp:
   # Samples to exclude from differential expression modeling (for example,
-  # outliers due to technical problems).
+  # outliers due to technical problems). List their `sample_name` column
+  # entry in the form of a YAML list.
   exclude:
+  # - sample_X
+  # - sample_Y
   # model for sleuth differential expression analysis
   # For an introduction to sleuth, see its online manual:
   # https://pachterlab.github.io/sleuth/about
@@ -136,13 +142,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).
@@ -329,20 +320,22 @@ def all_input(wildcards):
                 "results/tables/tpm-matrix/{model}.tpm-matrix.tsv",
                 "results/sleuth/{model}.samples.tsv",
                 "results/datavzrd-reports/diffexp-{model}",
-                "results/plots/diffexp-heatmap/{model}.diffexp-heatmap.{mode}.pdf",
+                "results/plots/diffexp-heatmap/{model}.diffexp-heatmap.{gene_list}.pdf",
             ],
             model=config["diffexp"]["models"],
-            mode=["topn"],
+            gene_list=["topn"],
         )
     )
     if config["diffexp"]["genes_of_interest"]["activate"]:
         wanted_input.extend(
             expand(
                 [
-                    "results/plots/diffexp-heatmap/{model}.diffexp-heatmap.{mode}.pdf",
+                    "results/plots/diffexp-heatmap/{model}.diffexp-heatmap.{gene_list}.pdf",
                 ],
                 model=config["diffexp"]["models"],
-                mode=["predefined"],
+                gene_list=lookup(
+                    within=config, dpath="diffexp/genes_of_interest/gene_lists"
+                ),
             )
         )
 

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"],
+        genes_of_interest=lookup(within=config, dpath="diffexp/genes_of_interest"),
     log:
         "logs/plots/bootstrap/{model}/{model}.plot_bootstrap.log",
     script:
@@ -258,18 +258,27 @@ rule tpm_matrix:
 rule plot_diffexp_heatmap:
     input:
         logcountmatrix_file="results/tables/logcount-matrix/{model}.logcount-matrix.tsv",
-        predef_genelist=input_genelist,
+        # default provides a dummy file path for the "topn" case (a file path
+        # that is always valid, but never loaded)
+        predef_gene_list=lookup(
+            within=config,
+            dpath="diffexp/genes_of_interest/gene_lists/{gene_list}",
+            default=lookup(within=config, dpath="samples"),
+        ),
     output:
         diffexp_heatmap=report(
-            "results/plots/diffexp-heatmap/{model}.diffexp-heatmap.{mode}.pdf",
+            "results/plots/diffexp-heatmap/{model}.diffexp-heatmap.{gene_list}.pdf",
             caption="../report/plot-heatmap.rst",
             category="Heatmaps",
-            labels={"model": "{model}-{mode}"},
+            labels={
+                "model": "{model}",
+                "gene list": "{gene_list}",
+            },
         ),
     params:
         model=get_model,
     log:
-        "logs/plots/diffexp-heatmap/{model}.diffexp-heatmap.{mode}.log",
+        "logs/plots/diffexp-heatmap/{model}.diffexp-heatmap.{gene_list}.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

@@ -19,15 +19,23 @@ top_genes <- results %>%
    distinct(ext_gene)
 
 
-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)
+if (snakemake@params[["genes_of_interest"]][["activate"]] == TRUE) {
+  genes = character()
+  for (gene_list in unname(snakemake@params[["genes_of_interest"]][["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())
+  # Create empty tibble when genes_of_interest is not activated
+  genes_of_interest <- tibble(ext_gene = character())
 }
 genes <- full_join(top_genes, genes_of_interest, by = "ext_gene") %>%
     add_row(ext_gene = "Custom") %>%

workflow/scripts/plot_diffexp_heatmap.R

@@ -12,7 +12,7 @@ sleuth_file <- read.csv(snakemake@input[["logcountmatrix_file"]],
 )
 # Check the config file if it contains pre-defined gene list
 
-if (snakemake@wildcards[["mode"]] == "topn") {
+if (snakemake@wildcards[["gene_list"]] == "topn") {
     # Replacing the rownames with transcript Id and gene name
     rownames(sleuth_file) <- paste(
         sleuth_file$transcript,
@@ -27,10 +27,10 @@ if (snakemake@wildcards[["mode"]] == "topn") {
     selectedgenes <- sleuth_file[row.names(sleuth_file)
     %in% names(vargenes[order(vargenes, decreasing = TRUE)][1:50]), ]
     # If the config file contains pre-defined gene list
-} else if (snakemake@wildcards[["mode"]] == "predefined") {
+} else {
     # Adding gene list to the variable
     predefine_genelist <-
-        read.table(snakemake@input[["predef_genelist"]],
+        read.table(snakemake@input[["predef_gene_list"]],
             sep = "\t"
         )
     selectedgenes <-