alignment-indexing

Name: alignment-indexing
Author: GPTomics

1.1kstars185forks

How to install alignment-indexing?

gptomics/bioskills/alignment-indexing

$npx -y skills add GPTomics/bioSkills --skill alignment-indexing

Installs into the current project.

Summary

Create and use BAI/CSI indices for BAM/CRAM files using samtools and pysam. Use when enabling random access to alignment files or fetching specific genomic regions.

Files · 1

View on GitHub

SKILL.md

1	## Version Compatibility
2
3	Reference examples tested with: pysam 0.22+, samtools 1.19+
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	- CLI: `<tool> --version` then `<tool> --help` to confirm flags
8
9	If code throws ImportError, AttributeError, or TypeError, introspect the installed
10	package and adapt the example to match the actual API rather than retrying.
11
12	# Alignment Indexing
13
14	Create indices for random access to alignment files using samtools and pysam.
15
16	"Index a BAM file" -> Create a .bai/.csi index enabling random access to genomic regions.
17	- CLI: `samtools index file.bam`
18	- Python: `pysam.index('file.bam')`
19
20	## Index Types
21
22	\| Index \| Extension \| Max contig \| Bin shift \| When required \|
23	\|-------\|-----------\|-----------\|-----------\|---------------\|
24	\| BAI \| `.bai` / `.bam.bai` \| 2^29 bp ≈ 537 Mbp \| fixed (16 kb) \| Default for human, mouse, fly, fish \|
25	\| CSI \| `.csi` / `.bam.csi` \| 2^(min_shift + depth3) \| configurable via `-m` \| Required* for any contig >537 Mbp \|
26	\| CRAI \| `.crai` / `.cram.crai` \| chunk-based \| n/a \| CRAM only \|
27	\| TBI \| `.tbi` \| 2^29-1 \| fixed \| tabix VCF/BED -- same limit as BAI \|
28
29	### Which Index for Which Genome
30
31	\| Genome \| Largest contig \| Index \|
32	\|--------\|---------------\|-------\|
33	\| GRCh38 / GRCh37 (human) \| 248 Mbp \| BAI \|
34	\| GRCm39 (mouse) \| 195 Mbp \| BAI \|
35	\| GRCz11 (zebrafish), TAIR10 (Arabidopsis) \| 78 Mbp / 30 Mbp \| BAI \|
36	\| Wheat IWGSC (Triticum aestivum) \| ~830 Mbp (chr3B) \| CSI \|
37	\| Pine, fir, axolotl, sugar pine \| multi-Gbp \| CSI with larger `-m` \|
38	\| Long-read assembly with very large contigs \| varies \| check `cut -f2 ref.fa.fai \\| sort -nr \\| head -1` \|
39
40	For polyploid plants and salamander-scale genomes, increase the bin shift:
41	```bash
42	# Default CSI matches BAI bin layout: 2^(14 + 5*3) = 2^29 bp ≈ 537 Mbp per contig
43	samtools index -c file.bam
44
45	# Larger min_shift for contigs >537 Mbp (wheat, axolotl, sugar pine)
46	samtools index -c -m 18 file.bam # 2^(18+15) = 2^33 = ~8.5 Gbp per contig
47	```
48
49	Index file precedence: htslib (and therefore the samtools CLI) tries `.csi` before `.bai` on auto-load, so when both exist the `.csi` is used -- a stale `.bai` left over after re-indexing to CSI is generally ignored. (Note: htsjdk/Java prefers `.bai`, the opposite order.) Removing the obsolete `.bai` still avoids confusion.
50
51	## samtools index
52
53	### Create BAI Index
54	```bash
55	samtools index input.bam
56	# Creates input.bam.bai
57	```
58
59	### Create CSI Index
60	```bash
61	samtools index -c input.bam
62	# Creates input.bam.csi
63	```
64
65	### Specify Output Name
66	```bash
67	samtools index input.bam output.bai
68	```
69
70	### Multi-threaded Indexing
71	```bash
72	samtools index -@ 4 input.bam
73	```
74
75	### Index CRAM
76	```bash
77	samtools index input.cram
78	# Creates input.cram.crai
79	```
80
81	## Index Requirements
82
83	Indexing requires coordinate-sorted files:
84	```bash
85	# Check sort order
86	samtools view -H input.bam \| grep "^@HD"
87	# Should show SO:coordinate
88
89	# Sort if needed, then index
90	samtools sort -o sorted.bam input.bam
91	samtools index sorted.bam
92	```
93
94	## Using Indices for Region Access
95
96	Goal: Extract reads overlapping specific genomic coordinates from an indexed BAM.
97
98	Approach: With the index present, `samtools view` or `pysam.fetch()` can jump directly to the relevant file offset instead of scanning the entire file.
99
100	### samtools view with Region
101	```bash
102	# Requires index file present
103	samtools view input.bam chr1:1000000-2000000
104	```
105
106	### Multiple Regions
107	```bash
108	samtools view input.bam chr1:1000-2000 chr2:3000-4000
109	```
110
111	### Regions from BED File
112	```bash
113	samtools view -L regions.bed input.bam
114	```
115
116	## pysam Python Alternative
117
118	### Create Index
119	```python
120	import pysam
121
122	pysam.index('input.bam')
123	# Creates input.bam.bai
124	```
125
126	### Create CSI Index
127	```python
128	# pysam.index passes through to samtools index; pass the -c flag for CSI.
129	pysam.index('-c', 'input.bam')
130	# Produces input.bam.csi.
131	```
132
133	### Fetch with Index
134	```python
135	with pysam.AlignmentFile('input.bam', 'rb') as bam:
136	# fetch() requires index
137	for read in bam.fetch('chr1', 1000000, 2000000):
138	print(read.query_name)
139	```
140
141	### Check if Indexed
142	```python
143	import pysam
144	from pathlib import Path
145
146	def is_indexed(bam_path):
147	bam_path = Path(bam_path)
148	return (bam_path.with_suffix('.bam.bai').exists() or
149	Path(str(bam_path) + '.bai').exists() or
150	bam_path.with_suffix('.bam.csi').exists())
151
152	if not is_indexed('input.bam'):
153	pysam.index('input.bam')
154	```
155
156	### Fetch Multiple Regions
157	```python
158	regions = [('chr1', 1000, 2000), ('chr1', 5000, 6000), ('chr2', 1000, 2000)]
159
160	with pysam.AlignmentFile('input.bam', 'rb') as bam:
161	for chrom, start, end in regions:
162	count = sum(1 for _ in bam.fetc

People who liked this also viewed

6 picks

Type

find-skillsHelps users discover and install agent skills when they ask questions like "how do I do X", "find a skill for X", "is there a skill that can...", or express interest in extending capabilities. This skSkillsJul 20262.6M27k
grill-meA relentless interview to sharpen a plan or design.SkillsJul 2026633k182k
vercel-react-best-practicesVercel's official collection of agent skillsSkillsJul 2026573k29k
improve-codebase-architectureScan a codebase for deepening opportunities, present them as a visual HTML report, then grill through whichever one you pick.SkillsJul 2026522k182k
setup-matt-pocock-skillsConfigure this repo for the engineering skills — set up its issue tracker, triage label vocabulary, and domain doc layout. Run once before first use of the other engineering skills.SkillsJul 2026454k182k
triageMove issues and external PRs through a state machine of triage roles — categorise, verify, grill if needed, and write agent-ready briefs.SkillsJul 2026424k182k

alignment-indexing

1.1kstars185forks

View on GitHub

How to install alignment-indexing?

gptomics/bioskills/alignment-indexing

$npx -y skills add GPTomics/bioSkills --skill alignment-indexing

Installs into the current project.

Summary

Create and use BAI/CSI indices for BAM/CRAM files using samtools and pysam. Use when enabling random access to alignment files or fetching specific genomic regions.

Files · 1

View on GitHub

SKILL.md

1	## Version Compatibility
2
3	Reference examples tested with: pysam 0.22+, samtools 1.19+
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	- CLI: `<tool> --version` then `<tool> --help` to confirm flags
8
9	If code throws ImportError, AttributeError, or TypeError, introspect the installed
10	package and adapt the example to match the actual API rather than retrying.
11
12	# Alignment Indexing
13
14	Create indices for random access to alignment files using samtools and pysam.
15
16	"Index a BAM file" -> Create a .bai/.csi index enabling random access to genomic regions.
17	- CLI: `samtools index file.bam`
18	- Python: `pysam.index('file.bam')`
19
20	## Index Types
21
22	\| Index \| Extension \| Max contig \| Bin shift \| When required \|
23	\|-------\|-----------\|-----------\|-----------\|---------------\|
24	\| BAI \| `.bai` / `.bam.bai` \| 2^29 bp ≈ 537 Mbp \| fixed (16 kb) \| Default for human, mouse, fly, fish \|
25	\| CSI \| `.csi` / `.bam.csi` \| 2^(min_shift + depth3) \| configurable via `-m` \| Required* for any contig >537 Mbp \|
26	\| CRAI \| `.crai` / `.cram.crai` \| chunk-based \| n/a \| CRAM only \|
27	\| TBI \| `.tbi` \| 2^29-1 \| fixed \| tabix VCF/BED -- same limit as BAI \|
28
29	### Which Index for Which Genome
30
31	\| Genome \| Largest contig \| Index \|
32	\|--------\|---------------\|-------\|
33	\| GRCh38 / GRCh37 (human) \| 248 Mbp \| BAI \|
34	\| GRCm39 (mouse) \| 195 Mbp \| BAI \|
35	\| GRCz11 (zebrafish), TAIR10 (Arabidopsis) \| 78 Mbp / 30 Mbp \| BAI \|
36	\| Wheat IWGSC (Triticum aestivum) \| ~830 Mbp (chr3B) \| CSI \|
37	\| Pine, fir, axolotl, sugar pine \| multi-Gbp \| CSI with larger `-m` \|
38	\| Long-read assembly with very large contigs \| varies \| check `cut -f2 ref.fa.fai \\| sort -nr \\| head -1` \|
39
40	For polyploid plants and salamander-scale genomes, increase the bin shift:
41	```bash
42	# Default CSI matches BAI bin layout: 2^(14 + 5*3) = 2^29 bp ≈ 537 Mbp per contig
43	samtools index -c file.bam
44
45	# Larger min_shift for contigs >537 Mbp (wheat, axolotl, sugar pine)
46	samtools index -c -m 18 file.bam # 2^(18+15) = 2^33 = ~8.5 Gbp per contig
47	```
48
49	Index file precedence: htslib (and therefore the samtools CLI) tries `.csi` before `.bai` on auto-load, so when both exist the `.csi` is used -- a stale `.bai` left over after re-indexing to CSI is generally ignored. (Note: htsjdk/Java prefers `.bai`, the opposite order.) Removing the obsolete `.bai` still avoids confusion.
50
51	## samtools index
52
53	### Create BAI Index
54	```bash
55	samtools index input.bam
56	# Creates input.bam.bai
57	```
58
59	### Create CSI Index
60	```bash
61	samtools index -c input.bam
62	# Creates input.bam.csi
63	```
64
65	### Specify Output Name
66	```bash
67	samtools index input.bam output.bai
68	```
69
70	### Multi-threaded Indexing
71	```bash
72	samtools index -@ 4 input.bam
73	```
74
75	### Index CRAM
76	```bash
77	samtools index input.cram
78	# Creates input.cram.crai
79	```
80
81	## Index Requirements
82
83	Indexing requires coordinate-sorted files:
84	```bash
85	# Check sort order
86	samtools view -H input.bam \| grep "^@HD"
87	# Should show SO:coordinate
88
89	# Sort if needed, then index
90	samtools sort -o sorted.bam input.bam
91	samtools index sorted.bam
92	```
93
94	## Using Indices for Region Access
95
96	Goal: Extract reads overlapping specific genomic coordinates from an indexed BAM.
97
98	Approach: With the index present, `samtools view` or `pysam.fetch()` can jump directly to the relevant file offset instead of scanning the entire file.
99
100	### samtools view with Region
101	```bash
102	# Requires index file present
103	samtools view input.bam chr1:1000000-2000000
104	```
105
106	### Multiple Regions
107	```bash
108	samtools view input.bam chr1:1000-2000 chr2:3000-4000
109	```
110
111	### Regions from BED File
112	```bash
113	samtools view -L regions.bed input.bam
114	```
115
116	## pysam Python Alternative
117
118	### Create Index
119	```python
120	import pysam
121
122	pysam.index('input.bam')
123	# Creates input.bam.bai
124	```
125
126	### Create CSI Index
127	```python
128	# pysam.index passes through to samtools index; pass the -c flag for CSI.
129	pysam.index('-c', 'input.bam')
130	# Produces input.bam.csi.
131	```
132
133	### Fetch with Index
134	```python
135	with pysam.AlignmentFile('input.bam', 'rb') as bam:
136	# fetch() requires index
137	for read in bam.fetch('chr1', 1000000, 2000000):
138	print(read.query_name)
139	```
140
141	### Check if Indexed
142	```python
143	import pysam
144	from pathlib import Path
145
146	def is_indexed(bam_path):
147	bam_path = Path(bam_path)
148	return (bam_path.with_suffix('.bam.bai').exists() or
149	Path(str(bam_path) + '.bai').exists() or
150	bam_path.with_suffix('.bam.csi').exists())
151
152	if not is_indexed('input.bam'):
153	pysam.index('input.bam')
154	```
155
156	### Fetch Multiple Regions
157	```python
158	regions = [('chr1', 1000, 2000), ('chr1', 5000, 6000), ('chr2', 1000, 2000)]
159
160	with pysam.AlignmentFile('input.bam', 'rb') as bam:
161	for chrom, start, end in regions:
162	count = sum(1 for _ in bam.fetc

People who liked this also viewed

6 picks

Type