DOCSEARCHbyAlgolia
Submit Docs

pysam Documentation

repository·master·Indexed 21 days ago

https://github.com/pysam-developers/pysam

A Python module for reading, manipulating, and writing genomic data, serving as a wrapper for the htslib C-API. It provides high-level and low-level interfaces for working with SAM, BAM, CRAM, VCF, BCF, BED, GFF, GTF, FASTA, and FASTQ formats, offering access to samtools and bcftools functionality.

Tokens
11K
Snippets
39
Records
56
Agent score
72%

What's inside pysam

  1. Overview of pysam capabilities

    master

    Pysam is a Python wrapper for the htslib C-API. It provides both a low-level wrapper and a high-level, pythonic API for genomic data processing.

    Supported File Formats:

    • SAM/BAM
    • VCF/BCF
    • BED
    • GFF
    • GTF
    • FASTA
    • FASTQ

    Key Features:

    • Access to samtools and bcftools command-line functionality.
    • Support for file compression.
    • Support for random access through indexing.

  2. How to use multiple iterators on the same AlignmentFile

    master

    By default, an AlignmentFile object maintains a single file position. If you attempt to create multiple iterators (e.g., by calling .fetch() multiple times) on the same file object, the iterators will interfere with each other because they share that single file position.

    To use multiple iterators simultaneously, use the multiple_iterators=True argument in the .fetch() method. This causes pysam to return an iterator backed by a newly opened file handle, preventing interference. Note that this incurs a performance penalty due to re-opening the file.

    Incorrect usage (interfering iterators):

    samfile = pysam.AlignmentFile("pysam_ex1.bam", "rb")

iter1 = samfile.fetch("chr1")
print(next(iter1).reference_id)
iter2 = samfile.fetch("chr2")
print(next(iter2).reference_id)
print(next(iter1).reference_id) # This will fail or return unexpected results

    Correct usage:

    samfile = pysam.AlignmentFile("pysam_ex1.bam", "rb")

iter1 = samfile.fetch("chr1", multiple_iterators=True)
print(next(iter1).reference_id)
iter2 = samfile.fetch("chr2")
print(next(iter2).reference_id)
print(next(iter1).reference_id) # This works correctly
    samfile = pysam.AlignmentFile("pysam_ex1.bam", "rb")

iter1 = samfile.fetch("chr1", multiple_iterators=True)
print(next(iter1).reference_id)
iter2 = samfile.fetch("chr2")
print(next(iter2).reference_id)
print(next(iter1).reference_id)

  3. How to handle pysam coordinate conventions

    master

    Pysam uses 0-based coordinates and half-open notation for ranges (consistent with Python). While pysam automatically translates coordinates to/from file formats (like the 1-based SAM format), there is one important exception:

    The region string passed to AlignmentFile.fetch() and AlignmentFile.pileup() follows the samtools command line convention (which is 1-based). Any coordinates passed directly to samtools utilities, such as pysam.mpileup, also follow this convention.

  4. Understand genomic region coordinate systems

    master

    When working with genomic regions in pysam, it is critical to distinguish between internal coordinate systems and external string notations:

    Internal pysam coordinates

    pysam uses 0-based half-open intervals.

    • The first base of a reference sequence is index 0.
    • The start position is inclusive (part of the interval).
    • The end position is exclusive (not part of the interval).

    samtools-compatible string notation

    When a region is provided as a single string (e.g., 'chr1:15001-20000'), the coordinates follow a 1-based closed interval convention. In this format, both the start and end positions are considered part of the interval.

  5. Understand CIGAR string representation in the Python API

    master

    In pysam, a CIGAR (Compact Idiosyncratic Gapped Alignment Report) string is represented as a list of tuples in the format (operation, length). Each tuple describes a specific alignment operation and the number of bases it spans.

    Example: [(0, 3), (1, 5), (0, 2)] represents an alignment with 3 matches, followed by 5 insertions, followed by 2 matches.

    [(0, 3), (1, 5), (0, 2)]

  6. Python language requirements and constraints

    master

    Pysam requires Python 3.9 as a minimum.

    While modern features like f-strings, the walrus operator (:=), and string methods like str.removeprefix are available, you must not use features introduced in Python 3.10 or later in the source code or infrastructure scripts, such as:

    • Union types using the pipe operator (e.g., type | None)
    • Grouping parentheses in with statements

  8. How to fetch from a BAM file without an index

    master

    The .fetch() method normally requires a BAM/CRAM index to perform random access. If you need to iterate over a file that lacks an index, use the until_eof=True flag to perform a linear scan.

    bf = pysam.AlignmentFile(fname, "rb")
for r in bf.fetch(until_eof=True):
    print(r)
    bf = pysam.AlignmentFile(fname, "rb")
for r in bf.fetch(until_eof=True):
    print(r)

  11. Use samtools and bcftools commands within Python

    master

    Pysam provides direct Python wrappers for samtools and bcftools commands. These are available as function calls in the pysam.samtools and pysam.bcftools namespaces, and samtools commands are also imported into the main pysam namespace. Command-line options are passed as arguments to the functions.

    Handling Output

    • Capture stdout: By default, stdout is captured and returned as the function's return value. Use catch_stdout=True (default).
    • Save to file: Use the save_stdout keyword argument with a filename, or pass "-o", "filename" as an argument and set catch_stdout=False to prevent pysam from overriding the redirection.
    • Discard stdout: Set catch_stdout=False. This is useful in environments like IPython notebooks to prevent large outputs from cluttering the interface. The function will return None.

    Error Handling and Debugging

    • Errors: Argument errors or command failures raise a pysam.SamtoolsError.
    • Stderr: Standard error messages are always captured and can be retrieved using the .get_messages() method on the command object.
    • Usage: To see the help/usage information for a specific command, call .usage().
    import pysam.samtools
# Equivalent to: samtools sort -o output.bam ex1.bam
pysam.samtools.sort("-o", "output.bam", "ex1.bam", catch_stdout=False)

import pysam.bcftools
# Equivalent to: bcftools index --csi ex2.vcf.gz
pysam.bcftools.index("--csi", "ex2.vcf.gz")

# Using samtools via main namespace
pysam.sort("-m", "1000000", "-o", "output.bam", "ex1.bam", catch_stdout=False)

# Getting usage information
print(pysam.sort.usage())

# Retrieving stderr messages
pysam.sort.get_messages()

  12. Read and iterate over SAM/BAM/CRAM files

    master

    Use pysam.AlignmentFile to open mapped short read sequence data. You can iterate over reads in a specific genomic region using the .fetch() method. Each iteration returns an AlignedSegment object representing a single read.

    Note: Pysam uses 0-based coordinates (Python convention), whereas SAM text files use 1-based coordinates.

    import pysam

# Open a BAM file for reading
samfile = pysam.AlignmentFile("ex1.bam", "rb")

# Fetch and iterate over reads in a specific region
for read in samfile.fetch('chr1', 100, 120):
    print(read)

samfile.close()