Pure-Rust rewrite of the TrackCluster long-read isoform clustering/counting pipeline.
Goals:
- No runtime dependency on
bedtools(native sort/intersect/cluster primitives). - CLI parity with the legacy Python
trackrun.pysurface (in-progress).
Source checkouts pin Rust 1.90.0 via rust-toolchain.toml to avoid a known
EXDEV artifact-write failure seen with newer toolchains in this environment.
Implemented subcommands:
flow: one-command end-to-end pipeline (recommended)preparedir: split reads into per-gene folders (and write<prefix>_gene.txt,<prefix>_dedup.bed,<prefix>_novel.bed)clusterj: junction-chain clustering (fast mode; SL-aware 5' merge controls; optimized truncation collapsing for large loci)cluster: overlap-based clustering (slower, more permissive)count: isoform expression countingcount-multi: per-sample (and optional per-group) isoform usage from pooled isoformsdesc: novel isoform description/classification vs referenceaddgene: assign gene names to reads by overlap with referencevalidate-bed: strict BED12/bigGenePred input validation, with explicit lenient repair reportsbam2bigg: convert genome-aligned BAM records to TrackCluster bigGenePred-compatible BED12+8gff2bigg: convert GFF3 or GTF exon annotations to a TrackCluster reference catalogexport: write transcript catalogs as GTF, GFF3, or a SQANTI3 input-audit tablemod-import-m6anet: normalize m6Anet RNA002 per-read probabilities to genomic sitesmod-import-dorado: normalize one MM/ML modification code from a genome-aligned modBAMmod-subsample: generate synchronized technical coverage partitions and ready-to-run manifestsmod-aggregate: join normalized read-site observations to unique isoform assignmentsmod-site-summary: reduce complete isoform/site tables to a per-site QC inventorymod-contrast: compute descriptive isoform, condition, and interaction effect sizes
Extra binary:
clusterj_batch: runclusterjper gene folder in parallel (manual junction-mode batched runner; overlap-mode batching is exposed throughtrackcluster flow --cluster-mode cluster)
Download a tarball for your platform from the latest GitHub release:
# Example: Linux x86_64
REPO=lrslab/trackcluster-rs
TAG="$(curl -fsSL "https://api.github.com/repos/${REPO}/releases/latest" | sed -n 's/.*"tag_name": "\([^"]*\)".*/\1/p' | head -n1)"
ARCHIVE="trackcluster-${TAG}-x86_64-unknown-linux-musl"
curl -fLO "https://github.com/${REPO}/releases/download/${TAG}/${ARCHIVE}.tar.gz"
curl -fLO "https://github.com/${REPO}/releases/download/${TAG}/SHA256SUMS"
grep -F " ${ARCHIVE}.tar.gz" SHA256SUMS > "${ARCHIVE}.sha256"
test -s "${ARCHIVE}.sha256"
if command -v sha256sum >/dev/null 2>&1; then
sha256sum -c "${ARCHIVE}.sha256"
else
shasum -a 256 -c "${ARCHIVE}.sha256"
fi
# Supply-chain verification when GitHub CLI is installed:
if command -v gh >/dev/null 2>&1 && gh auth status >/dev/null 2>&1; then
gh attestation verify "${ARCHIVE}.tar.gz" --repo "${REPO}"
fi
tar xzf "${ARCHIVE}.tar.gz"
# Current archives may be flat; newer self-contained archives use one top-level directory.
if [ -d "${ARCHIVE}" ]; then cd "${ARCHIVE}"; fi
./trackcluster --help
# Make the unpacked binaries available to the quickstart commands below.
export PATH="$PWD:$PATH"Available targets: Linux x86_64 (musl static), Linux ARM64 (glibc 2.31+), macOS Apple Silicon.
cargo build --release
./target/release/trackcluster --helpcargo install --path . --locked --bins
trackcluster --help
clusterj_batch --help# One-line flow: prepare per-gene inputs, run per-gene clustering, merge outputs, count, and desc
trackcluster flow -s examples/reads.bed -r examples/ref.bed -o out --prefix sample
# Genes above 5,000 reads are deterministically subsampled by default so loci such as
# mitochondrial cox1 cannot dominate runtime. Set `--max-reads-per-gene 0` to disable it.
# For a tighter targeted cap: `--downsample-gene cox1 --max-reads-per-gene 2000`.
# Independent per-gene downsampling is rejected when one molecule belongs to multiple genes;
# disable the cap or exclude every affected gene from downsampling in that case.
# Malformed or empty-ID read tracks are skipped individually by default and recorded in
# `<prefix>_rejected_reads.tsv` / `<gene-path-key>/rejected_reads.tsv`.
# Add `--invalid-read-policy fail` to restore strict read-track parsing.
# Gene-local failures are logged and excluded while verified genes continue through merge/count/desc.
# Add `--strict-gene-errors` to stop before downstream outputs when any gene fails.
# If per-gene clustering already finished, rerun only merge/count/desc outputs
trackcluster flow --count-only -r examples/ref.bed -o out --prefix sample
# Count from an existing output folder; unique assignment stays inside each gene folder
trackcluster count -r examples/ref.bed -o out --prefix sample
# Validate a BED12/bigGenePred file
trackcluster validate-bed -i examples/minimal.bed
# Convert your own genome-aligned BAM to TrackCluster BED12+8. The default MAPQ cutoff is 30.
# (A BAM is not bundled with the tiny text examples.)
# Output BED score is 0 because no SL evidence is imported. MAPQ only filters records.
trackcluster bam2bigg --bamfile alignments.bam --out reads.bed
# Convert the packaged GFF3 model to a deterministic reference BED12+8 catalog.
trackcluster gff2bigg --gff examples/annotation.gff3 --out reference.bed
# Junction-mode clustering (writes isoform.bed + mapping + unused)
trackcluster clusterj -s examples/reads.bed -r examples/ref.bed -o isoform.bed
# Platform presets:
# --platform-preset rna002 # junction offset 15; SL 5' offsets 20/25/20; 3' cluster offset 15
# --platform-preset rna004 # conservative defaults: junction offset 10; SL 5' offsets 15/25/15; 3' cluster offset 10
# Junction-mode defaults treat reads as no-SL (`--sw-score -1`). Pass
# `--sw-score 11` only when BED score is valid SL/SW 5' evidence.
# SL 5' merge behavior can be tuned with --sl-partial-5prime-offset,
# --sl-same-junction-5prime-offset, --sl-5prime-cluster-offset, and
# --sl-5prime-min-support.
# Same-junction 3' retention can be tuned with --same-junction-3prime-offset,
# --3prime-cluster-offset, and --3prime-min-support.
# SL evidence is optional. Reads without SL information use the normal junction
# correction and 5' truncation collapse path, but are not SL-protected isoforms.
# Supported same-junction 3' terminal clusters are retained as isoforms. On the
# minus strand the 3' end is tx_start; an early stop has a higher tx_start than
# the corresponding full-length isoform.
# Overlap-mode clustering (legacy-style two-round exon/intron overlap mode)
trackcluster cluster -s examples/reads.bed -r examples/ref.bed -o isoform.bed
# Full flow in overlap mode
trackcluster flow --cluster-mode cluster -s examples/reads.bed -r examples/ref.bed -o out --prefix sample
# Flow keeps its shared no-SL default (`--sw-score -1`) in either clustering mode.
# Pass `--sw-score 11` to opt into legacy score-based protection in overlap mode.
# Legacy low-level count from a standalone isoform BED. Default unique mode also
# writes isoform_count.provenance.tsv; fractional mode does not.
trackcluster count -s examples/reads.bed -r examples/ref.bed -i isoform.bed --read-to-isoform isoform.read_to_isoform.tsv --out isoform_count.csv
# Describe/classify isoforms vs reference (writes <prefix>_*.txt)
trackcluster desc --isoform isoform.bed --reference examples/ref.bed -o desc_outUse a manifest TSV to pool reads for clustering once, then quantify per-sample isoform usage.
Example manifest (samples.tsv):
sample group reads
S1 control /path/S1.reads.bed
S2 treated /path/S2.reads.bedRun full pooled flow:
trackcluster flow --manifest examples/samples.tsv -r examples/ref.bed -o out --prefix pooledAdd --emit-pooled-reads if you also want <prefix>_pooled_reads.bed written.
If clustering already completed and you only need to regenerate merged count/description outputs, use --count-only. Include --manifest when you want the multi-sample usage tables regenerated too:
trackcluster flow --count-only --manifest examples/samples.tsv -r examples/ref.bed -o out --prefix pooledOr run per-sample quantification from an existing pooled isoform BED:
trackcluster count-multi --manifest examples/samples.tsv -r examples/ref.bed -i out/pooled_isoform.bed -o out/pooledTip: with default --name2-mode coverage (or none), use --read-to-isoform out/pooled_read_to_isoform.tsv (or keep the TSV next to the isoform BED for auto-discovery).
For overlap-mode pooled clustering, add --cluster-mode cluster to the flow command above.
count-multi writes:
out/pooled.isoform_count.csvout/pooled.isoform_usage.long.tsvout/pooled.isoform_counts.matrix.tsvout/pooled.isoform_usage.group.tsv(when at least one sample has a non-emptygroup)out/pooled.unique_assignment.provenance.tsv(default unique mode)
In unique assignment mode, flow also writes <prefix>_read_to_isoform.unique.tsv, the exact read-to-isoform mapping used for final counts, plus <prefix>_unique_assignment.provenance.tsv with the effective --unique-assignment-junction-offset and one-to-one/no-collapse matching policy. The raw merged <prefix>_read_to_isoform.tsv remains the unselected mapping from per-gene clustering.
The aggregate out/pooled.isoform_count.csv is derived from the per-sample matrix: each isoform count is the sum of that isoform's sample columns. In flow --manifest, the main <prefix>_isoform_count.csv is synchronized from the same aggregate count, so total and per-sample counts use the same assignment result.
Modification callers are normalized before they are joined to TrackCluster's final unique read-to-isoform assignments. This keeps caller/model/chemistry provenance separate and prevents missing observations from being interpreted as unmodified calls.
# Normalize caller output. See docs/CLI.md for all required provenance options.
trackcluster mod-import-dorado --sample S1 --assay-id dorado_rna004_m6a \
--bam S1.aligned.bam --mod-code A+a --model-id rna004_m6a_model \
--candidate-rule all-target-canonical-bases \
--source-emission-threshold 0.05 --out S1.mod
# For a motif-restricted model, use --candidate-rule DRACH; only matching
# read-sequence candidates enter the denominator.
# Aggregate against the exact mapping used for unique-mode expression counts.
trackcluster mod-aggregate --manifest samples.tsv \
--isoforms out/pooled_isoform.bed \
--read-to-isoform out/pooled_read_to_isoform.unique.tsv \
--mod-manifest mod_samples.tsv \
--reference-fasta genome.fa \
--analysis-threshold dorado_rna004_m6a=0.5 \
--eligibility-profile strict --out out/pooled
# Optional technical coverage test: split one parent sample by molecule.
trackcluster mod-subsample --manifest samples.tsv \
--read-to-isoform out/pooled_read_to_isoform.unique.tsv \
--mod-manifest mod_samples.tsv --source-sample S1 \
--replicates 4 --reads-per-sample 5000 --mode disjoint \
--out-dir out/S1_low_inputsflow --manifest ... --mod-manifest ... can run the aggregation as an optional
post-processing step. V1 reports read-derived counts, fractions, Wilson
intervals, and explicit effect-only contrasts; it does not treat reads as
biological replicates, and contrast p/q values are NA. Exact schemas,
denominator rules, and m6Anet examples are documented in
docs/FORMATS.md and docs/CLI.md.
Pinned public-data and realistic-simulation checks are described in
docs/MODIFICATION_VALIDATION.md.
mod-subsample pseudo-samples are technical coverage partitions only; their
groups are intentionally blank and they must not be used as biological
replicates.
The default modification eligibility profile is exploratory. Use strict
for comparison-ready screening: it requires exact coverage BAMs, an indexed
reference FASTA, Dorado version/model/threshold provenance verified within one
coherent source @PG record, and configurable minimum covering/callable counts
plus candidate/covering and callable/covering rates.
Flow-integrated modification results are committed as hash-verified generations
and are current only while <prefix>.mod.current.json exists and validates.
New catalogs use deterministic tc_novel_v1: structural IDs for novel
isoforms and a percent-encoded tc_name2_v1: payload in --name2-mode full.
Count CSVs have columns gene,isoform_id,count and use standard CSV escaping.
Repeated read labels are treated as one abundance molecule; conflicting
structures for one label are rejected in unique-assignment mode. See
docs/FORMATS.md for the identity and migration contract.
Within the 0.2.0 format contract, rejected-read reporting does not otherwise
change BED, isoform, count, or description/classification schemas and rules.
Skipped reads do not contribute biological evidence, so result contents can
change when an input contains rejected tracks. I/O, reference, configuration,
and algorithm errors are not downgraded by --invalid-read-policy skip.
- Changelog
- Pipeline tutorial
- CLI reference
- File formats
- Modification validation
- Interchange formats
- Rust API policy
- Clustering behavior
- Description/classification behavior
These commands require the repository's source and test fixtures; they are not included as runnable inputs in pre-built binary archives.
cargo test --all --all-featuresJunction-cluster and count golden fixtures:
# Regenerate the clusterj and count goldens from the current Rust implementation
bash tests/generate_goldens.shLicensed under either of:
- MIT license (
LICENSE-MIT) - Apache License, Version 2.0 (
LICENSE-APACHE)
at your option.