Pysam[1]是一个 Python 模块，它打包了高通量测序库htslib[2]的 C-API，可用于读写基因组相关文件，如 Fasta/Fastq，SAM/BAM/CRAM，VCF 等。本文以 Fasta/Fastq 文件的读写为例，介绍 Pysam 的用法，详细教程请查看官网。
Install
pip install pysam
或者
conda install pysamFasta files
对于 Fasta 文件，可以实现随机访问，前提是要先创建 faidx 索引。
import pysam
# 构建FastaFile对象，随机访问需要先创建faidx，没有的话在这里会自动创建faidx
fa = pysam.FastaFile("ex1.fa")
# Fasta文件中序列的数量，结果是一个整数
print("number of reference sequences: %d" % fa.nreferences)
# Fasta文件中序列的名称，结果是一个列表
print("names of reference sequences: " + ",".join(fa.references))
# Fasta文件中序列的长度，结果是一个列表
print("lengths of reference sequences: " + ",".join([str(i) for i in fa.lengths]))
# 这里是关键，用fetch函数随机读取序列
# 1. 提取整条序列
chr2 = fa.fetch("chr2")
print("Random fetch chr2 sequence:\n%s" % chr2)
# 2. Python风格半开区间：提取chr2位置11-20之间的碱基
# 半开区间碱基位置编号从0开始,（10, 20），其中包含位置10，不包含位置20
front1 = fa.fetch("chr2", 10, 20)
print("Python style region(chr2, 10, 20): %s" % front1)
# 3. Samtools风格闭区间：提取chr2位置11-20之间的碱基，碱基位置编号从1开始
front2 = fa.fetch(region="chr2:11-20")
print("samtools style region(chr2:11-20): %s" % front2)结果显示：
number of reference sequences: 2
names of reference sequences: chr1,chr2
lengths of reference sequences: 1575,1584
Random fetch chr2 sequence:
TTCAAATGAACTTCTGTAATTGAAAAATTCATTTAAGAAATTACAAAATATAGTTGAAAG
CTCTAACAATAGACTAAACCAAGCAGAAGAAAGAGGTTCAGAACTTGAAGACAAGTCTCT
...
Python style region(chr2, 10, 20): CTTCTGTAAT
samtools style region(chr2:11-20): CTTCTGTAATFastx files
顺序访问 Fasta/Fastq 文件。
import pysam
with pysam.FastxFile("ex1.fa") as fh:
for record in fh:
print(record.name)
print(record.sequence)
print(record.comment)
print(record.quality)
with pysam.FastxFile("ex1.fa") as fin, open("out.fa", 'w') as fout:
for record in fin:
fout.write(str(record) + "\n")结果如下：
>chr1
CACTAGTGGCTCATTGTAAATGTGTGGTTTAACTCGTCCATGGCCCAGCATTAGGGAGCT
GTGGACCCTGCAGCCTGGCTGTGGGGGCCGCAGTGGCTGAGGGGTGCAGAGCCGAGTCAC
...
>chr2
TTCAAATGAACTTCTGTAATTGAAAAATTCATTTAAGAAATTACAAAATATAGTTGAAAG
CTCTAACAATAGACTAAACCAAGCAGAAGAAAGAGGTTCAGAACTTGAAGACAAGTCTCT
...SAM/BAM/CRAM files
序列比对文件操作一般用 samtools 软件。
VCF files
变异文件操作一般用 bcftools 软件。
Tabix files
对于 TAB 键分隔的基因组位置文件（BED, SAM, GFF, VCF），可用tabix软件创建索引，然后随机访问。
写在后面
Pysam 作为一个轮子读写基因组相关文件很好用，可以替代 Biopython 的这部分功能。。
但其实现方式是通过 Cython，Python 代码中混合 C 语言代码，说实话这种代码看着非常头大，我宁愿单独用 C/C++写好相关程序，然后通过 Python 来调用。
参考资料
[1]
Pysam: https://pysam.readthedocs.io/en/latest/index.html
[2]htslib: http://www.htslib.org/
关于简说基因
生信平台
Galaxy中国（UseGalaxy.cn）致力于打造中国人的云上生物信息基础设施。大量在线工具免费使用。无需安装，用完即走。活跃的用户社区，随时交流使用心得。
生信培训
简说基因的生信培训班，荣获学员的一致好评。如果你也对生物信息学感兴趣，欢迎来跟简说基因，学真生信。
生信分析
我们能够承接所有 NGS 组学数据分析业务，包括但不限于 WGS / WES / RNA-seq 等。基因组组装、注释，以及各种重测序业务都可以与简说基因合作。
