如何根据标识符列表提取参考文本的特定部分？

Question

我有一个参考文件 (.fasta) 和一个基因 ID 列表。对于基因ID列表中的每个ID，我需要将对应的序列放入一个文本文件中。如何实现自动化？

到目前为止我尝试过的事情：

1. sed

sed -n -e '/{GENEID1}/,/>/p' referencefile.fasta | sed $d >> seqs.txt

'>' 是我希望 sed 停止的字符。我需要第二个 sed 来删除最后一行，它也抓住了下一个序列的第一行。 如果我只运行一次，这将有效，但如果我尝试

cat geneID.txt | xargs sed -n -e '/{}/,/>/p' referencefile.fasta >> seqs.txt

然后我只得到一个 ID 列表，没有序列。它也需要很长时间，所以我假设 sed 正在读取参考文件，但我不明白为什么它不会抓取序列？

1. grep

grep -o -P '(?={GENEID}).*(?=>)

在这里我遇到了同样的问题 - 单独工作，但不适用于 xargs 或循环。

1. 使用 for 循环

    for LINE in $(cat geneIDs.txt); do
    echo $LINE >> seqs.txt
    sed -n -e '/$LINE/,/>/p' referencefile.fasta | sed $d >> seqs.txt
    done
   

我也愿意在 python 中尝试一些东西，尽管我还不是很精通它。我的初步尝试是基于this question here。我有一个 10 行的测试 ID 列表，我尝试这样运行：

t = open('test.txt', 'r')
test = t.readlines()
test = test.split()
t.close()

with open('referencefile.fasta', 'r') as ref:
    for line in ref:
        for i in test:
            if i in line:
                print(line)

这个，我什至无法从参考文件中得到一个序列，不管循环。

你们能发现问题吗？为什么这些都不会给我序列？

提前致谢！

编辑添加：

示例参考：

>000000F
ctatcttcgaggttgccacctgtatcgaggagttggcgtctagatcacgaacatgtattttagctatcgtgagctcacacctgacggatccagctttcgaggtcacatcctcaagtctcg


>000001F
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

>000002F

TGCGTGAGGTGCTAGGGATGACAATTGAAAAGAGGACATTGATCGATCACTTGACTCATTTCAGAAAGGAGTTTGGGTTGTCCAACAAGTTGAGGGGGATGATCATCAGGCATCCTGAGT

测试 ID： 000000F, 000001F

理想结果：

000000F ctatcttcgaggttgccacctgtatcgaggagttggcgtctagatcacgaacatgtattttagctatcgtgagctcacacctgacggatccagctttcgaggtcacatcctcaagtctcg

000001F NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

当前结果：

000000F 000001F

Answer 1

如果您的 fasta 文件中的一个geneId 后面总是有一行，这将有所帮助：

grep -A1 -Fwf geneIds.txt input.fasta

检查这个例子：

$  head -n 20 *
==> ids.txt <==
000000F
000001F

==> input.fasta <==
>000000F
Yes I want it!


>000001F
Yes I want it too!

>000002F
skip

>00000XYZ
skip

kent$  grep -A1 -Fwf ids.txt input.fasta
>000000F
Yes I want it!
--
>000001F
Yes I want it too!

Answer 2

取决于大小和访问模式以及您可以使用的其他序列，因为它可能是最简单的构建一个 BLAST 数据库，然后输入您的标识符，它会准确返回您要求的内容（格式正确的 FASTA 除外） .

优点是它设计精良、经过测试且速度快

缺点是它可能对你的任务来说太过分了

（但如果您继续在这个领域工作，仍然非常有用）

https://duckduckgo.com/?q=build+a+blast+database&ia=web

Answer 3

给定：

$ cat file
>000000F
ctatcttcgaggttgccacctgtatcgaggagttggcgtctagatcacgaacatgtattttagctatcgtgagctcacacctgacggatccagctttcgaggtcacatcctcaagtctcg


>000001F
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

使用awk，您可以读取paragraph mode 中由两个或多个\n 分隔的数据。这使您可以轻松地为该格式的文件构建关联数据库。

例如，按确切的字符串搜索：

awk -v RS= -v FS="\n" -v q=">000000F" '$1==q{print $2}' file
ctatcttcgaggttgccacctgtatcgaggagttggcgtctagatcacgaacatgtattttagctatcgtgagctcacacctgacggatccagctttcgaggtcacatcctcaagtctcg

或通过正则表达式搜索：

awk -v RS= -v FS="\n" -v q="[01]F$" '$1~q {print $2}' file
ctatcttcgaggttgccacctgtatcgaggagttggcgtctagatcacgaacatgtattttagctatcgtgagctcacacctgacggatccagctttcgaggtcacatcctcaagtctcg
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

或者，构建一个关联数组：

awk -v RS= -v FS="\n"   '{arr[$1]=$2} END{ "do something with the data in arr" }' file

您可以使用它从具有 id 列表的文件中打印：

cat ids
>000001F
>000000F

awk -v RS= -v FS="\n"  'FNR==NR{for(i=1; i<=NF; i++) ids[$i]; next}
$1 in ids{print $2}' ids file
ctatcttcgaggttgccacctgtatcgaggagttggcgtctagatcacgaacatgtattttagctatcgtgagctcacacctgacggatccagctttcgaggtcacatcctcaagtctcg
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN