Filter singleton-insertion sequences using phyz util alnqc

[matsen]

Filter singleton-insertion sequences using phyz util alnqc

Motivation

The flu scaling survey (matsengrp/phyz#663) showed that a handful of
sequences per dataset create disproportionate alignment bloat via
singleton insertions — columns where only 1 sequence has a residue and
all others are gapped. For example, in H3 (86K seqs), just 5 sequences
create 156 of 186 singleton columns (84%), inflating gap fraction from
~3% to 12.4%.

These singleton-insertion sequences are problematic for tree building:

• They inflate VCF size and slow UShER/matOptimize
• They create long terminal branches that distort parsimony scores
• They add noise to host-specific subtree extraction

The existing curate_and_extract_coding_seqs.py filters catch gross
quality problems (terminal gaps, high gap fraction, ambiguous bases,
CDS validation) but do not detect singleton insertions because those
sequences have normal ungapped lengths and pass all per-sequence checks.

Prerequisites

phyz upstream: matsengrp/phyz#725 added --aligned mode to
phyz util filter (merged in PR #730). Build phyz from main to
get this feature.

Approach

Use phyz util alnqc (matsengrp/phyz#668, merged in PR #669) to
identify and remove singleton-insertion outliers. The phyz toolchain
provides three composable utilities:

1. phyz util alnqc aligned.fa > report.tsv
   — per-sequence diagnostic (singleton_cols is column 3 of the TSV)
2. Filter the TSV with pandas to extract names exceeding the threshold
   — avoids reading the alignment a second time
3. phyz util filter aligned.fa --names remove.txt --invert --aligned
   — remove flagged sequences, keep the rest (preserving alignment)
4. phyz util trim-cols -
   — strip resulting all-gap columns (reads from stdin with -)

Note: phyz reads .gz (gzip) transparently but not .xz. The
Snakemake rule below decompresses xz before calling phyz. If the
pipeline switched from xz to gzip compression (.fasta.gz instead
of .fasta.xz), the decompress/cleanup steps could be eliminated
and phyz could read the files directly.

The full pipeline is:

# Decompress curated MSA
xz -dk results/HA/H3/curated_msa.fasta.xz

# Identify outliers (sequences with >3 singleton columns)
phyz util alnqc results/HA/H3/curated_msa.fasta \
    --remove-threshold 3 > remove.txt

# Remove outliers and strip empty columns
phyz util filter results/HA/H3/curated_msa.fasta \
    --names remove.txt --invert --aligned | \
  phyz util trim-cols - > cleaned.fasta

# Recompress
xz cleaned.fasta

No re-alignment is needed: singleton columns by definition contain
residues from only the removed sequence(s), so dropping those columns
preserves all pairwise relationships among remaining sequences.

Implementation plan

New Snakemake rule: filter_singleton_insertions

Insert between curate_and_extract_coding_seqs and randomize_alignment:

The rule uses a shell: block for the phyz commands and a helper
Python script (scripts/filter_alnqc_report.py) to extract the
remove list from the TSV report. This avoids awk and keeps the
filtering logic readable and consistent with the rest of the pipeline.

scripts/filter_alnqc_report.py:

"""Extract sequence names exceeding singleton column threshold from alnqc report."""
import pandas as pd
import argparse

parser = argparse.ArgumentParser()
parser.add_argument("--report", required=True, help="alnqc_report.tsv from phyz util alnqc")
parser.add_argument("--threshold", type=int, required=True, help="Remove seqs with singleton_cols > this")
parser.add_argument("--output", required=True, help="Output names file (one per line)")
args = parser.parse_args()

df = pd.read_csv(args.report, sep="\t")
flagged = df.loc[df["singleton_cols"] > args.threshold, "seq_name"]
flagged.to_csv(args.output, index=False, header=False)

rule filter_singleton_insertions:
    conda: "envs/python.yaml"
    input:
        curated_msa="results/{segment}/{subtype}/curated_msa.fasta.xz"
    output:
        filtered_msa="results/{segment}/{subtype}/filtered_msa.fasta.xz",
        remove_list="results/{segment}/{subtype}/singleton_remove.txt",
        alnqc_report="results/{segment}/{subtype}/alnqc_report.tsv"
    params:
        remove_threshold=config.get("singleton_remove_threshold", 3)
    log:
        "logs/{segment}/{subtype}/filter_singletons.log"
    shell:
        """
        # Decompress (phyz cannot read xz natively, only gzip)
        xz -dk {input.curated_msa}
        MSA="${{input.curated_msa%.xz}}"
        trap 'rm -f "$MSA"' EXIT

        # Diagnostic report (always produced)
        # Columns: seq_name, ungapped_len, singleton_cols, low_occ_cols,
        #          total_block_width, unique_ins_frac
        phyz util alnqc "$MSA" > {output.alnqc_report} 2>> {log}

        # Generate remove list from report via pandas
        python scripts/filter_alnqc_report.py \
            --report {output.alnqc_report} \
            --threshold {params.remove_threshold} \
            --output {output.remove_list}

        # Filter and trim (or pass through if nothing to remove)
        if [ -s {output.remove_list} ]; then
            phyz util filter "$MSA" \
                --names {output.remove_list} --invert --aligned | \
              phyz util trim-cols - | \
              xz > {output.filtered_msa}
        else
            cp {input.curated_msa} {output.filtered_msa}
        fi
        """

Update downstream rules

Three rules currently consume curated_msa.fasta.xz and should switch
to filtered_msa.fasta.xz:

1. randomize_alignment (Snakefile line 182) — alignment input for
   tree building
2. create_root_fasta (Snakefile line 419) — root sequence inference
3. create_host_samples_file (Snakefile line 463) — host group
   extraction

All three should use filtered_msa.fasta.xz so that removed sequences
are excluded consistently from tree building, rooting, and host analysis.

rule randomize_alignment:
    input:
        curated_msa="results/{segment}/{subtype}/filtered_msa.fasta.xz"
    ...

(Same pattern for the other two rules.)

Config additions

# Singleton insertion filtering (phyz util alnqc)
# Sequences with > this many singleton columns are removed.
# Set to a very large value (e.g., 999999) to disable filtering.
singleton_remove_threshold: 3

Environment: phyz binary

The rule needs phyz on PATH. phyz has no external dependencies —
just a Zig compiler.

Building phyz from source:

# 1. Install Zig 0.15+ (https://ziglang.org/download/)
#    On Linux, download and extract the tarball, then add to PATH:
wget https://ziglang.org/builds/zig-linux-x86_64-0.15.0.tar.xz
tar xf zig-linux-x86_64-0.15.0.tar.xz
export PATH="$PWD/zig-linux-x86_64-0.15.0:$PATH"

# 2. Clone and build phyz
git clone git@github.com:matsengrp/phyz.git
cd phyz
zig build -Doptimize=ReleaseFast

# 3. The binary is at zig-out/bin/phyz — copy or symlink to PATH
cp zig-out/bin/phyz ~/.local/bin/

No conda, pip, or system libraries required. The build takes ~1 minute.

Expected impact (from flu scaling survey data)

These are projections from phyz's greedy-removal analysis on alignments
produced by phyz aln, not re-measured values. Actual impact on
Nextclade-aligned flu-usher data may differ.

Dataset	Total seqs	Seqs removed	% removed	Gap% before	Gap% after (est.)
H3	86,227	~430	0.50%	12.4%	~2.9%
H1	89,896	~450	0.50%	15.5%	~4.8%
H7	3,198	~16	0.50%	15.5%	~10.0%
H5	21,040	~106	0.50%	8.3%	~2.6%

Removing <0.5% of sequences reduces gap fraction by 3-10 percentage
points and shrinks alignment width by 5-11%.

Validation

Alignment integrity check (must pass before using filtered data)

The claim that no re-alignment is needed rests on the fact that
removed columns have occupancy 1 (only the removed sequence has a
residue there). This must be verified empirically:

Add a validation script (scripts/validate_filtered_alignment.py)
that takes the original and filtered alignments and checks:

1. Row preservation: For every sequence in the filtered alignment,
   extract its row from the original alignment, drop all-gap columns
   (the same columns trim-cols removes), and confirm the result is
   identical character-by-character to the filtered row.
2. No column reordering: The relative order of retained columns is
   unchanged.

If either check fails, the filtering is not safe and we need to
re-align. Add this as a Snakemake rule that runs after filtering and
before randomize_alignment.

Visual inspection

Produce phyz alnview HTML renderings of a representative subset
(e.g., H3 and one smaller dataset like N8) before and after filtering.
Specifically, view a region around a known singleton block (e.g., the
57-column block owned by EPI_ISL_361993 in H3) to confirm:

• The singleton columns are gone in the filtered alignment
• Flanking columns are unchanged
• No new gaps or shifts appeared

Quantitative checks

• Run phyz util alnqc on each curated MSA and inspect the report TSV.
  The singleton_cols column (column 3) identifies per-sequence burden.
• Verify known H3 outliers appear (e.g., EPI_ISL_361993 with 57
  singleton columns in the phyz scaling survey — may differ in
  Nextclade alignment)
• Gap fraction of each filtered alignment should be below the
  corresponding curated alignment's gap fraction (strictly less)
• Filtered alignment width should be <= curated alignment width
• Confirm downstream tree building (UShER + matOptimize) completes
  without errors on filtered alignments
• Parsimony score on filtered alignment should not exceed parsimony
  score on unfiltered alignment (removing noisy columns should not
  hurt parsimony)

Non-goals

• Modifying the existing Python curation script
• Tree-aware filtering (future work using edge stats)

[matsen]

cc @Haddox let me know if this seems reasonable!