cmuts hmm

Purpose

Computing reactivity rates from alignment files via the pair-HMM.

Requires

• One or more coordinate-sorted alignment files of single-end or merged reads. SAM, BAM, and CRAM formats are all supported

• The FASTA library

Alignments must be single-end. Paired-end data is not supported; please merge upstream before passing data to cmuts hmm.

Model Parameters

The pair HMM is configured by five internal parameters. --dump-params writes them in the form --params reads.

cmuts hmm --dump-params > params.txt

You may specify a subset of the parameters to modify only them.

Error Checking

cmuts hmm verifies the FASTA against the alignment header, by comparing each sequence’s name and length, and its MD5 checksum when present. Any mismatch between the header and the FASTA ends the run early. This behavior is configurable via the --verify flag.

A read carrying the paired flag also ends the run. Two mates read one molecule, so counting them as separate reads would count their overlap twice. Merge the mates before aligning, which cmuts align does for paired-end input.

Note

The length check is required to avoid buffer overflows and cannot be disabled.

Output

coverage

Shape (n, l) · Type float32 · Fill zero

The number of reads in which this base was present.

reactivity

Shape (n, l) · Type float32 · Fill NaN

The mutation rate at this base, weighted by PHRED scores and in accordance with the HMM parameters.

error

Shape (n, l) · Type float32 · Fill NaN

Standard error of the reactivity values. Purely the statistical error introduced by finite read depths; does not account for experimental or systemic errors.

reads/lengths

Shape (n, 2l) · Type uint64 · Fill zero

The number of reads passing all filters, binned by length.

reads/counted

Shape (n,) · Type uint64 · Fill zero

The number of reads passing all filters.

reads/rejected

Shape (n,) · Type uint64 · Fill zero

The number of reads rejected by at least one filter, or which couldn’t be modelled by the HMM

reads/unmapped

Shape () · Type uint64 · Fill zero

The number of reads not aligned to any reference.

pairwise/correlation

Shape (n, l, l) · Type float32 · Fill NaN · Written with --pairwise correlation

The Pearson correlation of mutations between this pair of bases.

pairwise/conditional

Shape (n, l, l) · Type float32 · Fill NaN · Written with --pairwise conditional

The probability that the base on the first axis was mutated in a read, given that the base on the second axis was.

pairwise/coverage

Shape (n, l, l) · Type float32 · Fill zero · Written with --pairwise

The number of reads in which this pair of bases was present.

CLI Options

Arguments

Argument	Description
BAM...	coordinate-sorted alignments

Input and output

Option	Description
-f, --fasta FASTA	reference sequences (required)
-o, --output HDF5	write results to this file (required)
--overwrite	replace the output file if it already exists
--verify CHECKS	identity checks to make against the FASTA (name, checksum, none; default name,checksum)

Filtering

Option	Description
-q, --min-mapq N	discard alignments below this mapping quality (0 to 254; default 20)
--min-length N	discard reads shorter than this (default: no limit)
--max-length N	discard reads longer than this (default: no limit)
-s, --strand STRANDS	keep alignments on these strands (forward, reverse; default forward,reverse)

Counting

Option	Description
--band N	reference positions the marginal may look either side of the CIGAR (default 2)
--pairwise STATS	write these statistics of how often two positions are modified together (correlation, conditional, none; default none)
--min-depth D	evidence a position needs before its rate is written (default 1)
--params FILE	read the pair HMM’s rates from this file (default: built in)
--substitution-weight W	what a substitution counts towards the mutation total (0 to 1; default 1)
--deletion-weight W	what a deletion counts towards the mutation total (0 to 1; default 1)
--insertion-weight W	what an insertion counts towards the mutation total (0 to 1; default 0)

Performance

Option	Description
-j, --workers N	threads running the processing step (default 1)
--decode-threads N	htslib threads for BGZF decompression (default 0)

Information

Option	Description
-h, --help	show this help and exit
-V, --version	show the version and exit
--dump-params	write the rates in the form –params reads and exit

Advanced

Accepted, and left out of --help.

Option	Description
--queue-capacity N	reads in transit at once (default 4096)
--batch N	reads transferred per queue operation (default 64)
--live-refs N	references in flight (default 64)
--dump-options	describe every argument as JSON and exit
--dump-layout	describe the output format as JSON and exit