简介
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)