快速、节省内存、pythonic(和命令行)访问 fasta 序列文件
项目描述
- 电子邮件:
- 执照:
麻省理工学院
内容
执行
需要 Python >= 2.6。存储没有空格或标题的 fasta 文件的扁平版本,并使用 numpy 二进制格式的 mmap 或 fseek/fread,因此 永远不会将序列数据读入内存。保存 fasta 文件中每个标头的开始、停止(用于 fseek/mmap)位置的 pickle (.gdx) 以供内部使用。
用法
>>> from pyfasta import Fasta
>>> f = Fasta('tests/data/three_chrs.fasta')
>>> sorted(f.keys())
['chr1', 'chr2', 'chr3']
>>> f['chr1']
NpyFastaRecord(0..80)
切片
# get full the sequence:
>>> a = str(f['chr1'])
>>> b = f['chr1'][:]
>>> a == b
True
>>> f['chr1'][:10]
'ACTGACTGAC'
# get the 1st basepair in every codon (it's python yo)
>>> f['chr1'][::3]
'AGTCAGTCAGTCAGTCAGTCAGTCAGT'
# can query by a 'feature' dictionary (note this is one based coordinates)
>>> f.sequence({'chr': 'chr1', 'start': 2, 'stop': 9})
'CTGACTGA'
# same as:
>>> f['chr1'][1:9]
'CTGACTGA'
# use python, zero based coords
>>> f.sequence({'chr': 'chr1', 'start': 2, 'stop': 9}, one_based=False)
'TGACTGA'
# with reverse complement (automatic for - strand)
>>> f.sequence({'chr': 'chr1', 'start': 2, 'stop': 9, 'strand': '-'})
'TCAGTCAG'
按键功能
有时您的 fasta 会有一个长标题,例如:“AT1G51370.2 | 符号:| F-box 家族蛋白 | chr1:19045615-19046748 FORWARD”,当您只想关闭时:“AT1G51370.2”。在这种情况下,为构造函数指定 key_fn 参数:
>>> fkey = Fasta('tests/data/key.fasta', key_fn=lambda key: key.split()[0])
>>> sorted(fkey.keys())
['a', 'b', 'c']
麻木的
默认是使用 memmaped numpy 数组作为后端。在这种情况下,可以直接取回一个数组……
>>> f['chr1'].as_string = False >>> f['chr1'][:10] # doctest: +NORMALIZE_WHITESPACE memmap(['A', 'C', 'T', 'G', 'A', 'C', 'T', 'G', 'A', 'C'], dtype='|S1') >>> import numpy as np >>> a = np.array(f['chr2']) >>> a.shape[0] == len(f['chr2']) True >>> a[10:14] # doctest: +NORMALIZE_WHITESPACE array(['A', 'A', 'A', 'A'], dtype='|S1')
屏蔽子序列
>>> a[11:13] = np.array('N', dtype='S1')
>>> a[10:14].tostring()
'ANNA'
后端(记录类)
也可以指定另一个记录类作为切片和读取的基础工作。目前,只有默认值:
NpyFastaRecord 使用 numpy memmap
FastaRecord,使用 fseek/fread
MemoryRecord 将所有内容读入内存,并且每次都必须重新解析原始 fasta。
TRCecord 与 NpyFastaRecord 相同,除了它将索引保存在 TokyoCabinet 哈希数据库中,对于有足够记录将整个索引从 pickle 加载到内存中的情况是不明智的。(注意:无论哪种情况,序列都不会加载到内存中)。
可以将与record_class kwarg 一起使用的类指定给Fasta 构造函数:
>>> from pyfasta import FastaRecord # default is NpyFastaRecord
>>> f = Fasta('tests/data/three_chrs.fasta', record_class=FastaRecord)
>>> f['chr1']
FastaRecord('tests/data/three_chrs.fasta.flat', 0..80)
除了 repr,它的行为应该与 Npy 记录类后端完全相同
可以使用 FastaRecord 的子类创建自己的。有关详细信息,请参阅 pyfasta/records.py 中的源代码。
展平
为了有效地访问序列内容,pyfasta 保存一个单独的扁平文件,其中所有换行符和标题从序列中删除。对于大型 fasta 文件,可能不希望保存 2 个 5GG+ 文件的副本。在这种情况下,可以“就地”展平文件,保留所有标题,并保留 fasta 文件的有效性——唯一的变化是从每个序列中删除换行符。这可以通过flatten_inplace = True指定
>>> import os
>>> os.unlink('tests/data/three_chrs.fasta.gdx') # cleanup non-inplace idx
>>> f = Fasta('tests/data/three_chrs.fasta', flatten_inplace=True)
>>> f['chr1'] # note the difference in the output from above.
NpyFastaRecord(6..86)
# sequence from is same as when requested from non-flat file above.
>>> f['chr1'][1:9]
'CTGACTGA'
# the flattened file is kept as a place holder without the sequence data.
>>> open('tests/data/three_chrs.fasta.flat').read()
'@flattened@'
命令行界面
还有一个命令行界面来操作/查看 fasta 文件。pyfasta可执行文件是通过 setuptools 安装的,运行它将显示帮助文本。
将 fasta 文件拆分为 6 个大小相对均匀的新文件:
$ pyfasta split -n 6 original.fasta
将 fasta 文件拆分为每个标题的一个新文件,并在每个文件名中填写“%(seqid)s”。:
$ pyfasta split –header “%(seqid)s.fasta” original.fasta
创建 1 个新的 fasta 文件,将序列拆分为 10K-mers:
$ pyfasta split -n 1 -k 10000 original.fasta
2 个新的 fasta 文件,其序列分为 2K 重叠的 10K-mers:
$ pyfasta split -n 2 -k 10000 -o 2000 original.fasta
显示有关文件的一些信息(并显示 gc 内容):
$ pyfasta info –gc test/data/three_chrs.fasta
从文件中提取序列。使用标头标志创建一个新的 fasta 文件。args 是要提取的序列列表。
$ pyfasta extract –header –fasta test/data/three_chrs.fasta seqa seqb seqc
使用包含新文件中不需要的标头的文件从文件中提取序列:
$ pyfasta extract –header –fasta input.fasta –exclude –file seqids_to_exclude.txt
从具有复杂键的 fasta 文件中提取序列,我们只想根据空格之前的部分进行查找。
$ pyfasta extract –header –fasta input.with.keys.fasta –space –file seqids.txt
原地展平文件,以便 pyfasta 以后更快地使用,而无需创建另一个副本。(压扁)
$ pyfasta 展平 input.fasta
清理
(尽管对于实际使用,这些将保留以加快访问速度)
>>> os.unlink('tests/data/three_chrs.fasta.gdx')
>>> os.unlink('tests/data/three_chrs.fasta.flat')
测试
目前,这 2 个模块和所有包含的记录类的测试覆盖率 > 99%。运行测试:
$ python setup.py nosetests
变化
0.5.2
修复补码(@mruffalo)
0.5.0
python 3 兼容性感谢 mruffalo
0.4.5
pyfasta split 可以处理 > 52 个文件。(感谢德夫图利亚)
0.4.4
修复pyfasta提取物
0.4.3
在 sequence() 中添加基于 0 或 1 的间隔感谢@jamescasbon
0.4.2
更新最新的 numpy(无法关闭 memmap)
0.4.1
检查重复的标题。
0.4.0
将 key_fn kwarg 添加到构造函数
0.3.9
内存映射只需要'r'(不是r+)。
0.3.8
清理混合就地/非就地扁平文件的逻辑。如果就地可用,则始终使用它。
0.3.6/7
不要每次都重新展平文件!
在原始 fasta 中的标题前后允许空格。
0.3.5
更新 README.txt 中的文档以获取新的 CLI 内容。
允许就地展平。
摆脱 memmap(导致更快的解析)。
0.3.4
恢复python2.5兼容性。
CLI:添加从提取中排除序列的功能
CLI:允许基于标头进行拆分。
0.3.3
将此文件包含在 tar 球中(感谢 wen h。)
0.3.2
将后端分离到 records.py
使用鼻子测试(python setup.py 鼻子测试)
如果 tc 已(简单)安装,则为下一代测序添加 TRCecord 后端。
提高测试覆盖率。