我知道有很多类似的问题,我已经通读了很多。但我仍然无法使我的代码正常工作。有人能帮我指出这个问题吗?谢谢
(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.pep2
0:一次循环输入一行,默认情况下将其分配给$_
。在每次循环迭代后添加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] != ""
就是这么做的。
本例假设第一个文件只有两个字段(数据中没有空格)。如果有空格,我可以发布编辑。