我目前想按序列大小对 hudge fasta 文件(+10**8 行和序列)进行排序。 fasta 是生物学中用于存储序列(遗传或蛋白质)的明确定义的格式:
>id1
序列 1 # 可以位于多行
>id2
序列2
...
我运行了一个提供 tsv 格式的工具:
标识符、长度以及标识符的位置(以字节为单位)。
现在我正在做的是按长度列对这个文件进行排序,然后解析这个文件并使用eek来检索相应的序列,然后将其附加到一个新文件中。
# this fonction will get the sequence using seek
def get_seq(file, bites):
with open(file) as f_:
f_.seek(bites, 0) # go to the line of interest
line = f_.readline().strip() # this line is the begin of the
#sequence
to_return = "" # init the string which will contains the sequence
while not line.startswith('>') or not line: # while we do not
# encounter another identifiant
to_return += line
line = f_.readline().strip()
return to_return
# simply append to a file the id and the sequence
def write_seq(out_file, id_, sequence):
with open(out_file, 'a') as out_file:
out_file.write('>{}\n{}\n'.format(id_.strip(), sequence))
# main loop will parse the index file and call the function defined below
with open(args.fai) as ref:
indice = 0
for line in ref:
spt = line.split()
id_ = spt[0]
seq = get_seq(args.i, int(spt[2]))
write_seq(out_file=args.out, id_=id_, sequence=seq)
我的问题是以下速度真的很慢,这是否正常(需要几天)?我还有其他方法吗?我不是一个纯粹的信息学家,所以我可能会错过一些要点,但我相信索引文件并使用搜索是实现这一目标的最快方法,我错了吗?
似乎为每个序列打开两个文件可能会对运行时间产生很大影响。您可以将文件句柄传递给 get/write 函数而不是文件名,但我建议使用已建立的 fasta 解析器/索引器,如 biopython 或 samtools。这是使用 samtools 的(未经测试的)解决方案:
subprocess.call(["samtools", "faidx", args.i])
with open(args.fai) as ref:
for line in ref:
spt = line.split()
id_ = spt[0]
subprocess.call(["samtools", "faidx", args.i, id_, ">>", args.out], shell=True)
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)