MSAfix

MSAfix checks and cleans nucleotide or amino-acid multiple-sequence alignments. It can process one alignment or perform post-alignment QC across a manifest of independent loci.

Build

MSAfix requires a C++17 compiler and zlib.

make

Run the integration suite with:

make test

Clean one alignment

MSAfix -i input.fasta -o cleaned.fasta [options]

The input and output paths must differ.

Option	Description
-nt	Treat the alignment as nucleotide data (default).
-aa	Treat the alignment as amino-acid data.
-recoverTechnicalOffsets	Conservatively repair technical gap-placement offsets in a coding nucleotide MSA (-recoverFrameshifts is a compatibility alias).
-codingFrame 1|2|3	First coding base in the alignment (default 1).
-geneticCode INT	NCBI transl_table used for recovery validation (default 1).
-recoveryBand INT	Maximum displacement of each residue in alignment columns (default 3).
-recoveryReport FILE	Write per-sequence before/after evidence; also enables recovery.
-threads INT	Number of worker threads (default 1).
-rmGapColsGreater FLOAT	Remove columns whose gap fraction is greater than this value (0–1; default 1).
-maskBorderGap	Mask terminal gaps with N (nucleotide) or X (amino acid).
-maskLowID	Mask long regions with unusually low identity; this is a generic segment filter, not a frameshift test.
-minGoodPosFrac FLOAT	Remove sequences with a smaller fraction of usable positions (0–1; default 0).
-gapsNoPen	Do not penalize gaps in identity calculations.
-avgID FILE	Write each sequence's average identity to the other sequences.
-noFullDist FILE	Write the full pairwise distance matrix.

Example:

MSAfix -i input.fna -o cleaned.fna -nt -maskBorderGap \
-rmGapColsGreater 0.9 -minGoodPosFrac 0.5

Recover technical alignment-frame offsets

Recovery is for a technical alignment error: a row contains the right observed symbols, but a gap was placed at the wrong position and a block is therefore offset from homologous columns. It does not claim that the organism has a biological frameshift and it never invents a missing nucleotide.

MSAfix -i coding.fna -o recovered.fna -nt -recoverTechnicalOffsets \
-codingFrame 1 -geneticCode 1 -recoveryBand 3 \
-recoveryReport recovery.tsv -threads 4

For each row, MSAfix extracts the non-gap symbols and uses bounded dynamic programming to re-thread those same symbols through the existing alignment columns against a leave-one-out nucleotide profile. A proposal is accepted only when it has enough evidence and improves both nucleotide agreement and codon-level amino-acid agreement, without adding internal stops. All proposals are computed against the original immutable MSA and are applied only after the worker threads finish.

The safeguards are:

• The alignment width and the ungapped symbol string of every row are unchanged.
• At least three sequences are required.
• A row must improve by at least four nucleotide matches and two amino-acid matches, with at least 24 comparable nucleotides and eight comparable codons.
• Nucleotide identity must improve by at least 0.08 and amino-acid identity by at least 0.10; internal stops may not increase.
• The default band permits the expected one- or two-column technical offset.

Supported genetic-code IDs are the current NCBI tables: 1–6, 9–16, and 21–33. Undefined IDs (7, 8, and 17–20) are rejected.

Recovery assumes that most informative sequences are correctly aligned. It cannot identify a displacement shared by the whole MSA, repair a row with no available gap columns without changing the MSA width, or reconstruct an unobserved nucleotide. For true biological frameshifts or pseudogenes, use a frameshift-aware coding-sequence aligner rather than this recovery pass.

Post-alignment locus QC

Locus-QC mode is a native implementation of post_alignment_locus_qc.pl. It reads one alignment path per line from a manifest, writes a tab-separated report for every locus, and writes accepted paths to a keep-list:

MSAfix -manifest alignments.txt -report locus_qc.tsv -keep retained.txt

Blank lines and lines beginning with # are ignored. Alignments may be plain FASTA or gzip-compressed; if a listed path is absent, MSAfix also tries the same path with .gz appended.

Before any workers start, MSAfix resolves those fallbacks and rejects alignments that refer to the same file through duplicate text, relative/absolute aliases, symbolic links, hard links, or fallback/direct .gz names. The manifest, locus report, keep-list, and optional sequence-outlier report must be pairwise distinct, and every output must differ from every resolved alignment. A validation failure occurs before any alignment or output file is changed.

The QC calculation measures the fraction of valid alignment cells, creates an unambiguous column consensus, and measures each sufficiently comparable sequence's divergence from that consensus. Loci are rejected by permissive absolute divergence thresholds and, when enough structurally valid loci are available, by a cross-locus modified-Z outlier test. The report retains all loci and records each rejection reason.

The report also records called_cells, gc_cells, gc_fraction, and effective_sites. These are calculated after column-overlap filtering: a column is retained when at least ceil(sequence_count * minOverlapMSA) sequences are not missing (-, N, ?, or . for nucleotide data; -, X, ?, or . for amino-acid data). To match alignmentGCMetric, called_cells counts only A, C, G, and T, gc_cells counts G and C, and effective_sites is called_cells / sequence_count.

The original Perl flags and defaults are supported:

Option	Description	Default
-sequenceType nt|aa	Alignment alphabet.	nt
-maskSequenceOutliers	Remove conservatively detected isolated sequence outliers from each alignment before final locus QC.	Off
-sequenceOutlierReport FILE	Atomic audit report; required with -maskSequenceOutliers.	—
-sequenceOutlierExemptPrefix STRING	Exclude full headers beginning with this prefix from detection and removal; repeat for multiple prefixes.	None
-minSequences INT	Minimum comparable sequences.	3
-minOccupancy FLOAT	Minimum valid alignment-cell fraction.	0.35
-minOverlapMSA FLOAT	Minimum called-sequence fraction for GC/effective-site columns.	0.35
-minComparableFraction FLOAT	Minimum consensus-site fraction required per sequence.	0.25
-minComparableSites INT	Absolute comparable sites required per sequence.	NT 30; AA 10
-maxMedianDivergence FLOAT	Maximum median consensus divergence.	NT 0.18; AA 0.35
-maxP90Divergence FLOAT	Maximum 90th-percentile consensus divergence.	NT 0.30; AA 0.50
-relativeModifiedZ FLOAT	Cross-locus robust outlier threshold.	8.0
-minRelativeMedianDivergence FLOAT	Median-divergence floor for the relative filter.	NT 0.06; AA 0.12
-minRelativeP90Divergence FLOAT	P90-divergence floor for the relative filter.	NT 0.12; AA 0.25
-minLociForRelative INT	Structurally valid loci required for cross-locus QC.	8
-threads INT	Alignment files evaluated concurrently; report order remains the manifest order.	1
-deactivate	Disable locus QC and pass every manifest entry through to the keep-list. Report metrics are zero in this mode.	Off

Mask isolated sequence outliers

Sequence masking is opt-in and applies independently within each manifest alignment. Enable it with an audit destination:

MSAfix -manifest alignments.txt -report locus_qc.tsv -keep retained.txt \
-maskSequenceOutliers -sequenceOutlierReport sequence_outliers.tsv \
-sequenceOutlierExemptPrefix "outgroup|"

The prefix is compared with the complete FASTA header text after >. Exempt records neither contribute to the masking consensus nor qualify for removal, and the prefix option may be repeated.

The detector builds an unambiguous consensus from non-exempt records and requires at least eight evaluable non-exempt sequences. Evaluation requires at least max(100, ceil(0.25 * consensus sites)) comparable nucleotide sites or max(30, ceil(0.25 * consensus sites)) amino-acid sites. A candidate must have an upper modified Z score greater than 5.0 and divergence of at least 0.10 for nucleotides or 0.20 for amino acids.

At most max(1, floor(0.10 * evaluable sequences)) candidates are treated as isolated errors. If more sequences qualify, all are retained as a possible subpopulation or generally problematic locus and are reported as RETAINED_CLUSTER. Otherwise, each candidate FASTA record is removed and reported as MASKED.

Affected plain or gzip-compressed alignments are replaced atomically in place. The ordinary locus metrics and PASS/REJECT decision are then recalculated from the retained records. The audit is written atomically in manifest and FASTA order and records the comparable sites, differences, divergence, median, MAD, modified Z score, evaluable-sequence count, and isolation limit for every masked or cluster-retained candidate.

The ordinary locus report appends sequence_mask_status, sequence_mask_consensus_sites, sequence_mask_required_comparable_sites, sequence_mask_evaluable_sequences, sequence_mask_exempt_sequences, sequence_mask_candidate_sequences, sequence_mask_removed_sequences, and sequence_mask_policy_version. Status is DISABLED, INSUFFICIENT_EVIDENCE, NO_CANDIDATES, MASKED, or RETAINED_CLUSTER; the current policy version is 1. These fields distinguish a disabled pass, inadequate sequence/site evidence, a sufficiently evaluated clean locus, isolated removals, and a retained divergent subgroup.

Without -maskSequenceOutliers, alignments are never rewritten, no sequence-outlier audit is created, and the appended status is DISABLED.

Both one- and two-dash forms of locus-QC options are accepted (for example, -manifest and --manifest). -noLocusQC is accepted as an alias for -deactivate.

To temporarily deactivate filtering while keeping pipeline output files intact:

MSAfix -manifest alignments.txt -report locus_qc.tsv -keep retained.txt \
-deactivate

Use MSAfix -h for command help and MSAfix -v for version information.