Merge pull request #137 from sunbeam-labs/dev

Low-complexity read removal and misc bug fixes

Author: eclarke

Readme.md

@@ -9,7 +9,7 @@ Sunbeam is a pipeline written in [snakemake](http://snakemake.readthedocs.io)
 that simplifies and automates many of the steps in metagenomic sequencing
 analysis. It uses [conda](http://conda.io) to manage dependencies, so it
 doesn't have pre-existing dependencies or admin privileges, and can be deployed
-on most Linux and Mac workstations and clusters.
+on most Linux workstations and clusters.
 
 Sunbeam currently automates the following tasks:
 

Snakefile

@@ -4,7 +4,7 @@
 # Author: Erik Clarke <ecl@mail.med.upenn.edu>
 # Created: 2016-04-28
 #
-
+import os
 import re
 import sys
 import yaml
@@ -26,6 +26,14 @@ if not config:
         "No config file specified. Run `sunbeam init` to generate a "
         "config file, and specify with --configfile")
 
+sunbeam_dir = ""
+try:
+    sunbeam_dir = os.environ["SUNBEAM_DIR"]
+except KeyError:
+    raise SystemExit(
+        "$SUNBEAM_DIR environment variable not defined. Are you sure you're "
+        "running this from the Sunbeam conda env?")
+
 # Check for major version compatibility
 pkg_major, cfg_major = check_compatibility(config)
 if pkg_major > cfg_major:
@@ -40,7 +48,7 @@ elif pkg_major < cfg_major:
         "`sunbeam init` and update it using `sunbeam_mod_config`\n")
 
 # Load extensions
-sbxs = list(listfiles("extensions/{sbx_folder}/{sbx}.rules"))
+sbxs = list(listfiles(sunbeam_dir+"/extensions/{sbx_folder}/{sbx}.rules"))
 for sbx in sbxs:
     sys.stderr.write("Found extension {sbx} in folder {sbx_folder}\n".format(**sbx[1]))
 

docs/index.rst

@@ -18,7 +18,7 @@ Sunbeam is a pipeline written in `snakemake <http://snakemake.readthedocs.io>`_
 that simplifies and automates many of the steps in metagenomic sequencing
 analysis. It uses `conda <http://conda.io>`_ to manage dependencies, so it
 doesn't have pre-existing dependencies or admin privileges, and can be deployed
-on most Linux and Mac workstations and clusters.
+on most Linux workstations and clusters.
 
 Sunbeam currently automates the following tasks:
 

docs/quickstart.rst

@@ -10,11 +10,12 @@ Quickstart Guide
 Installation
 ************
 
-Download a copy of Sunbeam from our GitHub repository, and install.
+On a Linux machine, download a copy of Sunbeam from our GitHub repository, and
+install. We do not currently support non-Linux environments.
 
 .. code-block:: shell
 
-   git clone -b stable https://github.com/eclarke/sunbeam sunbeam-stable
+   git clone -b stable https://github.com/sunbeam-labs/sunbeam sunbeam-stable
    cd sunbeam-stable
    ./install.sh
    tests/run_tests.bash -e sunbeam
@@ -31,7 +32,7 @@ runs some tests to make sure everything was installed correctly.
    PATH=$PATH:$HOME/miniconda3/bin` > ~/.bashrc``
 
 If you see "Tests failed", check out our :ref:`troubleshooting` section or file an issue
-on our `GitHub <https://github.com/eclarke/sunbeam/issues>`_ page.
+on our `GitHub <https://github.com/sunbeam-labs/sunbeam/issues>`_ page.
 
 Setup
 *****

docs/usage.rst

@@ -7,6 +7,15 @@ User Guide
 .. contents::
    :depth: 3
 
+Requirements
+============
+
+- A relatively-recent Linux computer with more than 2Gb of RAM
+
+We do not currently support Windows or Mac. (You may be able to run this on
+Windows using the [WSL](https://docs.microsoft.com/en-us/windows/wsl/about), but
+it has not been tested.
+      
 .. _installation:
 Installation
 ============

rules/assembly/assembly.rules

@@ -2,14 +2,14 @@
 #
 # Contig building and other assembly rules
 #
-# Requires Megahit and CAP3.
+# Requires Megahit.
 
 rule all_assembly:
     """Build contigs for all samples."""
     input:
         TARGET_ASSEMBLY
 
-ruleorder: megahit_unpaired > megahit_paired
+ruleorder: megahit_paired > megahit_unpaired
 
 rule megahit_paired:
     input:

rules/qc/qc.rules

@@ -7,7 +7,7 @@ rule all_qc:
     input:
         TARGET_QC
 
-ruleorder: adapter_removal_unpaired > adapter_removal_paired
+ruleorder: adapter_removal_paired > adapter_removal_unpaired
 
 rule adapter_removal_unpaired:
     input:
@@ -79,7 +79,7 @@ rule adapter_removal_paired:
             ln -s {input.r1} {output.gr1} && ln -s {input.r2} {output.gr2}
             """)
 
-ruleorder: trimmomatic_unpaired > trimmomatic_paired
+ruleorder: trimmomatic_paired > trimmomatic_unpaired
 
 rule trimmomatic_unpaired:
     input:
@@ -149,21 +149,34 @@ rule fastqc:
         outdir = str(QC_FP/'reports')
     shell:
         "fastqc -o {params.outdir} {input.reads} -extract"
-           
+
+rule find_low_complexity:
+    input:
+        expand(
+            str(QC_FP/'02_trimmomatic'/'{{sample}}_{rp}.fastq.gz'),
+            rp=Pairs)
+    output:
+        str(QC_FP/'log'/'komplexity'/'{sample}.filtered_ids')
+    shell:
+        """
+        for rp in {input}; do
+          gzip -dc $rp | kz | \
+          awk '{{ if ($4<{Cfg[qc][kz_threshold]}) print $1 }}' >> {output}
+        done
+        """
+
 rule remove_low_complexity:
     input:
-        str(QC_FP/'02_trimmomatic'/'{sample}_{rp}.fastq.gz')
+        reads = str(QC_FP/'02_trimmomatic'/'{sample}_{rp}.fastq.gz'),
+        ids = str(QC_FP/'log'/'komplexity'/'{sample}.filtered_ids')
     output:
         str(QC_FP/'03_komplexity'/'{sample}_{rp}.fastq.gz')
-    run:
-        if Cfg['qc']['mask_low_complexity']:
-            shell("""
-            kz --mask -t {Cfg[qc][kz_threshold]} -w {Cfg[qc][kz_window]} \
-            < <(gzip -cd {input}) | gzip -c > {output}
-            """)
-        else:
-            shell("ln -s {input} {output}")
-            
+    shell:
+        """
+        gzip -dc {input.reads} | rbt fastq-filter {input.ids} |\
+        gzip > {output}
+        """
+
 rule qc_final:
     input:
         str(QC_FP/'03_komplexity'/'{sample}_{rp}.fastq.gz')

sunbeamlib/data/default_config.yml

@@ -39,9 +39,7 @@ qc:
   fwd_adapters: ['GTTTCCCAGTCACGATC', 'GTTTCCCAGTCACGATCNNNNNNNNNGTTTCCCAGTCACGATC']
   rev_adapters: ['GTTTCCCAGTCACGATC', 'GTTTCCCAGTCACGATCNNNNNNNNNGTTTCCCAGTCACGATC']
   # Komplexity
-  mask_low_complexity: true
   kz_threshold: 0.55
-  kz_window: 32
   # Decontam.py
   pct_id: 0.5
   frac: 0.6

tests/targets_singleend.txt

@@ -0,0 +1,28 @@
+annotation/blastn/bacteria/contig/dummybfragilis.xml
+annotation/blastn/bacteria/contig/dummyecoli.xml
+annotation/blastn/bacteria/contig/random.xml
+
+annotation/summary/dummybfragilis.tsv
+annotation/summary/dummyecoli.tsv
+annotation/summary/random.tsv
+
+assembly/contigs/dummybfragilis-contigs.fa
+assembly/contigs/dummyecoli-contigs.fa
+assembly/contigs/random-contigs.fa
+
+classify/kraken/all_samples.biom
+classify/kraken/all_samples.tsv
+classify/kraken/dummybfragilis-taxa.tsv
+classify/kraken/dummyecoli-taxa.tsv
+classify/kraken/random-taxa.tsv
+
+qc/decontam/dummybfragilis_1.fastq.gz
+qc/decontam/dummyecoli_1.fastq.gz
+qc/decontam/random_1.fastq.gz
+
+qc/cleaned/dummybfragilis_1.fastq.gz
+qc/cleaned/dummyecoli_1.fastq.gz
+qc/cleaned/random_1.fastq.gz
+
+mapping/human/coverage.csv
+mapping/phix174/coverage.csv

tests/test_suite.bash

@@ -43,7 +43,31 @@ function test_version_check {
     fi
 }
 
-# Test that we can 
+# Test that we detect and run extensions
 function test_extensions {
     sunbeam run --configfile $TEMPDIR/tmp_config.yml sbx_test | grep "SBX_TEST"
 }
+
+# Test that single-end sequencing configurations work
+function test_single_end {
+    rm -rf $TEMPDIR/sunbeam_output/qc
+    sunbeam config modify --str 'all: {paired_end: false}' \
+	    $TEMPDIR/tmp_config.yml > $TEMPDIR/single_end_config.yml
+    sunbeam run --configfile $TEMPDIR/single_end_config.yml
+    python tests/find_targets.py --prefix $TEMPDIR/sunbeam_output tests/targets_singleend.txt
+}
+
+# Fix for #131
+# Test that paired-end qc rules produce files with the same number of reads
+function test_pair_concordance {
+    rm -rf $TEMPDIR/sunbeam_output/qc
+    sunbeam run --configfile $TEMPDIR/tmp_config.yml all_decontam
+    for r1 in $TEMPDIR/sunbeam_output/qc/cleaned/*_1.fastq.gz; do
+	r1_lines=$(zcat $r1 | wc -l)
+	r2=${r1%_1.fastq.gz}_2.fastq.gz
+	r2_lines=$(zcat $r2 | wc -l)
+	if [ $r1_lines -ne $r2_lines ]; then
+	    exit 1
+	fi
+    done
+}