GPTomics/bio-sequence-statistics

Version Compatibility

Reference examples tested with: BioPython 1.83+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show <package> then help(module.function) to check signatures

If code throws ImportError, AttributeError, or TypeError, introspect the installed package and adapt the example to match the actual API rather than retrying.

Sequence Statistics

"Calculate N50 and other assembly statistics" -> Compute sequence count, length distribution, N50/L50, auN, GC content, and nucleotide composition for FASTA datasets.

• Python: SeqIO.parse(), gc_fraction() (BioPython)

Calculate comprehensive statistics for sequence datasets using Biopython.

The Governing Principle

N50 measures CONTIGUITY, not correctness. A misassembled scaffold that wrongly joins distant regions can post a large N50 while being biologically wrong; N50 says nothing about base accuracy or join correctness. Treat contiguity metrics as one axis of assembly quality, alongside completeness (BUSCO) and correctness (read-backed validation, reference alignment).

Two further traps shape every reported number:

• N50 is a discontinuous, threshold-based statistic. Near the 50% crossing, contig lengths can differ by megabases, so a tiny change can jump N50 by megabases. Prefer auN (smooth, threshold-free) for robust comparison.
• N50 uses ASSEMBLY size as the denominator, so it is not comparable across assemblies of the same genome. Use NG50 (GENOME size denominator) to compare assemblies on a common baseline.

Required Imports

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction
import statistics

N50, L50, and Nx Statistics

Goal: Report the length at which half the assembled bases reside in equal-or-longer contigs (N50), and how many contigs that takes (L50).

Approach: N50 is the minimal length x such that contigs of length >= x together cover >= 50% of total assembly length. Sort lengths DESCENDING, take the cumulative sum, and return the length at which the cumulative sum first reaches or crosses 50%. L50 is the COUNT of contigs at that crossing. Three details silently break naive implementations: the sort must be descending, the crossing test must be >= (not strict >), and the denominator must be the assembly total (not a genome estimate).

Reference (BioPython 1.83+):

def n50_l50(lengths):
    '''Return (N50 length, L50 count) for a list of contig lengths.'''
    sorted_lengths = sorted(lengths, reverse=True)
    half = sum(sorted_lengths) / 2
    cumsum = 0
    for count, length in enumerate(sorted_lengths, start=1):
        cumsum += length
        if cumsum >= half:
            return length, count
    return 0, 0

lengths = [len(r.seq) for r in SeqIO.parse('assembly.fasta', 'fasta')]
n50, l50 = n50_l50(lengths)
print(f'N50: {n50:,} bp  L50: {l50} contigs')

Any Nx Value (N75, N90)

def calculate_nx(lengths, x):
    '''Nx where x is a percentage (50 for N50, 90 for N90).'''
    sorted_lengths = sorted(lengths, reverse=True)
    threshold = sum(sorted_lengths) * (x / 100)
    cumsum = 0
    for length in sorted_lengths:
        cumsum += length
        if cumsum >= threshold:
            return length
    return 0

lengths = [len(r.seq) for r in SeqIO.parse('assembly.fasta', 'fasta')]
print(f'N50: {calculate_nx(lengths, 50):,} bp')
print(f'N90: {calculate_nx(lengths, 90):,} bp')

auN: Robust Contiguity (Preferred for Comparison)

Goal: Replace the discontinuous N50 with a smooth, threshold-free contiguity score that responds to every join.

Approach: auN (Heng Li, 2020) is the area under the Nx curve, equivalently a length-weighted average length: each contig contributes its own length weighted by the fraction of the assembly it represents. Connecting any two contigs always raises auN, even when N50 stays unchanged (joining two contigs both above, or both below, the N50 contig leaves N50 fixed). No single straddling contig arbitrarily sets the score.

Formula: auN = sum_i(L_i^2) / sum_j(L_j)

Reference (BioPython 1.83+):

def calculate_aun(lengths):
    '''auN = sum(L_i^2) / sum(L_j); a length-weighted mean length.'''
    total = sum(lengths)
    return sum(length * length for length in lengths) / total if total else 0

lengths = [len(r.seq) for r in SeqIO.parse('assembly.fasta', 'fasta')]
print(f'auN: {calculate_aun(lengths):,.0f} bp')

NG50, NGx, NA50, NGA50: Cross-Assembly and Misassembly-Aware

NG50/NGx use the GENOME size as the denominator instead of the assembly size, so two assemblies of the same genome share one baseline and become directly comparable. They require a known or estimated genome size (QUAST --est-ref-size or a reference). NA50/NGA50 are computed on alignment blocks broken at misassembly breakpoints; NGA50 markedly below NG50 signals misassemblies.

Reference (BioPython 1.83+):

def calculate_ngx(lengths, genome_size, x=50):
    '''NGx uses genome_size (not assembly size) as the denominator.'''
    sorted_lengths = sorted(lengths, reverse=True)
    threshold = genome_size * (x / 100)
    cumsum = 0
    for length in sorted_lengths:
        cumsum += length
        if cumsum >= threshold:
            return length
    return 0  # assembly never covers x% of the genome

lengths = [len(r.seq) for r in SeqIO.parse('assembly.fasta', 'fasta')]
print(f'NG50: {calculate_ngx(lengths, genome_size=3_100_000_000):,} bp')

Length Distribution

Median contig length is near-useless for assemblies: it is dominated by the many tiny contigs and sits among fragments, ignoring where the sequence mass lives. N50 and auN are mass-weighted precisely to answer "in contigs of what size does the bulk of the genome reside?" Report median for read-length QC, not for assembly contiguity.

lengths = [len(r.seq) for r in SeqIO.parse('sequences.fasta', 'fasta')]
print(f'Count: {len(lengths)}  Total: {sum(lengths):,} bp')
print(f'Min: {min(lengths):,}  Max: {max(lengths):,}  Mean: {statistics.mean(lengths):,.1f} bp')

Length Histogram Data

from collections import Counter

lengths = [len(r.seq) for r in SeqIO.parse('sequences.fasta', 'fasta')]
bin_size = 100  # 100-bp length bins
histogram = Counter((l // bin_size) * bin_size for l in lengths)
for length_bin in sorted(histogram):
    print(f'{length_bin}-{length_bin + bin_size}: {histogram[length_bin]}')

GC Content: Choose the Ambiguity Mode Explicitly

gc_fraction(seq, ambiguous=...) returns a FRACTION 0-1 (the old Bio.SeqUtils.GC() returned a PERCENT 0-100 and was REMOVED in 1.82; swapping names without rescaling is a silent 100x error). The ambiguous= argument changes the answer, so set it on purpose:

| Mode | Numerator | Denominator | GCGCNNNN | |------|-----------|-------------|------------| | remove (default) | G, C, S | only unambiguous + S/W (N excluded) | 1.0 | | ignore | G, C, S | full length (N dilutes GC) | 0.5 | | weighted | G, C, S + each ambiguous code x its expected GC | full length | 0.75 |

A naive (G + C) / len silently equals the ignore mode, under-reporting GC whenever N is present. remove reports GC among called bases; ignore reports GC over the full sequence including gaps/Ns; weighted apportions each IUPAC code its expected GC.

from Bio.Seq import Seq
from Bio.SeqUtils import gc_fraction

seq = Seq('GCGCNNNN')
gc_fraction(seq, ambiguous='remove')    # 1.0  - N dropped from both
gc_fraction(seq, ambiguous='ignore')    # 0.5  - N counted in denominator
gc_fraction(seq, ambiguous='weighted')  # 0.75 - N contributes 0.5 each

Per-Sequence GC Distribution

gc_values = [gc_fraction(r.seq, ambiguous='remove') for r in SeqIO.parse('sequences.fasta', 'fasta')]
print(f'Mean GC: {statistics.mean(gc_values):.1%}')
print(f'Median GC: {statistics.median(gc_values):.1%}')
print(f'Range: {min(gc_values):.1%} - {max(gc_values):.1%}')

Comprehensive Summary Report

Goal: Generate a complete QC summary (counts, lengths, N50/L50, auN, GC) for any FASTA file in one pass.

Approach: Load all records once, compute length and GC arrays, derive N50/L50 from the cumulative sorted lengths and auN from the squared-length sum, and package into a dictionary.

Reference (BioPython 1.83+):

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction
import statistics

def sequence_summary(fasta_file):
    records = list(SeqIO.parse(fasta_file, 'fasta'))
    lengths = [len(r.seq) for r in records]
    gc_values = [gc_fraction(r.seq, ambiguous='remove') for r in records]

    sorted_lengths = sorted(lengths, reverse=True)
    total_bp = sum(lengths)
    half = total_bp / 2

    cumsum, n50, l50 = 0, 0, 0
    for count, length in enumerate(sorted_lengths, start=1):
        cumsum += length
        if cumsum >= half:
            n50, l50 = length, count
            break

    aun = sum(length * length for length in lengths) / total_bp if total_bp else 0

    return {
        'file': fasta_file, 'sequences': len(records), 'total_bp': total_bp,
        'min_length': min(lengths), 'max_length': max(lengths),
        'mean_length': statistics.mean(lengths), 'median_length': statistics.median(lengths),
        'n50': n50, 'l50': l50, 'aun': aun,
        'gc_mean': statistics.mean(gc_values),
        'gc_std': statistics.stdev(gc_values) if len(gc_values) > 1 else 0,
    }

stats = sequence_summary('assembly.fasta')
print(f'Sequences: {stats["sequences"]:,}  Total: {stats["total_bp"]:,} bp')
print(f'N50: {stats["n50"]:,} bp (L50: {stats["l50"]})  auN: {stats["aun"]:,.0f} bp')
print(f'GC: {stats["gc_mean"]:.1%} (+/- {stats["gc_std"]:.1%})')

Compare Multiple Assemblies

Goal: Build a side-by-side table of key metrics across assembly files.

Approach: Run sequence_summary on each file and format the results into an aligned table; auN is the most reliable single column for ranking contiguity.

Reference (BioPython 1.83+):

from pathlib import Path

files = sorted(Path('assemblies/').glob('*.fasta'))
print(f'{"File":<30} {"Seqs":>8} {"Total bp":>15} {"N50":>12} {"auN":>12}')
print('-' * 80)
for fasta_file in files:
    s = sequence_summary(str(fasta_file))
    print(f'{fasta_file.name:<30} {s["sequences"]:>8,} {s["total_bp"]:>15,} {s["n50"]:>12,} {s["aun"]:>12,.0f}')

Nucleotide Composition

from collections import Counter

def nucleotide_composition(fasta_file):
    counts = Counter()
    for record in SeqIO.parse(fasta_file, 'fasta'):
        counts.update(str(record.seq).upper())
    total = sum(counts.values())
    return {base: count / total for base, count in counts.items()}

comp = nucleotide_composition('sequences.fasta')
for base in ['A', 'T', 'G', 'C', 'N']:
    if base in comp:
        print(f'{base}: {comp[base]:.2%}')

Common Errors

| Symptom | Cause | Fix | |---------|-------|-----| | N50 looks far too small | Sorted ascending instead of descending | Sort lengths with reverse=True before the cumulative sum | | N50 off by one contig near the crossing | Strict > test misses the exact-50% case | Use >= 50% (minimal length that reaches/exceeds half) | | N50 not comparable between assemblies | Used assembly size as denominator | Use NG50 with the genome size for cross-assembly comparison | | GC values off by 100x | Treated gc_fraction (0-1) like old GC() (0-100) | Multiply by 100 only for display; never mix the two | | GC silently low when Ns present | Default remove vs naive (G+C)/len (= ignore) | Pass ambiguous= explicitly to match intent | | Huge N50 on a wrong assembly | N50 measures contiguity, not correctness | Pair with BUSCO completeness and read-backed/reference validation; prefer auN |

References

• Li H (2020). "auN: a new metric to measure assembly contiguity." Technical blog post, https://lh3.github.io/2020/04/08/a-new-metric-on-assembly-contiguity (auN = area under the Nx curve; smooth, threshold-free contiguity).

Related Skills

• read-sequences - Parse sequences for statistics calculation
• batch-processing - Calculate stats across multiple files
• fastq-quality - Quality score statistics for FASTQ files
• sequence-manipulation/sequence-properties - Per-sequence GC content and properties
• alignment-files/bam-statistics - samtools stats/flagstat for alignment statistics