speconsense 0.7.6
High-quality clustering and consensus generation for Oxford Nanopore amplicon reads
License: BSD-3-Clause
Requires Python: >=3.8
Keywords: bioinformatics, nanopore, consensus, clustering, amplicon

# Speconsense

A tool for high-quality clustering and consensus generation from Oxford Nanopore amplicon reads.

## Overview

Speconsense is a specialized tool for generating high-quality consensus sequences from Oxford Nanopore Technology (ONT) amplicon data. It is specifically designed as an experimental alternative to NGSpeciesID in the fungal DNA barcoding pipeline.

The key features of Speconsense include:
- Robust clustering of amplicon reads using either Markov Clustering (MCL) graph-based or greedy algorithms
- **Automatic variant phasing**: Detects variant positions within clusters and splits reads into separate haplotypes
- Automatic merging of clusters with identical or similar consensus sequences
- High-quality consensus generation using SPOA
- Primer trimming for clean consensus sequences
- Read identity metrics for quality assessment
- IUPAC ambiguity codes for unphased heterozygous positions
- Optimized for fungal amplicon datasets but suitable for any amplicon sequencing application

## Installation

### Requirements

- Python 3.8 or higher
- External dependencies:
  - [SPOA (SIMD POA)](https://github.com/rvaser/spoa) - Required (install via conda)
  - [MCL](https://micans.org/mcl/) - Optional but recommended for graph-based clustering (install via conda)
  - [vsearch](https://github.com/torognes/vsearch) - Optional, required for scalability mode with large datasets (install via conda)

### Install from GitHub (Recommended)

The easiest way to install speconsense is directly from GitHub using pip. We recommend using a virtual environment to avoid dependency conflicts:

```bash
# Create and activate a virtual environment (recommended)
python -m venv speconsense-env
source speconsense-env/bin/activate  # On Windows: speconsense-env\Scripts\activate

# Install directly from GitHub
pip install git+https://github.com/joshuaowalker/speconsense.git

# External dependencies need to be installed separately
# SPOA: conda install bioconda::spoa
# MCL: conda install bioconda::mcl (optional but recommended)
```

After installation, the tools will be available as command-line programs:
- `speconsense` - Main clustering and consensus tool
- `speconsense-summarize` - Post-processing and summary tool
- `speconsense-synth` - Synthetic read generator for testing

To deactivate the virtual environment when done:
```bash
deactivate
```

### Installing External Dependencies

**MCL (Markov Clustering) - Recommended:**
```bash
# Via conda (easiest)
conda install bioconda::mcl

# Or from source (more complex)
# See https://micans.org/mcl/ for source installation
```

**SPOA (SIMD POA) - Required:**
```bash
# Via conda (easiest)
conda install bioconda::spoa

# Or install from GitHub releases or build from source
# See https://github.com/rvaser/spoa for source installation instructions
```

**vsearch - Optional (for scalability mode):**
```bash
# Via conda (easiest)
conda install bioconda::vsearch

# Or download from https://github.com/torognes/vsearch/releases
```

**Note:** If the mcl tool is not available, speconsense will automatically fall back to the greedy clustering algorithm.

## Usage

### Basic Usage

```bash
speconsense input.fastq
```

By default, this will:
1. Cluster reads using graph-based Markov Clustering (MCL) algorithm
2. Merge clusters with identical consensus sequences
3. Generate a consensus sequence for each cluster
4. Output FASTA files containing consensus sequences

### Common Options

```bash
# Use greedy clustering algorithm instead of Markov Clustering
speconsense input.fastq --algorithm greedy

# Set minimum cluster size
speconsense input.fastq --min-size 10

# Set minimum identity threshold for clustering (default: 0.9)
speconsense input.fastq --min-identity 0.85

# Control the maximum sample size for consensus generation (default: 100)
speconsense input.fastq --max-sample-size 200

# Specify output directory (default: clusters)
speconsense input.fastq --output-dir results/

# Disable automatic variant phasing
speconsense input.fastq --disable-position-phasing

# Using short form for output directory
speconsense input.fastq -O my_results/
```

### Using Profiles

Profiles save parameter configurations for different workflows:

```bash
# List available profiles
speconsense --list-profiles

# Use a profile (CLI arguments override profile values)
speconsense input.fastq -p herbarium
speconsense input.fastq -p herbarium --min-size 10
```

**Bundled profiles:**
- `compressed` — Compress variants into minimal IUPAC consensus sequences (aggressive merging with indels, 20% thresholds, full consensus, 20% selection size ratio)
- `herbarium` — High-recall for degraded DNA/type specimens
- `largedata` — Experimental settings for large input files
- `nostalgia` — Simulate older bioinformatics pipelines
- `strict` — High-precision for confident results

On first use, an `example.yaml` template is created in `~/.config/speconsense/profiles/` — copy and edit it to create custom profiles.

### Two-Phase Processing Philosophy

**speconsense** extracts distinct biological sequences from raw reads within a single specimen. It separates true biological variation (alleles, loci, contaminants) from sequencing noise by clustering reads and generating consensus sequences. Philosophy: *split rather than conflate* — better to produce separate clusters that can be merged downstream than to lose real biological distinctions.

**speconsense-summarize** curates and consolidates consensus sequences for downstream analysis. It handles residual ambiguity that consensus generation cannot resolve — primarily minor SNP differences that may represent the same biological entity — groups related variants, selects representatives, and can analyze patterns across a run. Philosophy: *conservative simplification* — merge only what is clearly equivalent, preserve distinctions when uncertain.

| Aspect | speconsense | speconsense-summarize |
|--------|-------------|----------------------|
| Input | Raw reads | Consensus sequences |
| Scope | Single specimen | Single or multiple specimens |
| Goal | Discovery/enumeration | Curation/simplification |
| Error model | Read-level noise | Residual ambiguity |
| Bias | Prefer over-splitting | Conservative merging |

### Post-processing with Speconsense-Summarize

After running speconsense, use the summarize tool to process and refine outputs:

```bash
# Generate summary with default settings
speconsense-summarize

# Custom minimum RiC (Reads in Consensus) threshold
speconsense-summarize --min-ric 5

# Process specific source directory and custom output
speconsense-summarize --source /path/to/speconsense/output --summary-dir MyResults
```

**Variant Detection and Haplotype Phasing:**

Speconsense can detect and isolate sequence variants within specimens (aka "haplotype phasing"). The graph-based clustering algorithm excels at discriminating between variants, and `speconsense-summarize` provides sophisticated tools for managing multiple variants per specimen, including:
- MSA-based variant merging with IUPAC ambiguity codes and size-weighted consensus
- Hierarchical variant grouping to separate contaminants from primary targets
- Size-based or diversity-based variant selection strategies

For detailed information on variant handling options, see the [Advanced Post-Processing](#advanced-post-processing) section below.

## Quick Start: Complete Workflow

Speconsense is designed to replace the NGSpeciesID step in the [ONT DNA Barcoding Fungal Amplicons protocol](https://www.protocols.io/view/primary-data-analysis-basecalling-demultiplexing-a-dm6gpbm88lzp/v4). Here's a complete workflow:

**1. Demultiplex with Specimux**

Follow the protocol through the "Demultiplex the Reads with Specimux" step.

**2. Run Speconsense**

Replace the NGSpeciesID step with speconsense:

```bash
# Instead of NGSpeciesID:
# ls *.fastq | parallel NGSpeciesID --ont --consensus --t 1 --abundance_ratio 0.2 ...

# Use Speconsense:
ls *.fastq | parallel speconsense {}
```

**3. Post-process with Speconsense-Summarize**

Process the outputs to prepare for downstream analysis:

```bash
speconsense-summarize
```

This will create a `__Summary__/` directory with:
- `summary.fasta` - All final consensus sequences (merged variants only)
- Individual FASTA files per specimen
- `summary.txt` - Statistics and metrics
- `quality_report.txt` - Prioritized list of sequences with potential quality concerns
- `FASTQ Files/` - Reads contributing to each final consensus
- `variants/` - Pre-merge variant FASTA and FASTQ files (for merged sequences)

## Output Files

### Speconsense Core Output

For each specimen, Speconsense generates:

1. **Main consensus FASTA**: `{sample_name}-all.fasta`
   - Contains all consensus sequences for the specimen (one per cluster)
   - Uses original cluster numbering: `{sample_name}-c{cluster_num}-RiC{size}`

2. **Debug directory** (`cluster_debug/`):
   - `{sample_name}-c{cluster_num}-RiC{size}-reads.fastq`: Original reads in each cluster
   - `{sample_name}-c{cluster_num}-RiC{size}-sampled.fastq`: Sampled reads used for consensus generation
   - `{sample_name}-c{cluster_num}-RiC{size}-untrimmed.fasta`: Untrimmed consensus sequences

### Speconsense-Summarize Output

When using `speconsense-summarize` for post-processing, creates `__Summary__/` directory with:

#### **Main Output Files** (summarization numbering: `-1.v1`, `-1.v2`, `-2.v1`):
- **Individual FASTA files**: `{sample_name}-{group}.v{variant}-RiC{reads}.fasta` (all final consensus variants)
- **Combined file**: `summary.fasta` - all final consensus sequences (merged variants only, excludes .raw files)
- **Statistics**: `summary.txt` - sequence counts and metrics
- **Quality report**: `quality_report.txt` - highlights sequences with potential quality concerns
- **Log file**: `summarize_log.txt` - complete processing log

#### **FASTQ Files/** (aggregated reads for final consensus):
- `{sample_name}-{group}.v{variant}-RiC{reads}.fastq` - all reads contributing to final variant consensus

#### **variants/** (pre-merge variant files for traceability):
- `{sample_name}-{group}.v{variant}.raw{N}-RiC{reads}.fasta` - FASTA for pre-merge variant N (when variants were merged)
- **FASTQ Files/**:
  - `{sample_name}-{group}.v{variant}.raw{N}-RiC{reads}.fastq` - reads for pre-merge variant N

### Naming Convention Summary

**Two distinct namespaces maintain traceability:**

| **Namespace** | **Used In** | **Format** | **Purpose** |
|---------------|-------------|------------|-------------|
| **Original** | Source `cluster_debug/` | `-c1`, `-c2`, `-c3` | Preserves speconsense clustering results |
| **Summarization** | `__Summary__/`, `FASTQ Files/`, `variants/` | `-1.v1`, `-1.v2`, `-2.v1`, `.raw1` | Post-processing groups and variants |
| **Full consensus** | `__Summary__/` | `-1.full` | IUPAC consensus from pre-merge components of surviving variants |

### Example Directory Structure
```
__Summary__/
├── sample-1.v1-RiC45.fasta                  # Primary variant (group 1, merged)
├── sample-1.v2-RiC23.fasta                  # Additional variant (not merged)
├── sample-1.full-RiC68.fasta                # Full IUPAC consensus for group 1 (surviving variants' components)
├── sample-2.v1-RiC30.fasta                  # Second organism group, primary variant
├── summary.fasta                            # All final consensus sequences (excludes .raw)
├── summary.txt                              # Statistics
├── quality_report.txt                       # Quality assessment report
├── summarize_log.txt                        # Processing log
├── FASTQ Files/                             # Reads for final consensus
│   ├── sample-1.v1-RiC45.fastq              # All reads for merged consensus
│   ├── sample-1.v2-RiC23.fastq              # Reads for additional variant
│   └── sample-2.v1-RiC30.fastq              # All reads for second group
└── variants/                                # Pre-merge variant files (for traceability)
    ├── sample-1.v1.raw1-RiC30.fasta         # Pre-merge variant 1 FASTA
    ├── sample-1.v1.raw2-RiC15.fasta         # Pre-merge variant 2 FASTA
    └── FASTQ Files/                         # Pre-merge variant FASTQ files
        ├── sample-1.v1.raw1-RiC30.fastq     # Reads for pre-merge variant 1
        └── sample-1.v1.raw2-RiC15.fastq     # Reads for pre-merge variant 2
```

### FASTA Header Metadata

Consensus sequence headers contain metadata fields separated by spaces:

**Core Fields (Always Present):**
- `size=N` - Total number of raw sequence reads in the cluster
- `ric=N` - **Reads in Consensus** - Number of reads actually used for consensus generation (may be less than size due to sampling limits)

**Optional Fields:**
- `rawric=N+N+...` - RiC values of .raw source variants (pre-merge, largest-first, only present in merged variants from speconsense-summarize)
- `rawlen=N+N+...` - Original sequence lengths before overlap merging (largest-first, only present in overlap-merged variants)
- `snp=N` - Number of SNP positions from IUPAC merging (only present in merged variants from speconsense-summarize)
- `ambig=N` - Count of IUPAC ambiguity codes in consensus (Y, R, W, S, K, M, etc.)
- `length=N` - Sequence length in bases (available via --fasta-fields option)
- `primers=list` - Comma-separated list of detected primer names (e.g., `primers=ITS1F,ITS4`)
- `rid=X.X` - Mean read identity percentage (0-100%, available via --fasta-fields qc preset)
- `rid_min=X.X` - Minimum read identity percentage (worst-case read, available via --fasta-fields qc preset)
- `group=N` - Variant group number (available via --fasta-fields option)
- `variant=vN` - Variant identifier within group (available via --fasta-fields option, only for variants)

**Example Headers:**
```
# Simple consensus from speconsense (debug files):
>sample-c1 size=50 ric=45 rid=97.2 rid_min=94.1 primers=ITS1F,ITS4

# Merged variant from speconsense-summarize (default fields):
>sample-1.v1 size=250 ric=250 rawric=100+89+61 snp=2 ambig=2 primers=ITS1F,ITS4

# Overlap-merged variant (different-length sequences merged):
>sample-1.v1 size=471 ric=471 rawric=248+223 rawlen=630+361 primers=ITS1F,ITS4

# With QC preset (--fasta-fields qc):
>sample-1.v1 size=250 ric=250 length=589 rid=98.5 ambig=1
```

**Notes:**
- Use `--fasta-fields` option to customize which fields appear in output headers (see [Customizing FASTA Header Fields](#customizing-fasta-header-fields))
- Read identity metrics (`rid`, `rid_min`) reflect homopolymer-normalized sequence identity
- Field name changes from earlier versions: `merged_ric` → `rawric`, `p50diff`/`p95diff` → `rid`/`rid_min`
- Variant merging only occurs between sequences with identical primer sets
- `snp` counts positions where IUPAC codes were introduced during merging; `ambig` counts total ambiguity codes in the final sequence

## Scalability Features

For large datasets (thousands of sequences), speconsense uses external tools to accelerate pairwise comparison operations. **Scalability mode is enabled by default** for datasets with 1000+ sequences.

### Configuration

```bash
# Default: scalability enabled for datasets with 1000+ sequences
speconsense input.fastq

# Custom threshold: enable scalability only for datasets with 500+ sequences
speconsense input.fastq --scale-threshold 500

# Disable scalability entirely (use brute-force for all sizes)
speconsense input.fastq --scale-threshold 0

# Same options work for speconsense-summarize
speconsense-summarize --scale-threshold 500 --source clusters/
```

The threshold parameter allows you to run a single command across multiple specimens of varying sizes. Smaller specimens will use brute-force (which has less overhead for small inputs), while larger specimens will use vsearch acceleration.

### Requirements

Scalability mode requires **vsearch** to be installed:

```bash
conda install bioconda::vsearch
```

If vsearch is not available, speconsense will automatically fall back to brute-force computation.

### How It Works

Scalability mode uses a two-stage approach to reduce the O(n^2) pairwise comparison cost:

1. **Fast candidate finding**: vsearch quickly identifies approximate sequence matches using heuristic search
2. **Exact scoring**: Only candidate pairs are scored using the full alignment algorithm

This reduces the computational complexity from O(n^2) to approximately O(n x k), where k is the number of candidates per sequence.

### Performance

Based on benchmarking, the break-even point is approximately 1000 sequences:
- Below 1000 sequences: brute-force is faster (less overhead)
- Above 1000 sequences: vsearch acceleration provides significant speedup
- At 10,000 sequences: ~10x faster than brute-force

### Implementation Notes

- The scalability layer is designed to be backend-agnostic; future versions may support alternative tools (BLAST, minimap2, etc.)
- For summarize's HAC clustering, scalability is most effective when there are many consensus sequences to cluster
- Scalability mode produces identical or near-identical results to brute-force mode

### Concurrency Control

The `--threads` option controls internal parallelism (vsearch, SPOA consensus generation). Use `0` for auto-detect.

**Default behavior differs between tools:**
- `speconsense`: defaults to `--threads 1` (safe for GNU parallel workflows)
- `speconsense-summarize`: defaults to `--threads 0` (auto-detect, since it runs once across all specimens)

```bash
# Many speconsense jobs via parallel - single-threaded by default (safe)
ls *.fastq | parallel speconsense {}

# Single large speconsense job - enable internal parallelism
speconsense large_dataset.fastq --threads 0

# speconsense-summarize auto-detects threads by default
speconsense-summarize --source clusters/
```

## Early Filtering

Speconsense can optionally apply size filtering early in the pipeline (after pre-phasing merge) to skip expensive variant phasing on small clusters that will ultimately be filtered out. This can significantly improve performance for large datasets.

### Options

- `--enable-early-filter`: Apply min-size and min-cluster-ratio filtering after pre-phasing merge. Small clusters skip variant phasing, improving performance.
- `--collect-discards`: Write all discarded reads (outliers + filtered clusters) to `cluster_debug/{sample}-discards.fastq` for inspection.

### Performance Optimization

For large datasets where performance is critical:

```bash
speconsense input.fastq --enable-early-filter
```

This skips variant phasing for clusters that will be filtered out anyway, while still using scalability optimizations for O(n²) comparisons.

## Algorithm Details

### Clustering Methods

Speconsense offers two clustering approaches with different characteristics:

#### **Graph-based clustering with Markov Clustering (MCL)** - Default and recommended
- Constructs a similarity graph between reads and applies the Markov Clustering algorithm to identify clusters
- **Tends to produce more clusters** by discriminating between sequence variants within a specimen
- Suitable when you want to identify multiple variants per specimen and are willing to interpret complex results
- Excellent for detecting subtle sequence differences and biological variation
- Relatively fast with high-quality results for most datasets

**Why MCL?** Markov Clustering was chosen for its ability to detect natural communities in densely connected graphs without requiring a priori specification of cluster numbers. The algorithm simulates flow in the graph, causing flow to spread within natural clusters and evaporate between different clusters. By varying a single parameter (inflation), MCL can detect clusters at different scales of granularity. MCL has proven highly effective in bioinformatics applications, particularly for protein family detection in large-scale sequence databases, making it well-suited for identifying sequence variants in amplicon data.

#### **Greedy clustering** - Fast and simple alternative
- Uses greedy star clustering that iteratively selects the sequence with the most similar neighbors as cluster centers
- **Tends to produce fewer, simpler clusters** by focusing on well-separated sequences
- Excellent at discriminating between distinct targets (e.g., primary sequence vs. contaminant)
- Usually does not separate variants within the same biological sequence
- Suitable when you want a fast algorithm with minimal complexity in interpreting results
- Ideal for applications where one consensus per specimen is sufficient

#### **Choosing the Right Algorithm**

**Use `--algorithm greedy` when:**
- You want one clear consensus sequence per specimen
- Speed is prioritized over detailed variant detection
- You prefer simpler output interpretation
- Your dataset has well-separated sequences (targets vs. contaminants)

**Use `--algorithm graph` (default) when:**
- You want to detect sequence variants within specimens
- You're willing to evaluate multiple consensus sequences per specimen
- You need fine-grained discrimination between similar sequences
- You want the most comprehensive analysis of sequence diversity
- **Note**: Use `speconsense-summarize` after clustering to manage multiple variants per specimen. Key options include `--merge-position-count` for merging variants differing by few SNPs/indels, `--group-identity` for grouping similar variants, and `--select-max-variants`/`--select-strategy` for controlling which variants to output

### Cluster Merging

After initial clustering, Speconsense automatically merges clusters with identical or homopolymer-equivalent consensus sequences:

**Default behavior (homopolymer-aware merging):**
- Clusters with identical consensus sequences are automatically merged
- Clusters with homopolymer-equivalent sequences are also merged (e.g., "AAA" vs "AAAAA" treated as identical)
- Homopolymer length differences are ignored as they typically represent ONT sequencing artifacts
- This helps eliminate redundant clusters that represent the same biological sequence but differ only in homopolymer lengths

**Strict identity merging (`--disable-homopolymer-equivalence`):**
- Only clusters with perfectly identical consensus sequences are merged
- Use this flag when homopolymer differences may have biological significance
- Results in more clusters but ensures no information loss from homopolymer variation

This automatic merging step helps consolidate redundant clusters that were separated during initial clustering. The adjusted identity scoring with homopolymer normalization provides more accurate assessment of sequence similarity for merging decisions, especially for nanopore data where homopolymer length calling can be inconsistent.

**Variant merging with IUPAC codes:** For more aggressive merging based on SNP and indel thresholds (which creates IUPAC consensus sequences with ambiguity codes), use the `--merge-position-count` and `--merge-indel-length` options in `speconsense-summarize` during post-processing. See the [Advanced Post-Processing](#advanced-post-processing) section.

### Cluster Size Filtering

Speconsense provides two complementary filters to control which clusters are output:

**Absolute size filtering (`--min-size`, default: 5):**
- Filters clusters by absolute number of sequences
- Applied **after merging** identical/homopolymer-equivalent clusters
- Set to 0 to disable and output all clusters regardless of size

**Relative size filtering (`--min-cluster-ratio`, default: 0.01):**
- Filters clusters based on size relative to the largest cluster
- Also applied **after merging** identical/homopolymer-equivalent clusters (post-merge sizes)
- **Based on original cluster sizes**, not sampled sizes from `--max-sample-size`
- Set to 0 to disable and keep all clusters that pass `--min-size`

**Processing order:**
1. Initial clustering produces raw clusters
2. Merge identical/homopolymer-equivalent clusters
3. Filter by `--min-size` (absolute threshold, using post-merge sizes)
4. Filter by `--min-cluster-ratio` (relative threshold, using post-merge sizes)
5. Sample sequences for consensus generation if cluster > `--max-sample-size`

This order ensures that small clusters with identical/homopolymer-equivalent consensus sequences can merge before size filtering is applied. This allows, for example, two clusters of size 3 with identical consensus to merge into a size-6 cluster that passes the `--min-size=5` threshold, rather than being discarded prematurely.

**Deferred filtering strategy:**
For maximum flexibility in detecting rare variants and contaminants, disable filtering in speconsense (`--min-size 0 --min-cluster-ratio 0`) and apply final quality thresholds using `--min-ric` in speconsense-summarize. This allows you to run expensive clustering once and experiment with different quality thresholds during post-processing. However, be aware that permissive filtering may allow more bioinformatic contamination through the pipeline. When using this approach, consider stricter filtering during upstream demultiplexing or perform careful manual review of low-abundance clusters.

### Consensus Generation

Consensus sequences are generated using SPOA (SIMD Partial Order Alignment) which efficiently handles the error profile of nanopore reads. For larger clusters, a random subset of reads (controlled by `--max-sample-size`) is used to generate the consensus.

### Read Identity and Variance Metrics

To evaluate the reliability of consensus sequences, Speconsense calculates read identity metrics by:
1. Aligning all reads in a cluster to the consensus sequence using SPOA multiple sequence alignment
2. Computing per-read identity scores with homopolymer normalization
3. Reporting mean read identity (`rid`) and minimum read identity (`rid_min`)

**Key metrics:**
- `rid` - Mean read identity: The average identity of all reads to the consensus (0-100%). Higher values indicate more homogeneous clusters.
- `rid_min` - Minimum read identity: The worst-case read identity. Low values may indicate outliers or mixed clusters.

**Homopolymer-normalized identity**: Speconsense uses the adjusted-identity algorithm with homopolymer normalization for more accurate sequence comparisons. This means that differences in homopolymer run lengths (e.g., AAA vs AAAAA) are treated as identical, which is particularly important for nanopore sequencing data where homopolymer length calling can be inconsistent. The identity metrics exclude homopolymer length differences, so low `rid` values indicate true substitutions or structural indels rather than sequencing artifacts.

**Automatic variant detection**: By default, Speconsense analyzes positional variation within clusters to detect true biological variants vs sequencing errors. Positions where variant alleles exceed the minimum frequency threshold are flagged as variant positions, and reads are automatically separated into distinct haplotypes (see [Automatic Variant Phasing](#automatic-variant-phasing)).

### Primer Trimming

When a primers file is provided via `--primers`, Speconsense will identify and trim primer sequences from the 5' and 3' ends of consensus sequences, producing clean amplicon sequences for downstream analysis.

**Automatic primer detection**: If `--primers` is not specified, Speconsense will automatically look for `primers.fasta` in the same directory as the input FASTQ file. If found, primer trimming will be enabled automatically.

### Automatic Variant Phasing

By default, Speconsense automatically detects and separates biological variants within clusters. This feature is particularly useful for heterozygous samples or mixed-species amplicons.

**How variant phasing works:**

1. **Variant detection**: After initial clustering, Speconsense analyzes positional variation using multiple sequence alignment. Positions where the minor allele frequency exceeds the threshold (default 10%) and meets minimum read count requirements are identified as variant positions.

2. **Position selection**: When multiple variant positions are detected, Speconsense selects the single best position for splitting (minimizing within-cluster error), then recursively regenerates MSA for each subcluster to discover additional variant positions. This hierarchical approach prevents over-fragmentation while allowing deep phasing when supported by the data.

3. **Haplotype separation**: Reads are grouped by their allele combinations at selected variant positions. Each unique combination becomes a separate sub-cluster (haplotype).

4. **Haplotype filtering**: Small haplotypes that don't meet minimum thresholds are reassigned to the nearest qualifying haplotype, preventing read loss.

5. **IUPAC ambiguity calling**: When only one haplotype qualifies (insufficient support for phasing), variant positions are encoded using IUPAC ambiguity codes (e.g., `Y` for C/T, `R` for A/G) rather than forcing a potentially incorrect consensus call. Ambiguity calling uses a lower count threshold (3 reads vs 5 for splitting), capturing variation that doesn't have enough support to form a viable separate cluster.

**Key parameters for variant phasing:**
- `--min-variant-frequency` - Minimum minor allele frequency to trigger cluster splitting (default: 0.10 = 10%)
- `--min-variant-count` - Minimum read count for minor allele to trigger splitting (default: 5)
- `--disable-position-phasing` - Disable variant phasing entirely

**Key parameters for IUPAC ambiguity calling:**
- `--min-ambiguity-frequency` - Minimum minor allele frequency for IUPAC codes (default: 0.10 = 10%)
- `--min-ambiguity-count` - Minimum read count for minor allele for IUPAC codes (default: 3)
- `--disable-ambiguity-calling` - Disable IUPAC codes for unphased variants

The frequency thresholds are unified at 10%, but the count thresholds differ: splitting requires 5 reads (matching `--min-size`) to ensure viable clusters, while ambiguity calling uses 3 reads since it doesn't create new clusters.

**Example:**
```bash
# Default behavior: variant phasing enabled
speconsense input.fastq

# More permissive variant detection (lower frequency threshold)
speconsense input.fastq --min-variant-frequency 0.05

# More aggressive IUPAC ambiguity calling
speconsense input.fastq --min-ambiguity-frequency 0.05 --min-ambiguity-count 2

# Disable variant phasing but keep ambiguity calling
speconsense input.fastq --disable-position-phasing

# Disable both phasing and ambiguity calling
speconsense input.fastq --disable-position-phasing --disable-ambiguity-calling
```

**When to use variant phasing (default):**
- Heterozygous specimens where you want to separate alleles
- Samples with potential mixed-species amplification
- Quality control to identify clusters with internal heterogeneity

**When to disable variant phasing:**
- Homozygous samples where variation indicates sequencing errors only
- When you prefer merged IUPAC consensus over separate haplotypes
- For faster processing when variant separation is not needed

## Advanced Post-Processing

The `speconsense-summarize` tool provides sophisticated options for managing multiple variants per specimen. This section covers advanced variant handling - for basic usage, see the [Usage](#usage) section above.

### Variant Grouping and Selection

When multiple variants exist per specimen, `speconsense-summarize` first groups similar variants together, then applies selection strategies within each group:

**Variant Grouping:**
```bash
speconsense-summarize --group-identity 0.9
```
- Uses **Hierarchical Agglomerative Clustering (HAC)** to group variants with sequence identity ≥ threshold
- Default threshold is 0.9 (90% identity)
- Variants within each group are considered similar enough to represent the same biological entity
- **Occurs before merging** to separate dissimilar sequences (e.g., contaminants from primary targets)
- Each group will contribute one primary variant plus additional variants based on selection strategy

**Variant Selection (within each group):**

When multiple variants exist per specimen, `speconsense-summarize` offers two strategies for selecting which variants to output from each group:

**Size-based selection (default):**
```bash
speconsense-summarize --select-strategy size --select-max-variants 2
```
- Selects variants by cluster size (largest first)
- Primary variant is always the largest cluster
- Additional variants are selected in order of decreasing size
- Best for identifying the most abundant sequence variants
- Suitable when read count reflects biological abundance
- **Default**: `--select-max-variants=-1` (outputs all variants, no limit)

**Diversity-based selection:**
```bash
speconsense-summarize --select-strategy diversity --select-max-variants 2
```
- Uses a **maximum distance algorithm** to select variants that are most genetically different from each other
- Primary variant is still the largest cluster
- Additional variants are selected to maximize sequence diversity in the output
- Iteratively selects the variant with the **maximum minimum distance** to all previously selected variants
- Best for capturing the full range of genetic variation in your sample
- Suitable when you want to detect distinct sequence types regardless of their abundance

**Algorithm Details for Diversity Selection (within each group):**
1. Primary variant = largest cluster within the group (by read count)
2. For each additional variant slot in the group:
   - Calculate the minimum sequence distance from each remaining candidate to all already-selected variants in this group
   - Select the candidate with the largest minimum distance (farthest from all selected in this group)
   - Repeat until select_max_variants reached for this group

**Overall Process:**
1. Group variants by sequence identity using HAC clustering
2. For each group independently:
   - Apply MSA-based merging to find largest compatible subsets
   - Apply variant selection strategy (size or diversity)
   - Output up to select_max_variants per group
3. Final output contains representatives from all groups, ensuring both biological diversity (between groups) and appropriate sampling within each biological entity (within groups)

**Selection Size Ratio Filtering:**
```bash
speconsense-summarize --select-min-size-ratio 0.2
```
- Filters out post-merge variants whose size is too small relative to the largest variant in their group
- Ratio calculated as `variant_size / largest_size` — must be ≥ threshold to keep
- Example: `--select-min-size-ratio 0.2` means a variant must have ≥20% the reads of the largest variant in its group
- Default is 0 (disabled) — all post-merge variants pass through to selection
- Applied after merging but before variant selection
- Useful for suppressing noise variants that survived merging but are too small to be meaningful
- Set to 0.2 in the `compressed` profile to match the 20% calling threshold theme

This two-stage process ensures that distinct biological sequences are preserved as separate groups, while providing control over variant complexity within each group.

### Customizing FASTA Header Fields

Control which metadata fields appear in FASTA headers using the `--fasta-fields` option:

**Presets:**
```bash
# Default preset (current behavior)
speconsense-summarize --fasta-fields default
# Output: >sample-1 size=638 ric=638 rawric=333+305 snp=1 ambig=1 primers=5'-ITS1F,3'-ITS4_RC

# Minimal headers - just the essentials
speconsense-summarize --fasta-fields minimal
# Output: >sample-1 size=638 ric=638

# QC preset - includes read identity and length
speconsense-summarize --fasta-fields qc
# Output: >sample-1 size=638 ric=638 length=589 rid=98.5 ambig=1

# Full metadata (all available fields)
speconsense-summarize --fasta-fields full
# Output: >sample-1 size=638 ric=638 length=589 rawric=333+305 snp=1 ambig=1 rid=98.5 primers=...

# ID only (no metadata)
speconsense-summarize --fasta-fields id-only
# Output: >sample-1
```

**Custom field selection:**
```bash
# Specify individual fields (comma-separated)
speconsense-summarize --fasta-fields size,ric,primers

# Combine presets and fields
speconsense-summarize --fasta-fields minimal,rid

# Combine multiple presets
speconsense-summarize --fasta-fields minimal,qc
```

**Available fields:**
- `size` - Total reads across merged variants
- `ric` - Reads in consensus
- `length` - Sequence length in bases
- `rawric` - RiC values of .raw source variants (only when merged)
- `rawlen` - Original sequence lengths before overlap merging (only when overlap-merged)
- `snp` - Number of IUPAC positions from merging (only when >0)
- `ambig` - Count of IUPAC ambiguity codes in consensus
- `rid` - Mean read identity percentage (when available)
- `rid_min` - Minimum read identity percentage (when available)
- `primers` - Detected primer names (when detected)
- `group` - Variant group number
- `variant` - Variant identifier within group (only for variants)

**Use cases:**
- **Downstream tool compatibility**: Use `minimal` or `id-only` for tools expecting simple headers
- **Quality control**: Use `qc` preset to include read identity metrics for assessing consensus quality
- **File size optimization**: Use `minimal` to reduce file size for large datasets
- **Custom workflows**: Combine presets and fields for workflow-specific needs

### Quality Assessment and Reporting

Speconsense-summarize automatically generates a `quality_report.txt` file to help prioritize manual review of sequences with potential quality concerns. This is particularly valuable for high-throughput workflows where human time for manual inspection is limited.

**Report Generation:**
- Created automatically in the summary output directory
- No configuration required - generated regardless of `--fasta-fields` settings
- Focuses exclusively on sequences that may need review (no "excellent" sequences listed)

**What the Report Includes:**

**1. Statistical Outliers (Low Read Identity):**
- Sequences with mean read identity (`rid`) below the statistical threshold
- Threshold is calculated as: mean - 2×standard deviation across all sequences
- These represent the ~2.5% of sequences with lowest internal consistency
- May indicate mixed clusters, contamination, or problematic consensus

**2. Sequences with IUPAC Ambiguity Codes:**
- Sequences containing IUPAC ambiguity codes (Y, R, W, S, K, M, etc.)
- May result from merging or from unphased heterozygous positions
- Count shown in `ambig` field

**3. Overlap Merge Analysis:**
- Lists specimens where partial overlap merging occurred (different-length sequences with shared overlap region)
- Shows merge iterations, overlap percentages, and prefix/suffix extensions
- Flags edge cases (overlap near threshold, large length ratios >3:1)
- Only includes partial overlaps (prefix or suffix extension >0); excludes full containment merges

**Understanding Read Identity Metrics:**

Read identity assessment works by:
1. Aligning all reads in a cluster to the consensus using SPOA multiple sequence alignment
2. Computing per-read identity with homopolymer normalization
3. Reporting mean identity (`rid`) and minimum identity (`rid_min`)

**Key points about identity calculation:**
- Uses homopolymer normalization - differences like AAA vs AAAAA don't reduce identity
- Therefore, low identity indicates **true substitutions or structural indels**
- `rid` typically ranges from 95-99% for good quality clusters
- `rid_min` shows the worst-case read - very low values may indicate outlier reads

**Example Report Entry:**
```
Type    Sequence                                              RiC    rid      Notes
------------------------------------------------------------------------------------------------
STAT    specimen-1                                            450    93.2     Statistical outlier (rid be