使用awk通过文件中的ID从multifasta文件中提取序列
extract sequences from multifasta file by ID in file using awk
我想从 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 and 中找到的 awk 方法,但是那里给出的代码对于我给出的示例仍然不能完美运行。原因如下:
(1)
awk -v seq="7P58X:01332:11636" -v RS='>' ' == seq {print RS [=13=]}' seq.fasta
此代码适用于多行序列,但 ID 必须单独插入到代码中。
(2)
awk 'NR==FNR{n[">"[=14=]];next} f{print f ORS [=14=];f=""} [=14=] in n{f=[=14=]}' id.txt seq.fasta
此代码可以从 id.txt 文件中获取 ID,但 returns 只能获取多行序列的第一行。
我想最好的办法是修改代码 (2) 中的 RS 变量,但到目前为止我的所有尝试都失败了。请问有人可以帮我吗?
以下 awk
可能会对您有所帮助。
awk 'FNR==NR{a[[=10=]];next} /^>/{val=[=10=];sub(/^>/,"",val);flag=val in a?1:0} flag' ids.txt fasta_file
$ awk -F'>' 'NR==FNR{ids[[=10=]]; next} NF>1{f=( in ids)} f' id.txt seq.fasta
>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC
我遇到了类似的问题。我的 multi-fasta 文件大小约为 25G。
我使用 sed 而不是 awk,尽管我的解决方案是一个丑陋的 hack。
首先,我将每个序列标题的行号提取到一个数据文件中。
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=
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)是第42520行的行号:>DM_0000016115 在 multi-fasta.idx.
42520(idx1)是行号>DM_0000016115 在 multi-fasta.fa.
idx2 是所需序列标题正下方的行号(>DM_0000016115).
最后,使用sed,我们可以提取idx1和idx2 minus 1,即标题和序列,此时可以使用grep -A .
这个ugly-hack的优点是它不需要multi-fasta文件中每个序列的特定行数。
困扰我的是这个过程很慢。对于我的 25G multi-fasta 文件,这样的提取需要几十秒。但是,它比使用 samtools faidx 快得多。
我想从 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 and
(1)
awk -v seq="7P58X:01332:11636" -v RS='>' ' == seq {print RS [=13=]}' seq.fasta
此代码适用于多行序列,但 ID 必须单独插入到代码中。
(2)
awk 'NR==FNR{n[">"[=14=]];next} f{print f ORS [=14=];f=""} [=14=] in n{f=[=14=]}' id.txt seq.fasta
此代码可以从 id.txt 文件中获取 ID,但 returns 只能获取多行序列的第一行。
我想最好的办法是修改代码 (2) 中的 RS 变量,但到目前为止我的所有尝试都失败了。请问有人可以帮我吗?
以下 awk
可能会对您有所帮助。
awk 'FNR==NR{a[[=10=]];next} /^>/{val=[=10=];sub(/^>/,"",val);flag=val in a?1:0} flag' ids.txt fasta_file
$ awk -F'>' 'NR==FNR{ids[[=10=]]; next} NF>1{f=( in ids)} f' id.txt seq.fasta
>7P58X:01332:11636
TTCAGCAAGCCGAGTCCTGCGTCGTTACTTCGCTT
CAAGTCCCTGTTCGGGCGCC
>7P58X:01334:11613
ACGAGTGCGTCAGACCCTTTTAGTCAGTGTGGAAAC
我遇到了类似的问题。我的 multi-fasta 文件大小约为 25G。 我使用 sed 而不是 awk,尽管我的解决方案是一个丑陋的 hack。
首先,我将每个序列标题的行号提取到一个数据文件中。
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=
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)是第42520行的行号:>DM_0000016115 在 multi-fasta.idx.
42520(idx1)是行号>DM_0000016115 在 multi-fasta.fa.
idx2 是所需序列标题正下方的行号(>DM_0000016115).
最后,使用sed,我们可以提取idx1和idx2 minus 1,即标题和序列,此时可以使用grep -A .
这个ugly-hack的优点是它不需要multi-fasta文件中每个序列的特定行数。
困扰我的是这个过程很慢。对于我的 25G multi-fasta 文件,这样的提取需要几十秒。但是,它比使用 samtools faidx 快得多。