简介

Pysam 是一个用于读取、操作和编写基因组数据集的 python 模块。它是 htslib C-API 的轻量级包装器，提供读取和写入 SAM / BAM / VCF / BCF / BED / GFF / GTF / FASTA / FASTQ 文件以及访问 samtools 和 bcftools 包的命令行功能的工具.

主要功能模块

读取 sam/bam/cram

import pysam
samfile = pysam.AlignmentFile("ex1.bam", "rb")
for read in samfile.fetch('chr1', 100, 120):
    print(read)

samfile.close()

使用 samtools 工具接口

命令行中 samtools 的功能模块都可以直接使用：
pysam.sort("-o", "output.bam", "ex1.bam")
等同于命令行的命令：
samtools sort -o output.bam ex1.bam

读取 fasta 文件

使用 fasta 文件获取 SRY 基因（Y:2654895-2655723）序列:

fa = pysam.FastaFile('b37.fasta')
SRY_seq = fa.fetch(region='Y:2654895-2655723')

读取 bcf/vcf 文件

过滤掉一些无用的染色体变异，另存一份干净的变异文件：

from pysam import VariantFile
bcf_in  = VariantFile("test_in.vcf", "r")
bcf_out = VariantFile("test_out.vcf", "w", header=bcf_in.header)
chroms = [str(i) for i in range(1,23)] + ['X', 'Y']
for rec in bcf_in.fetch():
    if rec.chrom not in chroms:
        continue
    bcf_out.write(rec)