pysam

Pysam - Python Genomic Data Processing Tool

Skills Overview

Pysam is a Python module for reading, manipulating, and writing genomic datasets. It provides read/write interfaces for SAM/BAM/CRAM alignment files, VCF/BCF variant files, and FASTA/FASTQ sequence files. It supports querying genomic regions, stacked coverage analysis, and executing samtools/bcftools commands.

Use Cases

1. NGS Sequencing Data Analysis

Suitable for processing high-throughput sequencing data, including reading alignment result files (BAM/CRAM), analyzing sequencing coverage, extracting reads from specific gene regions, and performing data quality control. It supports batch processing of large-scale sequencing datasets, enabling fast localization and analysis of particular genomic regions.

2. Detection and Annotation of Genomic Variants

Suitable for analyzing genetic variant data (VCF/BCF). It supports querying variant loci, extracting genotypes, quality filtering, and functional annotation. It can be used for variant analysis on single or multiple samples, and can verify and filter variants in combination with coverage information.

3. Sequence Extraction and Genomic Data Processing

Suitable for extracting gene sequences from a reference genome, validating variant sites, and processing raw FASTQ sequencing data. It supports random access to genomic regions via index files, and can be used to build bioinformatics analysis workflows, validate sequences, and compute coverage statistics.

Core Features

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

Provides an AlignmentFile class to handle sequencing alignments. It supports opening and reading BAM/SAM/CRAM files, retrieving reads by genomic region, filtering reads based on mapping quality and tags, computing coverage statistics, performing stacked analysis, and accessing read sequences and quality scores. It also supports writing modified alignment files, making it suitable for alignment result analysis, coverage calculation, and data quality control.

2. Variant File Operations (VCF/BCF)

Provides a VariantFile class to handle genetic variant data. It supports reading and writing VCF/BCF files, querying variants in specific regions, accessing variant information (position, alleles, quality), extracting genotype data for samples, filtering variants by quality or allele frequency, adding annotation information, and selecting specific samples or regions. It is suitable for variant analysis, filtering, annotation, and population genetics studies.

3. Sequence File Operations (FASTA/FASTQ)

Provides FastaFile for random access to reference sequences and FastxFile for reading raw sequencing data. It supports querying reference sequences by genomic coordinates, extracting sequences for genes or regions of interest, reading FASTQ files with quality scores, validating reference alleles for variants, computing sequence statistics, filtering reads by quality or length, and converting between FASTA and FASTQ formats. It is suitable for extracting gene sequences, validating variants, or processing raw reads.

Frequently Asked Questions

How does Pysam read a specific region in a BAM file?

After opening the BAM file with AlignmentFile, specify the chromosome and coordinate range using the fetch() method. Note that Pysam uses a 0-based half-open coordinate system, but the string arguments for fetch() follow samtools’ 1-based convention. Random access requires creating a .bai index file in advance, which can be generated using pysam.index().

How can I use Pysam to compute NGS data coverage?

Pysam provides two main methods: using pileup() for per-position stacked analysis, or using count() to count reads in a specific region. pileup() is suitable for scenarios that require detailed coverage information, as it can provide coverage depth, base distributions, and quality information for each position. For large-scale datasets, it is recommended to use index files and process regions in parallel to improve performance.

What are the basic ways to process VCF variant data with Pysam?

Use the VariantFile class to open a VCF or BCF file. You can iterate through all variants or use fetch() to query variants in a specific region. Access chromosome, position, reference and alternative alleles, quality scores, and other information via the VariantRecord object. For multi-sample VCFs, you can extract genotype data for each sample. Pysam also supports filtering variants by conditions such as quality, depth, and allele frequency, and allows adding custom annotation information or creating a new VCF file.