Step 3 Filtering and Quality Control

Step 3: Filtering & Quality Control

What is this?

Every step that can drop a sequence, or a whole specimen, between the sequencing run and the published tree: the exact rule each step applies, where the evidence is written, and how to change a filter safely.

Why does it matter?

Low-quality or mislabelled sequences distort trees, dates and clade frequencies. Filters keep them out, and you need to know which filter removed a sample when someone asks.

What do I do?

Nothing during a normal run: filtering is automatic. Use this page to explain a missing sample (Where did my sample go?), to read the report's QC section (Finding removals in the report), or before proposing a change (Changing a filter).

flowchart TD
    A["Run FASTAs from the sequencing team<br/>(already passed assembly QC)"] --> B["① prepare_source.py<br/>bad records, controls,<br/>specimen in two runs"]
    B --> C["② flusort + fludb upload<br/>untyped, no sequence,<br/>bad vaccine headers"]
    C --> D["③ download from fludb<br/>subtype not mapped to a build<br/>(e.g. xxxx)<br/>genome: not all 8 segments"]
    D --> E["④ Nextclade<br/>clade + QC per segment"]
    E --> F{"⑤ augur filter<br/>QC + coverage, length,<br/>exclude.tsv"}
    F -->|removed| X["results/…/filter_log.tsv"]
    F -->|kept| G{"⑥ augur refine<br/>clock filter<br/>(segment builds only)"}
    G -->|outlier| Y["Pruned from the tree"]
    G -->|kept| H["In the build"]

Loading

Sequences can be dropped at six steps. Steps ① and ② happen before the surveillance report's "submitted" count (① entirely, ② partly), so their evidence is in the prep report and run logs. Steps ③ to ⑥ are counted in the report's QC section.

All filters at a glance

#	Step (rule or script)	Applies to	Drops	Evidence
①a	scripts/prepare_source.py	JHH	Malformed records: bad name, segment not 1-8, empty, non-DNA characters	source_prep/output/prep_report.txt
①b	scripts/prepare_source.py	JHH	Negative/positive controls	prep_report.txt
①c	scripts/prepare_source.py	JHH	The other run's copy of a specimen sequenced in more than one run	prep_report.txt (NOTE)
②a	flusort → upload_genomes	JHH	Specimens with no BLAST hit on any segment (untyped)	.run/logs/upload_genomes/upload_genomes.log (counts only); report: No sequence data
②b	upload_genomes (--require-sequence)	JHH	Metadata rows with no stored segment	same
②c	upload_vaccines	Vaccines	Headers it cannot parse or map	.run/logs/upload_vaccines/upload_vaccines.log
③a	download_segments	All	Subtype not listed in config/subtype_assignment.yaml (e.g. xxxx, mixed, Yamagata)	source/flusort_JHH_metadata.tsv; report funnel
③b	download_genomes	Genome builds	Fewer than 8 segments stored	Report funnel
④	merge_quality_metrics	Segment builds	Nothing directly, but a specimen with no HA gets no QC values and fails ⑤ in every segment build	.run/logs/merge_quality_metrics/
⑤a	augur_filter	Segment builds	Nextclade QC bad or missing, coverage < 0.9	results/{subtype}/{segment}/filter_log.tsv
⑤b	augur_filter	Segment builds	Shorter than the segment minimum	same
⑤c	augur_filter, genome_filter	All builds	IDs in config/exclude.tsv	same
⑤d	genome_filter	Genome builds	Under 13,000 nt; no matching sequence (all vaccine strains, see below)	results/{subtype}/genome/filter_log.tsv
⑥	refine	Segment builds	Clock outliers (> 4 IQD)	Not in filter_log.tsv; report: Clock-filter outlier

Stage 1: preparing the source files (prepare_source.py)

Everything here happens on your computer before the pipeline runs (Step 1). Every decision is printed and saved in source_prep/output/prep_report.txt, which is overwritten each time the script runs.

Records skipped (WARNING), the rest of the file is used:

A record is skipped when...	Rule in the code
Its name is not <ID>_<segment number>	ID may contain only letters, digits, . and - (no _). Text after the first space or | is ignored; lower-case jh is fixed to JH
The segment number is not 1-8	SEGMENTS
It is empty after removing gaps (-, .) and spaces
It contains characters that are not IUPAC bases	Allowed: ACGTRYSWKMBDHVN

A whole run file stops the script (RESULT: STOPPED, nothing written) if it is empty, binary, not a FASTA, or more than half its records are unusable.

Decisions the script makes for you:

• Same segment twice in one run file. Identical copies are merged. Different copies: the one with more non-N bases is kept (the first one if tied), with a WARNING.
• Controls are never written. A name is a control if it is not JH + digits and contains NTC, NEG, BLANK, WATER, H2O, CTRL, CONTROL, PTC or POS (any case). A control segment of 200+ bases raises POSSIBLE CONTAMINATION. A control given a JH number is not recognised; list it in config/exclude.tsv.
• Same specimen ID in more than one run. One run is kept: the run with the most segments, or the higher-numbered run if tied. The other run's segments are discarded; segments are never mixed between runs. This is a NOTE, not a warning.

Warning

Specimens are matched by name only. If two different specimens share an ID (a typo, a reused name), the script treats them as one specimen sequenced twice and silently keeps only one. Check the "sequenced in more than one run" NOTE against the run sheets. Non-JH IDs also get a WARNING.

Not dropped here: short or high-N segments (only counted; dropped at ⑤) and IDs in config/exclude.tsv (only listed; dropped at ⑤c). source/vaccine.fasta is not checked by this script (Tutorial 3).

Stage 2: during ingest (ingest.smk)

Typing with flusort

flusort BLASTs every segment (blastn, e-value ≤ 20, best bitscore per segment) against scripts/flusort/blast_database/pyflute_ha_database, then types each specimen:

• Type: Influenza A if any segment hit an A segment, otherwise B. If no segment has a hit, the specimen has no type.
• Subtype: the H and N found among the hits; a missing one becomes xx (e.g. H3xx, xxN2). More than one H or N gives a mixed value such as H1, H3N2.
• Influenza B is always labelled Victoria. flusort does not tell B lineages apart, so a Yamagata specimen would go into the B/Victoria builds; look for it at review.

Upload to fludb

Script	Skips	Notes
upload_jhh.py (rule upload_genomes)	Records whose specimen has no type (②a); names that do not split into exactly <ID>_<segment>; records with no metadata row. Then metadata rows with no stored segment (--require-sequence, ②b)	The log prints counts, not IDs. Such specimens show in the report as No sequence data
upload_vaccine.py (rule upload_vaccines)	Headers with fewer than 7 | fields, a segment not 1-8, or a type/lineage it cannot map (②c). Mapped: A/H1 with lineage pdm09 → H1N1; A/H3 → H3N2; B Victoria; B Yamagata	Yamagata is stored but has no build. Every skipped header is printed in the log

Caution

fludb keys every record on sequence_ID. Duplicate IDs are merged, not dropped: each segment takes the last record's sequence, and the metadata takes the last non-empty value. prepare_source.py guarantees unique IDs; hand-edited source/ files do not, and can silently create one "genome" from two specimens. Always build source/ with the script.

Download from fludb

Rule	Keeps	Effect
download_segments	subtype listed for the build in config/subtype_assignment.yaml (③a): h1n1 = H1N1, H1xx, xxN1; h3n2 = H3N2, H3xx, xxN2; vic = Victoria	xxxx, mixed, untyped and Yamagata specimens are in no build. Incomplete genomes are kept: each segment build gets every specimen with that segment
download_genomes	Same subtypes, and all 8 segments stored (--complete-genomes, ③b)	Partial genomes are in segment builds only

Both rules run scripts/fludb_download_builds.py for their {subtype} (--segments or --genome).

Assumed subtypes. A specimen with only HA or only NA called (H1xx, xxN1, H3xx, xxN2) is assumed to be H1N1 or H3N2, because these are the only influenza A subtypes circulating seasonally in humans. It goes into that subtype's builds (without the segment that was not called), and the report labels it subtype assumed and explains the assumption in its Methods. flusort and fludb keep the original value (e.g. xxN1), and so does the build metadata. A specimen with neither called (xxxx) is not built.

Warning

If a reassortant (e.g. H1N2 or H3N1) starts circulating, a lone HA or NA call no longer identifies the subtype. Remove the affected values from assumed in config/subtype_assignment.yaml, or add a build (Tier 2).

Stage 3: segment builds (segments.smk)

④ Clades and QC with Nextclade

1. Datasets. Each run downloads the latest Nextclade HA and NA datasets (rule fetch_hana_datasets; git-ignored, their tags go into the run's provenance file); the other six segments use the datasets stored in the repo under nextclade/flu/{subtype}/{segment}/.
2. nextclade runs on every segment.
3. assign_clades copies the HA clade and subclade onto every segment's metadata (matched on sample_ID).
4. merge_quality_metrics adds that segment's own coverage, qc_overallScore and qc_overallStatus.

Important

Step 4 finds each segment's QC row through the specimen's HA Nextclade row. A specimen with no HA sequence gets no QC values for any segment, so all its segments fail ⑤a. The report shows this as No HA segment, or per sequence as Failed Nextclade QC: no Nextclade result. The same reason appears when Nextclade could not align a segment.

Why Nextclade instead of augur clades?

This differs from the standard Nextstrain workflow and is open for discussion. Nextclade offers:

• Comprehensive QC metrics (details): missing data, problematic sites, private mutations, mutation clusters, frameshifts, stop codons. The overall score is labelled good, mediocre or bad from these.
• Metadata-centric workflow: clades are stored in the metadata, so no separate augur clades step is needed.
• Built-in glycosylation prediction.
• Flexible filtering: QC metrics can be filtered with augur filter and inform judgement during analysis.

See the augur clades documentation for the alternative.

⑤ augur_filter

A sequence is kept only if all of these hold:

Filter	Value	Set in	filter_log.tsv name
Coverage	coverage >= 0.9	--query in segments.smk	filter_by_query
QC status	qc_overallStatus is good or mediocre. bad and missing fail	--query in segments.smk	filter_by_query
Minimum length	PB2 2000, PB1 2000, PA 1800, HA 1400, NP 1200, NA 1200, MP 700, NS 700. augur counts only A, C, G and T, so Ns do not count	min_lengths in the Snakefile	filter_by_min_length
Not excluded	ID not in config/exclude.tsv, matched on sample_ID	config/exclude.tsv	filter_by_exclude
Has a sequence	Metadata row has a matching FASTA record	built into augur	filter_by_sequence_ids

Each removed strain is logged once, under the first filter that removed it. The log columns are strain, filter and kwargs (the threshold used).

⑥ refine (clock filter)

augur refine --clock-filter-iqd 4 prunes tips whose root-to-tip divergence is more than 4 interquartile distances from the molecular-clock expectation. It uses fixed rates per subtype and segment (clock_rate() in segments.smk; 0.001 if a pair is missing; standard deviation = rate / 5). These removals are not in filter_log.tsv: a tip in results/{subtype}/{segment}/filtered.tsv but missing from tree.nwk was clock-filtered. The report works this out for you (Clock-filter outlier).

After refine, the remaining steps (traits, ancestral, translate, export, frequencies) remove nothing. A tip in the tree but missing from the Auspice JSON is reported as Dropped at export and should never happen.

Stage 3: genome builds (genomes.smk)

• genome_clades joins the HA Nextclade calls to the genome metadata. It removes nothing.
• genome_filter keeps genomes with at least 13,000 A/C/G/T across the 8 concatenated segments (min_genome_size in genomes.smk) that are not in config/exclude.tsv. No Nextclade QC filter is applied to genomes.
• genome_refine sets no clock filter.

Note

Vaccine strains never reach a genome tree. The genome FASTA is labelled with sequence_ID (the GISAID EPI_ISL_… ID for vaccines), but augur matches the metadata on sample_ID (the isolate name). JHH specimens are unaffected because their two IDs are identical. Vaccines are logged as filter_by_sequence_ids (9 in the 2026-10-08 production snapshot) and counted as No sequence data. Fixing it changes the download step, which is Tier 3 (Best Practices).

Finding removals in the report

Open reports/surveillance_report.html (confidential; see Step 4). Each question below is answered by a section of the report. Every table is also saved in reports/tables/ and reports/report.xlsx.

Question	Report section → item	Table file
How many specimens were lost at each step, and why?	Quality control: run-level funnel → Run-level funnel: in, removed, and in the final build (tabs By subtype / lineage and By type; run menu under the chart) and the table below it	qc_funnel_table.tsv
How did each run do?	Quality control: per sequencing run → Where each run's specimens ended up, QC summary per sequencing run, Run detail: specimens retained in each build	qc_runs.tsv
Which build loses most, and at which step?	Quality control: per build → Retention in each Nextstrain build, Where sequences were removed, by build, Stage counts per build	qc_builds.tsv
Is QC or coverage the problem?	Quality control: Nextclade metrics → Nextclade QC status by segment, HA coverage vs Nextclade QC score (dashed line = coverage threshold)	–
Why was this sequence removed?	Quality control: detailed tables → Every sequence removed from a build, with the reason (search by specimen, run or reason)	qc_removed.tsv
What happened to this specimen in every build?	Quality control: detailed tables → Specimen journey: status in each build	qc_specimens.tsv

The funnel's steps are: Submitted → Sequence data in fludb → Subtype/lineage assigned → HA in build input → HA passed augur filter → HA in final build, then Whole genome sequenced (8 segments) → Genome in final build.

What each reason means:

Report label (short code in Specimen journey)	Meaning	Step
No sequence data (no seq)	In the input metadata but not in fludb; or, in a build, no matching sequence	②a/②b, ⑤d
Subtype/lineage not assigned	flusort subtype not mapped to a build in config/subtype_assignment.yaml (e.g. xxxx)	③a
No HA segment	No HA sequence, so no HA build and no QC for other segments	③a, ④
Failed Nextclade QC / coverage (QC)	The reason says which: coverage 0.85 < 0.9, QC status 'bad', or no Nextclade result	⑤a
Shorter than minimum length / HA too short (short)	Below the segment or genome minimum	⑤b, ⑤d
Manually excluded (exclude.tsv) (excl)	Listed in config/exclude.tsv	⑤c
Removed by augur filter (other) (filter)	Any other augur filter	⑤
Clock-filter outlier (clock)	Passed the filter, pruned by refine	⑥
Dropped at export (export)	In the tree, missing from the Auspice JSON	export
Segment not sequenced (–)	Not a removal: the segment never existed	–

Not in the report:

• Stage 1 removals: read source_prep/output/prep_report.txt.
• Controls: hidden using control_pattern in config/report.yaml.
• Skipped vaccine headers: read .run/logs/upload_vaccines/upload_vaccines.log.

The report explains QC failures using qc_min_coverage and qc_pass_status in config/report.yaml. If those differ from the --query in segments.smk, the stated reasons will be wrong.

Where did my sample go?

1. Is it in source/JHH_metadata.txt? If not, read prep_report.txt (control? other run kept? record skipped?) and Step 1.
2. Search the report's Specimen journey table for its ID. The Outcome column and one status per build usually answer the question. No report? Continue below.
3. What subtype did flusort give it (source/flusort_JHH_metadata.tsv)? Not exactly H1N1/H3N2/Victoria means no build; an empty type means no BLAST hit.
4. Is it in data/{subtype}/{segment}/metadata.tsv? (Genome build: all 8 segments needed.) Does it have an HA? Without HA it fails QC in every segment build.
5. Search results/{subtype}/{segment}/filter_log.tsv for its ID: the filter column says why.
6. In filtered.tsv but not in the tree: clock-filtered at refine.

Changing a filter

Filters decide what the lab sees, so changes are controlled. Read Best Practices → Change tiers first:

• Tier 2 (maintainer approves before you start): thresholds, config/exclude.tsv, clock filter and rates, Nextclade datasets.
• Tier 3 (as Tier 2, plus the maintainer reviews the code): the rules in prepare_source.py, flusort, the fludb scripts, which subtypes are downloaded.

Where each setting lives

To change	Edit	Keep in sync	Tier
Coverage or QC status	--query in rule augur_filter, workflow/snakemake_rules/segments.smk	qc_min_coverage, qc_pass_status in config/report.yaml	2
Segment minimum length	min_lengths in Snakefile	MIN_LENGTH in scripts/prepare_source.py (reporting only)	2
Genome minimum length	min_genome_size in genomes.smk	–	2
Exclude a sequence	config/exclude.tsv (rules below)	–	2
Clock filter strictness	clock_filter_iqd in rule refine, segments.smk	–	2
Clock rates	clock_rate() in segments.smk; genome rates in rule genome_refine	–	2
Nextclade datasets	rule fetch_hana_datasets (HA/NA); nextclade/flu/{subtype}/{segment}/ for the rest	–	2
Control names	CONTROL_NAME in scripts/prepare_source.py	control_pattern in config/report.yaml (the two differ today)	3
Record checks, re-sequenced specimens	scripts/prepare_source.py	Step 1	3
Which flusort calls go into which build (incl. assumed subtypes)	config/subtype_assignment.yaml (read by scripts/fludb_download_builds.py and scripts/report.py)	Assumed subtypes	2
Typing, which subtypes get builds	scripts/flusort/, config/subtype_assignment.yaml, wildcards and rule all	Checklist: adding code	3

config/exclude.tsv rules

• One ID per line, exactly the sample_ID: the JH number for JHH, the full isolate name with passage suffix for vaccines (e.g. A/Darwin/9/2021-E3). One entry removes the ID from every build of every subtype.
• Put a comment line above new entries with the date, your initials and the reason: # 2026-10-07 EA: outside clock bounds in HA and NA builds. Never put a comment on the same line as an ID: the report would no longer recognise that ID.
• Never remove an entry without the maintainer.

Step by step

1. Agree the change and the reason with the maintainer, and note the tier.

2. Branch (Branches and pull requests): git switch main && git pull && git switch -c yourname/short-description.

3. Record a baseline before editing. Run the tutorial (Tutorial 1) and save the filter logs and report tables outside the repo:

   snakemake --profile profiles/default --cores 8 --configfile config/tutorial.yaml
   mkdir -p ../filter-before
   rsync -a --include='*/' --include='filter_log.tsv' --include='filtered.tsv' --exclude='*' results/ ../filter-before/results/
   cp -R reports/tables ../filter-before/report-tables
   snakemake --cores 8 snapshot_clean

4. Make the smallest change in one place, update its sync partner from the table above, then dry-run with and without the tutorial config: snakemake -n --cores 1 and snakemake -n --cores 1 --configfile config/tutorial.yaml.

5. Run the tutorial again and compare:

   snakemake --profile profiles/default --cores 8 --configfile config/tutorial.yaml
   for f in results/*/*/filter_log.tsv; do
     echo "$f"
     echo "  before: $(tail -n +2 ../filter-before/$f | cut -f2 | sort | uniq -c | tr -s ' ' | tr '\n' ';')"
     echo "  after:  $(tail -n +2 $f | cut -f2 | sort | uniq -c | tr -s ' ' | tr '\n' ';')"
   done
   diff ../filter-before/report-tables/qc_builds.tsv reports/tables/qc_builds.tsv

   Then open the trees (Step 5 review checklist): a looser filter can let outliers in. The tutorial is small, so for threshold changes the maintainer also compares a production run.

6. Update the docs in the same pull request: this page, plus Step 1 for prepare_source.py changes (Keep the docs in step).

7. Commit with a descriptive message: old value → new value, why, and the effect (Commit messages). Stage files by name. Never commit or paste the results/ or reports/ files: they contain JHH IDs and GISAID-derived data.

8. Open a pull request with the tier, the reason, who approved it, before/after counts (no specimen IDs), and the docs you updated. Agent-written changes also follow the Required review process. The maintainer reviews and merges Tier 2-3 changes.

9. After the merge, check the next production report's Quality control: per build section matches what you expected, and record the change with the upload (Step 5).

Troubleshooting / FAQ

• A whole segment build lost most sequences. Check Nextclade QC for that segment (results/{subtype}/{segment}/nextclade.tsv, or Nextclade QC status by segment in the report); a dataset change can shift scores.
• A specimen is missing from every segment build, with "no Nextclade result". It has no HA sequence (or HA failed alignment); see ④.
• Every clade is "unassigned". The Nextclade HA download failed or changed format; check .run/logs/fetch_hana_datasets/.
• Vaccine strains are missing from the genome trees. Expected for now; see the note under genome builds.
• An ID in exclude.tsv is still in the tree. It must match sample_ID exactly (vaccines include the passage suffix) and sit on its own line.
• Two different specimens share an ID. prepare_source.py keeps only one. Ask the sequencing team to rename one at the source; do not edit source/ by hand.
• Should a long-branch outlier be excluded? Possibly; it is a Tier 2 decision for the maintainer (Step 5).

Next step → Step 4: Reports & Outputs