我曾尝试针对以下问题制定解决方案:我有一个 .gff3 文件,我想将基因标头替换为简化名称。原始基因标题和新基因名称都在一个单独的文件中给出,原始名称在第 1 列中,新名称在第 2 列中。 如何使用 sed(我认为 sed 最适合这里)来替换所有出现在 .gff3 文件的第二列中使用新的缩短名称?
示例行 .gff3 文件:
tulip_contig_65_pilon_pilon . contig 1 93354 . . . ID=tulip_contig_65_pilon_pilon;Name=tulip_contig_65_pilon_pilon
tulip_contig_65_pilon_pilon maker gene 19497 23038 . + . ID=maker-tulip_contig_65_pilon_pilon-augustus-gene-0.4;Name=maker-tulip_contig_65_pilon_pilon-augustus-gene-0.4
tulip_contig_65_pilon_pilon maker mRNA 19497 23038 . + . ID=maker-tulip_contig_65_pilon_pilon-augustus-gene-0.4-mRNA-1;Parent=maker-tulip_contig_65_pilon_pilon-augustus-gene-0.4;Name=maker-tulip_contig_65_pilon_pilon-augustus-gene-0.4-mRNA-1;_AED=0.00;_eAED=0.00;_QI=418|1|1|1|0|0|3|2100|206
Run Code Online (Sandbox Code Playgroud)
示例行替换文件:
Run Code Online (Sandbox Code Playgroud)augustus_masked-tulip_contig_306_pilon_pilon-processed-gene-0.1 gene1 maker-tulip_contig_306_pilon_pilon-augustus-gene-0.12 gene2 maker-tulip_contig_65_pilon_pilon-augustus-gene-0.4 gene3
预期结果:
Run Code Online (Sandbox Code Playgroud)tulip_contig_65_pilon_pilon . contig 1 93354 . . . ID=tulip_contig_65_pilon_pilon;Name=tulip_contig_65_pilon_pilon tulip_contig_65_pilon_pilon maker gene 19497 23038 . + . ID=gene3;Name=gene3 tulip_contig_65_pilon_pilon maker mRNA 19497 23038 . + . ID=gene3-mRNA-1;Parent=gene3;Name=gene3-mRNA-1;_AED=0.00;_eAED=0.00;_QI=418|1|1|1|0|0|3|2100|206
我曾尝试使用:
while read -r pattern replacement; do sed -i "s/$pattern/$replacement/" file.gff3 ; done < rename.txt
但是根据下面的答案,我现在使用 AWK。我使用脚本(与 Ed Morton 给出的缩进完全相同,但在此处复制它会稍微改变它):
Run Code Online (Sandbox Code Playgroud)NR==FNR { map[$1] = $2 next } { for (old in map) { gsub(old,map[old]) } print }
要运行我使用:
awk -f tst.awk rename.txt original.gff3 > new.gff3
Run Code Online (Sandbox Code Playgroud)
但是,此脚本适用于部分正则表达式匹配,而它应该是完全匹配的。如何更改此 awk 脚本以使其完全匹配?
gff 文件有 7369803 行。rename.txt 文件有 18477 行。
这里欢迎任何建议!
这会在 .gff3 的每一行上从之后=到结尾进行完整的字符串匹配,-gene=<number>并且应该比我们之前所做的更快、更健壮几个数量级,因为它只替换了在每一行中实际找到的 1-3 个字符串original.gff3 文件,而不是尝试替换 rename.txt 文件中存在的所有 18,000 多个字符串:
$ cat tst.awk
NR==FNR {
map[$1] = $2
next
}
{
head = ""
tail = $0
while ( match(tail,/((ID|Parent|Name)=)([^;]+-gene-[0-9]+\.[0-9]+)(.*)/,a) ) {
old = a[3]
head = head substr(tail,1,RSTART-1) a[1] (old in map ? map[old] : old)
tail = a[4]
}
print head tail
}
Run Code Online (Sandbox Code Playgroud)
.
$ awk -f tst.awk rename.txt original.gff3
tulip_contig_65_pilon_pilon . contig 1 93354 . . . ID=tulip_contig_65_pilon_pilon;Name=tulip_contig_65_pilon_pilon
tulip_contig_65_pilon_pilon maker gene 19497 23038 . + . ID=gene3;Name=gene3
tulip_contig_65_pilon_pilon maker mRNA 19497 23038 . + . ID=gene3-mRNA-1;Parent=gene3;Name=gene3-mRNA-1;_AED=0.00;_eAED=0.00;_QI=418|1|1|1|0|0|3|2100|206
Run Code Online (Sandbox Code Playgroud)
它使用 GNU awk 作为第三个参数 match() - 我假设您有 GNU awk 可用(或可以安装它),因为您使用的是 GNU sed。
因此,match()被隔离的字符串(即随后被存储在old从当前行)original.gff3这可能是rename.txt(存储在map[]第一挡),然后old in map是测试如果该字符串实际上是rename.txt或不是,如果是这样,替换old为相应的新值map[]。while只要match()不断寻找新的字符串作为当前行上要替换的候选字符串,这一切都在循环中。
因此,而不是下面的原始 awk 脚本(以及您问题中的 sed 脚本)为 18,000 多行中的每一行rename.txt循环一次,上面只为当前行中original.gff3可能需要替换的每个字符串循环一次,根据您发布的示例输入,最多只有 3 次。
原始答案仅基于加速调用 sed 的 shell 循环:
像这样的东西是你需要的:
$ cat tst.awk
NR==FNR {
map[$1] = $2
next
}
{
for (old in map) {
gsub(old,map[old])
}
print
}
Run Code Online (Sandbox Code Playgroud)
.
$ awk -f tst.awk repl.txt foo.gff3
tulip_contig_65_pilon_pilon . contig 1 93354 . . . ID=tulip_contig_65_pilon_pilon;Name=tulip_contig_65_pilon_pilon
tulip_contig_65_pilon_pilon maker gene 19497 23038 . + . ID=gene3;Name=gene3
tulip_contig_65_pilon_pilon maker mRNA 19497 23038 . + . ID=gene3-mRNA-1;Parent=gene3;Name=gene3-mRNA-1;_AED=0.00;_eAED=0.00;_QI=418|1|1|1|0|0|3|2100|206
Run Code Online (Sandbox Code Playgroud)
有一些关于字符串与正则表达式匹配以及完全匹配与部分匹配的决定也适用于您的 shell+sed 循环,因此请考虑您的完整要求并提供示例输入/输出以进行测试,然后我们可以对其进行调整以适应它做你想做的。现在它正在做部分正则表达式匹配,就像你问题中的 sed 命令一样。
| 归档时间: |
|
| 查看次数: |
207 次 |
| 最近记录: |