Generate consensus sequences and manage reference files using samtools. Use when creating consensus from alignments, indexing references, or creating sequence dictionaries.
Standardmäßig ist der Prompt ausgewählt, der zuerst die Quelle prüft. Sie können zu einem direkten Befehl wechseln oder eine lokale Kopie herunterladen.
Quelldateien prüfen
Lesen Sie SKILL.md und alle von SkillsMP angezeigten Begleitdateien, bevor Sie sich für eine Installation entscheiden.
Mit Codex oder Claude installieren Kopieren Sie diesen Prompt, fügen Sie ihn in Codex, Claude oder einen anderen Assistant ein und lassen Sie die Skill-Seite prüfen und installieren.
Ein direkter Befehl überspringt den Prüf-Prompt. Prüfen Sie die Quelle, bevor Sie ihn ausführen.
Generate consensus sequences and manage reference files using samtools. Use when creating consensus from alignments, indexing references, or creating sequence dictionaries.
Before using code patterns, verify installed versions match. If versions differ:
Python: pip show <package> then help(module.function) to check signatures
CLI: <tool> --version then <tool> --help to confirm flags
If code throws ImportError, AttributeError, or TypeError, introspect the installed
package and adapt the example to match the actual API rather than retrying.
Reference Operations
Generate consensus sequences and manage reference files using samtools.
"Prepare a reference genome" -> Index the FASTA and create a sequence dictionary for downstream tools.
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.
GRCh38 Is Not One Reference
Reference flavor
ALT
Decoy
EBV
HLA
Use case
GRCh38 no-alt
no
no
no
no
Conservative analyses
GRCh38 + decoy + EBV (1000G analysis set)
no
yes
yes
no
Cohort projects
GRCh38 ALT + decoy + EBV + HLA (Broad / hs38DH)
yes
yes
yes
yes
GATK Best Practices
T2T-CHM13 v2.0
n/a
n/a
n/a
n/a
Distinct coordinates -- NOT interchangeable
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.
Contig Naming: The Silent Killer
Convention
Source
chr1
mitochondrion
UCSC (hg19, hg38)
UCSC Genome Browser
chr1
chrM
Ensembl (GRCh37, GRCh38)
Ensembl, ENA
1
MT
NCBI RefSeq (recent)
NCBI
chr1
chrM
1000G analysis sets
1000G GRCh38 analysis set
chr1
chrM
A BAM with @SQ SN:chr1 cannot be analyzed against a 1-named reference (and vice versa). Detect:
samtools view -H sample.bam | grep '^@SQ' | head -3
samtools dict ref.fa | head -3
Convert: bcftools annotate --rename-chrs for VCF; for BAM there is no clean conversion -- re-align.
# Minimum depth to call base
samtools consensus -d 5 input.bam -o consensus.fa
# Call all positions (including low coverage)
samtools consensus -a input.bam -o consensus.fa
IUPAC Ambiguity for Heterozygotes
# Emit IUPAC codes (R, Y, S, W, K, M, B, D, H, V, N) for heterozygous columns# --ambig is REQUIRED -- without it, output is restricted to A,C,G,T,N,*
samtools consensus --ambig --het-fract 0.2 --call-fract 0.5 input.bam -o consensus.fa
--het-fract controls the fraction of the second-most-common base relative to the most common required to call a heterozygote (verify the default for the installed version with samtools consensus --help; the manpage documents none). Without --ambig, columns where the second base passes --het-fract resolve to N rather than the IUPAC code. --show-ins / --show-del control insertion / deletion display, not ambiguity.
Platform-Aware Consensus
# Default: Bayesian algorithm (no --config needed)
samtools consensus -f fasta input.bam -o consensus.fa
# Platform-specific profiles (samtools 1.17+; verify via samtools consensus --help for installed version)
samtools consensus --config hifi input.bam -o consensus.fa # PacBio HiFi
samtools consensus --config r10.4_sup input.bam -o consensus.fa # ONT R10.4+ (r10.4_dup for duplex)
samtools consensus --config ultima input.bam -o consensus.fa # Ultima Genomics
samtools consensus --config hiseq input.bam -o consensus.fa # Illumina# Report ref base where consensus unavailable (low coverage; -T added in samtools 1.22)
samtools consensus -T ref.fa input.bam -o consensus.fa
samtools consensus vs bcftools consensus
Different operations -- conflating them produces nonsense:
Tool
Input
Output
Use case
samtools consensus
BAM
Consensus FASTA derived from reads (Bayesian)
Viral, de novo / amplicon, low-coverage species
bcftools consensus
reference + VCF
Reference with VCF variants applied
Apply called variants (haplotype reconstruction, custom ref for re-mapping)
For bacterial / phage assembly polishing, prefer Pilon (short-read) or medaka (ONT); samtools consensus is not iterative.
pysam Python Alternative
Fetch from Indexed FASTA
import pysam
with pysam.FastaFile('reference.fa') as ref:
seq = ref.fetch('chr1', 999, 2000) # 0-basedprint(seq)
Get Reference Lengths
with pysam.FastaFile('reference.fa') as ref:
for name in ref.references:
length = ref.get_reference_length(name)
print(f'{name}: {length:,} bp')
Fetch All Chromosomes
with pysam.FastaFile('reference.fa') as ref:
for chrom in ref.references:
seq = ref.fetch(chrom)
print(f'>{chrom}')
print(seq[:100] + '...')
Generate Simple Consensus
import pysam
from collections import Counter
defconsensus_at_position(bam, chrom, pos):
bases = Counter()
for pileup in bam.pileup(chrom, pos, pos + 1, truncate=True):
if pileup.pos == pos:
for read in pileup.pileups:
ifnot read.is_del andnot read.is_refskip:
bases[read.alignment.query_sequence[read.query_position]] += 1if bases:
return bases.most_common(1)[0][0]
return'N'with pysam.AlignmentFile('input.bam', 'rb') as bam:
consensus = consensus_at_position(bam, 'chr1', 1000000)
print(f'Consensus at chr1:1000000 = {consensus}')
Build Consensus Sequence (Pedagogical Only)
The Python majority-vote consensus below is illustrative, NOT production. samtools consensus is Bayesian, quality-aware, and platform-aware; majority vote ignores base qualities and produces wrong calls on low-coverage / low-quality regions. Use for teaching pileup iteration mechanics; use samtools consensus for any real consensus.
import pysam
from collections import Counter
defbuild_consensus(bam_path, chrom, start, end, min_depth=3):
consensus = []
with pysam.AlignmentFile(bam_path, 'rb') as bam:
for pileup in bam.pileup(chrom, start, end, truncate=True):
bases = Counter()
for read in pileup.pileups:
ifnot read.is_del andnot read.is_refskip:
base = read.alignment.query_sequence[read.query_position]
bases[base] += 1ifsum(bases.values()) >= min_depth:
consensus.append(bases.most_common(1)[0][0])
else:
consensus.append('N')
return''.join(consensus)
Create Dictionary Header
import pysam
defcreate_dict_header(fasta_path):
header = {'HD': {'VN': '1.6', 'SO': 'unsorted'}, 'SQ': []}
with pysam.FastaFile(fasta_path) as ref:
for name in ref.references:
length = ref.get_reference_length(name)
header['SQ'].append({'SN': name, 'LN': length})
return header
header = create_dict_header('reference.fa')
for sq in header['SQ'][:5]:
print(f'{sq["SN"]}: {sq["LN"]:,} bp')
Reference Preparation Workflow
Goal: Set up a reference genome with all indices needed by common analysis tools.
Approach: Create FASTA index (.fai), sequence dictionary (.dict), and aligner-specific indices in sequence.
Prepare Reference for Analysis
# 1. Index FASTA for samtools/pysam
samtools faidx reference.fa
# 2. Create sequence dictionary for GATK/Picard
samtools dict reference.fa -o reference.dict
# 3. Pre-populate CRAM REF_CACHE (for offline HPC nodes)
seq_cache_populate.pl -root $REF_CACHE_DIR reference.fa
For aligner-specific indices (BWA, Bowtie2, STAR, minimap2, Salmon), see read-alignment.
Check Reference Setup
# Verify FAI existsls -la reference.fa.fai
# Verify dict existshead reference.dict
# Test fetch
samtools faidx reference.fa chr1:1-100
Common Operations
Extract Chromosome
samtools faidx reference.fa chr1 > chr1.fa
samtools faidx chr1.fa # Index the subset
