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]}')