Pysam是一个Python模块，它打包了高通量测序库htslib的C-API，可用于读写基因组相关文件，如Fasta/Fastq，SAM/BAM/CRAM，VCF等。本文以Fasta/Fastq文件的读写为例，介绍Pysam的用法，详细教程请查看官网。

Install

1
2
3

pip install pysam
或者
conda install pysam

Fasta files

对于Fasta文件，可以实现随机访问，前提是要先创建faidx索引。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27

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)

结果显示：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10

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): CTTCTGTAAT

Fastx files

顺序访问Fasta/Fastq文件。

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12

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")

结果如下：

1
2
3
4
5
6
7
8

>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来调用。