Skip to main content
 首页 » 编程设计

python-3.x之按序列大小对 fasta 进行排序

2025年01月19日24zhwl

我目前想按序列大小对 hudge fasta 文件(+10**8 行和序列)进行排序。 fasta 是生物学中用于存储序列(遗传或蛋白质)的明确定义的格式:

>id1

sequence 1 # 可以在几行

>id2

序列 2

...

我运行了一个以 tsv 格式提供给我的工具:

标识符、长度和标识符的字节位置。

现在我正在做的是按长度列对这个文件进行排序,然后我解析这个文件并使用 seek 来检索相应的序列,然后将它附加到一个新文件。

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

我的问题是以下真的很慢这正常吗(需要几天)?我有另一种方法吗?我不是一个纯粹的信息学家,所以我可能会漏掉一些要点,但我相信索引文件并使用搜索是实现这一目标的最佳方式,我错了吗?

请您参考如下方法:

似乎为每个序列打开两个文件可能对运行时间有很大影响。您可以将文件句柄而不是文件名传递给您的获取/写入函数,但我建议使用已建立的 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)