如何从 fasta 文件中获取索引列表序列？

【问题标题】：How to get the indexed list sequences from a fasta file?如何从 fasta 文件中获取索引列表序列？
【发布时间】：2016-12-03 02:47:41
【问题描述】：

我比较了两个 fasta 文件（具有不同长度的序列和名称）并将共享的序列名称放入列表中。我试图用列表中的名称获取序列。

文件1：

SRR3350720.1

atccaaccactaaagcagtggtatcaacgcagagtacatggggacattcagtgattatggcatgcactgggtc

SRR3350720.3

caggtgcaggtggtgcagtctggggctgaggtgaagaagcctggctcctcggtgaagatctcatgcaaggctt

SRR3350720.5

caggtccagctggtacagtctggggctgaggtgaagaagcctggggcctcagtgaaggtctcctgcaaggttg

SRR3350720.6

caggtgcagttccagccgtggggcgcaggactgttgaagccttcggagaccctgtccctcacctgcgctgtct

list = ['SRR3350720.1','SRR3350720.5']

我在 python 中尝试了脚本。

import HTSeq
fasta_file = HTSeq.FastaReader('file1.fasta', 'r')
for line in fasta_file:
    for ls in list:
        if str(line.name) == ls:
            print str(line)

但是对于每个 NGS 测序，我在列表中有百万个序列和 1 万个序列 ID。如何改进脚本并高效处理数据。

【问题讨论】：

标签： python fasta

【解决方案1】：

您的进程主要是I/O bound，即它不会帮助您在多个 CPU 上并行化代码。加速它的最佳方法是将fasta复制到SSD或直接复制到内存。

关于您的代码，假设您的序列 ID 存储在名为 seq_ids 的列表中。

import HTSeq
seq_ids = ['SRR3350720.1','SRR3350720.5']
fasta_file = HTSeq.FastaReader('file1.fasta', 'r')
for read in fasta_file:
    if str(read.name) in seq_ids:
        print str(read)

解释：

str(read.name) in seq_ids

检查read.name 是否在您的列表中，然后才打印读取本身。

在您的代码中，您正在遍历搜索列表以进行条目读取，即使一个读取匹配您仍在继续您的循环。

如果您只需要标题和一行序列，请尝试使用grep

grep -A1 -w -f list.txt file1.fasta

-w 打印匹配的行

-A1 匹配后打印行

-f 使用 list.txt 中的模式

file1.fasta 被搜索的文件

【讨论】：