Skip to content
wdecosterPublic

About

No description, website, or topics provided.

Resources

Stars

263 stars

Watchers

6 watching

Forks

Latest commit

 

History

193 Commits

Folders and files

Repository files navigation

chopper

Rust implementation of NanoFilt+NanoLyse, both originally written in Python. This tool, intended for long read sequencing such as PacBio or ONT, filters and trims a fastq or unaligned BAM (uBAM) file.

Filtering is based on average read quality, minimum or maximum read length, and GC content percentage.

On the other hand, trimming is performed using one of four approaches:

  • Fixed cropping (head crop at the start of the read and tail crop at the end),

  • Quality-based trimming using a threshold,

  • Extracting the highest-quality sub-read based on a quality score,

  • Splitting reads by low-quality segments and outputting high-quality parts.

Optionally, poly(A) tails and poly(T) heads are removed from RNA and cDNA reads, see Poly(A) tails.

Reads that pass the filters are printed to standard output (STDOUT), in the same format as the input unless --output-format says otherwise. See Unaligned BAM for how BAM tags and headers are handled.

Base modification tags (MM/ML) carried over from a modified-basecalled BAM can be kept in step with the trimming with --update-mods, see Base modification tags.

Compared to the Python implementation the scope is to deliver the same results, almost the same functionality, at much faster execution times. At the moment this tool does not support filtering using a sequencing_summary file. If those features are of interest then please reach out.

Installation

Preferably, for most users, download a ready-to-use binary for your system to add directory on your $PATH from the releases.
You may have to change the file permissions to execute it with chmod +x chopper

Alternatively, use conda to install
conda install -c bioconda chopper

Usage

Reads on stdin and writes to stdout. The input format (FASTQ or unaligned BAM) is detected from the contents.

Usage: chopper [OPTIONS]

Options:
  -h, --help
          Print help (see a summary with '-h')

  -V, --version
          Print version

Filtering Options:
  -q, --quality <MINQUAL>
          Sets a minimum Phred average quality score
          
          [default: 0]

      --maxqual <MAXQUAL>
          Sets a maximum Phred average quality score
          
          [default: 1000]

  -l, --minlength <MINLENGTH>
          Sets a minimum read length
          
          [default: 1]

      --maxlength <MAXLENGTH>
          Sets a maximum read length
          
          [default: INF]

      --mingc <MINGC>
          Filter min GC content

      --maxgc <MAXGC>
          Filter max GC content

  -c, --contam <CONTAM>
          Filter contaminants against a fasta

Trimming Options:
      --trim-approach <TRIM_APPROACH>
          Select the trimming strategy to apply to the reads

          Possible values:
          - fixed-crop:           Remove a fixed number of bases from the start and/or end of the read. At least one of --headcrop and --tailcrop must be greater than 0
          - trim-by-quality:      Trim low-quality bases from the ends of the read until reaching a base with quality ≥ --cutoff
          - best-read-segment:    Extract the highest-quality read segment based on --cutoff, trimming low-quality bases from both ends
          - split-by-low-quality: Split reads by low-quality segments and output high-quality parts on the left and right, provided they pass the length filter

      --cutoff <CUTOFF>
          Set the minimum quality score (Q-score) threshold for trimming low-quality bases from read ends. Required when using the `trim-by-quality`, `best-read-segment`, or `split-by-low-quality` trimming approaches

      --split-window <SPLIT_WINDOW>
          Minimum number of consecutive bases below --cutoff required to split a read. Shorter low-quality runs are tolerated and kept inside the segment. Only used with the `split-by-low-quality` trimming approach
          
          [default: 1]

      --headcrop <HEADCROP>
          Trim N bases from the start of each read. Used only with the `fixed-crop` trimming approach, which requires at least one of --headcrop and --tailcrop to be greater than 0
          
          [default: 0]

      --tailcrop <TAILCROP>
          Trim N bases from the end of each read. Used only with the `fixed-crop` trimming approach, which requires at least one of --headcrop and --tailcrop to be greater than 0
          
          [default: 0]

      --update-mods
          Recompute the base modification tags (MM, ML and MN), from unaligned BAM or as written by `samtools fastq -T MM,ML,MN`, so that they match the trimmed read. Only those tags are corrected; other position- or quality-dependent tags (qs, ns, ts, du, ...) are passed through unchanged and go stale. Reads whose tags are malformed or do not describe their sequence are reported as an error rather than written out with wrong coordinates

      --trim-poly-tails [<MIN_LENGTH>]
          Remove a poly(T) head from the start and a poly(A) tail from the end of each read (U counts as T), allowing up to 20% other bases in them. Runs shorter than MIN_LENGTH [default: 10] are left alone, since reads end in a few A's by chance. Tails must reach the end of the read (a few other bases past them are removed along with them, within that 20%), so remove adapters and primers first. Applied after --trim-approach, if any

Setup Options:
  -t, --threads <THREADS>
          Use N parallel threads
          
          [default: 4]

  -i, --input <INPUT>
          Input filename [default: read from stdin]. FASTQ (optionally gzip, bzip2 or xz compressed) or unaligned BAM; the format is detected from the contents

      --output-format <OUTPUT_FORMAT>
          Output format [default: same as the input]. FASTQ output from BAM writes every BAM tag as a tab-separated SAM field after the read name, like `samtools fastq -T '*'`. BAM output keeps the input BAM header and tags, adding a @PG line

          Possible values:
          - fastq
          - bam:   Unaligned BAM

      --inverse
          Output the opposite of the normal results

Examples

gunzip -c reads.fastq.gz | chopper -q 10 -l 500 | gzip > filtered_reads.fastq.gz
chopper -q 10 -l 500 -i reads.fastq > filtered_reads.fastq
chopper -q 10 -l 500 -i reads.fastq.gz | gzip > filtered_reads.fastq.gz

# Trim low-quality bases from read ends
chopper --trim-approach trim-by-quality --cutoff 15 -i reads.fastq > trimmed_reads.fastq

# Extract the highest-quality segment from each read
chopper --trim-approach best-read-segment --cutoff 15 -l 100 -i reads.fastq > best_segments.fastq

# Split reads by low-quality segments and output high-quality parts
chopper --trim-approach split-by-low-quality --cutoff 15 -l 50 -i reads.fastq > split_reads.fastq

# Only split when at least 5 consecutive bases fall below the cutoff (tolerate shorter dips)
chopper --trim-approach split-by-low-quality --cutoff 15 --split-window 5 -l 50 -i reads.fastq > split_reads.fastq

# Keep base modification tags in step with the trimming
samtools fastq -T MM,ML,MN reads.bam \
  | chopper --trim-approach fixed-crop --headcrop 20 --tailcrop 20 --update-mods \
  > trimmed_reads.fastq

# Unaligned BAM straight from dorado: filter and trim, write uBAM
chopper -q 10 -l 500 -i calls.bam > filtered.bam
chopper --trim-approach fixed-crop --headcrop 20 --tailcrop 20 --update-mods -i calls.bam > trimmed.bam

# uBAM in, FASTQ out (every BAM tag kept as a SAM field in the header)
chopper -q 10 -i calls.bam --output-format fastq > filtered.fastq

Unaligned BAM

chopper reads unaligned BAM (as written by dorado and other ONT basecallers) as well as FASTQ, from a file or on stdin, and detects which one it got. Output is in the same format as the input; --output-format fastq or --output-format bam converts. Filtering and trimming are exactly the same for both formats.

What happens to the header and tags:

Input → output Header Tags
uBAM → uBAM Passed through unchanged (@HD, @RG, @PG, @CO, ...), with a @PG line for chopper chained after the last program. A sort order (SO, GO, SS) is kept only when the output keeps the input order: with --threads 1 and without split-by-low-quality. Otherwise it becomes SO:unsorted Every tag is kept with its exact type and in its original order. Only the read name (_segment_N suffix when split), sequence and qualities change, plus MM/ML/MN under --update-mods. Flags are kept
uBAM → FASTQ Not written (FASTQ has none) Every tag is written as a tab-separated SAM field after the read name, as samtools fastq -T '*' does
FASTQ → uBAM Built: @HD VN:1.6 SO:unknown and a @PG line for chopper Tab-separated SAM fields after the read name become typed BAM tags; integers get the smallest fitting type, as htslib does. Free text after the read name (e.g. runid=... ch=...) is kept in a CO:Z: tag. Records get flag 4 (unmapped)
FASTQ → FASTQ — Unchanged from earlier releases, byte for byte

No tag is ever dropped. The usual caveats still apply: tags that trimming invalidates other than MM/ML/MN (qs, ns, ts, du, mv, ...) are kept as they are and go stale, see Base modification tags. A FASTQ round trip through BAM does not keep integer widths, since SAM text has a single integer type, and an RG:Z: tag from FASTQ does not get a matching @RG header line.

chopper stops with an error rather than write something wrong when the BAM holds an aligned read (trimming would invalidate its position and CIGAR), a secondary or supplementary record (the read would be written twice), a record flagged as reverse complemented (MM counts bases on the original strand), one of a read pair (filtering each mate on its own would leave orphans still flagged as paired), or a read without base qualities. FASTQ to BAM also stops on header fields after a tab that are not valid SAM tags, and on read names BAM cannot hold (empty, longer than 254 characters, or containing @ or spaces).

BAM input is decompressed on a single thread. BAM output is compressed on a quarter of --threads (at least one) extra threads.

Poly(A) tails

--trim-poly-tails removes a poly(T) head from the start and a poly(A) tail from the end of each read, as left by RNA and cDNA protocols (cDNA reads come in both orientations). It is off by default.

chopper --trim-poly-tails -i cdna_reads.fastq > trimmed.fastq
# with quality trimming first, so that low-quality bases past a tail do not hide it
chopper --trim-approach trim-by-quality --cutoff 10 --trim-poly-tails -i cdna_reads.fastq > trimmed.fastq

Tails are found the way cutadapt's --poly-a finds them: scanning in from the end of the read, each A (or T at the start) scores +1 and any other base -2, and the tail ends where the score peaks, provided at most 20% of it are other bases. That tolerates sequencing errors inside the tail without running on into the transcript. Runs shorter than 10 bases are left alone, since a read ends in a few A's by chance; give another minimum as --trim-poly-tails 15. U is treated as T.

Two things to keep in mind:

  • Tails must reach the end of the read. Adapters and primers have to be removed first, for example by the basecaller or dorado trim. A few leftover bases past a tail are removed along with it, as long as the tail as a whole stays within the 20%.
  • A body that happens to end in A or start with T loses those bases with the tail, as they cannot be told apart from it.

It is applied after --trim-approach, to the outer ends of the read only, and works with --update-mods and BAM output. As with the trimming approaches, filters see the untrimmed read and only --minlength is checked again afterwards.

On simulated reads with ONT-like errors (half deletions, a quarter substitutions and a quarter insertions), at most 1.5 tail bases per read remained on average up to a 15% error rate, and 0.1% of reads without a tail were trimmed.

Base modification tags

Modified basecalls from dorado are carried in the unaligned BAM, which chopper reads directly, or into fastq as SAM tags:

samtools fastq -T MM,ML,MN reads.bam > reads.fastq

MM stores the position of each modified base as a count of canonical bases to skip, so trimming or splitting a read silently invalidates it: the coordinates now point at the wrong bases, or past the end of the read entirely. Passing --update-mods makes chopper recompute MM, ML and MN for whatever part of the read it keeps, for every trimming approach. Calls that fall outside the kept range are dropped along with their ML probabilities. It is off by default, and without it the tags are written through untouched, as before.

Only MM, ML and MN are corrected. Other tags that trimming also invalidates are passed through unchanged and will be stale, including the mean quality qs, the signal offsets ns and ts, and the duration du. If those matter downstream, recompute or drop them yourself.

Reads whose modification tags cannot be rewritten are a fatal error rather than a warning, because writing them through would attach modification coordinates to a sequence they no longer describe. This happens when the tags are malformed, when MM and ML disagree on how many calls there are, or when the tags describe a different sequence than the read carries (for example because the read was already trimmed by another tool without updating them). Rerun without --update-mods to write such reads through unchanged.

Note that --update-mods only fixes the tags. It does not make chopper aware of modification calls when deciding what to trim.

Performance

chopper's filtering is fast, and in practice the runtime is often dominated by (de)compressing gzipped FASTQ rather than by the filtering itself. Since chopper reads from stdin and writes to stdout, the recommended way to run it is to pipe data through pigz:

pigz -dc reads.fastq.gz | chopper -q 10 -l 500 | pigz > filtered_reads.fastq.gz

This helps in two ways:

  • Output compression is genuinely parallelised. pigz compresses the filtered output across multiple cores, which is usually much faster than a single-threaded gzip.
  • Decompression runs in its own process. Note that gzip decompression is inherently single-threaded (a gzip stream must be inflated sequentially, so pigz cannot parallelise it either). The benefit of piping is that the decompression happens in a separate process, overlapping with chopper's filtering on other cores. chopper itself already uses the fast zlib-rs backend for decompression, so reading a .gz file directly with -i is also efficient.

Citation

If you use this tool, please consider citing our publication.

About

No description, website, or topics provided.

Resources

Stars

263 stars

Watchers

6 watching

Forks

Releases

Packages

Used by

Contributors

Languages