Pysam

Dependencies

• Python: 3.10+. Repository baseline for current packaged skills.
• Third-party packages: not explicitly version-pinned in this skill package. Add pinned versions if this skill needs stricter environment control.

Example Usage

Skill directory: 20260316/scientific-skills/Data Analytics/pysam
No packaged executable script was detected.
Use the documented workflow in SKILL.md together with the references/assets in this folder.

Example run plan:

1. Read the skill instructions and collect the required inputs.
2. Follow the documented workflow exactly.
3. Use packaged references/assets from this folder when the task needs templates or rules.
4. Return a structured result tied to the requested deliverable.

Implementation Details

See ## Overview above for related details.

• Execution model: validate the request, choose the packaged workflow, and produce a bounded deliverable.
• Input controls: confirm the source files, scope limits, output format, and acceptance criteria before running any script.
• Primary implementation surface: instruction-only workflow in SKILL.md.
• Reference guidance: references/ contains supporting rules, prompts, or checklists.
• Parameters to clarify first: input path, output path, scope filters, thresholds, and any domain-specific constraints.
• Output discipline: keep results reproducible, identify assumptions explicitly, and avoid undocumented side effects.

Overview

Pysam is a Python module for reading, manipulating, and writing genomic datasets. It provides a Python-style interface to htslib, supporting reading/writing SAM/BAM/CRAM alignment files, VCF/BCF variant files, and FASTA/FASTQ sequences. It can also query tabix-indexed files, perform pileup analysis for coverage calculation, and execute samtools/bcftools commands.

When to Use This Skill

Use this skill in the following scenarios:

• Processing sequencing alignment files (BAM/CRAM)
• Analyzing genetic variants (VCF/BCF)
• Extracting reference sequences or gene regions
• Processing raw sequencing data (FASTQ)
• Calculating coverage or sequencing depth
• Implementing bioinformatics analysis pipelines
• Quality control of sequencing data
• Variant calling and annotation workflows

Quick Start

Installation

uv pip install pysam

Basic Examples

Reading alignment files:

import pysam

# Open BAM file and fetch reads in specified region
samfile = pysam.AlignmentFile("example.bam", "rb")
for read in samfile.fetch("chr1", 1000, 2000):
    print(f"{read.query_name}: {read.reference_start}")
samfile.close()

Reading variant files:

# Open VCF file and iterate through variant sites
vcf = pysam.VariantFile("variants.vcf")
for variant in vcf:
    print(f"{variant.chrom}:{variant.pos} {variant.ref}>{variant.alts}")
vcf.close()

Querying reference sequences:

# Open FASTA and extract sequence
fasta = pysam.FastaFile("reference.fasta")
sequence = fasta.fetch("chr1", 1000, 2000)
print(sequence)
fasta.close()

Core Capabilities

1. Alignment File Operations (SAM/BAM/CRAM)

Use the AlignmentFile class to work with aligned sequencing reads. This is suitable for analyzing alignment results, calculating coverage, extracting reads, or quality control.

Common operations:

• Open and read BAM/SAM/CRAM files
• Fetch reads from specific genomic regions
• Filter reads by alignment quality (mapping quality), flag, or other criteria
• Write filtered or modified alignment data
• Calculate coverage statistics
• Perform pileup analysis (per-base coverage)
• Access read sequences, quality values, and alignment information

Reference: For detailed documentation, see references/alignment_files.md:

• Opening and reading alignment files
• AlignedSegment attributes and methods
• Region-based extraction using fetch()
• Pileup analysis for coverage analysis
• Writing and creating BAM files
• Coordinate systems and indexing
• Performance optimization tips

2. Variant File Operations (VCF/BCF)

Use the VariantFile class to work with genetic variants from variant calling pipelines. This is suitable for variant analysis, filtering, annotation, or population genetics studies.

Common operations:

• Read/write VCF/BCF files
• Query variants in specific regions
• Access variant information (position, alleles, quality scores)
• Extract genotype data for samples
• Filter variants by quality, allele frequency, or other criteria
• Annotate variants with additional information
• Subset samples or regions

Reference: For detailed documentation, see references/variant_files.md:

• Opening and reading variant files
• VariantRecord attributes and methods
• Accessing INFO and FORMAT fields
• Handling genotypes and samples
• Creating and writing VCF files
• Filtering and extracting variant subsets
• Multi-sample VCF operations

3. Sequence File Operations (FASTA/FASTQ)

Use FastaFile for random access to reference sequences, and FastxFile for reading raw sequencing data. This is suitable for extracting gene sequences, validating variants against reference, or processing raw reads.

Common operations:

• Query reference sequences by genomic coordinates
• Extract sequences for genes or regions of interest
• Read FASTQ files with quality values
• Validate reference alleles for variants
• Calculate sequence statistics
• Filter reads by quality or length
• Convert between FASTA and FASTQ formats

Reference: For detailed documentation, see references/sequence_files.md:

• FASTA file access and indexing
• Extracting sequences by region
• Handling reverse complement sequences for genes
• Sequential reading of FASTQ files
• Quality score conversion and filtering
• Processing tabix-indexed files (BED, GTF, GFF)
• Common sequence processing patterns

4. Integrated Bioinformatics Workflows

Pysam excels at integrating multiple file types for comprehensive genomic analysis. Common workflows combine alignment files, variant files, and reference sequences.

Common workflows:

• Calculate coverage statistics for specific regions
• Verify variants using aligned reads
• Annotate variants with coverage information
• Extract sequences around variant positions
• Filter alignments or variants based on multiple criteria
• Generate coverage tracks for visualization
• Quality control across multiple data types

Reference: For detailed examples, see references/common_workflows.md:

• Quality control workflows (BAM statistics, reference consistency)
• Coverage analysis (per-base coverage, low coverage detection)
• Variant analysis (annotation, read support-based filtering)
• Sequence extraction (variant context, gene sequences)
• Read filtering and subset extraction
• Integration patterns (BAM+VCF, VCF+BED, etc.)
• Performance optimization for complex workflows

Key Concepts

Coordinate Systems

Critical: Pysam uses 0-based, left-closed right-open coordinates (Python convention):

• Start positions begin at 0 (first base is at position 0)
• End positions are exclusive (not included in range)
• Region 1000-2000 contains bases 1000-1999 (1000 bases total)

Exception: Region strings in fetch() follow samtools convention (1-based):

samfile.fetch("chr1", 999, 2000)      # 0-based: positions 999-1999
samfile.fetch("chr1:1000-2000")       # 1-based string: positions 1000-2000

VCF files: Use 1-based coordinates in file format, but VariantRecord.start is 0-based.

Indexing Requirements

Random access to specific genomic regions requires index files:

• BAM files: Require .bai index (created with pysam.index())
• CRAM files: Require .crai index
• FASTA files: Require .fai index (created with pysam.faidx())
• VCF.gz files: Require .tbi tabix index (created with pysam.tabix_index())
• BCF files: Require .csi index

If no index is available, use fetch(until_eof=True) for sequential reading.

File Modes

Specify format when opening files:

• "rb" - Read BAM (binary)
• "r" - Read SAM (text)
• "rc" - Read CRAM
• "wb" - Write BAM
• "w" - Write SAM
• "wc" - Write CRAM

Performance Considerations

1. Always use index files for random access operations
2. Use pileup() for column-wise analysis instead of repeated fetch operations
3. Use count() for counting instead of manual iteration counting
4. Process regions in parallel when analyzing independent genomic regions
5. Explicitly close files to release resources
6. Use until_eof=True for sequential processing without indexing
7. Avoid multiple iterators unless necessary (use multiple_iterators=True if needed)

Common Pitfalls

1. Coordinate confusion: Remember the difference between 0-based and 1-based systems in different contexts.
2. Missing indexes: Many operations require index files - create them first.
3. Partial overlap: fetch() returns reads overlapping region boundaries, not just reads fully contained within.
4. Iterator scope: Keep pileup iterator references alive to avoid "PileupProxy accessed after iterator finished" error.
5. Quality value editing: Cannot modify query_qualities in-place after modifying query_sequence - must create a copy first.
6. Streaming limitations: Only stdin/stdout streams are supported, not arbitrary Python file objects.
7. Thread safety: While GIL is released during I/O, full thread safety is not completely verified.

Command-Line Tools

Pysam provides access to samtools and bcftools commands:

# Sort BAM file
pysam.samtools.sort("-o", "sorted.bam", "input.bam")

# Index BAM
pysam.samtools.index("sorted.bam")

# View specific region
pysam.samtools.view("-b", "-o", "region.bam", "input.bam", "chr1:1000-2000")

# BCF tools
pysam.bcftools.view("-O", "z", "-o", "output.vcf.gz", "input.vcf")

Error handling: