Use cases and analysis guidelines#

What to run, and what the numbers mean afterwards. Each section is a real question with a configuration that answers it and a read-out that does not overstate what was measured.

Bulk RNA-seq: extract a repertoire from a transcriptome #

The primary case. A bulk library is 0.02–3 % receptor, so almost all the work is proving that reads are not receptor reads — which is what the mode’s --prefilter removes (1.99×, at a measured ~0.15 % of mapped reads; --exact turns it off). --ec-mode rnaseq is this mode’s denoising default.

arda rnaseq --r1 R1.fq.gz --r2 R2.fq.gz -d out -p sample --threads 16

Read-outs, and what each is worth:

clones.tsv: One row per clonotype. duplicate_count is reads encompassing the junction (coverage, the true expression estimate); consensus_count is distinct fragments. Use consensus_count when you need molecules, duplicate_count when you need expression.
arda.json: mapped_reads / total_reads is the receptor fraction. On bulk that is normally 0.02–3 %; a value far above it on a non-lymphoid tissue is a reason to look at the input, not a result.

Warning

Depth is not a nicety here. A 200,000-read subsample of a bulk library yields a handful of clonotypes — measured: 1 for GM12878, 3 for Jurkat — and no diversity, overlap or expansion statistic computed on that means anything. Full depth on the same GM12878 library gives a real repertoire. If you must subsample, say so beside every number.

Targeted amplicon / RepSeq #

arda amplicon --r1 R1.fq.gz --r2 R2.fq.gz -d out -p sample --threads 16

The mode carries --two-pass --fast-segments --v-only-on-segment and --ec-mode amplicon. Those are not interchangeable with the bulk ones and do not compose, which is why the regime is the command name; see Running on a cluster (SLURM). On an IGH RepSeq amplicon the amplicon mode is 4.15× faster than --exact at 2.7× less memory.

Monoclonal QC: is my cell line what it says it is?#

A monoclonal line has one productive rearrangement per expressed locus, so the read-out is purity — the dominant clonotype’s share of its locus’s reads — and, if the clone is published, how many reads land on it.

arda rnaseq --r1 R1.fq.gz -d out -p Jurkat \
    --ec-mode amplicon --clonotype-key junction

Measured on Jurkat (14,531 junction-bearing reads):

configuration	clonotypes	reads	TRB purity	reads on the 2 published clones
`fast` (default)	90	14,531	.98963	14,177 (.9756)
`amplicon` + `junction`	10	14,531	.99990	14,313 (.9850)

Important

A clonotype count falling is the job; a read count falling is a bug. Every arda denoising mode moves reads onto a parent and never discards them, so duplicate_count summed over the table is invariant. If it drops when you change --ec-mode, that is a defect — report it.

Warning

Do not reach for “drop clonotypes with 1 read” instead. Measured on the same library it loses 82 reads outright and plateaus at .99398 purity, because the largest error class by reads is call splits carrying 33 and 30 reads — no abundance rule separates those. Raising the threshold from 1 to 2 buys 0.00000 purity and loses 2 more reads.

Negative controls #

K562 and HepG2 carry no V(D)J recombinase, so any clonotype is a false positive by construction. Run them exactly as your real samples and expect zero.

Note

Passing them validates cross-lineage specificity and nothing else. In this project four different reference configurations — including one with a --min-score low enough to take precision from 94.3 % to 65.5 % — all yielded zero clonotypes on both. A negative control that passes is not evidence that a threshold is right.

Low-frequency variants: spike-ins, MRD, minor clones #

arda correct -i s1.airr.tsv -o clones.tsv --ec-mode accurate --error-rate 1e-5

Two things decide whether a rare real variant survives:

--error-rate: At the shipped 1e-3 the abundance model erases both published MIGEC spike-in variants. At 1e-5 both are recovered exactly. ⛔ This is a per-library calibration, not a default change: the right value depends on the library’s error-cloud abundance ratio.
--ec-mode accurate: Adds a Phred gate on the base that discriminates a clonotype from its parent — evidence the abundance model does not have. A real low-frequency variant is a good read of a rare molecule: the MIGEC variants sit at median mean junction Phred 37.2–37.5 against the parent’s 37.8.

Warning

At signal-to-noise ≈ 1 no abundance threshold can separate a real variant from PCR error. On the MIGEC data the published V2 variant is less abundant than the worst 2-substitution PCR error (V2/Err2 = 0.28 on raw reads). That is not a tool limitation; it is why UMI consensus exists. If your variant is at that level, the answer is chemistry, not parameters.

Somatic hypermutation #

v_mutations / j_mutations record SHM in germline coordinates; see Somatic hypermutation.

Warning

Local zero-loss is not zero-loss. Three times in this project a clean result turned out to be a property of the test library’s SHM level. Validate anything SHM-sensitive on a hypermutated library (median V identity ~92 %), not on a naive one — and report both, because a change that helps a hypermutated repertoire can cost a naive one.

Comparing arda against another tool #

Important

Name the stage. TRUST4’s read and false-positive counts are scored at candidate extraction; arda’s are its final post---min-score clonotype table. Those are different quantities and differencing them produces a number that means nothing.

Benchmark every tool at its best configuration. MiXCR’s shipped rna-seq preset discards reads it aligned (minSumScore=200), understating its recall roughly ninefold. Measured on a 200 k Jurkat subsample: 154 aligned reads at the default preset, 5,136 with -OallowNoCDR3PartAlignments=true -OminSumScore=40, at the same wall time.

Compare at gene level across tools. MiXCR suffixes every allele *00, so an allele-level comparison scores it zero by construction. And treat an ambiguous call as an abstention: comparing only the first element of a comma-joined tie list scores two tools that agree the answer is ambiguous between the same genes as a disagreement.

Give each call metric its own denominator. j_gene agreement divided by all truth reads measures how many reads reach the J — a property of the library. On one IGH library that read .4678; over the reads where both the truth and the tool named a J it is .9418.

The one test that needs no external truth #

If you have the same material sequenced two ways — a deep targeted amplicon and bulk RNA-seq — the strongest validation available is that they agree:

what fraction of the amplicon repertoire is recovered from RNA-seq, weighted by amplicon abundance (did the tool find the clones that are actually there?);
rank correlation of abundances on the shared clones — Spearman, not Pearson, because the two protocols have completely different amplification biases;
the RNA-seq-only rate, as a candidate false-positive rate and never as an FP count: the amplicon covers only its target loci and misses low expressers.