biopython

What this skill does

# Biopython - Bioinformatics Library

Industry-standard Python library for computational biology and bioinformatics workflows.

## When to Use

- Parsing and manipulating biological sequences (DNA, RNA, protein)
- Reading and writing sequence files (FASTA, FASTQ, GenBank, EMBL, SwissProt)
- Performing sequence alignments (pairwise and multiple)
- Running and parsing BLAST searches
- Analyzing protein structures from PDB files
- Calculating sequence statistics and molecular weights
- Translating DNA to protein sequences
- Finding restriction enzyme sites
- Building and analyzing phylogenetic trees
- Accessing NCBI databases (Entrez, PubMed)
- Computing sequence motifs and patterns
- Analyzing next-generation sequencing data

## Reference Documentation

**Official docs**: https://biopython.org/  
**Tutorial**: https://biopython.org/DIST/docs/tutorial/Tutorial.html  
**Search patterns**: `SeqIO.parse`, `Seq`, `AlignIO`, `NCBIWWW.qblast`, `PDBParser`

## Core Principles

### Use Biopython For

| Task | Module | Example |
|------|--------|---------|
| Create sequences | `Seq` | `Seq("ATCG")` |
| Read sequence files | `SeqIO` | `SeqIO.parse("file.fasta", "fasta")` |
| Pairwise alignment | `pairwise2` | `pairwise2.align.globalxx(s1, s2)` |
| Multiple alignment | `AlignIO` | `AlignIO.read("align.fasta", "fasta")` |
| BLAST searches | `NCBIWWW` | `NCBIWWW.qblast("blastn", "nr", seq)` |
| PDB structures | `PDB.PDBParser` | `PDBParser().get_structure()` |
| Phylogenetic trees | `Phylo` | `Phylo.read("tree.xml", "phyloxml")` |
| NCBI databases | `Entrez` | `Entrez.esearch(db="nucleotide")` |

### Do NOT Use For

- High-performance genome assembly (use SPAdes, Canu)
- Variant calling from BAM files (use GATK, BCFtools)
- RNA-seq differential expression (use DESeq2, edgeR)
- Protein structure prediction (use AlphaFold, RoseTTAFold)
- Large-scale metagenomics (use specialized pipelines)

## Quick Reference

### Installation

```bash
# pip (recommended)
pip install biopython

# With optional dependencies
pip install biopython[extra]

# conda
conda install -c conda-forge biopython

# Development version
pip install git+https://github.com/biopython/biopython.git
```

### Standard Imports

```python
# Core sequence handling
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio import SeqIO, AlignIO

# Sequence alignment
from Bio import pairwise2
from Bio.Align import MultipleSeqAlignment
from Bio.Align.Applications import ClustalwCommandline

# BLAST
from Bio.Blast import NCBIWWW, NCBIXML

# Structure analysis
from Bio.PDB import PDBParser, PDBIO, Select
from Bio.PDB.DSSP import DSSP
from Bio.PDB.Polypeptide import PPBuilder

# Phylogenetics
from Bio import Phylo
from Bio.Phylo.TreeConstruction import DistanceCalculator, DistanceTreeConstructor

# NCBI Entrez
from Bio import Entrez

# Additional tools
from Bio.SeqUtils import GC, molecular_weight
from Bio.Restriction import *
```

### Basic Pattern - Sequence Creation

```python
from Bio.Seq import Seq

# Create DNA sequence
dna = Seq("ATGGCCATTGTAATGGGCCGC")

# Transcribe to RNA
rna = dna.transcribe()

# Translate to protein
protein = dna.translate()

# Reverse complement
rev_comp = dna.reverse_complement()

print(f"DNA:     {dna}")
print(f"RNA:     {rna}")
print(f"Protein: {protein}")
print(f"RevComp: {rev_comp}")
```

### Basic Pattern - File Reading

```python
from Bio import SeqIO

# Read FASTA file (iterator - memory efficient)
for record in SeqIO.parse("sequences.fasta", "fasta"):
    print(f"ID: {record.id}")
    print(f"Length: {len(record.seq)}")
    print(f"Sequence: {record.seq[:50]}...")

# Read single sequence
record = SeqIO.read("single.fasta", "fasta")
```

### Basic Pattern - Sequence Alignment

```python
from Bio import pairwise2
from Bio.Seq import Seq

seq1 = Seq("ACCGT")
seq2 = Seq("ACGT")

# Global alignment
alignments = pairwise2.align.globalxx(seq1, seq2)

# Print best alignment
best = alignments[0]
print(pairwise2.format_alignment(*best))
```

## Critical Rules

### ✅ DO

- **Use iterators for large files** - `SeqIO.parse()` not `SeqIO.to_dict()`
- **Validate sequences** - Check alphabet and length before operations
- **Handle file formats correctly** - Match parser to actual file format
- **Check alignment quality** - Verify gaps and identity percentages
- **Use appropriate genetic code** - Specify table for translation
- **Close file handles** - Use context managers or explicit close
- **Filter FASTQ by quality** - Don't trust all reads equally
- **Verify PDB structure** - Check for missing atoms/residues
- **Set Entrez email** - Required for NCBI API usage
- **Handle translation frames** - Consider all three reading frames

### ❌ DON'T

- **Load entire FASTQ into memory** - Use streaming
- **Ignore sequence type** - DNA, RNA, and protein need different handling
- **Skip quality filtering** - FASTQ quality scores matter
- **Use wrong genetic code** - Different organisms use different tables
- **Forget stop codons** - Handle them explicitly in translation
- **Mix alphabets** - Don't compare DNA with protein sequences
- **Trust all BLAST hits** - Filter by e-value and identity
- **Ignore chain breaks** - PDB structures may have gaps
- **Hammer NCBI servers** - Use rate limiting (3 requests/sec without API key)
- **Compare raw sequences** - Align first, then compare

## Anti-Patterns (NEVER)

```python
# ❌ BAD: Loading entire file into memory
records = list(SeqIO.parse("huge.fastq", "fastq"))
for record in records:
    process(record)  # OOM for large files!

# ✅ GOOD: Stream processing
for record in SeqIO.parse("huge.fastq", "fastq"):
    if min(record.letter_annotations["phred_quality"]) >= 20:
        process(record)

# ❌ BAD: Translating without checking frame
protein = dna_seq.translate()
# May include stop codons or wrong frame!

# ✅ GOOD: Specify frame and stop codon handling
protein = dna_seq.translate(to_stop=True, table=1)

# ❌ BAD: No email for Entrez
Entrez.esearch(db="nucleotide", term="human")
# NCBI will block you!

# ✅ GOOD: Always set email
Entrez.email = "[email protected]"
handle = Entrez.esearch(db="nucleotide", term="human")

# ❌ BAD: Comparing sequences directly
if str(seq1) == str(seq2):
    print("Same")  # Ignores gaps, misalignments!

# ✅ GOOD: Align then compare
from Bio import pairwise2
alignments = pairwise2.align.globalxx(seq1, seq2)
score = alignments[0][2]  # Alignment score

# ❌ BAD: Wrong file format
records = SeqIO.parse("file.gb", "fasta")  # Wrong!

# ✅ GOOD: Match format to file
records = SeqIO.parse("file.gb", "genbank")
```

## Sequence Objects

### Creating and Manipulating Sequences

```python
from Bio.Seq import Seq
from Bio.Alphabet import IUPAC

# DNA sequence
dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")

# Basic operations
length = len(dna)
gc_count = dna.count("G") + dna.count("C")
gc_content = (gc_count / length) * 100

# Find subsequence
position = dna.find("ATG")  # Returns index or -1

# Slicing
first_codon = dna[:3]
last_10 = dna[-10:]

print(f"Length: {length}")
print(f"GC content: {gc_content:.2f}%")
print(f"ATG at position: {position}")
```

### Transcription and Translation

```python
from Bio.Seq import Seq

# DNA to RNA transcription
dna = Seq("ATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG")
rna = dna.transcribe()

# RNA to protein translation
protein = rna.translate()

# DNA to protein (direct)
protein_direct = dna.translate()

# Back-transcription
dna_back = rna.back_transcribe()

# Reverse complement
rev_comp = dna.reverse_complement()

print(f"DNA:        {dna}")
print(f"RNA:        {rna}")
print(f"Protein:    {protein}")
print(f"Rev Comp:   {rev_comp}")
```

### Translation with Different Genetic Codes

```python
from Bio.Seq import Seq

dna = Seq("ATGGGCTAG")

# Standard genetic code (table 1)
protein_standard = dna.translate(table=1)

# Mitochondrial genetic code (table 2)
protein_mito = dna.translate(table=2)

# Stop at first stop codon
protein_to_stop = dna.translate(to_stop=True)

# All three reading fra