Recipes#
Runnable snippets for the things people actually do after (and around) a run. Every command and
every script on this page was executed against tests/data/rnaseq_real before it was written
down; the printed numbers are that fixture’s, not illustrations.
Usage is the reference for what each flag means, Samples split across files for multi-file samples, Use cases and analysis guidelines for choosing a mode and a denoising preset.
FASTQ to clonotypes#
One command per library type. The mode name picks the speed configuration — the two speed levers do not compose, and each is a loss in the other’s regime, so do not hand-assemble them.
# bulk RNA-seq / WTS: a few receptor reads in a lot of transcriptome
arda rnaseq --r1 R1.fq.gz --r2 R2.fq.gz --out-prefix PT01 -d results/ --threads 16
# targeted RepSeq / 5'RACE amplicon: nearly every read is a receptor
arda amplicon --r1 R1.fq.gz --r2 R2.fq.gz --out-prefix PT01 -d results/ --threads 16
Four files come out, and the paths are echoed on stdout (one per line) while progress goes to
stderr — so $(arda rnaseq ... | head -1) is a usable idiom in a shell script:
results/PT01.airr.tsv one AIRR Rearrangement row per mapped read
results/PT01.clones.tsv one row per clonotype, with duplicate_count
results/PT01.arda.json the run report: what was read, what mapped, wall time, peak RSS
results/PT01.stats.tsv run QC in long format (scope / key / metric / value)
On the 1,320-read test fixture that is 453/1,320 reads mapped (34.32 %) and 49 clonotypes.
Single cell is a different entry point, because it starts from UMI consensus per molecule rather than raw reads — demultiplexing, barcode correction and UMI collapse belong upstream:
arda cells asm/PBMC.consensus.fq.gz -p results/PBMC --cells ref/PBMC.cells.tsv --plot svg
The upstream migec commands, what each output table holds and how to read the QC plots:
Single cell.
Annotate sequences you already have#
No reads, just sequences — assembled contigs, Sanger, a synthesised panel, someone else’s
consensus. arda annotate streams FASTA/FASTQ to AIRR and is memory-bounded:
arda annotate -i contigs.fasta -o contigs.airr.tsv --organism human
cut -f1,3,4,7,55 contigs.airr.tsv | head -3
sequence_id locus v_call j_call junction_aa
PZ235980.1 IGH IGHV3-9*01 IGHJ6*03 CARDIGAGGFGDNFYFFYYMDVW
PV083657.1 IGK IGKV1-33*01,IGKV1D-33*01 IGKJ5*01 CQQYDSLPYTF
Bare (junction_aa, V, J) records — a VDJdb-style table, a published supplement — go through
arda markup, which marks up the junction and repairs an anchor the record lost:
arda markup -i junctions.tsv -o fixed.tsv --id-col id
# 7 records -> fixed.tsv (5 repaired, 1 failed)
FLVGPQGSSASKIIF comes back as CLVGPQGSSASKIIF (the Cys104 anchor restored) and
CAIRDDKII as CAIRDDKIIF. Add --vdjdb to read VDJdb’s own column names.
From Python#
The library entry point takes sequences and returns AIRR record dicts — no files, no subprocess:
from arda import annotate_sequences
rows = annotate_sequences(
[("q1", "TGTGCCAGCAGCTTAGCGGGAGGGAACACCGGGGAGCTGTTTTTTGGA")], seqtype="nt"
)
print({k: rows[0][k] for k in ("locus", "v_call", "j_call", "junction_aa")})
# {'locus': 'TRB', 'v_call': 'TRBV7-2*04', 'j_call': 'TRBJ2-2*01',
# 'junction_aa': 'CASSLAGGNTGELFF'}
seqtype="aa" annotates amino-acid input. For a whole FASTQ use the CLI or
arda.rnaseq.pipeline.run() — annotate_sequences holds its batch in memory.
A cohort from one sheet, and samples split across lanes#
Lanes and in-house chunks are read groups of one sample: one sample still gives one clonotype
table. Repeat --r1/--r2 and name the sample with --id, or list them in a sheet whose
columns are nf-core’s, so an existing nf-core samplesheet works unmodified:
cat > sheet.tsv <<'SHEET'
sample fastq_1 fastq_2
PT01 PT01_S1_L001_R1_001.fastq.gz PT01_S1_L001_R2_001.fastq.gz
PT01 PT01_S1_L002_R1_001.fastq.gz PT01_S1_L002_R2_001.fastq.gz
PT02 PT02_S2_L001_R1_001.fastq.gz PT02_S2_L001_R2_001.fastq.gz
SHEET
arda rnaseq --samples sheet.tsv -d results/ --threads 16 # one box, samples in turn
arda cluster submit-samples --samples sheet.tsv --work-dir work/ -d results/ \
--regime rnaseq --threads 16 --partition medium --submit # SLURM, one task per read group
Repeated sample values merge in row order. Do not cat the lanes first: arda maps each
read group and concatenates after Stage 1, which is byte-identical to the same reads in one file
and skips a full copy of the data. Samples split across files has the command-line form, the ordering rule and
the one-worker-per-read-group recipe.
Analysing the clonotype table#
Important
Read arda’s TSVs with quoting off. Every arda writer uses quote_style="never" and every
reader quote_char=None. Read one with polars’ default quoting and a blank
chimera_parents (written by arda correct --flag-chimeras) becomes the two-character value
"" — every clonotype then reads as chimeric.
import polars as pl
READ = dict(separator="\t", quote_char=None)
clones = pl.read_csv("results/PT01.clones.tsv", **READ)
Top clonotypes within a locus, with their clonal fraction. The fraction is per locus, because a TRB frequency computed over a table that also holds IGK is not a TRB frequency:
top = (
clones.filter(pl.col("locus") == "TRB")
.with_columns(
(pl.col("duplicate_count") / pl.col("duplicate_count").sum()).alias("frequency")
)
.sort("duplicate_count", descending=True)
.select("junction_aa", "v_call", "j_call", "duplicate_count", "frequency")
.head(5)
)
┌─────────────────┬─────────────┬────────────┬─────────────────┬───────────┐
│ junction_aa ┆ v_call ┆ j_call ┆ duplicate_count ┆ frequency │
╞═════════════════╪═════════════╪════════════╪═════════════════╪═══════════╡
│ CSQSGGFGADTQYF ┆ TRBV29-1*03 ┆ TRBJ2-3*01 ┆ 2 ┆ 0.181818 │
│ CSATPPDSWTGELFF ┆ TRBV20-1*07 ┆ TRBJ2-2*01 ┆ 2 ┆ 0.181818 │
│ CASSLRGSYEQYF ┆ TRBV7-8*01 ┆ TRBJ2-7*01 ┆ 2 ┆ 0.181818 │
└─────────────────┴─────────────┴────────────┴─────────────────┴───────────┘
Which chains the library actually carries, in clonotypes and in reads:
chains = (
clones.group_by("locus")
.agg(pl.len().alias("clonotypes"), pl.col("duplicate_count").sum().alias("reads"))
.sort("reads", descending=True)
)
# IGL 17/23, IGK 14/19, TRB 8/11, TRA 6/8, IGH 4/6 on the test fixture
V-gene usage as a frequency, gene not allele. Split on * rather than trusting a gene
column: v_call is allele-level by default and can be a comma-separated tie:
usage = (
clones.filter(pl.col("locus") == "IGH")
.with_columns(pl.col("v_call").str.split("*").list.first().alias("v_gene"))
.group_by("v_gene")
.agg(pl.col("duplicate_count").sum().alias("reads"))
.with_columns((pl.col("reads") / pl.col("reads").sum()).alias("frequency"))
.sort("reads", descending=True)
)
Run arda rnaseq --call-level gene if you want the collapse done upstream, at clonotype-key
level, instead.
Overlap between two repertoires, on the amino-acid junction within one locus:
a = pl.read_csv("results/PT01.clones.tsv", **READ).filter(pl.col("locus") == "TRB")
b = pl.read_csv("results/PT02.clones.tsv", **READ).filter(pl.col("locus") == "TRB")
shared = a.join(b.select("junction_aa"), on="junction_aa")
f = shared.height / min(a.height, b.height) # overlap as a fraction of the smaller set
Note
junction_aa includes the Cys104 and [FW]118 anchors; IMGT cdr3_aa excludes them and is
two residues shorter. Joining one table’s junction_aa to another’s cdr3_aa silently
finds nothing. MiXCR’s CDR3 column is AIRR junction.
Somatic hypermutation per read, from the AIRR table — v_mutations is a comma-separated
list of the mutations found outside the junction, so its length is the count:
airr = pl.read_csv("results/PT01.airr.tsv", **READ)
shm = (
airr.filter(pl.col("locus").str.starts_with("IG"))
.with_columns(
pl.when(pl.col("v_mutations").is_null()).then(0)
.otherwise(pl.col("v_mutations").str.split(",").list.len())
.alias("v_mutation_count")
)
.group_by("locus")
.agg(pl.col("v_mutation_count").mean().round(2).alias("mean_v_mutations"))
)
Somatic hypermutation explains what is and is not counted (the junction is excluded on purpose) and how to
recount an existing AIRR file with arda shm.
Gate a run on its own report#
Do not parse the log. <prefix>.arda.json carries the same numbers as structured data, and a
pipeline should assert on them:
import json
rep = json.loads(open("results/PT01.arda.json").read())
m = rep["map"]
print(f"mapped {m['mapped_fraction']:.3f} of {m['total_reads']} reads, "
f"{rep['correct']['clonotypes_out']} clonotypes, {rep['wall_seconds']} s")
assert m["mapped_fraction"] > 0.01, "nothing mapped — wrong organism or wrong reference"
assert rep["correct"]["clonotypes_out"] > 0, "no clonotypes"
mapped 0.343 of 1320 reads, 49 clonotypes, 2.47 s
A sample that arrived as several read groups adds read_groups and reports
wall_seconds_max / wall_seconds_sum instead of a single wall_seconds for the map
stage: summing forty array tasks’ wall time and calling it “wall seconds” would be a lie.
reads_per_second, peak_rss_mb-style memory and every count are present either way.
<prefix>.stats.tsv is the same run in long format, plus per-chain and per-gene QC —
scope/key/metric/value, which pivots:
qc = pl.read_csv("results/PT01.stats.tsv", **READ)
per_chain = qc.filter(pl.col("scope") == "chain").pivot(
on="metric", index="key", values="value"
)
Quality control lists every scope and metric it writes.
Cohort QC in one table, and one page#
Every run already wrote its own QC table. arda qc batch joins them — reading only those
tables, never the AIRR or the clonotypes, so this stays cheap on a cohort of any size:
arda qc batch -d results/ -o results/cohort --samples sheet.tsv
[arda] qc batch: 6 samples, 2205 rows, 41 flagged
results/cohort.qc.tsv
results/cohort.qc.wide.tsv
results/cohort.qc.json
--samples is read only for the sheet’s optional project and batch columns: they name
the group each sample is compared within. cohort.qc.wide.tsv is one row per sample:
$ cut -f1,2,3 results/cohort.qc.wide.tsv | head -4
sample project batch
S0 TRIAL9 RUN1
S1 TRIAL9 RUN1
S2 TRIAL9 RUN1
Nothing is filtered and no threshold is shipped. What the long table adds is each metric’s group median, its MAD and a robust z, so the question “is this sample like its batch” has an answer:
import polars as pl
long = pl.read_csv("results/cohort.qc.tsv", separator="\t",
infer_schema_length=0, quote_char=None)
flagged = long.filter((pl.col("outlier") == "1") & (pl.col("scope") == "sample"))
print(flagged.select("sample", "metric", "value", "median", "z"))
┌────────┬───────────────────────────┬───────────┬───────────┬────────────┐
│ sample ┆ metric ┆ value ┆ median ┆ z │
╞════════╪═══════════════════════════╪═══════════╪═══════════╪════════════╡
│ S0 ┆ clonotype_junction_aa_max ┆ 21 ┆ 16.000000 ┆ 6.745000 │
│ S4 ┆ clonotype_junction_aa_max ┆ 21 ┆ 16.000000 ┆ 6.745000 │
│ S2 ┆ j_gene_coverage_reads ┆ 0.0485437 ┆ 0.121359 ┆ -10.117264 │
│ S3 ┆ j_gene_coverage_reads ┆ 0.165049 ┆ 0.121359 ┆ 6.070361 │
└────────┴───────────────────────────┴───────────┴───────────┴────────────┘
S2’s J-gene coverage is a third of its cohort’s, which is where to look first: an order of
magnitude below the batch usually means the wrong organism, the wrong regime, or a library that
did not work. Note these six “samples” are contiguous cuts of one small fixture, so the spread
is the fixture’s, not a real cohort’s.
Finally, the whole thing as one file you can send someone:
arda qc report -i results/cohort.qc.json -o results/cohort.qc.html
It inlines its own data and fetches nothing, so it opens on an air-gapped login node or as an
email attachment long after results/ is gone. See Quality control.
Export the reference for a genome browser#
The markup arda transfers is itself exportable — every in-frame V·J scaffold with IgBLAST-quality FR1–4 / CDR1–3 coordinates, in three kinds × four formats:
arda export-ref --kind scaffolds --locus TRB --format gff3 -o trb.gff3
arda export-ref --kind alleles --locus IGH --format tsv -o igh_alleles.tsv
Coordinates are 1-based closed, as AIRR and GFF3 both are, so they load unshifted. Exporting the reference covers the kinds, the formats and the round-trip guarantee.