使用awk通过文件中的ID从multifasta文件中提取序列

Question

我想从 multifasta 文件中提取与单独的 ID 列表给出的 ID 相匹配的序列。

FASTA 文件 seq.fasta:

>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11605
TTCAGCAAGCCGAGTCCTGCGTCGAGAGTTCAAGTC
CCTGTTCGGGCGCCACTGCTAG
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC
>7P58X:01334:11635
TTCAGCAAGCCGAGTCCTGCGTCGAGAGATCGCTTT
CAAGTCCCTGTTCGGGCGCCACTGCGGGTCTGTGTC
GAGCG
>7P58X:01336:11621
ACGCTCGACACAGACCTTTAGTCAGTGTGGAAATCT
CTAGCAGTAGAGGAGATCTCCTCGACGCAGGACT

ID文件id.txt：

7P58X:01332:11636
7P58X:01334:11613

我想获取仅包含与 id.txt 文件中的 ID 匹配的序列的 fasta 文件：

>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC

我真的很喜欢我在答案here和here中找到的awk方法，但是那里给出的代码对于我给出的示例仍然不能完美工作。原因如下：

(1)

awk -v seq="7P58X:01332:11636" -v RS='>' '$1 == seq {print RS $0}' seq.fasta

此代码适用于多行序列，但必须将 ID 单独插入到代码中。

(2)

awk 'NR==FNR{n[">"$0];next} f{print f ORS $0;f=""} $0 in n{f=$0}' id.txt seq.fasta

此代码可以从 id.txt 文件中获取 ID，但仅返回多行序列的第一行。

我想最好的办法是修改代码（2）中的 RS 变量，但到目前为止我所有的尝试都失败了。请问有人可以帮我吗？

Answer 1

$ awk -F'>' 'NR==FNR{ids[$0]; next} NF>1{f=($2 in ids)} f' id.txt seq.fasta
>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC

Answer 2

以下

awk

可能会对您有所帮助。

awk 'FNR==NR{a[$0];next} /^>/{val=$0;sub(/^>/,"",val);flag=val in a?1:0} flag' ids.txt  fasta_file

Answer 3

我也面临着类似的问题。我的 multi-fasta 文件大小约为 25G。 我使用 sed 而不是 awk，尽管我的解决方案是一个丑陋的黑客。

首先，我将每个序列的标题的行号提取到数据文件中。

grep -n ">" multi-fasta.fa > multi-fasta.idx 

我得到的是这样的：

1:>DM_0000000004
5:>DM_0000000005
11:>DM_0000000007
19:>DM_0000000008
23:>DM_0000000009

然后，我通过标题提取了想要的序列，例如。 DM_0000000004，使用下面的脚本。

seqnm=$1
idx0_idx1=`grep -n $seqnm multi-fasta.idx`

idx0=`echo $idx0_idx1 | cut -d ":" -f 1`
idx0plus1=`expr $idx0 + 1`

idx1=`echo $idx0_idx1 | cut -d ":" -f 2`

idx2=`head -n $idx0plus1 multi-fasta.idx | tail -1 | cut -d ":" -f 1`
idx2minus1=`expr $idx2 - 1`

sed ''"$idx1"','"$idx2minus1"'!d' multi-fasta.fa > ${seqnm}.fasta

例如我想提取DM_0000016115的序列。 idx0_idx1变量给了我：

7507:42520:>DM_0000016115

7507 (idx0) 是 multi-fasta.idx 中第 42520:>DM_0000016115 的行号。

42520 (idx1) 是multi-fasta.fa中>DM_0000016115行的行号。

idx2 是所需标题正下方的序列标题的行号 (>DM_0000016115)。

最后，使用sed，我们可以提取idx1和idx2减去1之间的行，即标题和序列，在这种情况下你可以使用grep -A。

这个丑陋的 hack 的优点是它不需要多 fasta 文件中每个序列的特定行数。

令我困扰的是这个过程很慢。对于我的25G multi-fasta 文件来说，这样的提取需要几十秒。不过，它比使用 samtools faidx 快得多。

Answer 4

这里只是为仍然遇到问题的人提供有关解决方案的说明：

请务必首先使用

dos2unix id.txt

删除每行末尾隐藏的

^M

Windows 插入内容。

此外，此解决方案要求

id.txt

上的名称在每行开头不包含

>

。