前言

基因组序列的提取，有不少强大的工具像samtools,bedtools，之前也提到 pybedtools 提取序列。 不过pybedtools是对bedtools提供一个Python接口，除了安装pybedtools外， 还需要bedtools。最近想要一款纯Python实现，不依赖于其他的软件，能随机索引获取基因组序列的库，了解到了pyfaidx。

pyfaidx

快速进行随机索引的关键是index，有了index，就能以较少的内存消耗，快速获取任意序列。像samtools, bedtools都会构建index。 而pyfaidx 也利用纯Python实现索引算法,可以FASTA的序列index, 序列获取，原地修改操作。

使用

下载

下载用一行命令就行：

pip install pyfaidx 

用法

用法也比较简单

from pyfaidx import Fasta

# GRCh38.primary_assembly.genome.fa 换成自己fasta 序列路径
genes = Fasta('GRCh38.primary_assembly.genome.fa')
genes 
# Fasta("GRCh38.primary_assembly.genome.fa")

查看多少个染色体序列

>>> genes.keys()
odict_keys(['chr1', 'chr2', 'chr3', 'chr4', 'chr5', 'chr6', 'chr7', 'chr8', 'chr9', 'chr10', 'chr11', 'chr12', 'chr13', 'chr14', 'chr15', 'chr16', 'chr17', 'chr18', 'chr19', 'chr20', 'chr21', 'chr22', 'chrX', 'chrY', 'chrM'])

获取指点区间的序列

>>> genes["chrM"][0:5]
>chrM:1-5
GATCA
>>> genes["chrM"][0:5].seq
GATCA
>>> genes["chrM"][0:5].name
'chrM'
>>> genes["chrM"][0:5].name
'chrM'
>>> genes["chrM"][0:5].start
1
>>> genes["chrM"][0:5].end
5

互补，反向序列的获取

>>> genes["chrM"][0:5]
>chrM:1-5
GATCA
# 单纯反向序列
>>> genes["chrM"][0:5].reverse
>chrM:5-1
ACTAG
# 互补序列
>>> genes["chrM"][0:5].complement
>chrM:1-5 (complement)
CTAGT
# 反向互补序列
>>> genes["chrM"][0:5].complement.reverse
>chrM:5-1 (complement)
TGATC
# 反向互补序列，还可以这样
>>> -genes["chrM"][0:5]
>chrM:5-1 (complement)
TGATC

原文地址：http://www.cnblogs.com/huanping/p/16875533.html