atavide: Comprehensive Snakemake Workflow for Advanced Metagenomics Analysis and Annotation

public 1yr ago Version: Version 1 0 bookmarks

View Workflow

Help improve this workflow!

This workflow has been published but could be further improved with some additional meta data:

Keyword(s) in categories input, output

You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .

atavide is a complete workflow for metagenomics data analysis, including QC/QA, optional host removal, assembly and cross-assembly, and individual read based annotations. We have also built in some advanced analytics including tools to assign annotations from reads to contigs, and to generate metagenome-assembled genomes in several different ways, giving you the power to explore your data!

atavide is 100% snakemake and conda, so you only need to install the snakemake workflow, and then everything else will be installed with conda.

Steps:

QC/QA with prinseq++
optional host removal using bowtie2 and samtools, as described previously . To enable this, you need to provide a path to the host db and a host db.

Metagenome assembly

pairwise assembly of each sample using megahit
extraction of all reads that do not assemble using samtools flags
assembly of all unassembled reads using megahit
compilation of all contigs into a single unified set using Flye
comparison of reads -> contigs to generate coverage

MAG creation

metabat
concoct
Pairwise comparisons using turbocor followed by clustering

Read-based annotations

Want something else added to the suite? File an issue on github and we'll add it ASAP!

Installation

You will need to install

The NCBI taxonomy database somewhere
The superfocus databases somewhere, and set the SUPERFOCUS_DB environmental variable

Everything else should install automatically.

Code Snippets

shell:
    """
    python3 {params.sct} --file {input} --output {output.h5} --header --dataset contigs --indexfile {output.idx}
    """

SnakeMake From line 14 of rules/atavide_clusters.snakefile

shell:
    """
    turbocor compute {input} {output.tmp} --dataset contigs;
    turbocor topk 1000000 {output.tmp} > {output.pearson}
    """

SnakeMake turbocor From line 27 of rules/atavide_clusters.snakefile

shell:
    """
    python3 {params.sct} --file {input} --output {output} \
        --separator , --threshold {params.pearson_threshold}
    """

SnakeMake From line 43 of rules/atavide_clusters.snakefile

shell:
    """
    python3 {params.sct} --fasta {input.fa} --clusters {input.cl} \
    --directory {output.directory} --idx {input.idx}
    """

SnakeMake From line 60 of rules/atavide_clusters.snakefile

shell:
    """
    mkdir --parents {params.basename};
    jgi_summarize_bam_contig_depths --outputDepth {output.depth} {input.bam}
    """

SnakeMake From line 23 of rules/binning.snakefile

shell:
    """
    touch {output}
    metabat2 -i {input.contigs} -a {input.depth} -m 1500 -o {params.base} -t {threads}
    """

SnakeMake MetaBAT 2 From line 42 of rules/binning.snakefile

shell:
    """
    concoct --composition_file {input.contigs} --coverage_file {input.tsv} -b {params.od} -t {threads}
    """

SnakeMake CONCOCT From line 61 of rules/binning.snakefile

shell:
    """
    touch {output}
    extract_fasta_bins.py {input.contigs} {input.odir} --output_path {params.base}
    """      

SnakeMake From line 76 of rules/binning.snakefile

shell:
    """
    python {params.sct} --bam {input.bam} \
        --singlem {input.singlem} \
        --kraken {input.kraken} \
        --output {params.out} \
        --verbose
    """

SnakeMake From line 22 of rules/combine_read_annotations.snakefile

shell:
    """
    mkdir -p {output.d};
    focus -q {input.r1} -q {input.r2} -o {output.d} -t {threads}
    """

SnakeMake FOCUS From line 44 of rules/focus.snakefile

shell:
    """
    for F in {input}; do zip -j $F.zip $F; done
    """

SnakeMake From line 73 of rules/focus.snakefile

shell:
    """
    mkdir -p {params.d};
    seqtk fqchk {input.r1} > {output}
    """

SnakeMake seqtk From line 19 of rules/fqchk.snakefile

shell:
    """
    cut -f 1,8 {input} | tail -n +2 | grep -v ALL | sed -e 's/avgQ/{params.s}/' > {output}
    """

SnakeMake From line 33 of rules/fqchk.snakefile

shell:
    """
    perl {params.sct} -h {input} > {output}
    """

SnakeMake From line 46 of rules/fqchk.snakefile

shell:
    """
    seqtk sample {input.r1} {params.f} > {output}
    """

SnakeMake seqtk From line 33 of rules/kraken_rarefaction.snakefile

    shell:
        """
        kraken2 --report {output.rt} \
                --threads {threads} \
                {input.r1}
        """

"""

I think the next three rules will make this work, but need to double check
the rule. Essentially we just mix some bash in here too!

cd RBADIR
WD=$PWD; for SEQ in FAME0000*; do echo $SEQ; cd $SEQ/kraken; \
        echo -e "Fraction\t$SEQ" > $SEQ.kraken_rarefaction.tsv; \
        for FRX in 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9; \
        do awk -v F=$FRX '$4 == "S" {c++} END {print F,"\t",c}' $SEQ.report.$FRX.tsv; \
        done >> $SEQ.kraken_rarefaction.tsv; \
        awk '$4 == "S" {c++} END {print "1.0\t",c}' $SEQ.report.tsv >> $SEQ.kraken_rarefaction.tsv; \
        cd $WD; done

joinlists.pl -h */kraken/*kraken_rarefaction.tsv > ../statistics/kraken_rarefaction.tsv
"""

SnakeMake kraken2 From line 51 of rules/kraken_rarefaction.snakefile

shell:
    """
    touch {output}
    """

SnakeMake From line 82 of rules/kraken_rarefaction.snakefile

shell:
    """
    echo -e "Fraction\t{params.sample}" > {output};
    for FRX in {params.frx}; do awk -v F=$FRX '$4 == "S" {{c++}} END {{print F,"\t",c}}' {params.smp}.report.$FRX.tsv; done >> {output};
    awk '$4 == "S" {{c++}} END {{print "1.0\t",c}}' {params.smp}.report.tsv >> {output}
    """

SnakeMake From line 97 of rules/kraken_rarefaction.snakefile

shell:
    """
    perl {params.sct} -h {input} > {output}
    """

SnakeMake From line 113 of rules/kraken_rarefaction.snakefile

shell:
    """
    kraken2 --report {output.rt} \
            --output {output.ot} \
            --threads {threads} \
            {input.r1}
    """

SnakeMake kraken2 From line 18 of rules/kraken_taxonomy.snakefile

shell:
    """
    grep -v ^U {input} | \
    taxonkit lineage -j {threads} -i 3 --data-dir {params.t} | \
    taxonkit reformat -j {threads} --data-dir {params.t} -i 6 -f \
    "Root; d__{{k}}; p__{{p}}; c__{{c}}; o__{{o}}; f__{{f}}; g__{{g}}" \
    --fill-miss-rank | cut -f 2,3,7 > {output}
    """

SnakeMake TaxonKit From line 40 of rules/kraken_taxonomy.snakefile

shell:
    """
    echo "Family\t{wildcards.sample}" > {output};
    cat {input} | sed -e 's/Candidatus//' | \
    awk '{{if ($4 == "P") {{printf "%s\\t%f\\n", $6, $1}}; }}' >> {output}
    """

SnakeMake From line 55 of rules/kraken_taxonomy.snakefile

shell:
    """
    echo "Family\t{wildcards.sample}" > {output};
    cat {input} | sed -e 's/Candidatus//' | \
    awk '{{if ($4 == "F") {{printf "%s\\t%f\\n", $6, $1}}; }}' >> {output}
    """

SnakeMake From line 69 of rules/kraken_taxonomy.snakefile

shell:
    """
    echo "Family\t{wildcards.sample}" > {output};
    cat {input} | sed -e 's/Candidatus//' | \
    awk '{{if ($4 == "G") {{printf "%s\\t%f\\n", $6, $1}}; }}' >> {output}
    """

SnakeMake From line 83 of rules/kraken_taxonomy.snakefile

shell:
    """
    echo "Species\t{wildcards.sample}" > {output};
    cat {input} | sed -e 's/Candidatus//' | \
    awk '{{if ($4 == "S") {{for (i=6; i<=NF; ++i) {{printf "%s ", $i}} printf "\\t%f\\n", $1}}; }}' \
    >> {output}
    """

SnakeMake From line 96 of rules/kraken_taxonomy.snakefile

shell:
    """
    perl {params.sct} -z -h {input} > {output}
    """

SnakeMake From line 112 of rules/kraken_taxonomy.snakefile

shell:
    """
    perl {params.sct} -z -h {input} > {output}
    """

SnakeMake From line 125 of rules/kraken_taxonomy.snakefile

shell:
    """
    perl {params.sct} -z -h {input} > {output}
    """

SnakeMake From line 137 of rules/kraken_taxonomy.snakefile

shell:
    """
    perl {params.sct} -z -h {input} > {output}
    """

SnakeMake From line 149 of rules/kraken_taxonomy.snakefile

shell:
    """
        prinseq++ -min_len 60 -min_qual_mean 25 -ns_max_n 1 -derep 1 \
                -out_format 0 -trim_tail_left 5 -trim_tail_right 5 \
                -ns_max_n 5  -trim_qual_type min -trim_qual_left 30 \
                -trim_qual_right 30 -trim_qual_window 10 \
                -threads {threads} \
                -out_name {params.o} \
                -out_bad /dev/null \
                -out_bad2 /dev/null \
                -fastq {input.r1} \
                -fastq2 {input.r2};
    """

SnakeMake From line 18 of rules/qc_qa.snakefile

shell:
    """
    for F in {input}; do gzip -c $F > $F.gz; done
    """

SnakeMake From line 54 of rules/qc_qa.snakefile

shell:
    """
    rmdir {params.odir} ; 
    megahit -1 {input.r1} -2 {input.r2} -r {input.s1} -r {input.s2} \
            -o {params.odir} -t {threads} --use-gpu --mem-flag 2
    """

SnakeMake MEGAHIT From line 38 of rules/round1_assembly.snakefile

shell:
    """
    rmdir {params.odir} ; 
    megahit -1 {input.r1} -2 {input.r2} -r {input.s1} -r {input.s2} \
            -o {params.odir} -t {threads} --mem-flag 2
    """

SnakeMake MEGAHIT From line 69 of rules/round1_assembly.snakefile

script:
    "../scripts/renumber_merge_fasta_smk.py"

SnakeMake From line 90 of rules/round1_assembly.snakefile

    shell:
        """
        bowtie2 --mm -x {params.contigs} -1 {input.r1} -2 {input.r2} \
        -U {input.s1} -U {input.s2} --threads {threads} | \
        samtools view -bh | samtools sort -o {output} -
        """

"""
Using samtools to extract unmapped reads from the bam files
A samtools flag of -f 77 means the read is paird, neither the read nor mate is mapped, and the it is the first read in the pair, while a flag of -f 141 means the same except it is the second mate in the pair.
Then we use two flags: -f 4 (the read is unmapped) and -F 1 (the read is not paired) to find the single reads that are not mapped

See: https://edwards.sdsu.edu/research/command-line-deconseq/
"""

SnakeMake SAMtools Bowtie 2 From line 147 of rules/round1_assembly.snakefile

shell:
    """
    samtools view -@ {threads} -h {input} | 
            awk 'BEGIN {{FS="\t"; OFS="\t"}} 
            {{if (/^@/ && substr($2, 3, 1)==":") {{print}} 
            else if (and($2, 0x1) && and($2, 0x40) && 
            (and($2, 0x4) || and($2, 0x8))) {{print}}}}' \
            | samtools bam2fq -@ {threads} > {output}
    """

SnakeMake SAMtools From line 175 of rules/round1_assembly.snakefile

shell:
    """
    samtools view -@ {threads} -h {input} | 
            awk 'BEGIN {{FS="\t"; OFS="\t"}} 
            {{if (/^@/ && substr($2, 3, 1)==":") {{print}} 
            else if (and($2, 0x1) && and($2, 0x80) && 
            (and($2, 0x4) || and($2, 0x8))) {{print}}}}' \
            | samtools bam2fq -@ {threads} > {output}
    """

SnakeMake SAMtools From line 198 of rules/round1_assembly.snakefile

    shell:
        "samtools fastq -@ {threads} -f 4 -F 1 {input} > {output}"

"""
We concatanate the unassembled reads into separate R1/R2/s files so
we can assmble them all together.

We also do this in 3 separate threads to take advantage of parallelization
and to make the command easier
"""

SnakeMake SAMtools From line 221 of rules/round1_assembly.snakefile

shell:
    """
    for F in {input}; do gzip $F; done
    """

SnakeMake From line 242 of rules/round1_assembly.snakefile

shell:
    "cat {input} > {output}"

SnakeMake From line 257 of rules/round1_assembly.snakefile

shell:
    "cat {input} > {output}"

SnakeMake From line 268 of rules/round1_assembly.snakefile

shell:
    "cat {input} > {output}"

SnakeMake From line 279 of rules/round1_assembly.snakefile

shell:
    """
    rmdir {params.odir};
    megahit -1 {input.r1} -2 {input.r2} -r {input.s0} -o {params.odir} -t {threads}  --use-gpu --mem-flag 2
    """

SnakeMake MEGAHIT From line 35 of rules/round2_assembly.snakefile

    shell:
        """
        rmdir {params.odir};
        megahit -1 {input.r1} -2 {input.r2} -r {input.s0} -o {params.odir} -t {threads} --mem-flag 2
        """

"""
Combine all the contigs and use metaflye to merge the subassmblies
"""

SnakeMake MEGAHIT From line 64 of rules/round2_assembly.snakefile

script:
    "../scripts/renumber_merge_fasta_smk.py"

SnakeMake From line 90 of rules/round2_assembly.snakefile

shell:
    """
    flye --meta --subassemblies {input.contigs} -o {CCMO} --threads {threads}
    """

SnakeMake Flye From line 111 of rules/round2_assembly.snakefile

script:
    "../scripts/renumber_merge_fasta_smk.py"

SnakeMake From line 130 of rules/round2_assembly.snakefile

shell:
    """
    bowtie2 --mm -x {params.contigs} -1 {input.r1} -2 {input.r2} \
    -U {input.s1} -U {input.s2} --threads {threads} | \
    samtools view -bh | samtools sort -o {output} -
    """

SnakeMake SAMtools Bowtie 2 From line 187 of rules/round2_assembly.snakefile

shell:
    "samtools index {input}"

SnakeMake SAMtools From line 201 of rules/round2_assembly.snakefile

shell:
    """
    samtools idxstats -@ {threads} {input.bam} | cut -f 1,3 > {output}
    """

SnakeMake SAMtools From line 215 of rules/round2_assembly.snakefile

script:
    "../scripts/rpkm.py"

SnakeMake From line 229 of rules/round2_assembly.snakefile

shell:
    """
    mkdir --parents {output.d};
    singlem pipe --forward {input.r1} --reverse {input.r2} --otu_table {output.otu} --output_extras --threads {threads}
    """

SnakeMake singlem From line 14 of rules/singlem.snakefile

script:
    "../scripts/countfastq.py"

SnakeMake From line 14 of rules/statistics.snakefile

script:
    "../scripts/countfastq.py"

SnakeMake From line 25 of rules/statistics.snakefile

shell:
    """
    python3 {params.sct} -f {input} > {output}
    """

SnakeMake From line 35 of rules/statistics.snakefile

shell:
    """
    perl {params.sct} -t {input} > {output}
    """

SnakeMake From line 49 of rules/statistics.snakefile

shell:
    """
    perl {params.sct} -t {input} > {output}
    """

SnakeMake From line 63 of rules/statistics.snakefile

shell:
    """
    python3 {params.sct} -f {input} -ln > {output}
    """

SnakeMake From line 76 of rules/statistics.snakefile

shell:
    """
    perl {params.sct} -t {input} > {output}
    """

SnakeMake From line 90 of rules/statistics.snakefile

shell:
    """
    perl {params.sct} {input}  > {params.out} && gzip {params.out}
    """

SnakeMake From line 104 of rules/statistics.snakefile

    shell:
        """
        superfocus -q {input.r1} -q {input.r2} -dir {params.d} -a mmseqs2 -t {threads} -n 0 -tmp $(mktemp -d -p {params.TMPDIR})
        """


"""
In the rules below we use the m8 files from superfocus
to crate taxonomy tables.
"""

SnakeMake SUPER-FOCUS From line 36 of rules/superfocus.snakefile

shell:
    """
    perl -F"\\t" -lane 'if ($F[1] =~ /fig\|(\d+)\.\d+/ && $1 != 6666666) \
       {{print "$F[0]\\t$1"}}' {input.m8} | \
    taxonkit lineage -j {threads} -i 2 --data-dir {params.t} | \
    taxonkit reformat -j {threads} --data-dir {params.t} -i 3 \
      -f "Root; d__{{k}}; p__{{p}}; c__{{c}}; o__{{o}}; f__{{f}}; g__{{g}}" \
      --fill-miss-rank | \
    cut -f 1,2,4 > {output}
    """

SnakeMake TaxonKit From line 60 of rules/superfocus.snakefile

script:
    "../scripts/count_superfocus_taxonomy.py"

SnakeMake From line 77 of rules/superfocus.snakefile

shell:
    """
    for F in {input}; do zip -j $F.zip $F; done
    """

SnakeMake From line 92 of rules/superfocus.snakefile

shell:
    """
    for F in {input}; do zip -j $F.zip $F; done
    """

SnakeMake From line 105 of rules/superfocus.snakefile

import os
import sys
from atavide_lib import stream_fastq

__author__ = 'Rob Edwards'


# snakemake.input.fqdir is the directory of fastq files to count
# snakemake.output.stats is the output file to write

with open(snakemake.output.stats, 'w') as out:
    print("File\tNumber of Sequences\tTotal bp\tShortest length\tLongest length\tN50\tN75\tAuN", file=out)
    for f in os.listdir(snakemake.input.fqdir):
        if 'fastq' in f or 'fq' in f:
            print(f"Counting {f}", file=sys.stderr)
            lens = []
            for (sid, label, seq, qual) in stream_fastq(os.path.join(snakemake.input.fqdir, f)):
                lens.append(len(seq))
            lens.sort()
            length=sum(lens)

            len_so_far = 0
            n50 = ""
            n75 = ""
            auN = 0
            for i in lens:
                len_so_far += i
                if not n50 and len_so_far >= length * 0.5:
                    n50 = i
                if not n75 and len_so_far >= length * 0.75:
                    n75 = i
                auN += i**2
            if length > 0:
                auN /= length
            else:
                auN = 0

            print(f"{f}\t{len(lens):,}\t{length:,}\t{lens[0]:,}\t" \
                  + f"{lens[-1]:,}\t{n50:,}\t{n75:,}\t{int(auN):,}", file=out)
        else:
            print(f"Skipped {os.path.join(subdir, f)}. Does not appear to be fastq", file=sys.stderr)

Python From line 8 of scripts/countfastq.py

from atavide_lib import read_fasta, colours

__author__ = 'Rob Edwards'

# snakemake.input is an array of files to merge
# snakemake.output.contigs is the name of the final contigs file
# snakemake.output.ids is the name of the ids file to write3
# snakemake.params.sample_id is the sample_id that will be prepended to the contig names


idmap = open(snakemake.output.ids, 'w')
out = open(snakemake.output.contigs, 'w')

counter = 0

for f in snakemake.input:
    fa = read_fasta(f)
    for id in fa:
        counter += 1
        if hasattr(snakemake, 'params') and hasattr(snakemake.params, 'sample_id'):
            out.write(f">{snakemake.params.sample_id}_{counter}\n{fa[id]}\n")
            idmap.write(f"{f}\t{id}\t{snakemake.params.sample_id}_{counter}\n")
        else:
            out.write(f">{counter}\n{fa[id]}\n")
            idmap.write(f"{f}\t{id}\t{counter}\n")

idmap.close()
out.close()

Python From line 10 of scripts/renumber_merge_fasta_smk.py

import os
import sys
from atavide_lib import stream_fastq

__author__ = 'Rob Edwards'

if not os.path.exists(snakemake.input.r1):
    sys.stderr.write(f"ERROR: snakemake.input.r1 {snakemake.input.r1} not found\n")
if not os.path.exists(snakemake.input.contig_len):
    sys.stderr.write(f"ERROR: snakemake.input.contig_len {snakemake.input.contig_len} not found\n")
if not os.path.exists(snakemake.input.hits):
    sys.stderr.write(f"ERROR: snakemake.input.hits {snakemake.input.hits} not found\n")


numseqs = 0
for seqid, header, seq, qualscores in stream_fastq(snakemake.input.r1):
    numseqs+=1
numseqs /= 1e6

contigs = {}
with open(snakemake.input.contig_len, 'r') as f:
    for l in f:
        r = l.strip()
        if not r:
            continue
        p = r.split("\t")
        if len(p) < 2:
            sys.stderr.write(f"Skipped {r} in {f}\n")
            continue
        contigs[p[0]] = int(p[1])

with open(snakemake.input.hits, 'r') as f, open(snakemake.output.rpkm, 'w') as out:
    for l in f:
        p = l.strip().split("\t")
        if p[0] == '*':
            continue
        if p[0] not in contigs:
            sys.stderr.write(f"ERROR: contig {p[0]} does not have a length.\n")
            contigs[p[0]]=1
        rpkm = int(p[1])/numseqs/contigs[p[0]]
        print(f"{p[0]}\t{rpkm}", file=out)