Benchmarking adapter and quality trimming tools

public 1yr ago 0 bookmarks

View Workflow

Benchmarking adapter and quality trimming tools

Original implementation by Brian Bushnell (2014): http://seqanswers.com/forums/showthread.php?t=42776

The following tools are compared:

bbduk
cutadapt
fastp
trimmomatic

Fake adapters, "gruseq"

The fake truseq adapters, "gruseq", provided by Brian Bushnell, downloaded from: http://seqanswers.com/forums/attachment.php?attachmentid=2993&d=1398383571

Test data

A single sample is downloaded from SRA. Feel free to replace it with whatever you want.

Running

Run the benchmarking workflow with snakemake --use-conda --jobs 10

Code Snippets

__author__ = "Fredrik Boulund"
__date__ = "2019"

from sys import argv, exit
from pathlib import Path
import argparse

import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt


parser = argparse.ArgumentParser(description=f"{__doc__}. {__author__} (c) {__date__}")
parser.add_argument("benchmarks", metavar="FILE", nargs="+", help="Benchmark output(s)")
parser.add_argument("--output", default="plot.pdf", 
        help="Plot output. Will produce a png variant as well")

if len(argv) < 2:
    parser.print_help()
    exit()

args = parser.parse_args()

dfs = []
for benchmark in args.benchmarks:
    tool = Path(benchmark).stem.split(".", maxsplit=1)[0]
    df = pd.read_csv(benchmark, sep="\t")
    df["Tool"] = tool
    dfs.append(df)

table = pd.concat(dfs)


fig, (ax1, ax2) = plt.subplots(1,2, figsize=(10,5))

ax1.set_title("Average time")
ax1.set_ylabel("seconds")
table\
    .groupby("Tool")\
    .mean()["s"]\
    .plot(kind="bar", ax=ax1)

ax2.set_title("Average max_vms")
ax2.set_ylabel("Megabytes")
table\
    .groupby("Tool")\
    .mean()["max_vms"]\
    .plot(kind="bar", ax=ax2)


fig.savefig(args.output, bbox_inches="tight")
fig.savefig(args.output.replace(".pdf", ".png"), bbox_inches="tight")

Python Pandas matplotlib From line 3 of scripts/plot_benchmarks.py

shell:
    """
    wget ftp://ftp.sra.ebi.ac.uk/vol1/fastq/SRR921/004/SRR9218144/SRR9218144.fastq.gz
    """

SnakeMake From line 27 of master/Snakefile

shell:
    """
    addadapters.sh \
        in={input.fastq} \
        out={output} \
        qout=33 \
        ref={input.adapters} \
        right \
        int=f \
        2> {log}
    """

SnakeMake From line 40 of master/Snakefile

shell:
    """
    cutadapt \
        --cores {threads} \
        --minimum-length 10 \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATATGATACTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATAACTGCGTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATAGGTCCATGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATAGCTAATTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATATATCGCTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATACAATTGTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATAATCTGATGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATATAGGCTTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATACTGATCTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATAGTCAGGTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATACCAGTATGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATAAGGCGTTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATATCGATTATTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATATCGGAACGTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATATGCGATCTTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATAAACGAAACTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATACGAACATATGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATACGCTTTACTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATACGCCAAGGTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATACGGGACCTTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATAACGTACGTTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATACTCGCCTGTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATATAGCTGTGTGAGACGTGCAACGAGGAGCAGGC \
        --anywhere CTGACCTTCTCATATACGAGCTTAGAATCGATATGGAAGGGTGAGACGTGCAACGAGGAGCAGGC \
        {input.fastq} \
        > {output} \
        2> {log}
    """

SnakeMake Cutadapt From line 63 of master/Snakefile

shell:
    """
    trimmomatic \
        SE \
        -phred33 \
        -threads {threads} \
        {input.fastq} \
        {output} \
        ILLUMINACLIP:gruseq.fa:2:28:10 \
        MINLEN:10 \
        2> {log}
    """

SnakeMake Trimmomatic From line 108 of master/Snakefile

shell:
    """
    bbduk.sh \
        in={input.fastq} \
        out={output} \
        ref={input.adapters} \
        ktrim=r \
        mink=12 \
        hdist=1 \
        minlen=10 \
        threads={threads} \
        2> {log}
    """

SnakeMake BBMap From line 133 of master/Snakefile

shell:
    """
    fastp \
        --in1 {input.fastq} \
        --out1 {output.fq} \
        --adapter_fasta {input.adapters} \
        --thread {threads} \
        --html {output.html} \
        --json {output.json} \
        --length_required 10 \
        2> {log}
    """

SnakeMake fastp From line 160 of master/Snakefile

shell:
    """
    addadapters.sh \
        in={input} \
        grade \
        2> {output}
    """

SnakeMake From line 179 of master/Snakefile

shell:
    """
    scripts/plot_benchmarks.py \
        --output {output.benchmarks} \
        {input.benchmarks}
    scripts/plot_grades.py \
        --output {output.grades} \
        {input.grades}
    """