Bio Pileup Generation

Version Compatibility

Reference examples tested with: bcftools 1.19+, pysam 0.22+, samtools 1.19+

Before using code patterns, verify installed versions match. If versions differ:

• Python: pip show then help(module.function) to check signatures
• CLI: --version then --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.

Pileup Generation

Generate pileup data for variant calling and position-level analysis.

"Generate pileup from BAM" -> Produce per-position read summaries showing depth, bases, and qualities.

• CLI: samtools mpileup -f ref.fa input.bam
• Python: bam.pileup(chrom, start, end) (pysam)

"Count alleles at a position" -> Extract per-base read support at a specific genomic coordinate.

• Python: iterate pileup_column.pileups and count bases (pysam)

What is Pileup?

Pileup shows all reads covering each position in the reference, used for:

• Variant calling (with bcftools)
• Coverage analysis
• Allele frequency calculation
• SNP/indel detection

samtools mpileup vs bcftools mpileup (Deprecation)

samtools mpileup -g/-u (BCF output for variant calling) has been deprecated since samtools 1.9 -- the genotype-likelihood code now lives in bcftools mpileup, which keeps mpileup logic versioned alongside bcftools call and avoids version-skew bugs.

| Use case | Recommended tool | |----------|------------------| | Quick allele counts at known sites | samtools mpileup or pysam pileup | | Germline variant calling (small genomes, simple cohorts) | bcftools mpileup -> bcftools call | | Germline WGS / WES production | DeepVariant or HaplotypeCaller (not mpileup) | | Somatic SNV/indel | Mutect2 / VarDict / VarScan2 (direct from BAM) | | Long-read small variants | clair3 / DeepVariant ONT (direct from BAM) | | Long-read SV | Sniffles / cuteSV (direct from BAM) | | Ultra-low-frequency (ctDNA / MRD) | fgbio consensus -> bcftools call or hot-spot Mutect2 | | Per-position allele counts (custom) | pysam pileup |

samtools mpileup (without -g) is still the standard tool for human-readable per-position read summaries.

Basic Pileup

samtools mpileup -f reference.fa input.bam > pileup.txt

Pileup Specific Region

samtools mpileup -f reference.fa -r chr1:1000000-2000000 input.bam

Regions from BED

samtools mpileup -f reference.fa -l targets.bed input.bam

Multiple BAM Files

samtools mpileup -f reference.fa sample1.bam sample2.bam sample3.bam > pileup.txt

Output Format

Text pileup format (6 columns per sample):

chr1    1000    A    15    ...............    FFFFFFFFFFF
chr1    1001    T    12    ............      FFFFFFFFFFFF

| Column | Description | |--------|-------------| | 1 | Chromosome | | 2 | Position (1-based) | | 3 | Reference base | | 4 | Read depth | | 5 | Read bases | | 6 | Base qualities |

Read Bases Encoding

| Symbol | Meaning | |--------|---------| | . | Match on forward strand | | , | Match on reverse strand | | ACGT | Mismatch (uppercase = forward) | | acgt | Mismatch (lowercase = reverse) | | ^Q | Start of read (Q = MAPQ as ASCII) | | $ | End of read | | +NNN | Insertion of N bases | | -NNN | Deletion of N bases | | * | Deleted base | | > / `') else: qpos = pileupread.queryposition base = pileupread.alignment.querysequence[qpos] if base.upper() == refbase.upper(): bases.append('.' if not pileupread.alignment.isreverse else ',') else: bases.append(base.upper() if not pileupread.alignment.is_reverse else base.lower())

print(f'{chrom}\t{pos+1}\t{ref_base}\t{depth}\t{"".join(bases)}')

pileup_text('input.bam', 'reference.fa', 'chr1', 1000000, 1000100)

## Pileup Options Summary

| Option | Description | Common pitfall |
|--------|-------------|----------------|
| `-f FILE` | Reference FASTA | Triggers BAQ ON by default |
| `-r REGION` | Restrict to region | |
| `-l FILE` | BED file of regions | |
| `-q INT` | Min mapping quality | Aligner-dependent semantics |
| `-Q INT` | Min base quality | `-Q 0` with default overlap detection has subtle behavior |
| `-d INT` | Max depth | **Default 8000 silently truncates**; bcftools mpileup default is 250 |
| `-B` | Disable BAQ | Often correct for long reads, SV, viral, consensus |
| `-A` | Count anomalous pairs | Required for amplicon |
| `-aa` | Output all positions | Required for consensus generation |
| `--ignore-overlaps` | Disable mate-overlap correction | Rarely correct |
| `--max-BQ INT` | Cap BQ (default 60) | Useful for ONT (Q values inflated) |
| `-g` (DEPRECATED) | Old BCF output | Use `bcftools mpileup` instead |

## Quick Reference

| Task | Command |
|------|---------|
| Basic pileup | `samtools mpileup -f ref.fa in.bam` |
| Quality filter | `samtools mpileup -f ref.fa -q 20 -Q 20 in.bam` |
| Region | `samtools mpileup -f ref.fa -r chr1:1-1000 in.bam` |
| To bcftools | `bcftools mpileup -f ref.fa -d 1000000 in.bam \| bcftools call -mv` |

## Common Errors

| Error | Cause | Solution |
|-------|-------|----------|
| `No FASTA reference` | Missing -f option | Add `-f reference.fa` |
| `Reference mismatch` | Wrong reference | Use same reference as alignment |
| Out of memory | High coverage region | Use `-d` to cap depth |

## Related Skills

- alignment-filtering - Filter BAM before pileup
- reference-operations - Index reference for pileup; M5 cross-check
- bam-statistics - mosdepth, depth tool selection
- variant-calling/variant-calling - Full variant calling workflows
- variant-calling/vcf-basics - VCF/BCF I/O
- variant-calling/joint-calling - Multi-sample joint calling

## Source & license

This open-source skill is cataloged on AgentStack and links to its original source — we do not rehost the code.

- **Author:** [GPTomics](https://github.com/GPTomics)
- **Source:** [GPTomics/bioSkills](https://github.com/GPTomics/bioSkills)
- **License:** MIT

Install and usage instructions live in the source repository linked above.