awk合并来自两个文件的信息(fasta文件头)



我知道有很多类似的问题,我已经通读了很多。但我仍然无法使我的代码正常工作。有人能帮我指出这个问题吗?谢谢

(base) $ head Sample.pep2
>M00000032072 gene=G00000025773 seq_id=ChrM type=cds
MFKQNPSPGWKECPPSSDKEGTTPERLDEGREMRRGKEKAFGDREISFLLHRKRRRPRIA
YGACYLKGARFFDRGAMIAGASPRSARWPIGIAACGLCLPIRIIIKNSGSARESAGNNRK
EGVHVAAAPAPLLSQWGSRFGIY*
>M00000032073 gene=G00000025774 seq_id=ChrM type=cds
MFKQNPSPGWKECPPSSDKEGTTPERLDEGREMRRGKEKAFGDREISFLLHRKRRRPRIA
YGACYLKGARFFDRGAMIAGASPRSARWPIGIAACGLCLPIRIIIKNSGSARESAGNNRK
EGVHVAAAPAPLLSQWGSSIASMILGALAAMAQTKVKRPLAHSSIGHVGYIRTGFSCGTI
EGIQSLLIGIFIYALMTMDAFAIVSALRQTRVKYIADLGALAKTNPISAITFSITMFSYA
GIPPLAGFCSKFYLFFAALGCGAYFLAPVGVVTSVIGRWAAGRLPRISKFGGPKAVLRAP
$ head -n 3 mRNA.function
M00000032074 locus=g17091;makerName=TCONS_00021197.p2
M00000032073 Dbxref=MobiDBLite:mobidb-lite;locus=g17092;makerName=TCONS_00021198.p3
M00000032072 Dbxref=MobiDBLite:mobidb-lite;locus=g17093;makerName=TCONS_00021199.p1

我想要输出

>M00000032072 gene=G00000025773 seq_id=ChrM type=cds Dbxref=MobiDBLite:mobidb-lite;locus=g17093;makerName=TCONS_00021199.p1
MFKQNPSPGWKECPPSSDKEGTTPERLDEGREMRRGKEKAFGDREISFLLHRKRRRPRIA
YGACYLKGARFFDRGAMIAGASPRSARWPIGIAACGLCLPIRIIIKNSGSARESAGNNRK
EGVHVAAAPAPLLSQWGSRFGIY*
>M00000032073 gene=G00000025774 seq_id=ChrM type=cds Dbxref=MobiDBLite:mobidb-lite;locus=g17092;makerName=TCONS_00021198.p3
MFKQNPSPGWKECPPSSDKEGTTPERLDEGREMRRGKEKAFGDREISFLLHRKRRRPRIA
YGACYLKGARFFDRGAMIAGASPRSARWPIGIAACGLCLPIRIIIKNSGSARESAGNNRK
EGVHVAAAPAPLLSQWGSSIASMILGALAAMAQTKVKRPLAHSSIGHVGYIRTGFSCGTI
EGIQSLLIGIFIYALMTMDAFAIVSALRQTRVKYIADLGALAKTNPISAITFSITMFSYA
GIPPLAGFCSKFYLFFAALGCGAYFLAPVGVVTSVIGRWAAGRLPRISKFGGPKAVLRAP

我的命令是CCD_ 1。但它不起作用。。。我不知道哪里错了。。。

这里有一个perl解决方案。

perl -lpe 'BEGIN { %id_to_function = map { /^(S+)s+(.*)/ } `cat mRNA.function`; } s{^>(S+)(.*)}{>$1$2 $id_to_function{$1}};' sample.pep2

在读取fasta文件之前,代码执行BEGIN { ... }块。在那里,具有id和函数的文件被读取到散列%id_to_function
在代码的主体中,替换运算符s{...}{...}使用哈希查找$id_to_function{$1}将对应id的函数附加到fasta标头
$1$2分别是第一个和第二个捕获组,它们在前面的regex:^>(S+)(.*)中用括号捕获。

Perl单行使用以下命令行标志:
-e:neneneba告诉Perl在线查找代码,而不是在文件中查找代码
awk 'NR==FNR{id[$1]=$2; next} /^>/ {print $0=$0,id[$1]}' mRNA.function Sample.pep20:一次循环输入一行,默认情况下将其分配给$_。在每次循环迭代后添加print $_
-l:在线执行代码之前,剥去输入行分隔符(默认为*NIX上的"n"),并在打印时附加。

另请参阅:
perldoc perlrun:如何执行Perl解释器:命令行开关
perldoc perlrequick:Perl正则表达式快速启动

当前代码很接近,只需要a)将第一个字段(Sample.pep2)与相应的数组条目(如果存在)匹配,b)确保打印所有输入行:

awk 'NR==FNR{id[$1]=$2; next} /^>/ {key=substr($1,2); if (key in id) $0=$0 OFS id[key]} 1' mRNA.function Sample.pep2

这将生成:

>M00000032072 gene=G00000025773 seq_id=ChrM type=cds Dbxref=MobiDBLite:mobidb-lite;locus=g17093;makerName=TCONS_00021199.p1
MFKQNPSPGWKECPPSSDKEGTTPERLDEGREMRRGKEKAFGDREISFLLHRKRRRPRIA
YGACYLKGARFFDRGAMIAGASPRSARWPIGIAACGLCLPIRIIIKNSGSARESAGNNRK
EGVHVAAAPAPLLSQWGSRFGIY*
>M00000032073 gene=G00000025774 seq_id=ChrM type=cds Dbxref=MobiDBLite:mobidb-lite;locus=g17092;makerName=TCONS_00021198.p3
MFKQNPSPGWKECPPSSDKEGTTPERLDEGREMRRGKEKAFGDREISFLLHRKRRRPRIA
YGACYLKGARFFDRGAMIAGASPRSARWPIGIAACGLCLPIRIIIKNSGSARESAGNNRK
EGVHVAAAPAPLLSQWGSSIASMILGALAAMAQTKVKRPLAHSSIGHVGYIRTGFSCGTI
EGIQSLLIGIFIYALMTMDAFAIVSALRQTRVKYIADLGALAKTNPISAITFSITMFSYA
GIPPLAGFCSKFYLFFAALGCGAYFLAPVGVVTSVIGRWAAGRLPRISKFGGPKAVLRAP
awk '
NR==FNR {a[">"$1]=$2}
NR!=FNR && a[$1] != "" {$0=$0" "a[$1]}
NR!=FNR' mRNA.function Sample.pep2

你的阵列走在了正确的轨道上。在文件2中,在附加数据之前,您需要检查该行是否以文件1中的匹配名称开头。a[$1] != ""就是这么做的。

本例假设第一个文件只有两个字段(数据中没有空格)。如果有空格,我可以发布编辑。

最新更新