Merge pull request #217 from ctmrbio/develop

boulund · web-flow · commit 39f1cd0de4b4 · 2023-06-01T13:21:45.000+02:00
StaG v6.1
diff --git a/CHANGELOG.md b/CHANGELOG.md
@@ -14,16 +14,39 @@ committed to the master branch that does not trigger any of the aforementioned
 situations.
 
 
-## [0.6.1] Unreleased
+## [0.6.1] 2023-06-01
 ### Added
+- BBMap: now outputs sorted BAM file, added options `keep_sam` and `keep_bam`.
+- Bowtie2: added option `keep_bam`.
 
 ### Fixed
+- KrakenUniq: environment variable `LC_ALL` has been added to Singularity image
+  to prevent unnecessary warning messages related to it being undefined.
+- KrakenUniq: now able to run when host removal is skipped, solved by adding
+  `krakenuniq_merge_reads` rule to create a temporary merged fasta file with
+  input data for KrakenUniq to avoid giving KrakenUniq symlinks as input.
 
 ### Changed
+- MetaPhlAn: Updated to v4.0.6
+- HUMAnN3: Updated to v3.7
+- HUMAnN3: Changed the way the temporary directory is resolved, now using
+  Snakemake's built-in `resources.tmpdir`. This should prevent HUMAnN from
+  creating large temporary directories outside of Slurm job folders so that
+  they cannot be automatically cleaned up if the Slurm job times out or fails
+  before HUMAnN can clean up after itself.
+- KrakenUniq: Concatenate reads with BBMap's `fuse.sh` with a padding of one
+  `N` instead of interleaving the paired inputs into a single FASTA to avoid
+  KrakenUniq treating paired reads independently.
 
 ### Deprecated
+- Kaiju, Kraken2, MetaPhlAn: area plot removed due to repeatedly leading to
+  failed runs in cached Singularity containers. The script still works as
+  intended in newer matplotlib versions and will remain in the scripts folder
+  for potential manual use if desired.
 
 ### Removed
+- Groot: Removed settings related to read length window as that feature was
+  removed in a previous StaG release.
 
 
 ## [0.6.0] 2023-04-17
diff --git a/config/config.yaml b/config/config.yaml
@@ -17,8 +17,7 @@ input_fn_pattern: "{sample}_{readpair}.fq.gz"
 samplesheet: ""           # Three-column samplesheet with sample_id,fastq_1,fastq_2 columns. Used instead of inputdir
 outdir: "output_dir"
 logdir: "output_dir/logs"
-tmpdir: "/tmp"            # Specifying a tmpdir is not normally necessary but some tools like HUMAnN requires it. If required please enter path to a temp directory e.g. /scratch
-dbdir: "databases"        # Databases will be downloaded to this dir, if requested
+dbdir: "databases"        # Databases will be downloaded to this dir
 report: "StaG_report-"    # Filename prefix for report file ("-{datetime}.html" automatically appended)
 email: ""                 # Email to send status message after completed/failed run.
 
@@ -161,8 +160,6 @@ humann:
 #########################
 groot:
     index_dir: ""            # [Required] Path to groot indexDir
-    minlength: 110           # Minlength for groot index
-    maxlength: 125           # Maxlength for groot index
     covcutoff: 0.97          # Coverage cutoff for groot report
     lowcov: False            # Report ARGs with no 5' or 3' coverage. Overrides covcutoff.
 
@@ -185,6 +182,8 @@ bbmap:
     - db_name: ""              # [Required] Custom name for BBMap database
       db_path: ""              # [Required] Path to BBMap database (folder should contain a 'ref' folder)
       min_id: 0.76             # Minimum id for read alignment, BBMap default is 0.76
+      keep_sam: False          # Set to True to keep intermediary SAM file
+      keep_bam: True           # Set to False to remove bam files after counting annotations
       extra: ""                # Extra BBMap command line parameters
       counts_table:
           annotations: ""      # Tab-separated annotation file with headers, first column is full FASTA header of reference sequences
@@ -196,6 +195,7 @@ bbmap:
           extra: ""            # Extra featureCount command line parameters
 bowtie2:
     - db_prefix: ""            # [Required] Full path to Bowtie2 index (not including file extension)
+      keep_bam: True           # Set to False to remove bam files after counting annotations
       extra: ""                # Extra bowtie2 commandline parameters
       counts_table:
           annotations: ""      # Tab-separated annotation file with headers, first column is full FASTA header of reference sequences
diff --git a/docs/source/conf.py b/docs/source/conf.py
@@ -56,9 +56,9 @@
 # built documents.
 #
 # The short X.Y version.
-version = '0.6.0'
+version = '0.6.1'
 # The full version, including alpha/beta/rc tags.
-release = '0.6.0-dev'
+release = '0.6.1'
 
 # reStructuredText prolog contains a string of reStructuredText that will be
 # included at the beginning of every source file that is read.
diff --git a/docs/source/modules.rst b/docs/source/modules.rst
@@ -12,7 +12,7 @@
 .. _featureCounts: https://subread.sourceforge.net/featureCounts.html
 .. _HUMAnN: https://github.com/biobakery/biobakery/wiki/humann3
 .. _GTF format: https://genome.ucsc.edu/FAQ/FAQformat.html#format4
-.. _SAF format: http://bioinf.wehi.edu.au/featureCounts/
+.. _SAF format: https://subread.sourceforge.net/featureCounts.html
 .. _MEGAHIT: https://github.com/voutcn/megahit
 .. _MultiQC: https://multiqc.info/
 .. _amrplusplus: https://megares.meglab.org/amrplusplus/latest/html/what_AMR++_produces.html
@@ -259,13 +259,13 @@ the default unit is counts per million (cpm).
 
    HUMAnN uses the taxonomic profiles produced by MetaPhlAn as input,
    so all MetaPhlAn-associated steps are run regardless of whether it is actually
-   enabled in ``config.yaml`` or not. It is important to use a MetaPhlAn database
-   compatible with HUMAnN3, e.g. mpa_v30_CHOCOPhlAn_201901 (run the metaphlan
-   step with the extra ``--mpa3`` flag in the StaG config file).
+   enabled in ``config.yaml`` or not.
 
-Due to temporary disk space issues with running HUMAnN it is now a requirement
-to specify a $TMPDIR in ``config.yaml``, e.g. ``/scratch`` or ``/tmp`` depending 
-on your system's configuration.
+HUMAnN requires large amounts of temporary disk space when processing a sample
+and will automatically use a suitable temporary directory from system
+environment variable ``$TMPDIR``, using Snakemake's resources feature to
+evaluate the variable at runtime (which means it can utilize node-local
+temporary disk if executing on a compute cluster).
 
 
 Antibiotic resistance
@@ -296,9 +296,8 @@ mapped graphs of all detected antibiotic resistance genes.
    The feature may reintroduced in future versions of GROOT but is not
    available in StaG now.
 
-The read lengths input to `groot`_ must conform to the settings used during
-`groot`_ database construction. The length window can be configured in the
-config file.
+The read lengths used with `groot`_ should preferably conform to the settings
+used during `groot`_ database construction.
 
 AMRPlusPlus_v2
 -------
@@ -353,10 +352,11 @@ BBMap
 :Tool: `BBMap`_
 :Output folder: ``bbmap/<database_name>``
 
-This module maps read using `BBMap`_. The output is in gzipped SAM format. It
-is possible to configure the mapping settings almost entirely according to
-preference, with the exception of changing the output format from gzipped SAM.
-Use the configuration parameter ``bbmap:extra`` to add any standard BBMap
+This module maps read using `BBMap`_. The output is in sorted and indexed BAM
+format (with an option to keep the intermediary SAM file used to create the
+BAM). It is possible to configure the mapping settings almost entirely
+according to preference, with the exception of changing the output format. Use
+the configuration parameter ``bbmap:extra`` to add any standard BBMap
 commandline parameter you want.
 
 Bowtie2
@@ -381,6 +381,8 @@ signifies a list)::
         - db_name: ""
           db_path: ""
           min_id: 0.76
+          keep_sam: False
+          keep_bam: True
           extra: ""
           counts_table:
               annotations: ""
@@ -400,6 +402,8 @@ configuration options, but with different settings. For example, to map against
         - db_name: "db1"
           db_path: "/path/to/db1"
           min_id: 0.76
+          keep_sam: False
+          keep_bam: True
           extra: ""
           counts_table:
               annotations: ""
@@ -412,6 +416,8 @@ configuration options, but with different settings. For example, to map against
         - db_name: "db2"
           db_path: "/path/to/db2"
           min_id: 0.76
+          keep_sam: False
+          keep_bam: True
           extra: ""
           counts_table:
               annotations: "/path/to/db2/annotations.txt"
diff --git a/workflow/Snakefile b/workflow/Snakefile
@@ -21,7 +21,7 @@ from scripts.common import UserMessages, SampleSheet
 
 user_messages = UserMessages()
 
-stag_version = "0.6.0"
+stag_version = "0.6.1"
 singularity_branch_tag = "-master"  # Replace with "-master" before publishing new version
 
 configfile: "config/config.yaml"
@@ -31,7 +31,6 @@ citations = {publications["StaG"], publications["Snakemake"]}
 INPUTDIR = Path(config["inputdir"])
 OUTDIR = Path(config["outdir"])
 LOGDIR = Path(config["logdir"])
-TMPDIR = Path(config["tmpdir"])
 DBDIR = Path(config["dbdir"])
 all_outputs = []
 
diff --git a/workflow/envs/Singularity.biobakery b/workflow/envs/Singularity.biobakery
@@ -11,7 +11,7 @@ From: mambaorg/micromamba:0.17.0
 		-c bioconda \
 		-c conda-forge \
 		-c biobakery \
-		python=3.7 metaphlan=4.0.3 humann=3.6 krona=2.8.1
+		python=3.7 metaphlan=4.0.6 humann=3.7 krona=2.8.1
 	micromamba clean --yes --all
 
 	humann_databases --download utility_mapping full /opt/humann --update-config yes
diff --git a/workflow/envs/Singularity.krakenuniq b/workflow/envs/Singularity.krakenuniq
@@ -8,17 +8,20 @@ From: debian:11
     AUTHOR boulund
     VERSION 1.0
 
+%environment
+    export LC_ALL="C"
+
 %post
     apt-get update && \
-	apt-get install -y \
+    apt-get install -y \
       bash \
       perl \
       make \
       g++ \
       libbz2-dev \
       zlib1g-dev \
-	  file \
-	  wget
+      file \
+      wget
     
     mkdir -pv /opt/krakenuniq /opt/krakenuniq-src
     
@@ -28,9 +31,9 @@ From: debian:11
     tar -xf v1.0.3.tar.gz
     cd krakenuniq-1.0.3
 
-	./install_krakenuniq.sh -l /usr/local/bin /opt/krakenuniq
+    ./install_krakenuniq.sh -l /usr/local/bin /opt/krakenuniq
 
-	rm -fv /opt/krakenuniq-src/v1.0.3.tar.gz
+    rm -fv /opt/krakenuniq-src/v1.0.3.tar.gz
 
 
 %runscript
diff --git a/workflow/envs/humann.yaml b/workflow/envs/humann.yaml
@@ -5,4 +5,4 @@ channels:
     - bioconda
     - defaults
 dependencies:
-    - humann=3.6
+    - humann=3.7
diff --git a/workflow/envs/metaphlan.yaml b/workflow/envs/metaphlan.yaml
@@ -3,7 +3,6 @@ channels:
     - conda-forge
     - bioconda
     - defaults
-    - biobakery
 dependencies:
-    - metaphlan =4.0.3
+    - metaphlan =4.0.6
     - krona =2.8.1
diff --git a/workflow/rules/antibiotic_resistance/groot.smk b/workflow/rules/antibiotic_resistance/groot.smk
@@ -45,8 +45,6 @@ rule groot_align:
     threads: 8
     params:
         index_dir=groot_config["index_dir"],
-        minlength=groot_config["minlength"],
-        maxlength=groot_config["maxlength"],
     shell:
         """
         groot align \
diff --git a/workflow/rules/functional_profiling/humann.smk b/workflow/rules/functional_profiling/humann.smk
@@ -21,10 +21,6 @@ if config["functional_profile"]["humann"]:
         err_message += "If you do not want to run HUMAnN for functional profiling, set functional_profile:humann: False in config.yaml"
         err_message += "If you want to download the HUMAnN databases please see https://github.com/biobakery/humann"
         raise WorkflowError(err_message)
-    if not Path(config["tmpdir"]).exists():
-        err_message = "Please specify a tmpdir in config.yaml, if specified tmpdir does not exist, create it.\n"
-        err_message += "Do not specify a subdirectory of \"/scratch\" as tmpdir, only specifying \"/scratch\" is sufficient."
-        raise WorkflowError(err_message)
 
     merged_humann_tables = expand(f"{OUTDIR}/humann/all_samples.humann_{{output_type}}.txt",
             output_type=("genefamilies", "pathcoverage", "pathabundance"))
@@ -56,19 +52,18 @@ rule humann:
     threads: 20
     params:
         outdir=f"{OUTDIR}/humann/",
-        tmpdir=f"{TMPDIR}/{{sample}}",
         nucleotide_db=h_config["nucleotide_db"],
         protein_db=h_config["protein_db"],
         extra=h_config["extra"],
     shell:
         """
-        mkdir -pv {params.tmpdir} >> {log.stdout}
+        mkdir -pv {resources.tmpdir} >> {log.stdout}
         mkdir -pv {params.outdir} >> {log.stdout}
-        cat {input.read1} {input.read2} > {params.tmpdir}/concat_input_reads.fq.gz
+        cat {input.read1} {input.read2} > {resources.tmpdir}/concat_input_reads.fq.gz
 
         humann \
-            --input {params.tmpdir}/concat_input_reads.fq.gz \
-            --output {params.tmpdir} \
+            --input {resources.tmpdir}/concat_input_reads.fq.gz \
+            --output {resources.tmpdir} \
             --nucleotide-database {params.nucleotide_db} \
             --protein-database {params.protein_db} \
             --output-basename {wildcards.sample} \
@@ -78,9 +73,9 @@ rule humann:
             >> {log.stdout} \
             2> {log.stderr}
 
-        cp {params.tmpdir}/{wildcards.sample}*.tsv {params.outdir}
-        cp {params.tmpdir}/{wildcards.sample}_humann_temp/{wildcards.sample}.log {log.log}
-        rm -rfv {params.tmpdir} >> {log.stdout}
+        cp {resources.tmpdir}/{wildcards.sample}*.tsv {params.outdir}
+        cp {resources.tmpdir}/{wildcards.sample}_humann_temp/{wildcards.sample}.log {log.log}
+        rm -rfv {resources.tmpdir} >> {log.stdout}
         """
 
 
diff --git a/workflow/rules/mappers/bbmap.smk b/workflow/rules/mappers/bbmap.smk
@@ -23,7 +23,7 @@ for bbmap_config in config["bbmap"]:
         bbmap_alignments = expand(str(OUTDIR/"bbmap/{db_name}/{sample}.{output_type}"),
                 db_name=db_name,
                 sample=SAMPLES,
-                output_type=("sam.gz", "covstats.txt", "rpkm.txt"))
+                output_type=("bam", "covstats.txt", "rpkm.txt"))
         counts_table = expand(str(OUTDIR/"bbmap/{db_name}/counts.{column}.tsv"),
                 db_name=db_name,
                 column=map(str.strip, bbmap_config["counts_table"]["columns"].split(",")))
@@ -58,9 +58,11 @@ for bbmap_config in config["bbmap"]:
             read1=OUTDIR/"host_removal/{sample}_1.fq.gz",
             read2=OUTDIR/"host_removal/{sample}_2.fq.gz",
         output:
-            sam=bbmap_output_folder/"{sample}.sam.gz",
+            sam=bbmap_output_folder/"{sample}.sam.gz" if bbmap_config["keep_sam"] else temp(bbmap_output_folder/"{sample}.sam.gz"),
             covstats=bbmap_output_folder/"{sample}.covstats.txt",
             rpkm=bbmap_output_folder/"{sample}.rpkm.txt",
+            bamscript=temp(bbmap_output_folder/"{sample}.bamscript.sh"),
+            bamfile=bbmap_output_folder/"{sample}.bam" if bbmap_config["keep_bam"] else temp(bbmap_output_folder/"{sample}.bam"),
         log:
             stdout=str(bbmap_logdir/"{sample}.bbmap.stdout.log"),
             stderr=str(bbmap_logdir/"{sample}.bbmap.statsfile.txt"),
@@ -88,9 +90,14 @@ for bbmap_config in config["bbmap"]:
                 out={output.sam} \
                 covstats={output.covstats} \
                 rpkm={output.rpkm} \
+                bamscript={output.bamscript} \
                 {params.extra} \
                 > {log.stdout} \
                 2> {log.stderr}
+
+            sed -i 's/_sorted//g' {output.bamscript}
+
+            ./{output.bamscript} 2>> {log.stderr} >> {log.stdout}
             """
 
 
@@ -140,7 +147,7 @@ for bbmap_config in config["bbmap"]:
         f"""Summarize feature counts for {db_name}"""
         name: f"bbmap_feature_counts_{db_name}"
         input:
-            bams=expand(str(OUTDIR/"bbmap/{db_name}/{sample}.sam.gz"),
+            bams=expand(str(OUTDIR/"bbmap/{db_name}/{sample}.bam"),
                     db_name=db_name,
                     sample=SAMPLES)
         output:
diff --git a/workflow/rules/mappers/bowtie2.smk b/workflow/rules/mappers/bowtie2.smk
@@ -59,7 +59,7 @@ for bt2_config in config["bowtie2"]:
             sample=[OUTDIR/"host_removal/{sample}_1.fq.gz",
                     OUTDIR/"host_removal/{sample}_2.fq.gz"]
         output:
-            OUTDIR/"bowtie2/{db_name}/{{sample}}.bam".format(db_name=bt2_db_name)
+            OUTDIR/"bowtie2/{db_name}/{{sample}}.bam".format(db_name=bt2_db_name) if bt2_config["keep_bam"] else temp(OUTDIR/"bowtie2/{db_name}/{{sample}}.bam".format(db_name=bt2_db_name)),
         log:
             str(LOGDIR/"bowtie2/{db_name}/{{sample}}.log".format(db_name=bt2_db_name))
         message:
diff --git a/workflow/rules/taxonomic_profiling/kaiju.smk b/workflow/rules/taxonomic_profiling/kaiju.smk
@@ -33,7 +33,7 @@ if config["taxonomic_profile"]["kaiju"]:
     all_outputs.extend(kaiju)
     all_outputs.extend(kaiju_reports)
     all_outputs.append(kaiju_joined_table)
-    #all_outputs.append(kaiju_area_plot)  # Buggy in stag 4.1
+    #all_outputs.append(kaiju_area_plot)  # Buggy in stag 4.1, deprecated in 6.1
 
     citations.add(publications["Kaiju"])
 
diff --git a/workflow/rules/taxonomic_profiling/kraken2.smk b/workflow/rules/taxonomic_profiling/kraken2.smk
@@ -48,11 +48,11 @@ if config["taxonomic_profile"]["kraken2"]:
     all_outputs.extend(kreports_mpa_style)
     all_outputs.append(joined_kreport_mpa_style)
     all_outputs.append(combined_kreport)
+    #all_outputs.append(kraken_area_plot)  # Deprecated in 6.1
     citations.add(publications["Kraken2"])
 
     if kraken2_config["run_krona"]:
         all_outputs.append(kraken_krona)
-        all_outputs.append(kraken_area_plot)
         citations.add(publications["Krona"])
     
 
@@ -265,7 +265,7 @@ if config["taxonomic_profile"]["kraken2"] and kraken2_config["bracken"]["kmer_di
     all_outputs.extend(brackens_mpa_style)
     all_outputs.append(bracken_krona)
     all_outputs.append(all_table_mpa)
-    all_outputs.append(bracken_area_plot)
+    #all_outputs.append(bracken_area_plot)  # Deprecated in 6.1
 
 
 rule bracken_kreport:
diff --git a/workflow/rules/taxonomic_profiling/krakenuniq.smk b/workflow/rules/taxonomic_profiling/krakenuniq.smk
diff --git a/workflow/rules/taxonomic_profiling/metaphlan.smk b/workflow/rules/taxonomic_profiling/metaphlan.smk
