Gimme maelstrom

CHANGELOG.md

@@ -10,6 +10,8 @@ All changed fall under either one of these types: `Added`, `Changed`, `Deprecate
 
 ### Added
 
+- Seq2science makes a biological replicates count table, which is the mean of the biological replicates.
+- Seq2science now supports differential motif analysis by gimme maelstrom!!!
 - configurable setting `niceness` which sets a niceness prefix to all shell commands.
 
 ### Fixed

docs/content/workflows/atac_seq.md

@@ -66,11 +66,18 @@ Seq2science currently supports four different normalisation methods: quantile no
 
 After these normalisations the counts are log normalised, and the base can be set with `logbase` and defaults to 2.
 As a final step the count tables are mean-centered.
-This final count table can be used for tools like [gimme maelstrom](https://gimmemotifs.readthedocs.io/en/master/reference.html#command-gimme-maelstrom) to scan for enriched transcription factor motifs.
 Note that this table contains **all** peaks, and no selection on differential peaks has been made.
 
+#### Differential motif analysis
+
+Seq2science supports differential motif analysis!
+This analysis is based on the tool [gimme maelstrom](https://gimmemotifs.readthedocs.io/en/master/reference.html#command-gimme-maelstrom).
+As input to gimme maelstrom is the count table of the biological reps (average between reps) in the consensus peakset.
+It then tries to solve a system of linear equations, where the output is the read counts in a peak, and the input is the motif score in the peak times the "motif activity".
+This motif activity can then be compared across biological replicates for differential motifs.
+
 #### Differential peak analysis
-Seq2science can optionally use the raw peak counts table to perform differential peak analysis.
+On top of that, Seq2science can also use the raw peak counts table to perform differential peak analysis.
 See the [Differential gene/peak analysis page](../DESeq2.html) for more information!
 
 #### Trackhub

docs/content/workflows/chip_seq.md

@@ -75,11 +75,18 @@ Seq2science currently supports four different normalisation methods: quantile no
 
 After these normalisations the counts are log normalised, and the base can be set with `logbase` and defaults to 2.
 As a final step the count tables are mean-centered.
-This final count table can be used for tools like [gimme maelstrom](https://gimmemotifs.readthedocs.io/en/master/reference.html#command-gimme-maelstrom) to scan for enriched transcription factor motifs.
 Note that this table contains **all** peaks, and no selection on differential peaks has been made.
 
+#### Differential motif analysis
+
+Seq2science supports differential motif analysis!
+This analysis is based on the tool [gimme maelstrom](https://gimmemotifs.readthedocs.io/en/master/reference.html#command-gimme-maelstrom).
+As input to gimme maelstrom is the count table of the biological reps (average between reps) in the consensus peakset.
+It then tries to solve a system of linear equations, where the output is the read counts in a peak, and the input is the motif score in the peak times the "motif activity".
+This motif activity can then be compared across biological replicates for differential motifs.
+
 #### Differential peak analysis
-Seq2science can optionally use the raw peak counts table to perform differential peak analysis.
+On top of that, Seq2science can also use the raw peak counts table to perform differential peak analysis.
 See the [Differential gene/peak analysis page](../DESeq2.html) for more information!
 
 #### Trackhub

seq2science/envs/gimme.yaml

@@ -4,7 +4,9 @@ channels:
   - bioconda
   - defaults
 dependencies:
-  - bioconda::gimmemotifs=0.15.1
-  - bioconda::biofluff=3.0.3
-  - conda-forge::pandas=1.1.5  # won't be necessary after https://github.com/vanheeringen-lab/gimmemotifs/issues/168 is fixed
+  - bioconda::gimmemotifs=0.17.1
+  - bioconda::biofluff=3.0.4
+  - bioconda::gffread=0.12.7
+  - bioconda::orthofinder=2.5.4
+  - conda-forge::cpulimit=0.2
   - conda-forge::conda-ecosystem-user-package-isolation=1.0

seq2science/rules/configuration_explain.smk

@@ -195,8 +195,9 @@ else:
         "fastp_PE": f"Paired-end reads were trimmed with {href_v('fastp')}{options('trimoptions')}.",
         "chipseeker": f"A peak feature distribution plot and peak localization plot relative to TSS were made with {hyperref('chipseeker')}.",  # TODO: replace with href_v
         "combine_peaks": f"A consensus set of summits was made with {href_v('gimmemotifs',text='gimmemotifs.combine_peaks',env='gimme')}.",
+        "gimme_maelstrom": f"Differential peaks analysis on the consensus peakset was performed with {href_v('gimmemotifs',text='gimme maelstrom',env='gimme')}.",
         "bed_slop": f"All summits were extended with {config.get('slop')} bp to get a consensus peakset.",
-        "coverage_table": f"Finally, a count table from the consensus peakset with gimmemotifs.",  # already cited in "combine_peaks"
+        "coverage_table": f"Finally, a count table from the consensus peakset was made with gimmemotifs.coverage_table.",  # already cited in "combine_peaks"
         "sce": f"{hyperref('sce')} S4 class was used to store scRNA-seq count tables and saved to RDATA format",
         "sctk": f"scRNA count post-processing was performed using the {hyperref('sctk')} toolkit.",
     }

seq2science/rules/configuration_generic.smk

@@ -199,7 +199,11 @@ def parse_samples():
         os._exit(1)  # noqa
 
     # remove unused columns, and populate empty cells in used columns
-    samples_df = dense_samples(samples_df, config)
+    samples_df = dense_samples(samples_df, 
+                               config.get("technical_replicates") == "keep", 
+                               config.get("biological_replicates") == "keep", 
+                               config.get("ignore_strandedness", True),
+                               config.get("create_trackhub", False))
 
     # check if replicate names are unique between assemblies
     if "technical_replicates" in samples_df:
@@ -335,15 +339,21 @@ if "assembly" in samples:
     def parse_assemblies():
         # list assemblies that are used in this workflow
         used_assemblies = list(set(samples["assembly"]))
-
+        # we make a temporary _used assemblies as gimme maelstrom might need more assemblies downloaded
+        # and those assemblies need to be saved in the local/remote assemblies variable
+        if "motif2factors_reference" in config and config["run_gimme_maelstrom"]:
+            _used_assemblies = used_assemblies + config["motif2factors_reference"] + config["motif2factors_database_references"]
+        else:
+            _used_assemblies = used_assemblies
+
         # dictionary with which providers to use per genome
         providers = PickleDict(os.path.join(CACHE_DIR, "providers.p"))
 
         # split into local and remote assemblies, depending on if all required files can be found
         annotation_required = "rna_seq" in WORKFLOW or config["aligner"] == "star"
         ftype = "annotation" if annotation_required else "genome"
-        local_assemblies = [a for a in used_assemblies if is_local(a, ftype, config)]
-        remote_assemblies = [a for a in used_assemblies if a not in local_assemblies]
+        local_assemblies = [a for a in _used_assemblies if is_local(a, ftype, config)]
+        remote_assemblies = [a for a in _used_assemblies if a not in local_assemblies]
         search_assemblies = [assembly for assembly in remote_assemblies if providers.get(assembly, {}).get(ftype) is None]
 
         # custom assemblies without provider (for local/remote annotations)
@@ -531,6 +541,9 @@ def any_given(*args, prefix="", suffix=""):
         if column_name == "sample":
             elements.extend(samples.index)
         elif column_name in samples:
+            if (column_name == "assembly") and ("motif2factors_reference" in config):
+                elements.extend(config["motif2factors_database_references"])
+                elements.extend(config["motif2factors_reference"])
             elements.extend(samples[column_name])
 
     elements = [prefix + element + suffix for element in elements if isinstance(element, str)]

seq2science/rules/configuration_workflows.smk

@@ -62,9 +62,15 @@ if config.get("peak_caller", False):
 
         if "--broad" in config["peak_caller"]["macs2"]:
             config["macs2_types"].extend(["peaks.broadPeak", "peaks.gappedPeak"])
+            config["run_gimme_maelstrom"] = False
         else:
             config["macs2_types"].extend(["summits.bed", "peaks.narrowPeak"])
 
+    if "run_gimme_maelstrom" in config and config["run_gimme_maelstrom"]:
+        assert len(breps.index) > 1, (
+            "To run gimme maelstrom you need more than one biological replicate!"
+        )
+
 # make sure that both maximum and minimum insert sizes are existing when one of them is used
 if config.get("min_template_length") and not config.get("max_template_length"):
     config["max_template_length"] = 1_000_000_000

seq2science/rules/deseq2.smk

@@ -20,7 +20,7 @@ def deseq_input(wildcards):
     elif "atac" in WORKFLOW or "chip" in WORKFLOW:
         # only uses a single peak caller ------------------------------------------v
         # TODO different peak callers can probably be supported with wildcard_constraint peak_caller (.*) <-- empty allowed
-        return (expand("{counts_dir}/{peak_caller}/{{assembly}}_raw.tsv", **config)[0],)
+        return (expand("{counts_dir}/{peak_caller}/{{assembly}}_raw_technical_reps.tsv", **config)[0],)
     else:
         logger.error(
             f"The workflow you are running ({WORKFLOW}) does not support deseq2. "

seq2science/rules/motif_scan.smk

@@ -0,0 +1,76 @@
+"""
+all rules/logic related to motif scanning of peaks is found here.
+"""
+def get_motif2factors_input_genomes(wildcards):
+    all_out = []
+    if str(wildcards.assembly)[0:4] in ["GRCh", "GRCm", "hg19", "hg38", "mm10", "mm39"]:
+        return all_out
+
+    for assembly in config.get("motif2factors_database_references", []) + config.get("motif2factors_reference", []) + [wildcards.assembly]:
+        all_out.append(expand(f"{{genome_dir}}/{assembly}/{assembly}.annotation.gtf", **config)[0])
+        all_out.append(expand(f"{{genome_dir}}/{assembly}/{assembly}.fa", **config)[0])
+    return all_out
+
+
+rule motif2factors:
+    """
+    Create a gimme pfm/motif2factor file with the gimme motif2factors command.
+    For human/mouse it just uses the default database, for other species it is based
+    on TF orthologs.
+    """
+    input:
+        genome=rules.get_genome.output,
+        all_genomes=get_motif2factors_input_genomes
+    output:
+        expand("{result_dir}/gimme/{{assembly}}.{{gimme_database}}.pfm", **config),
+    log:
+        expand("{log_dir}/gimme/motif2factors/{{assembly}}-{{gimme_database}}.log", **config),
+    benchmark:
+        expand("{benchmark_dir}/gimme/motif2factors/{{assembly}}-{{gimme_database}}.log", **config)[0],
+    params:
+        genomes_dir=config.get("genome_dir"),
+        database=config.get("gimme_maelstrom_database"),
+        motif2factors_reference=config.get("motif2factors_reference")
+    threads: 24
+    conda:
+        "../envs/gimme.yaml"
+    script:
+        f"{config['rule_dir']}/../scripts/motif2factors.py"
+
+
+rule gimme_maelstrom:
+    """
+    Gimme maelstrom is a method to infer differential motifs between two or more biological replicates.
+    """
+    input:
+        genome=rules.get_genome.output,
+        count_table=expand("{counts_dir}/{{peak_caller}}/{{assembly}}_log2_quantilenorm_biological_reps.tsv", **config),
+        pfm=rules.motif2factors.output
+    output:
+        output=directory(expand("{result_dir}/gimme/maelstrom/{{assembly}}-{{gimme_database}}-{{peak_caller}}", **config)),
+        cache=temp(directory(expand("{result_dir}/gimme/maelstrom/cache/{{assembly}}-{{gimme_database}}-{{peak_caller}}", **config))),
+    params: config.get("gimme_maelstrom_params", "")
+    log:
+        expand("{log_dir}/gimme_maelstrom/{{assembly}}-{{gimme_database}}-{{peak_caller}}.log", **config),
+    benchmark:
+        expand("{benchmark_dir}/gimme_maelstrom/{{assembly}}-{{gimme_database}}-{{peak_caller}}.log", **config)[0],
+    resources:
+        mem_gb=40,
+    message: EXPLAIN["gimme_maelstrom"]
+    conda:
+        "../envs/gimme.yaml"
+    threads: 24
+    shell:
+        """
+        NEW_CACHE={output.cache}
+        mkdir -p $NEW_CACHE
+        if [ -z $XDG_CACHE_HOME ]; then
+            XDG_CACHE_HOME=$HOME/.cache
+        fi
+        cp -r $XDG_CACHE_HOME/gimmemotifs $NEW_CACHE/
+        export XDG_CACHE_HOME=$NEW_CACHE
+        """ +
+        ("cpulimit --include-children -l {threads}00 --\\" if config.get("cpulimit", True) else "") +
+        """
+        gimme maelstrom {input.count_table} {input.genome} {output.output} --pfmfile {input.pfm} --nthreads {threads} {params} > {log} 2>&1
+        """

seq2science/rules/peak_counts.smk

@@ -14,28 +14,40 @@ def count_table_output():
     if ftype != "narrowPeak":
         return []
 
-    return expand(
+    expanddict = {
+        "assemblies": ALL_ASSEMBLIES,
+            "peak_caller": config["peak_caller"].keys(),
+            "normalization": [
+                "raw",
+                f"meancenter_log{config['logbase']}_quantilenorm",
+                f"meancenter_log{config['logbase']}_TMM",
+                f"meancenter_log{config['logbase']}_RLE",
+                f"meancenter_log{config['logbase']}_upperquartile",
+            ],
+            "mc": ["", "meancenter_"],
+        }
+
+    count_tables = expand(
         [
-            "{counts_dir}/{peak_caller}/{assemblies}_{normalization}.tsv",
+            "{counts_dir}/{peak_caller}/{assemblies}_{normalization}_technical_reps.tsv",
             "{counts_dir}/{peak_caller}/{assemblies}_onehotpeaks.tsv",
         ],
-        **{
-            **config,
-            **{
-                "assemblies": ALL_ASSEMBLIES,
-                "peak_caller": config["peak_caller"].keys(),
-                "normalization": [
-                    "raw",
-                    f"meancenter_log{config['logbase']}_quantilenorm",
-                    f"meancenter_log{config['logbase']}_TMM",
-                    f"meancenter_log{config['logbase']}_RLE",
-                    f"meancenter_log{config['logbase']}_upperquartile",
-                ],
-                "mc": ["", "meancenter_"],
-            },
-        },
+        **{**config, **expanddict}
     )
+    if breps is not treps:
+        expanddict["normalization"] = [
+                f"log{config['logbase']}_quantilenorm",
+                f"log{config['logbase']}_TMM",
+                f"log{config['logbase']}_RLE",
+        ]
+        count_tables.extend(
+            expand(
+                "{counts_dir}/{peak_caller}/{assemblies}_{normalization}_biological_reps.tsv",
+                **{**config, **expanddict}
+            )
+        )
 
+    return count_tables
 
 def get_peakfile_for_summit(wildcards):
     ftype = get_peak_ftype(wildcards.peak_caller)
@@ -167,7 +179,7 @@ rule coverage_table:
         replicates=get_coverage_table_replicates("bam"),
         replicate_bai=get_coverage_table_replicates("bam.bai"),
     output:
-        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_raw.tsv", **config),
+        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_raw_technical_reps.tsv", **config),
     log:
         expand("{log_dir}/coverage_table/{{assembly}}-{{peak_caller}}.log", **config),
     benchmark:
@@ -209,7 +221,7 @@ rule quantile_normalization:
     input:
         rules.coverage_table.output,
     output:
-        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_quantilenorm.tsv", **config),
+        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_quantilenorm_technical_reps.tsv", **config),
     log:
         expand("{log_dir}/quantile_normalization/{{assembly}}-{{peak_caller}}-quantilenorm.log", **config),
     benchmark:
@@ -234,7 +246,7 @@ rule edgeR_normalization:
     input:
         rules.coverage_table.output,
     output:
-        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_{{normalisation,(TMM|RLE|upperquartile)}}.tsv", **config),
+        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_{{normalisation,(TMM|RLE|upperquartile)}}_technical_reps.tsv", **config),
     log:
         expand("{log_dir}/edgeR_normalization/{{assembly}}-{{peak_caller}}-{{normalisation}}.log", **config),
     benchmark:
@@ -254,9 +266,9 @@ rule log_normalization:
     Log1p normalization of a count table.
     """
     input:
-        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_{{normalisation}}.tsv", **config),
+        rules.edgeR_normalization.output
     output:
-        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_log{{base}}_{{normalisation}}.tsv", **config),
+        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_log{{base}}_{{normalisation}}_technical_reps.tsv", **config),
     run:
         import pandas as pd
         import numpy as np
@@ -287,9 +299,9 @@ rule mean_center:
     Mean centering of a count table.
     """
     input:
-        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_log{{base}}_{{normalisation}}.tsv", **config),
+        rules.log_normalization.output
     output:
-        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_meancenter_log{{base}}_{{normalisation}}.tsv", **config),
+        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_meancenter_log{{base}}_{{normalisation}}_technical_reps.tsv", **config),
     run:
         import pandas as pd
 
@@ -311,6 +323,27 @@ rule mean_center:
         )
 
 
+rule combine_biological_reps:
+    """
+    Combine biological replicates by taking their mean!
+    """
+    input:
+        rules.log_normalization.output
+    output:
+        expand("{counts_dir}/{{peak_caller}}/{{assembly}}_log{{base}}_{{normalisation}}_biological_reps.tsv", **config),
+    log:
+        expand("{log_dir}/combine_biological_reps/{{peak_caller}}/{{assembly}}_log{{base}}_{{normalisation}}_biological_reps.log", **config),
+    benchmark:
+        expand("{benchmark_dir}/combine_biological_reps/{{peak_caller}}/{{assembly}}_log{{base}}_{{normalisation}}_biological_reps.benchmark.txt", **config)[0]
+    params:
+        samples=lambda wildcards: samples[samples["assembly"] == wildcards.assembly].replace(' ', '_', regex=True).to_string(index_names=False),
+        breps=lambda wildcards: breps[breps["assembly"] == wildcards.assembly].index.to_list(),
+        technical_replicates=config.get("technical_replicates") == "keep",
+        biological_replicates=config.get("biological_replicates") == "keep",
+    script:
+        f"{config['rule_dir']}/../scripts/combine_biological_reps.py"
+
+
 def get_all_narrowpeaks(wildcards):
     return [
         f"{config['result_dir']}/{{peak_caller}}/{{assembly}}-{replicate}_peaks.narrowPeak"

seq2science/rules/qc.smk

@@ -1036,4 +1036,12 @@ def get_peak_calling_qc(sample, wildcards):
             )
         )
 
+    if config["run_gimme_maelstrom"]:
+        output.extend(
+            expand(
+                "{qc_dir}/gimme/{{assembly}}-{gimme_maelstrom_database}-{peak_caller}_mqc.html",
+                **config
+            )
+        )
+
     return output

seq2science/rules/qc_peaks.smk

@@ -124,3 +124,17 @@ rule upset_plot_peaks:
         "../envs/upset.yaml"
     script:
         f"{config['rule_dir']}/../scripts/upset.py"
+
+
+rule maelstrom_report_preparation:
+    """
+    This rule injects the symlinked motif logos into the html so they get rendered in the multiqc report.
+    """
+    input:
+        rules.gimme_maelstrom.output
+    output:
+        expand("{qc_dir}/gimme/{{assembly}}-{{gimme_database}}-{{peak_caller}}_mqc.html", **config)
+    log:
+        expand("{log_dir}/maelstrom_report_preparation/{{assembly}}-{{gimme_database}}-{{peak_caller}}.log", **config),
+    script:
+        f"{config['rule_dir']}/../scripts/maelstrom_report.py"