msa-statistics

Name: msa-statistics
Author: GPTomics

1.1kstars185forks

How to install msa-statistics?

gptomics/bioskills/msa-statistics

$npx -y skills add GPTomics/bioSkills --skill msa-statistics

Installs into the current project.

Summary

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.

Files · 1

View on GitHub

SKILL.md

1	## Version Compatibility
2
3	Reference examples tested with: BioPython 1.83+, numpy 1.26+
4
5	Before using code patterns, verify installed versions match. If versions differ:
6	- Python: `pip show <package>` then `help(module.function)` to check signatures
7
8	If code throws ImportError, AttributeError, or TypeError, introspect the installed
9	package and adapt the example to match the actual API rather than retrying.
10
11	# MSA Statistics
12
13	Calculate sequence identity, conservation scores, substitution counts, and other alignment metrics.
14
15	## Required Import
16
17	Goal: Load modules for alignment I/O, substitution scoring, and statistical calculations.
18
19	Approach: Import AlignIO for reading alignments, Counter for column analysis, numpy for matrix operations, and math for entropy calculations.
20
21	```python
22	from Bio import AlignIO
23	from Bio.Align import substitution_matrices
24	from collections import Counter
25	import numpy as np
26	import math
27	```
28
29	## Pairwise Identity
30
31	"Calculate percent identity" -> Compute the fraction of identical aligned residues between sequence pairs.
32
33	Goal: Measure sequence similarity as percent identity for individual pairs or across all sequences in an alignment.
34
35	Approach: Count matching non-gap positions divided by total aligned positions; optionally compute a full N-by-N identity matrix.
36
37	### Percent Identity Definitions
38
39	There are four common denominators, producing up to 11.5% difference on the same alignment. Combined with different alignment algorithms, variation reaches 22%. Always report which method was used.
40
41	\| Method \| Denominator \| Code \|
42	\|--------\|-------------\|------\|
43	\| PID1 \| Aligned positions including internal gaps \| `sum(a != '-' or b != '-' for a, b in zip(s1, s2))` \|
44	\| PID2 \| Aligned residue pairs only (no gaps) \| `sum(a != '-' and b != '-' for a, b in zip(s1, s2))` \|
45	\| PID3 \| Shorter sequence length (ungapped) \| `min(len(s1.replace('-', '')), len(s2.replace('-', '')))` \|
46	\| PID4 \| Mean sequence length (ungapped) \| `(len(s1.replace('-', '')) + len(s2.replace('-', ''))) / 2` \|
47
48	PID2 always gives the highest value; PID4 correlates best with structural similarity (Raghava & Barton 2006 BMC Bioinf, r=0.86 with STAMP's Sc score) and is recommended for evolutionary analyses.
49
50	Length-asymmetry pathology: When sequences differ greatly in length, PID4 and PID2 diverge sharply. Example: 80 matches between a 500-residue protein and a 100-residue domain fragment yields PID2 ~84% (matches over aligned residue pairs) but PID4 ~27% (matches over mean ungapped length). Neither is wrong; they answer different questions:
51	- PID2 -> "how similar is the aligned region?" (motif/domain detection)
52	- PID4 -> "how similar are the full sequences?" (structural similarity benchmarks)
53
54	For ortholog identification at the protein level (full-length, similar size), PID4 is recommended. For domain detection or fragment-vs-genome alignment, PID2 with explicit length annotation is more interpretable. Always report alignment length alongside any percent identity to disambiguate.
55
56	### Calculate Identity Between Two Sequences
57	```python
58	def pairwise_identity(seq1, seq2, method='pid1'):
59	matches = sum(a == b and a != '-' for a, b in zip(seq1, seq2))
60	if method == 'pid1':
61	denom = sum(a != '-' or b != '-' for a, b in zip(seq1, seq2))
62	elif method == 'pid2':
63	denom = sum(a != '-' and b != '-' for a, b in zip(seq1, seq2))
64	elif method == 'pid3':
65	denom = min(len(seq1.replace('-', '')), len(seq2.replace('-', '')))
66	elif method == 'pid4':
67	denom = (len(seq1.replace('-', '')) + len(seq2.replace('-', ''))) / 2
68	return matches / denom if denom > 0 else 0
69
70	alignment = AlignIO.read('alignment.fasta', 'fasta')
71	seq1, seq2 = str(alignment[0].seq), str(alignment[1].seq)
72	for method in ['pid1', 'pid2', 'pid3', 'pid4']:
73	print(f'{method}: {pairwise_identity(seq1, seq2, method) * 100:.1f}%')
74	```
75
76	### Identity Matrix for All Sequences
77
78	The double-loop is O(N^2 * L) and fine for hundreds of sequences; for thousands, vectorize via numpy broadcasting:
79
80	```python
81	def identity_matrix_vectorized(alignment):
82	# Build N x L character array; for each row, broadcast equality and validity masks against all rows
83	...
84	```
85
86	Full implementation: `examples/identity_matrix.py`. For very large alignments (>10k sequences), switch to k-mer-based distance estimation (e.g. mash) -- exact pairwise identity becomes prohibitive.
87
88	## Conservation Scoring Methods
89
90	Pick a conservation score by what the downstream task needs:
91
92	\| Pick this \| When \|
93	\|-----------\|------\|
94	\| Majority fraction or Shannon entropy \| Quick screening; DNA/RNA logos; coarse column QC \|
95	\| Capra-Singh J

People who liked this also viewed

6 picks

Type

msa-statistics

1.1kstars185forks

View on GitHub

How to install msa-statistics?

gptomics/bioskills/msa-statistics

$npx -y skills add GPTomics/bioSkills --skill msa-statistics

Installs into the current project.

Summary

Files · 1

View on GitHub

SKILL.md

1	## Version Compatibility
2
3	Reference examples tested with: BioPython 1.83+, numpy 1.26+
4
5	Before using code patterns, verify installed versions match. If versions differ:
6	- Python: `pip show <package>` then `help(module.function)` to check signatures
7
8	If code throws ImportError, AttributeError, or TypeError, introspect the installed
9	package and adapt the example to match the actual API rather than retrying.
10
11	# MSA Statistics
12
13	Calculate sequence identity, conservation scores, substitution counts, and other alignment metrics.
14
15	## Required Import
16
17	Goal: Load modules for alignment I/O, substitution scoring, and statistical calculations.
18
19	Approach: Import AlignIO for reading alignments, Counter for column analysis, numpy for matrix operations, and math for entropy calculations.
20
21	```python
22	from Bio import AlignIO
23	from Bio.Align import substitution_matrices
24	from collections import Counter
25	import numpy as np
26	import math
27	```
28
29	## Pairwise Identity
30
31	"Calculate percent identity" -> Compute the fraction of identical aligned residues between sequence pairs.
32
33	Goal: Measure sequence similarity as percent identity for individual pairs or across all sequences in an alignment.
34
35	Approach: Count matching non-gap positions divided by total aligned positions; optionally compute a full N-by-N identity matrix.
36
37	### Percent Identity Definitions
38
39	There are four common denominators, producing up to 11.5% difference on the same alignment. Combined with different alignment algorithms, variation reaches 22%. Always report which method was used.
40
41	\| Method \| Denominator \| Code \|
42	\|--------\|-------------\|------\|
43	\| PID1 \| Aligned positions including internal gaps \| `sum(a != '-' or b != '-' for a, b in zip(s1, s2))` \|
44	\| PID2 \| Aligned residue pairs only (no gaps) \| `sum(a != '-' and b != '-' for a, b in zip(s1, s2))` \|
45	\| PID3 \| Shorter sequence length (ungapped) \| `min(len(s1.replace('-', '')), len(s2.replace('-', '')))` \|
46	\| PID4 \| Mean sequence length (ungapped) \| `(len(s1.replace('-', '')) + len(s2.replace('-', ''))) / 2` \|
47
48	PID2 always gives the highest value; PID4 correlates best with structural similarity (Raghava & Barton 2006 BMC Bioinf, r=0.86 with STAMP's Sc score) and is recommended for evolutionary analyses.
49
50	Length-asymmetry pathology: When sequences differ greatly in length, PID4 and PID2 diverge sharply. Example: 80 matches between a 500-residue protein and a 100-residue domain fragment yields PID2 ~84% (matches over aligned residue pairs) but PID4 ~27% (matches over mean ungapped length). Neither is wrong; they answer different questions:
51	- PID2 -> "how similar is the aligned region?" (motif/domain detection)
52	- PID4 -> "how similar are the full sequences?" (structural similarity benchmarks)
53
54	For ortholog identification at the protein level (full-length, similar size), PID4 is recommended. For domain detection or fragment-vs-genome alignment, PID2 with explicit length annotation is more interpretable. Always report alignment length alongside any percent identity to disambiguate.
55
56	### Calculate Identity Between Two Sequences
57	```python
58	def pairwise_identity(seq1, seq2, method='pid1'):
59	matches = sum(a == b and a != '-' for a, b in zip(seq1, seq2))
60	if method == 'pid1':
61	denom = sum(a != '-' or b != '-' for a, b in zip(seq1, seq2))
62	elif method == 'pid2':
63	denom = sum(a != '-' and b != '-' for a, b in zip(seq1, seq2))
64	elif method == 'pid3':
65	denom = min(len(seq1.replace('-', '')), len(seq2.replace('-', '')))
66	elif method == 'pid4':
67	denom = (len(seq1.replace('-', '')) + len(seq2.replace('-', ''))) / 2
68	return matches / denom if denom > 0 else 0
69
70	alignment = AlignIO.read('alignment.fasta', 'fasta')
71	seq1, seq2 = str(alignment[0].seq), str(alignment[1].seq)
72	for method in ['pid1', 'pid2', 'pid3', 'pid4']:
73	print(f'{method}: {pairwise_identity(seq1, seq2, method) * 100:.1f}%')
74	```
75
76	### Identity Matrix for All Sequences
77
78	The double-loop is O(N^2 * L) and fine for hundreds of sequences; for thousands, vectorize via numpy broadcasting:
79
80	```python
81	def identity_matrix_vectorized(alignment):
82	# Build N x L character array; for each row, broadcast equality and validity masks against all rows
83	...
84	```
85
86	Full implementation: `examples/identity_matrix.py`. For very large alignments (>10k sequences), switch to k-mer-based distance estimation (e.g. mash) -- exact pairwise identity becomes prohibitive.
87
88	## Conservation Scoring Methods
89
90	Pick a conservation score by what the downstream task needs:
91
92	\| Pick this \| When \|
93	\|-----------\|------\|
94	\| Majority fraction or Shannon entropy \| Quick screening; DNA/RNA logos; coarse column QC \|
95	\| Capra-Singh J

People who liked this also viewed

6 picks

Type