Calculate alignment statistics including sequence identity, conservation scores, substitution matrices, and similarity metrics. Use when comparing alignment quality, measuring sequence divergence, and analyzing evolutionary patterns.
Instalar com Codex ou Claude Copie este prompt, cole no Codex, Claude ou outro assistente e deixe que ele revise a página da skill e instale para você.
Um comando direto ignora o prompt de revisão. Verifique a origem antes de executá-lo.
Instruções da origem · Visualização somente leitura
name
bio-alignment-msa-statistics
description
Calculate alignment statistics including sequence identity, conservation scores, substitution matrices, and similarity metrics. Use when comparing alignment quality, measuring sequence divergence, and analyzing evolutionary patterns.
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.
MSA Statistics
Calculate sequence identity, conservation scores, substitution counts, and other alignment metrics.
Required Import
Goal: Load modules for alignment I/O, substitution scoring, and statistical calculations.
Approach: Import AlignIO for reading alignments, Counter for column analysis, numpy for matrix operations, and math for entropy calculations.
from Bio import AlignIO
from Bio.Align import substitution_matrices
collections Counter
numpy np
math
from
import
import
as
import
Pairwise Identity
"Calculate percent identity" → Compute the fraction of identical aligned residues between sequence pairs.
Goal: Measure sequence similarity as percent identity for individual pairs or across all sequences in an alignment.
Approach: Count matching non-gap positions divided by total aligned positions; optionally compute a full N-by-N identity matrix.
Calculate Identity Between Two Sequences
defpairwise_identity(seq1, seq2):
matches = sum(a == b and a != '-'for a, b inzip(seq1, seq2))
aligned_positions = sum(a != '-'or b != '-'for a, b inzip(seq1, seq2))
return matches / aligned_positions if aligned_positions > 0else0
alignment = AlignIO.read('alignment.fasta', 'fasta')
seq1, seq2 = str(alignment[0].seq), str(alignment[1].seq)
identity = pairwise_identity(seq1, seq2)
print(f'Identity: {identity * 100:.1f}%')
Identity Matrix for All Sequences
defidentity_matrix(alignment):
n = len(alignment)
matrix = np.zeros((n, n))
for i inrange(n):
for j inrange(i, n):
seq_i = str(alignment[i].seq)
seq_j = str(alignment[j].seq)
ident = pairwise_identity(seq_i, seq_j)
matrix[i, j] = matrix[j, i] = ident
return matrix
alignment = AlignIO.read('alignment.fasta', 'fasta')
mat = identity_matrix(alignment)
seq_ids = [r.idfor r in alignment]
print('Pairwise Identity Matrix:')
print(f'{"":>10}', ' '.join(f'{s[:8]:>8}'for s in seq_ids))
for i, row inenumerate(mat):
print(f'{seq_ids[i][:10]:>10}', ' '.join(f'{v*100:>7.1f}%'for v in row))
Conservation Score
Goal: Quantify per-column and overall alignment conservation to identify conserved and variable regions.
Approach: Calculate the fraction of the most common residue at each column, optionally ignoring gaps, and smooth with a sliding window.
defconservation_profile(alignment, window=10):
profile = []
for i inrange(alignment.get_alignment_length()):
start = max(0, i - window // 2)
end = min(alignment.get_alignment_length(), i + window // 2)
scores = [column_conservation(alignment, j) for j inrange(start, end)]
profile.append(sum(scores) / len(scores))
return profile
profile = conservation_profile(alignment, window=10)
Substitution Counts
Goal: Tabulate observed substitution frequencies from the alignment for evolutionary analysis or custom scoring matrices.
Approach: Enumerate all pairwise non-gap character comparisons at each column and tally substitution pairs.
Count Substitutions from Alignment
defsubstitution_counts(alignment):
from collections import defaultdict
counts = defaultdict(int)
for col_idx inrange(alignment.get_alignment_length()):
column = alignment[:, col_idx]
chars = [c for c in column if c != '-']
for i, c1 inenumerate(chars):
for c2 in chars[i+1:]:
if c1 != c2:
pair = tuple(sorted([c1, c2]))
counts[pair] += 1returndict(counts)
subs = substitution_counts(alignment)
print('Substitution counts:')
for pair, count insorted(subs.items(), key=lambda x: -x[1])[:10]:
print(f' {pair[0]}<->{pair[1]}: {count}')
Build Substitution Matrix from MSA
defbuild_substitution_matrix(alignment):
from collections import defaultdict
matrix = defaultdict(lambda: defaultdict(int))
for col_idx inrange(alignment.get_alignment_length()):
column = alignment[:, col_idx]
chars = [c for c in column if c != '-']
for c1 in chars:
for c2 in chars:
matrix[c1][c2] += 1return {k: dict(v) for k, v in matrix.items()}
sub_matrix = build_substitution_matrix(alignment)
Using Alignment.substitutions (Pairwise Alignments)
For pairwise alignments created with PairwiseAligner, use the built-in .substitutions property:
defgap_statistics(alignment):
num_seqs = len(alignment)
num_cols = alignment.get_alignment_length()
total_positions = num_seqs * num_cols
total_gaps = sum(str(r.seq).count('-') for r in alignment)
gaps_per_seq = [str(r.seq).count('-') for r in alignment]
gaps_per_col = [alignment[:, i].count('-') for i inrange(num_cols)]
return {
'total_gaps': total_gaps,
'gap_fraction': total_gaps / total_positions,
'gappiest_seq': max(range(num_seqs), key=lambda i: gaps_per_seq[i]),
'gappiest_col': max(range(num_cols), key=lambda i: gaps_per_col[i]),
'gap_free_cols': sum(1for g in gaps_per_col if g == 0),
}
stats = gap_statistics(alignment)
print(f"Total gaps: {stats['total_gaps']}")
print(f"Gap fraction: {stats['gap_fraction']*100:.1f}%")
print(f"Gap-free columns: {stats['gap_free_cols']}")
Alignment Quality Metrics
Goal: Score alignment quality using sum-of-pairs or simple match/mismatch/gap scoring across all columns.
Approach: For each column, score all pairwise residue comparisons and sum across the alignment.
Overall Alignment Score
defalignment_score(alignment, match=1, mismatch=-1, gap=-2):
total_score = 0for col_idx inrange(alignment.get_alignment_length()):
column = alignment[:, col_idx]
for i, c1 inenumerate(column):
for c2 in column[i+1:]:
if c1 == '-'or c2 == '-':
total_score += gap
elif c1 == c2:
total_score += matchelse:
total_score += mismatch
return total_score
score = alignment_score(alignment)
print(f'Alignment score: {score}')
Sum of Pairs Score
defsum_of_pairs(alignment, substitution_matrix=None):
if substitution_matrix isNone:
substitution_matrix = substitution_matrices.load('BLOSUM62')
total = 0for col_idx inrange(alignment.get_alignment_length()):
column = alignment[:, col_idx]
for i, c1 inenumerate(column):
for c2 in column[i+1:]:
if c1 != '-'and c2 != '-':
total += substitution_matrix.get((c1, c2), 0)
return total
Position-Specific Score Matrix (PSSM)
Goal: Build a position-specific score matrix (PSSM) from the alignment for motif analysis or sequence scoring.
Approach: Count non-gap character frequencies at each column, producing a list of per-position dictionaries.
defposition_specific_score_matrix(alignment):
pssm = []
for col_idx inrange(alignment.get_alignment_length()):
column = alignment[:, col_idx]
counts = Counter(column)
if'-'in counts:
del counts['-']
pssm.append(dict(counts))
return pssm
alignment = AlignIO.read('alignment.fasta', 'fasta')
pssm = position_specific_score_matrix(alignment)
for i, row inenumerate(pssm[:10]):
print(f'Position {i}: {row}')
Note on Bio.Align.AlignInfo
The AlignInfo.SummaryInfo class is deprecated in recent Biopython versions. Use the custom functions in this skill instead:
For PSSM: use position_specific_score_matrix() above
For information content: use information_content() function earlier in this skill
For consensus: see msa-parsing skill
Quick Reference: Metrics
Metric
Description
Range
Identity
Fraction of identical residues
0-1
Conservation
Most common residue frequency
0-1
Shannon Entropy
Variability measure
0 to log2(alphabet)
Information Content
Max entropy - observed entropy
0 to log2(alphabet)
Gap Fraction
Proportion of gaps
0-1
Common Errors
Error
Cause
Solution
ZeroDivisionError
Empty column after gap removal
Check for gap-only columns
KeyError
Character not in substitution matrix
Handle gaps separately
Negative IC
Wrong alphabet size
Use 4 for DNA, 20 for protein
Related Skills
msa-parsing - Parse and manipulate alignments
alignment-io - Read/write alignment files
pairwise-alignment - Create and score pairwise alignments
