Merge pull request #24 from NBISweden/develop

Merge develop

**Changelog**:
- Update sourmash and contigtax versions
- Update and make sourmash database link configurable
- Add sourmash kmer size to config
- Update normalization in R scripts to handle one-sample data

Author: johnne

.test/config/annotate.yaml

@@ -140,6 +140,12 @@ taxonomy:
   # hash fraction to use for sourmash when computing signatures for contigs
   # this is evaluated as 1/<sourmash_fraction>
   sourmash_fraction: 100
+  # kmer size to use when computing signatures
+  sourmash_kmer_size: 31
+  # sourmash database to use
+  # see https://sourmash.readthedocs.io/en/latest/databases.html
+  # and choose an LCA database link matching the sourmash_kmer_size
+  sourmash_database_url: https://osf.io/ypsjq/download
   # ranks to report taxonomy for
   ranks:
     - "superkingdom"

config/config.yaml

@@ -139,10 +139,18 @@ taxonomy:
   # modify assignments without having to rerun the search step
   search_params: "--evalue 0.01 --top 10"
   # parameters for contigtax assign
-  assign_params: "--evalue 0.001 --top 5"
+  # by default the lowest common ancestor is assigned from the list of top hits
+  # to instead use a voting strategy set --mode rank_vote
+  assign_params: "--evalue 0.001 --top 5 --mode rank_lca"
   # hash fraction to use for sourmash when computing signatures for contigs
   # this is evaluated as 1/<sourmash_fraction>
   sourmash_fraction: 100
+  # kmer size to use when computing signatures
+  sourmash_kmer_size: 31
+  # sourmash database to use
+  # see https://sourmash.readthedocs.io/en/latest/databases.html
+  # and choose an LCA database link matching the sourmash_kmer_size
+  sourmash_database_url: https://osf.io/ypsjq/download
   # ranks to report taxonomy for
   ranks:
     - "superkingdom"

workflow/envs/sourmash.yml

@@ -0,0 +1,8 @@
+name: sourmash
+channels:
+  - bioconda
+  - conda-forge
+  - defaults
+dependencies:
+  - python=3.9.6
+  - sourmash=4.2.2

workflow/envs/taxonomy.yml

@@ -1,9 +1,9 @@
+name: taxonomy
 channels:
   - conda-forge
   - bioconda
   - defaults
 dependencies:
   - python=3.7.6
-  - contigtax=0.5.9
+  - contigtax=0.5.10
   - diamond=0.9.24
-  - sourmash=3.3.0

workflow/rules/annotation.smk

@@ -301,7 +301,7 @@ if config["runOnUppMax"]:
         params:
             resource_dir=lambda wildcards, input: os.path.dirname(input[1]),
             tmpdir=temppath+"/{assembly}-eggnog",
-            out=temppath+"/annotation/{assembly}/{assembly}",
+            out=results+"/annotation/{assembly}/{assembly}",
             flags="--no_file_comments"
         log:
             results+"/annotation/{assembly}/{assembly}.emapper.annotations.log"
@@ -313,6 +313,7 @@ if config["runOnUppMax"]:
             runtime=lambda wildcards, attempt: attempt**2*60
         shell:
             """
+            if [ -z ${{SLURM_JOB_ID+x}} ]; then SLURM_JOB_ID="emapper_annotate_hits_uppmax"; fi
             #Copy eggnog.db
             mkdir -p /dev/shm/$SLURM_JOB_ID
             cp {params.resource_dir}/eggnog.db {params.resource_dir}/eggnog_proteins.dmnd /dev/shm/$SLURM_JOB_ID
@@ -351,6 +352,8 @@ rule parse_emapper:
     input:
         annotations = results+"/annotation/{assembly}/{assembly}.emapper.annotations",
         info = "resources/kegg/kegg_{db}.tsv"
+    wildcard_constraints:
+        db="enzymes|pathways|kos|modules"
     output:
         results+"/annotation/{assembly}/{db}.parsed.tsv"
     script:

workflow/rules/classification.smk

@@ -2,6 +2,7 @@ from scripts.common import classify_input, krona_input, metaphlan_krona_string
 
 localrules:
     classify,
+    krona_taxonomy,
     download_kraken_build,
     download_centrifuge_build,
     centrifuge_kreport,

workflow/rules/quantification.smk

@@ -13,8 +13,9 @@ localrules:
 
 rule quantify:
     input:
-        expand(results+"/annotation/{assembly}/gene_{counts_type}.tsv",
-               assembly=assemblies.keys(), counts_type=["counts", "rpkm"])
+        expand("{results}/annotation/{assembly}/gene_{counts_type}.tsv",
+            results=[config["paths"]["results"]], assembly=assemblies.keys(),
+            counts_type=["counts", "rpkm"])
 
 
 rule write_featurefile:

workflow/rules/taxonomy.smk

@@ -201,15 +201,17 @@ rule contigtax_assign:
 
 rule download_sourmash_db:
     output:
-        "resources/sourmash/genbank-k31.lca.json"
+        lca="resources/sourmash/sourmash_db.lca.json",
+        version="resources/sourmash/version.txt"
     log:
         "resources/sourmash/download.log"
     params:
-        url="https://osf.io/4f8n3/download"
+        url=config["taxonomy"]["sourmash_database_url"]
     shell:
         """
-        curl -L -o {output}.gz {params.url} > {log} 2>&1
-        gunzip {output}.gz
+        curl -L -v -o {output.lca}.gz {params.url} > {log} 2>&1
+        grep filename {log} | cut -f2 -d ';' > {output.version}
+        gunzip {output.lca}.gz
         """
 
 rule sourmash_compute:
@@ -220,19 +222,20 @@ rule sourmash_compute:
     log:
         results+"/assembly/{assembly}/sourmash_compute.log"
     conda:
-        "../envs/taxonomy.yml"
+        "../envs/sourmash.yml"
     params:
-        frac=config["taxonomy"]["sourmash_fraction"]
+        frac=config["taxonomy"]["sourmash_fraction"],
+        k=config["taxonomy"]["sourmash_kmer_size"]
     shell:
         """
         sourmash compute --singleton --scaled {params.frac} \
-            -k 31 -o {output} {input} > {log} 2>&1
+            -k {params.k} -o {output} {input} > {log} 2>&1
         """
 
 rule sourmash_classify:
     input:
         sig=results+"/assembly/{assembly}/final_contigs.fa.sig",
-        db="resources/sourmash/genbank-k31.lca.json"
+        db="resources/sourmash/sourmash_db.lca.json"
     output:
         csv=results+"/annotation/{assembly}/taxonomy/sourmash.taxonomy.csv"
     log:
@@ -243,7 +246,7 @@ rule sourmash_classify:
     resources:
         runtime=lambda wildcards, attempt: attempt**2*30
     conda:
-        "../envs/taxonomy.yml"
+        "../envs/sourmash.yml"
     shell:
         """
         sourmash lca classify --db {input.db} --scaled {params.frac} \