SeqFeature and Location

In one line: an annotated region of a sequence — and Biopython stores all locations 0-based half-open regardless of the source format.

The anatomy

for f in rec.features:
    f.type            # 'CDS', 'gene', 'exon', 'source', 'mRNA', ...
    f.location        # SimpleLocation or CompoundLocation
    f.qualifiers      # dict of str → LIST of str
    f.id

    f.location.start  # 0-based
    f.location.end    # exclusive
    f.location.strand # +1, -1, or None
    len(f)            # total length (sums the parts of a join)

qualifiers values are always lists, even for a single value:

f.qualifiers["gene"]              # ['TP53']
f.qualifiers.get("gene", [""])[0] # 'TP53'  — the safe access idiom
f.qualifiers.get("product", ["?"])[0]
f.qualifiers.get("db_xref", [])

Coordinates: 0-based, half-open, always

A GenBank record says CDS 100..200 (1-based, inclusive). Biopython stores [99:200].

print(f.location)          # [99:200](+)
f.location.start           # 99
f.location.end             # 200
len(f)                     # 101

Do not adjust by 1 yourself. Biopython already converted from the source format's convention. This applies to GenBank, EMBL and GFF alike — all 1-based on disk, all 0-based in memory. See Coordinate Systems.

The convenient consequence: rec.seq[f.location.start:f.location.end] is a valid Python slice with no arithmetic.

extract() — the method to use

sub = f.extract(rec.seq)                # or f.extract(rec) for a SeqRecord
protein = sub.translate(cds=True)

extract handles three things you would otherwise get wrong:

1. Joins. A multi-exon CDS is join(100..200,300..400). extract concatenates the exons in the correct order.
2. Strand. A minus-strand feature is reverse-complemented automatically.
3. Order. For a minus-strand join, the exons are assembled in transcription order, not coordinate order.

Manual slicing gets all three wrong in ways that produce a plausible-looking but incorrect protein. Always use extract.

Location types

from Bio.SeqFeature import SeqFeature, SimpleLocation, CompoundLocation

SimpleLocation(99, 200, strand=+1)                  # a contiguous span
CompoundLocation([SimpleLocation(99, 200, 1),        # a join
                  SimpleLocation(299, 400, 1)])

loc.parts                    # list of SimpleLocation
len(loc.parts)               # number of exons
loc.operator                 # 'join' or 'order'

SimpleLocation was called FeatureLocation in older versions; the old name still works as an alias.

Fuzzy positions

GenBank can express uncertainty, and Biopython models it:

from Bio.SeqFeature import ExactPosition, BeforePosition, AfterPosition, UnknownPosition

# GenBank '<100..200' → BeforePosition(99)  — "starts at or before here"
# GenBank '100..>200' → AfterPosition(200)  — "ends at or after here"

int(f.location.start)                       # coerce to a plain int
isinstance(f.location.start, ExactPosition)  # check before trusting it

< and > mark partial features — a CDS running off the end of a contig. If you translate one without checking, you get a truncated protein and no warning. Test isinstance(..., ExactPosition) in pipeline code.

Building features

feat = SeqFeature(
    SimpleLocation(99, 200, strand=+1),
    type="CDS",
    qualifiers={"gene": ["TP53"], "product": ["tumor protein p53"]},
)
rec.features.append(feat)

Searching features

cds = [f for f in rec.features if f.type == "CDS"]
genes = [f for f in rec.features if f.qualifiers.get("gene", [None])[0] == "TP53"]

# features overlapping a position (0-based)
overlapping = [f for f in rec.features if pos in f]         # __contains__ works!
overlapping = [f for f in rec.features
               if f.location.start <= pos < f.location.end]

# features within a window
in_window = [f for f in rec.features
             if f.location.start < end and f.location.end > start]

pos in feature is the neat one — it respects compound locations, so an intronic position correctly reports as not in a spliced CDS.

For many queries against many features, this linear scan is too slow. Build an interval tree (intervaltree, ncls) or use pyranges. See Applied - GFF and Genomic Intervals.

Bioinformatics examples

# extract every CDS as a protein
proteins = []
for rec in SeqIO.parse("genome.gb", "genbank"):
    for f in rec.features:
        if f.type != "CDS":
            continue
        try:
            prot = f.extract(rec.seq).translate(cds=True)
        except Exception as e:
            print(f"skipping {f.qualifiers.get('gene', ['?'])[0]}: {e}")
            continue
        proteins.append(SeqRecord(prot,
                                  id=f.qualifiers.get("locus_tag", ["?"])[0],
                                  description=f.qualifiers.get("product", [""])[0]))

# exon structure of a gene
for f in rec.features:
    if f.type == "mRNA" and f.qualifiers.get("gene", [""])[0] == "TP53":
        for i, part in enumerate(f.location.parts, 1):
            print(f"exon {i}: {part.start}-{part.end} ({len(part)} bp)")

# features → a DataFrame for pandas work
df = pd.DataFrame([
    {"type": f.type,
     "start": int(f.location.start), "end": int(f.location.end),
     "strand": f.location.strand,
     "gene": f.qualifiers.get("gene", [None])[0],
     "product": f.qualifiers.get("product", [None])[0],
     "n_exons": len(f.location.parts)}
    for f in rec.features
])

# which feature contains a variant? (VCF pos is 1-based → subtract 1)
hits = [f for f in rec.features if (vcf_pos - 1) in f and f.type == "CDS"]

# intron coordinates, derived from exons
parts = sorted(cds.location.parts, key=lambda p: p.start)
introns = [(a.end, b.start) for a, b in zip(parts, parts[1:])]

The (vcf_pos - 1) in that second-to-last example is the coordinate conversion that this whole note exists to make explicit. VCF is 1-based; Biopython is 0-based.

Common mistakes

• Adjusting coordinates by 1 when Biopython already did.
• Manual slicing instead of extract(), breaking on joins and minus strands.
• Treating qualifiers values as strings instead of lists.
• Ignoring fuzzy positions and translating a partial CDS.
• Forgetting slicing a record drops straddling features. See SeqRecord.
• Linear feature scans in a loop over millions of variants.
• Assuming exons are in coordinate order on the minus strand — in a CompoundLocation they are in transcription order.
• f.location.strand is None for features with no strand — check before comparing to -1.

See also

SeqRecord · Seq · SeqIO · Coordinate Systems · Applied - GFF and Genomic Intervals · Applied - Genomic File Formats