reference-operations
starsforks
How to install reference-operations?
gptomics/bioskills/reference-operations
$npx -y skills add GPTomics/bioSkills --skill reference-operationsGenerate consensus sequences and manage reference files using samtools. Use when creating consensus from alignments, indexing references, or creating sequence dictionaries.
SKILL.md
| 1 | ## Version Compatibility |
| 2 | |
| 3 | Reference examples tested with: GATK 4.5+, bcftools 1.19+, 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 | # Reference Operations |
| 13 | |
| 14 | Generate consensus sequences and manage reference files using samtools. |
| 15 | |
| 16 | **"Prepare a reference genome"** -> Index the FASTA and create a sequence dictionary for downstream tools. |
| 17 | - CLI: `samtools faidx ref.fa` + `samtools dict ref.fa -o ref.dict` |
| 18 | - Python: `pysam.FastaFile('ref.fa')` (auto-uses .fai index) |
| 19 | |
| 20 | **"Build a consensus from BAM"** -> Derive the most-supported base at each position from aligned reads. |
| 21 | - CLI: `samtools consensus input.bam -o consensus.fa` |
| 22 | - Python: iterate pileup columns and take majority base (pysam) |
| 23 | |
| 24 | ## samtools faidx - Index Reference FASTA |
| 25 | |
| 26 | Create index for random access to reference sequences. |
| 27 | |
| 28 | ### Create Index |
| 29 | ```bash |
| 30 | samtools faidx reference.fa |
| 31 | # Creates reference.fa.fai |
| 32 | ``` |
| 33 | |
| 34 | ### Fetch Region from Reference |
| 35 | ```bash |
| 36 | samtools faidx reference.fa chr1:1000-2000 |
| 37 | ``` |
| 38 | |
| 39 | ### Fetch Multiple Regions |
| 40 | ```bash |
| 41 | samtools faidx reference.fa chr1:1000-2000 chr2:3000-4000 |
| 42 | ``` |
| 43 | |
| 44 | ### Fetch Entire Chromosome |
| 45 | ```bash |
| 46 | samtools faidx reference.fa chr1 |
| 47 | ``` |
| 48 | |
| 49 | ### Output to File |
| 50 | ```bash |
| 51 | samtools faidx reference.fa chr1:1000-2000 > region.fa |
| 52 | ``` |
| 53 | |
| 54 | ### Reverse Complement |
| 55 | ```bash |
| 56 | samtools faidx -i reference.fa chr1:1000-2000 |
| 57 | ``` |
| 58 | |
| 59 | ### FAI File Format |
| 60 | ``` |
| 61 | chr1 248956422 6 60 61 |
| 62 | chr2 242193529 253105708 60 61 |
| 63 | ``` |
| 64 | Columns: name, length, offset, line bases, line width |
| 65 | |
| 66 | ## samtools dict - Create Sequence Dictionary |
| 67 | |
| 68 | Create SAM header dictionary for reference (used by GATK, Picard). |
| 69 | |
| 70 | ### Create Dictionary |
| 71 | ```bash |
| 72 | samtools dict reference.fa -o reference.dict |
| 73 | ``` |
| 74 | |
| 75 | ### With Assembly Info |
| 76 | ```bash |
| 77 | samtools dict -a GRCh38 -s "Homo sapiens" reference.fa -o reference.dict |
| 78 | ``` |
| 79 | |
| 80 | ### Dictionary Format |
| 81 | ``` |
| 82 | @HD VN:1.0 SO:unsorted |
| 83 | @SQ SN:chr1 LN:248956422 M5:6aef897c3d6ff0c78aff06ac189178dd UR:file:reference.fa |
| 84 | @SQ SN:chr2 LN:242193529 M5:f98db672eb0993dcfdabafe2a882905c UR:file:reference.fa |
| 85 | ``` |
| 86 | |
| 87 | The `M5:` (MD5) tag is the only definitive reference-identity check -- two references named "GRCh38" with different decoy/alt content have different M5s. CRAM enforces M5 match on read-back. See alignment-validation for BAM-vs-reference M5 cross-check. |
| 88 | |
| 89 | ### GRCh38 Is Not One Reference |
| 90 | |
| 91 | | Reference flavor | ALT | Decoy | EBV | HLA | Use case | |
| 92 | |------------------|-----|-------|-----|-----|----------| |
| 93 | | GRCh38 no-alt | no | no | no | no | Conservative analyses | |
| 94 | | GRCh38 + decoy + EBV (1000G analysis set) | no | yes | yes | no | Cohort projects | |
| 95 | | GRCh38 ALT + decoy + EBV + HLA (Broad / hs38DH) | yes | yes | yes | yes | GATK Best Practices | |
| 96 | | T2T-CHM13 v2.0 | n/a | n/a | n/a | n/a | Distinct coordinates -- NOT interchangeable | |
| 97 | |
| 98 | Mixing no-alt and ALT-aware BAMs in one cohort produces inconsistent multi-mapping behavior at HLA, KIR, and segmental-duplication regions. Standardize before joint calling. |
| 99 | |
| 100 | ### Contig Naming: The Silent Killer |
| 101 | |
| 102 | | Convention | Source | chr1 | mitochondrion | |
| 103 | |-----------|--------|------|---------------| |
| 104 | | UCSC (hg19, hg38) | UCSC Genome Browser | chr1 | chrM | |
| 105 | | Ensembl (GRCh37, GRCh38) | Ensembl, ENA | 1 | MT | |
| 106 | | NCBI RefSeq (recent) | NCBI | chr1 | chrM | |
| 107 | | 1000G analysis sets | 1000G GRCh38 analysis set | chr1 | chrM | |
| 108 | |
| 109 | A BAM with `@SQ SN:chr1` cannot be analyzed against a `1`-named reference (and vice versa). Detect: |
| 110 | ```bash |
| 111 | samtools view -H sample.bam | grep '^@SQ' | head -3 |
| 112 | samtools dict ref.fa | head -3 |
| 113 | ``` |
| 114 | |
| 115 | Convert: `bcftools annotate --rename-chrs` for VCF; for BAM there is no clean conversion -- re-align. |
| 116 | |
| 117 | ## samtools consensus - Generate Consensus |
| 118 | |
| 119 | Create consensus sequence from alignments. |
| 120 | |
| 121 | ### Basic Consensus |
| 122 | ```bash |
| 123 | samtools consensus input.bam -o consensus.fa |
| 124 | ``` |
| 125 | |
| 126 | ### From Specific Region |
| 127 | ```bash |
| 128 | samtools consensus -r chr1:1000-2000 input.bam -o region_consensus.fa |
| 129 | ``` |
| 130 | |
| 131 | ### Output Formats |
| 132 | ```bash |
| 133 | # FASTA (default) |
| 134 | samtools consensus -f fasta input.bam -o consensus.fa |
| 135 | |
| 136 | # FASTQ (includes quality) |
| 137 | samtools consensus -f fastq input.bam -o consensus.fq |
| 138 | ``` |
| 139 | |
| 140 | ### Quality Options |
| 141 | ```bash |
| 142 | # Minimum depth to call base |
| 143 | samtools consensus -d 5 input.bam -o consensus.fa |
| 144 | |
| 145 | # Call all positions (including low coverage) |
| 146 | samtools consensus -a input.bam -o consensus.fa |
| 147 | ``` |
| 148 | |
| 149 | ### IUPAC Ambiguity for Heterozygotes |
| 150 | ```bash |
| 151 | # Emit IUPAC codes (R, Y, S, W, K, M, B, D, H, V, N) for heterozygous columns |
| 152 | # --ambig is REQUIRED -- |
People who liked this also viewed
6 picksType
- 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 2026
- grill-meA relentless interview to sharpen a plan or design.SkillsJul 2026
- vercel-react-best-practicesVercel's official collection of agent skillsSkillsJul 2026
- improve-codebase-architectureScan a codebase for deepening opportunities, present them as a visual HTML report, then grill through whichever one you pick.SkillsJul 2026
- 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 2026
- triageMove issues and external PRs through a state machine of triage roles — categorise, verify, grill if needed, and write agent-ready briefs.SkillsJul 2026
