How to Filter Sequences by Length and GC Content in Python

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-11-11

By the end of this walkthrough you will be able to take any FASTA file, compute the length and GC content of every sequence in it, keep only the records that fall inside cutoffs you choose, write those records to a new file, and generate a report showing exactly why each record passed or failed. You will also have a reusable command-line template that turns the same logic into a tool for any project.

To follow along you need three things: a working Python 3 installation, the Biopython package (version 1.79 or later, which is when the gc_fraction function was introduced), and a sequence file in FASTA format to practice on. Everything here runs on an ordinary laptop with no other dependencies.

Filtering by length and GC content is a natural first script for anyone learning Python for bioinformatics, and it doubles as a compact Biopython tutorial. Python has a long track record in sequence analysis, with dedicated modules for segmentation and hidden Markov model modeling dating back to the mid-2000s [1]. The two measurements covered here are the standard first look at any sequence dataset, and published comparative work relies on them as well. One Biopython-based study comparing the SARS-CoV-2 genome with other virus genomes used sequence length and GC content among its primary descriptors [2].

Why Length and GC Content Come First

Length is the cheapest sanity check available. Assemblies arrive with adapter remnants, primer dimers, contamination, and chimeric or truncated contigs (a contig is one contiguous assembled sequence). A minimum length cutoff removes obvious fragments before you spend compute on alignment or annotation. Length also tests completeness against a reference. The published mitochondrial genome (the small, maternally inherited chromosome found in mitochondria) of the Indian rock python Python molurus molurus is a single 17,258 bp molecule carrying 37 genes, so an assembled snake mitogenome landing far below that size is incomplete or contaminated [3].

GC content, the percentage of bases that are guanine or cytosine, is a compositional fingerprint. Sharp shifts in GC content along a genome often flag contamination, and extreme local GC values cause practical failures downstream. Codon optimization (rewriting a gene to match a host organism's codon usage) for heterologous protein expression can leave local GC highs and lows that interfere with primer design and with synthesizing a whole gene from PCR fragments, which is why the Python program Visual codon exists specifically to inspect and adjust GC content after optimization [4]. Automated pipelines treat GC content as a core measurement too. The CGAS chloroplast genome pipeline computes GC content alongside gene content and structural features when batch processing dozens of chloroplast genomes [5].

Both numbers are deterministic, cheap to compute, and easy to audit, which makes them ideal as automated gates in a workflow.

The Two Numbers, Precisely Defined

Length is the number of characters in the sequence. In Python, len(sequence) returns it. The count includes every character: ambiguity codes, gap characters, and lowercase letters all count unless you remove them first.

GC content is expressed as a percentage:

GC percent = 100 × (number of G + number of C) / (A + T + G + C)

Worked by hand: the sequence ATGGCGTTC has 9 bases, 3 G and 2 C, so GC content = 100 × 5 / 9 = 55.56 percent.

The denominator is where implementations disagree. Some tools divide by every character in the string, including N and gaps. Others divide only by unambiguous A, T, G, and C. Both are defensible. What matters is that you know which rule your tool applies and that you apply the same rule everywhere in a project. Biopython exposes the choice through an argument covered below.

A related pair of numbers compares complementary bases on one strand. AT skew is (A - T) / (A + T), and the corresponding GC skew compares G with C, with sign conventions that vary between papers. These matter when composition is strand-asymmetric, which is common in mitochondrial genomes. The P. molurus molurus mitogenome reports positive AT and CG skews, meaning more A than T and more C than G on the positive strand [3].

Installing Python and Biopython

Check Python and install Biopython:

python --version
pip install biopython

Then confirm the two imports this guide depends on:

# check_biopython.py
import Bio
from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

print("Biopython", Bio.__version__, "is ready")

Running python check_biopython.py prints the installed version number, which varies by release. Any version from 1.79 onward includes gc_fraction. On older installations the import fails with:

ImportError: cannot import name 'gc_fraction' from 'Bio.SeqUtils'

Fix that with pip install -U biopython. If you cannot upgrade, the manual formula later in this guide produces the same number.

One version note that prevents silent errors: the legacy function Bio.SeqUtils.GC returned a fraction between 0 and 1, and many older scripts and forum answers still use it. The current gc_fraction returns a percentage between 0 and 100. Mixing the two conventions is the most common way people end up with GC thresholds that filter nothing at all.

Worked Example: Filter Sequences by Length and GC Content in Python

Create a small practice file called demo.fasta:

>contig_001
ATGGCGTTCGTA
>contig_002
GGCCGGCCNNNN
>contig_003
ATATATATAT
>contig_004
GCGCGCGCGCGCGCGCGC

A FASTA file stores each sequence as a header line starting with > followed by the sequence itself, wrapped over as many lines as needed. The four demo records are deliberately tiny so you can verify every number by hand. The cutoffs used below (length 10 to 20,000 bases, GC 40 to 100 percent) are illustrative for the demo, not recommendations for real data.

Step 1: Inspect the file before scripting

head -2 demo.fasta
grep -c ">" demo.fasta

Expected output:

>contig_001
ATGGCGTTCGTA
4

head shows the first header and the first sequence line. grep -c ">" counts header lines, which equals the record count in a well-formed FASTA file, where > appears only in headers. Never count lines with wc, because FASTA sequences wrap at arbitrary widths.

Step 2: Read the records with SeqIO

SeqIO is Biopython's sequence input and output module. SeqIO.parse returns an iterator (an object that yields one record at a time), where each record carries an identifier and a sequence:

from Bio import SeqIO

records = list(SeqIO.parse("demo.fasta", "fasta"))
print("records:", len(records))

for record in records:
    print(record.id, len(record.seq))

Expected output:

records: 4
contig_001 12
contig_002 12
contig_003 10
contig_004 18

record.id is the first word of the header line. len(record.seq) is the length in bases.

Step 3: Calculate GC content per record

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

for record in SeqIO.parse("demo.fasta", "fasta"):
    print(record.id, len(record.seq), f"{gc_fraction(record.seq):.2f}")

Expected output:

contig_001 12 50.00
contig_002 12 66.67
contig_003 10 0.00
contig_004 18 100.00

Check contig_001 by hand: 12 bases, 4 G and 2 C, so 6/12 = 50.00 percent. gc_fraction returns a percentage, and the format spec .2f rounds to two decimals.

Step 4: Apply the cutoffs

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

MIN_LEN, MAX_LEN = 10, 20000
MIN_GC, MAX_GC = 40.0, 100.0

kept = []
for record in SeqIO.parse("demo.fasta", "fasta"):
    length = len(record.seq)
    gc = gc_fraction(record.seq)
    if MIN_LEN <= length <= MAX_LEN and MIN_GC <= gc <= MAX_GC:
        kept.append(record)

print("kept:", len(kept))

Expected output:

kept: 3

contig_003 passes the length test (exactly 10 bases) but fails GC with 0.00 percent. The chained comparison MIN_LEN <= length <= MAX_LEN is standard Python and reads the way you would say it.

Step 5: Write the surviving records

count = SeqIO.write(kept, "filtered.fasta", "fasta")
print("wrote:", count)

Expected output:

wrote: 3

filtered.fasta now contains:

>contig_001
ATGGCGTTCGTA
>contig_002
GGCCGGCCNNNN
>contig_004
GCGCGCGCGCGCGCGCGC

SeqIO.write returns the number of records written, which is the quickest cross-check available. One behavior to know: if the selection is empty, SeqIO.write raises ValueError: No records found in handle rather than producing an empty file. Test len(kept) before writing whenever an empty result is plausible.

Step 6: Produce a pass/fail report

An audit trail matters more than the filtered file itself once a real dataset is involved:

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

MIN_LEN, MAX_LEN = 10, 20000
MIN_GC, MAX_GC = 40.0, 100.0

with open("report.tsv", "w") as out:
    out.write("id\tlength\tgc_percent\tstatus\n")
    for record in SeqIO.parse("demo.fasta", "fasta"):
        length = len(record.seq)
        gc = gc_fraction(record.seq)
        passed = MIN_LEN <= length <= MAX_LEN and MIN_GC <= gc <= MAX_GC
        status = "pass" if passed else "fail"
        out.write(f"{record.id}\t{length}\t{gc:.2f}\t{status}\n")
        print(record.id, length, f"{gc:.2f}", status, sep="\t")

Expected output:

contig_001	12	50.00	pass
contig_002	12	66.67	pass
contig_003	10	0.00	fail
contig_004	18	100.00	pass

The report.tsv file contains the same rows plus a header line. Keep these reports with your project notes. They answer the question that always arrives later: why is this sequence missing?

A Reusable Command-Line Template

The following script, filter_seqs.py, packages the walkthrough into a template you can drop into any project:

#!/usr/bin/env python3
"""Filter sequence records by length and GC content."""
import argparse
import gzip

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

def open_input(path):
    if path.endswith(".gz"):
        return gzip.open(path, "rt")
    return open(path, "rt")

def main():
    parser = argparse.ArgumentParser(description="Filter records by length and GC content.")
    parser.add_argument("input", help="input sequence file")
    parser.add_argument("output", help="output file for records that pass")
    parser.add_argument("--min-length", type=int, default=0, help="minimum length in bases")
    parser.add_argument("--max-length", type=int, default=None, help="maximum length in bases")
    parser.add_argument("--min-gc", type=float, default=0.0, help="minimum GC percent (0 to 100)")
    parser.add_argument("--max-gc", type=float, default=100.0, help="maximum GC percent (0 to 100)")
    parser.add_argument("--ambiguous", choices=["ignore", "remove", "weighted"], default="ignore",
                        help="how ambiguous bases enter the GC calculation")
    parser.add_argument("--format", default="fasta", help="SeqIO format for input and output")
    parser.add_argument("--report", help="optional tab-separated report path")
    args = parser.parse_args()

    kept, dropped = [], []

    with open_input(args.input) as handle:
        for record in SeqIO.parse(handle, args.format):
            length = len(record.seq)
            gc = gc_fraction(record.seq, ambiguous=args.ambiguous)
            entry = (record, length, gc)
            passes = (
                length >= args.min_length
                and (args.max_length is None or length <= args.max_length)
                and args.min_gc <= gc <= args.max_gc
            )
            (kept if passes else dropped).append(entry)

    count = SeqIO.write([entry for entry in kept], args.output, args.format)
    print(f"kept {count} of {count + len(dropped)} records")

    if args.report:
        with open(args.report, "w") as out:
            out.write("id\tlength\tgc_percent\tstatus\n")
            for record, length, gc in kept:
                out.write(f"{record.id}\t{length}\t{gc:.2f}\tpass\n")
            for record, length, gc in dropped:
                out.write(f"{record.id}\t{length}\t{gc:.2f}\tfail\n")

if __name__ == "__main__":
    main()

Run it on the demo file:

python filter_seqs.py demo.fasta filtered.fasta --min-length 10 --min-gc 40 --report report.tsv

Expected output:

kept 3 of 4 records

On a real assembly you would use something closer to:

python filter_seqs.py assembly.fasta contigs_1kb.fasta --min-length 1000 --report dropped.tsv

The -max-length default of None means no upper bound, so a minimum-only filter needs no extra argument.

Choosing the format string

The -format value is any name SeqIO understands, and the template uses it for both reading and writing. The common ones:

• fasta: plain DNA or protein records, the format used throughout this guide.
• fastq: sequencing reads with per-base quality scores. Legacy score encodings have their own names, fastq-solexa and fastq-illumina.
• genbank: annotated records with features. Reading genbank and writing fasta in one pass requires a small change to the template, since the same format name is used in both directions.
• clustal and phylip: aligned formats, useful when you want length and GC measured on an alignment.

Files whose names end in .gz are opened through Python's gzip module automatically, so compressed inputs need no preprocessing.

Error messages you will meet, and the fixes

• FileNotFoundError: [Errno 2] No such file or directory: 'demo.fasta'. The path is misspelled or your working directory is not where you think it is. Print Path.cwd() and the full path before opening.
• ImportError: cannot import name 'gc_fraction' from 'Bio.SeqUtils'. Biopython is too old. Run pip install -U biopython, or use the manual formula from the verification section.
• ValueError: Unknown format 'fa'. The format string is not a SeqIO name. Use fasta, not fa or .fasta.
• ValueError: No records found in handle. You passed an empty selection to SeqIO.write. Check the count first and emit a clear warning instead of a crash buried in a pipeline log.

Key Options and Formulas at a Glance

Ambiguous Bases, Lowercase Letters, and Gaps

Three details change the numbers, and all three are under your control.

Ambiguity codes. IUPAC ambiguity codes stand for more than one base: N means any of A, T, G, C, S means G or C, R means A or G, and so on. Sequencing gaps and low-confidence calls appear as runs of N. Biopython's gc_fraction takes an ambiguous argument with three modes:

• ignore (the default): ambiguous letters are not counted as G or C but stay in the denominator.
• remove: ambiguous letters are stripped from the sequence before either number is computed.
• weighted: each ambiguous letter contributes the fraction of its possible bases that are G or C, so S counts fully, N counts as 0.5, and W counts as zero.

The demo record contig_002 (GGCCGGCCNNNN) shows the spread:

None of these is wrong. A contig with many Ns is less informative, and remove says so by reporting GC over usable bases only. Pick one mode per project and record the choice in your methods.

Lowercase letters. Repeat-masked FASTA files report repeats in lowercase, a convention called soft-masking. Biopython's functions handle mixed case, but plain string counting does not:

seq = "atggcgttcgta"
seq.count("G") + seq.count("C")   # 0, wrong
gc_fraction(seq)                  # 50.0, correct

If you compute GC manually, call seq.upper() first.

Gap characters. Alignments carry - columns, and len counts them:

aligned = "ATGG--CGTTC"
len(aligned)                    # 11
len(aligned.replace("-", ""))   # 9

GC content on the gapped string is 5/11 = 45.45 percent, and on the ungapped string 5/9 = 55.56 percent. For unaligned genomes the question never arises, but for alignments decide explicitly whether gaps count as sequence.

Two edge cases round this out. An empty sequence returns 0.0 from gc_fraction, while a manual division raises ZeroDivisionError: division by zero, so guard the manual version with a length check. And GC content has no meaning for protein sequences, where the letters are amino acids. Filter proteins by length alone.

Scaling Up: Large Files and Batch Jobs

SeqIO.parse returns a generator, which yields one record at a time so only one record is resident in memory. The walkthrough wrapped it in list() for clarity. On a real assembly with tens of thousands of contigs, or a read set with millions of records, stream instead:

def passing_records(handle, min_len, min_gc):
    for record in SeqIO.parse(handle, "fasta"):
        if len(record.seq) >= min_len and gc_fraction(record.seq) >= min_gc:
            yield record

kept = passing_records(open("assembly.fasta"), 1000, 30.0)
SeqIO.write(kept, "contigs_1kb.fasta", "fasta")

Memory use stays flat. The tradeoff is that a generator can be consumed once, so compute the report in the same pass or re-read the file.

To process a directory of files, loop over paths and rerun the same logic:

from pathlib import Path
from Bio import SeqIO

for path in Path("genomes").glob("*.fasta"):
    lengths = [len(r.seq) for r in SeqIO.parse(path, "fasta")]
    print(path.name, len(lengths), "records")

This is also the natural place to say that hand-rolled loops are not the ceiling. Automated Python pipelines now run entire comparative analyses, including GC content, across batches of dozens of genomes in a single reproducible run [5]. Performance-focused packages push further still, with recent work demonstrating orders-of-magnitude runtime improvements through GPU acceleration and optimized core operations [6]. For the filtering task itself, a laptop handles genome-scale FASTA files in seconds to minutes, and dedicated command-line tools become attractive only at read-set scale.

How to Check That Your Filter Worked

Three checks catch nearly every filtering bug.

First, recompute the criteria on the output file. The output should satisfy the rules by construction, and if it does not, something is wrong with the write step:

from Bio import SeqIO
from Bio.SeqUtils import gc_fraction

kept = list(SeqIO.parse("filtered.fasta", "fasta"))
assert len(kept) == 3

for record in kept:
    assert len(record.seq) >= 10
    assert gc_fraction(record.seq) >= 40.0

print("all checks passed")

Expected output:

all checks passed

Second, confirm the totals. Records kept plus records dropped must equal records read. The template prints this as kept 3 of 4 records. If the arithmetic fails, records are being skipped or duplicated.

Third, spot-check one record by hand. contig_002 has 12 characters: 4 G, 4 C, and 4 N. With the default mode, GC = 8/12 = 66.67 percent, matching the report. A one-record hand check takes a minute and validates the logic of the entire pipeline.

As a final cross-check, compute GC two ways and compare. Differences between the Biopython value and the manual formula below point directly at ambiguity content:

def gc_percent_manual(seq):
    seq = seq.upper()
    gc = seq.count("G") + seq.count("C")
    total = sum(seq.count(base) for base in "ATGC")
    return 100.0 * gc / total if total else 0.0

for record in SeqIO.parse("demo.fasta", "fasta"):
    print(record.id,
          f"{gc_fraction(record.seq):.2f}",
          f"{gc_percent_manual(record.seq):.2f}")

Expected output:

contig_001 50.00 50.00
contig_002 66.67 100.00
contig_003 0.00 0.00
contig_004 100.00 100.00

The two methods agree everywhere except contig_002, whose N bases sit in the denominator of one calculation and not the other. Agreement everywhere else confirms both implementations.

Common Mistakes and Limitations

1. Assuming gc_fraction returns a fraction. It returns a percentage. The legacy GC function returned a fraction between 0 and 1, and old snippets still circulate. A threshold written for one convention silently passes or rejects everything when applied to the other. Print one value before trusting a run.

1. Case-sensitive manual counting. seq.count("G") misses lowercase g in soft-masked files, reporting 0 percent GC on a perfectly normal sequence. Uppercase first or use gc_fraction.

1. Forgetting that N stays in the denominator by default. A draft assembly full of N runs looks artificially GC-poor under the ignore mode. Choose ignore, remove, or weighted deliberately and state the choice in your methods.

1. Counting gap characters in length. In aligned FASTA, len includes every - column. Strip gaps before measuring unless you genuinely want alignment width.

1. Writing an empty selection. SeqIO.write raises ValueError: No records found in handle when the record list is empty. Check the count first and fail loudly with a clear message.

1. Overwriting the input file. Give the output a different name. A filter that destroys its own input cannot be debugged or rerun.

1. Reading GC content as a verdict on its own. Composition varies across genome types and even between strands, so a single GC percentage is a summary, not a diagnosis. Treat an out-of-range value as a flag that triggers closer inspection with alignment and taxonomic tools.

1. Choosing thresholds blind. Cutoffs copied from another project can delete good data. Tabulate or plot the length and GC distributions of your own file first, then place cutoffs in the gaps you can actually see.

1. Parsing FASTA with line-based assumptions. Sequences wrap at arbitrary widths. A parser that assumes one sequence per line corrupts data silently. Use SeqIO, or if you must parse manually, accumulate lines until the next > header.

1. Applying contig rules to reads. For FASTQ reads, per-base quality usually matters more than GC, and read length distributions differ from contig length distributions by orders of magnitude. Match the filter to the data type.

One limitation sits above all of these: filtering removes records that fail your rules, nothing more. It does not repair misassemblies, identify contaminants, or certify completeness. Those questions need alignment, annotation, and taxonomic checks downstream.

Where to Go Next

Sliding-window GC. Whole-genome GC averages hide local extremes. Compute GC in fixed-size windows (for example 100 bases) to find regions that would block primer design or gene synthesis, the exact problem Visual codon addresses after codon optimization [4].

Strand skews. AT and GC skews are two lines of Python on top of the counting you already have, and they reveal strand-asymmetric composition that a single GC percentage hides [3].

Segmentation. When the boundaries between compositional domains matter, hidden Markov models and partitioning algorithms assign segments formally. The Sarment package provides Python modules built for exactly this kind of sequence segmentation [1].

Variant-level work. Once sequences pass your filters, variant calling and reporting follow. The hgvs library parses, validates, maps, and formats variants under HGVS nomenclature and is a standard Python tool for that layer [7].

Pipeline thinking. Wrap your filter, its thresholds, and its report into one script with a fixed interface, as the template above does. Reproducibility comes from the recorded parameters and the report, not from memory.

Frequently Asked Questions

How do I filter sequences by length and GC content in Python?

Read the file with Bio.SeqIO.parse, test each record with len(record.seq) and Bio.SeqUtils.gc_fraction against your cutoffs, and write the passing records with SeqIO.write. The six-step walkthrough above shows every line of code with its expected output.

Does Biopython's gc_fraction return a fraction or a percentage?

It returns a percentage between 0 and 100. The older Bio.SeqUtils.GC function returned a fraction between 0 and 1, which is why older scripts multiply by 100.

How are ambiguous bases like N handled in GC calculations?

By default, gc_fraction leaves N and other ambiguity codes in the denominator without counting them as G or C. Pass ambiguous="remove" to strip them from both numerator and denominator, or ambiguous="weighted" to score them fractionally.

Can I use the same code to filter FASTQ files?

Yes. Change the format string from fasta to fastq and the same length and GC logic applies, though for sequencing reads you normally also filter on per-base quality scores.

Why does my GC content differ from another tool's result?

The difference almost always comes from how ambiguous bases, lowercase soft-masked letters, or gap characters are counted in the denominator. Match the ambiguity mode and the gap handling, and the numbers will agree.

How do I save the sequences that failed the filter?

Append failing records to a second list while you loop, then write that list to its own file with SeqIO.write. The template in this guide also writes a report listing every dropped record with its measurements.

Can I filter sequences without installing Biopython?

Yes. Plain Python can read a FASTA file and count characters, and the manual formula in this guide needs no libraries. You take on the work of handling multiline records, mixed case, and format quirks yourself.

What length and GC thresholds should I use?

Derive them from a reference for your organism or marker, then inspect the length and GC distributions of your own file and place cutoffs in visible gaps. A threshold that suits one genome type can be wrong for another.

Sources

1. Sarment: Python modules for HMM analysis and partitioning of sequences.
2. Bio Python Application for Comparative Analysis of COVID-19 Virus Genome with Other Calamitous Genome
3. Complete mitochondrial genome sequence from an endangered Indian snake, Python molurus molurus (Serpentes, Pythonidae).
4. Visual codon: a user-friendly Python program for viewing and optimizing gene GC content.
5. CGAS (Chloroplast Genome Analysis Suite): An automated python pipeline for comprehensive comparative chloroplast genomics
6. pyVIPER: a fast and scalable Python package for protein activity estimation and master regulator analysis of single-cell RNA sequencing data.
7. A Python package for parsing, validating, mapping and formatting sequence variants using HGVS nomenclature.