Bio Sequence Slicing

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 Slicing

Extract sub-regions, splice non-contiguous regions, and concatenate sequences and annotated records.

"Extract a subsequence" -> Slice a Seq with 0-based half-open coordinates.

• Python: seq[start:end] (Bio.Seq)

"Pull out a sub-region but keep its quality scores and features" -> Slice the SeqRecord, not the bare Seq.

• Python: record[start:end] (Bio.SeqRecord)

"Splice exons into a transcript" -> Extract each region and concatenate.

• Python: sum((seq[s:e] for s, e in coords), Seq('')) or the + operator

The Governing Principle

Slicing a bare Seq is pure string math: seq[start:end] returns a new Seq, half-open, with no metadata to lose. Slicing a SeqRecord carries metadata, and the rule for WHAT survives is the single most error-prone part of this skill:

record[start:end] (verified against Bio/SeqRecord.__getitem__):

• PRESERVES id, name, description, and molecule_type.
• AUTO-SLICES letter_annotations (per-letter data such as PHRED phred_quality) to match the new coordinates -- this is why a FASTQ slice keeps the right per-base qualities for free.
• KEEPS only features FULLY CONTAINED in [start:end]; their locations are recalculated relative to the new start.
• SILENTLY DROPS the annotations dict (organism, taxonomy, references, comments), the dbxrefs list, and any feature that STRADDLES the slice boundary (dropped whole, never truncated). A non-trivial stride (record[::2]) drops features entirely.

Nothing warns when annotations vanish. The GenBank source feature spans the whole record, so it straddles almost any slice and disappears along with organism/taxonomy. To carry that metadata across, copy it explicitly:

sub = record[start:end]
sub.annotations = record.annotations.copy()

and re-add any boundary-straddling feature manually (with a clamped, recalculated location) if a truncated copy is needed.

Required Imports

from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO

Coordinate Systems: the 0-based vs 1-based trap

Python and Biopython slicing is 0-based and half-open: seq[start:end] includes start, excludes end, and returns end - start letters. File formats disagree, and mixing them is a SILENT off-by-one (no error, just the wrong bases):

The asymmetry is the catch: convert a 1-based inclusive interval by subtracting 1 from the START only; the end already lands correctly because Python's exclusive end cancels the inclusive end. Reading a coordinate straight off a GFF and slicing seq[start:end] without the -1 silently shifts everything one base left.

Bio.SeqFeature locations sidestep this entirely: they store a 0-based start and a Python-style end, so int(feature.location.start):int(feature.location.end) slices the parent directly, and feature.extract(record.seq) does the same automatically (handling strand and compound/joined locations).

def extract_1based(seq, start, end):
    '''Extract a 1-based inclusive interval (GenBank/GFF style).'''
    return seq[start - 1:end]

Slicing a Bare Seq

Slicing returns a Seq (not a string); negative indices and strides behave exactly like str (Seq has behaved like str since BioPython 1.78).

seq = Seq('ATGCGATCGATCG')
seq[0]       # 'A'  single base, 0-indexed -> returns a str
seq[-1]      # 'G'  last base
seq[0:3]     # Seq('ATG')   first 3 bases
seq[-5:]     # Seq('GATCG')  last 5
seq[::2]     # Seq('AGGTGTG')  every 2nd base (stride)
seq[::-1]    # Seq('GCTAGCTAGCGTA')  reversed (not the reverse complement)

str(record.seq) returns the raw string, but raises UndefinedSequenceError when the record's sequence content is undefined (e.g. Seq(None, length=n) from a header-only FASTA or a pysam-backed record). Guard with len() (always defined) before forcing the content to a string.

Code Patterns

Splice Non-Contiguous Regions (Exons -> Transcript)

Goal: Join several separated regions of a genomic sequence into one continuous sequence.

Approach: Extract each region with half-open coordinates and concatenate. sum() needs an explicit Seq('') start value because the default 0 cannot be added to a Seq.

def extract_regions(seq, regions):
    '''Concatenate multiple [start, end) regions in order.'''
    return sum((seq[start:end] for start, end in regions), Seq(''))

exon_coords = [(0, 50), (100, 150), (200, 250)]
mrna = extract_regions(genomic_seq, exon_coords)

For a real annotated transcript, let the feature do the work -- feature.extract honors strand and joined exon locations:

for feature in record.features:
    if feature.type == 'mRNA':
        transcript = feature.extract(record.seq)