bio-alignment-pairwise

FreedomIntelligence/OpenClaw-Medical-SkillsBiology Medicine and BioinformaticsAudit passedSnapshot b8f55bee95f3

Installs: 0

Summary

Perform pairwise sequence alignment using Biopython Bio.Align.PairwiseAligner. Use when comparing two sequences, finding optimal alignments, scoring similarity, and identifying local or global matches between DNA, RNA, or protein sequences.

SKILL.md

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.

Pairwise Sequence Alignment

"Align two sequences" → Compute an optimal alignment between a pair of sequences using dynamic programming.

Python: PairwiseAligner() (BioPython Bio.Align)
CLI: needle (global) or water (local) from EMBOSS
R: pairwiseAlignment() (Biostrings)

Align two sequences using dynamic programming algorithms (Needleman-Wunsch for global, Smith-Waterman for local).

Required Import

Goal: Load modules needed for pairwise alignment operations.

Approach: Import the PairwiseAligner class along with sequence and I/O utilities from Biopython.

from Bio.Align import PairwiseAligner
from Bio.Seq import Seq
from Bio import SeqIO

Core Concepts

Mode	Algorithm	Use Case
`global`	Needleman-Wunsch	Full-length alignment, similar-length sequences
`local`	Smith-Waterman	Find best matching regions, different-length sequences

Creating an Aligner

Configure a PairwiseAligner with appropriate scoring for the sequence type.

# Basic aligner with defaults
aligner = PairwiseAligner()

# Configure mode and scoring
aligner = PairwiseAligner(mode='global', match_score=2, mismatch_score=-1, open_gap_score=-10, extend_gap_score=-0.5)

# For protein alignment with substitution matrix
from Bio.Align import substitution_matrices
aligner = PairwiseAligner(mode='global', substitution_matrix=substitution_matrices.load('BLOSUM62'))

seq1 = Seq('ACCGGTAACGTAG')
seq2 = Seq('ACCGTTAACGAAG')

# Get all optimal alignments
alignments = aligner.align(seq1, seq2)
print(f'Found {len(alignments)} optimal alignments')
print(alignments[0])  # Print first alignment

# Get score only (faster for large sequences)
score = aligner.score(seq1, seq2)

target            0 ACCGGTAACGTAG 13
                  0 |||||.||||.|| 13
query             0 ACCGTTAACGAAG 13

alignment = alignments[0]

# Basic properties
print(alignment.score)                    # Alignment score
print(alignment.shape)                    # (num_seqs, alignment_length)
print(len(alignment))                     # Alignment length

# Get aligned sequences with gaps
target_aligned = alignment[0, :]          # First sequence (target) with gaps
query_aligned = alignment[1, :]           # Second sequence (query) with gaps

# Get coordinate mapping
print(alignment.aligned)                  # Array of aligned segment coordinates
print(alignment.coordinates)              # Full coordinate array

alignment = alignments[0]
counts = alignment.counts()

print(f'Identities: {counts.identities}')
print(f'Mismatches: {counts.mismatches}')
print(f'Gaps: {counts.gaps}')

# Calculate percent identity
total_aligned = counts.identities + counts.mismatches
percent_identity = counts.identities / total_aligned * 100
print(f'Percent identity: {percent_identity:.1f}%')

aligner = PairwiseAligner(mode='global', match_score=2, mismatch_score=-1, open_gap_score=-10, extend_gap_score=-0.5)

from Bio.Align import substitution_matrices
blosum62 = substitution_matrices.load('BLOSUM62')
aligner = PairwiseAligner(mode='global', substitution_matrix=blosum62, open_gap_score=-11, extend_gap_score=-1)

aligner = PairwiseAligner(mode='local', match_score=2, mismatch_score=-1, open_gap_score=-10, extend_gap_score=-0.5)

# Allow free end gaps on query (useful for primer alignment)
aligner = PairwiseAligner(mode='global')
aligner.query_left_open_gap_score = 0
aligner.query_left_extend_gap_score = 0
aligner.query_right_open_gap_score = 0
aligner.query_right_extend_gap_score = 0

from Bio.Align import substitution_matrices
print(substitution_matrices.load())  # List all available matrices

# Common matrices
blosum62 = substitution_matrices.load('BLOSUM62')  # General protein
blosum80 = substitution_matrices.load('BLOSUM80')  # Closely related proteins
pam250 = substitution_matrices.load('PAM250')      # Distantly related proteins

from Bio import SeqIO

records = list(SeqIO.parse('sequences.fasta', 'fasta'))
seq1, seq2 = records[0].seq, records[1].seq

aligner = PairwiseAligner(mode='global', match_score=1, mismatch_score=-1)
alignments = aligner.align(seq1, seq2)

# Limit number of alignments returned (memory efficient)
aligner.max_alignments = 100

for i, alignment in enumerate(alignments):
    print(f'Alignment {i+1}: score={alignment.score}')
    if i >= 4:
        break

alignment = alignments[0]
substitutions = alignment.substitutions

# View as array (rows=target, cols=query)
print(substitutions)

# Access specific substitution counts
# substitutions['A', 'T'] gives count of A aligned to T

alignment = alignments[0]

# Various output formats
print(format(alignment, 'fasta'))     # FASTA format
print(format(alignment, 'clustal'))   # Clustal format
print(format(alignment, 'psl'))       # PSL format (BLAT)
print(format(alignment, 'sam'))       # SAM format

Parameter	Description	Typical DNA	Typical Protein
`match_score`	Score for identical bases	1-2	Use matrix
`mismatch_score`	Penalty for mismatches	-1 to -3	Use matrix
`open_gap_score`	Cost to start a gap	-5 to -15	-10 to -12
`extend_gap_score`	Cost per gap extension	-0.5 to -2	-0.5 to -1
`substitution_matrix`	Scoring matrix	N/A	BLOSUM62

Error	Cause	Solution
`OverflowError`	Too many optimal alignments	Set `aligner.max_alignments`
Low scores	Wrong scoring scheme	Use substitution matrix for proteins
No alignments in local mode	Scores all negative	Ensure `match_score` > 0

Need full-length comparison?
├── Yes → Use mode='global'
│   └── Sequences similar length?
│       ├── Yes → Standard global
│       └── No → Consider semiglobal (free end gaps)
└── No → Use mode='local'
    └── Find best matching regions only

bio-alignment-pairwise

Summary

SKILL.md

Version Compatibility

Pairwise Sequence Alignment

Required Import

Core Concepts

Creating an Aligner

Performing Alignments

Alignment Output Format

Accessing Alignment Data

Alignment Counts (Identities, Mismatches, Gaps)

Common Scoring Configurations

DNA/RNA Alignment

Protein Alignment

Local Alignment (Find Best Region)

Semiglobal (Overlap/Extension)

Available Substitution Matrices

Working with SeqRecord Objects

Iterating Over Multiple Alignments

Substitution Matrix from Alignment

Export Alignment to Different Formats

Quick Reference: Scoring Parameters

Common Errors

Decision Tree: Choosing Alignment Mode

Related Skills

Related skills