sed使用while循环非常慢

Question

我有gff文件，内容如下(制表符分隔):

# start gene 1Chr.g1
1Chr    AUGUSTUS        gene    3636    5916    0.1     +       .       ID=1Chr.g1
1Chr    AUGUSTUS        transcript      3636    5916    0.1     +       .       ID=1Chr.g1.t1;Parent=1Chr.g1
1Chr    AUGUSTUS        transcription_start_site        3636    3636    .       +       .       Parent=1Chr.g1.t1
1Chr    AUGUSTUS        exon    3636    3913    .       +       .       Parent=1Chr.g1.t1
1Chr    AUGUSTUS        start_codon     3760    3762    .       +       0       Parent=1Chr.g1.t1
1Chr    AUGUSTUS        intron  3914    3995    1       +       .       
1Chr    AUGUSTUS        CDS     3760    3913    1       +       0       ID=1Chr.g1.t1.cds;Parent=1Chr.g1.t1
1Chr    AUGUSTUS        stop_codon      5628    5630    .       +       0       Parent=1Chr.g1.t1
1Chr    AUGUSTUS        transcription_end_site  5916    5916    .       +       .       Parent=1Chr.g1.t1
# start gene 1Chr.g2
1Chr    AUGUSTUS        gene    5938    8761    0.17    -       .       ID=1Chr.g2
1Chr    AUGUSTUS        transcript      5938    8761    0.17    -       .       ID=1Chr.g2.t1;Parent=1Chr.g2
1Chr    AUGUSTUS        transcription_end_site  5938    5938    .       -       .       Parent=1Chr.g2.t1
1Chr    AUGUSTUS        exon    5938    6594    .       -       .       Parent=1Chr.g2.t1
1Chr    AUGUSTUS        stop_codon      6428    6430    .       -       0       Parent=1Chr.g2.t1
1Chr    AUGUSTUS        intron  6595    7156    0.8     -       .       Parent=1Chr.g2.t1
1Chr    AUGUSTUS        CDS     6428    6594    0.89    -       2       ID=1Chr.g2.t1.cds;Parent=1Chr.g2.t1
# start gene 2Chr.g1
2Chr    AUGUSTUS        gene    11612   13481   0.09    -       .       ID=2Chr.g1
2Chr    AUGUSTUS        transcript      11612   13481   0.09    -       .       ID=2Chr.g1.t1;Parent=2Chr.g1
2Chr    AUGUSTUS        transcription_end_site  11612   11612   .       -       .       Parent=2Chr.g1.t1
2Chr    AUGUSTUS        exon    11612   13481   .       -       .       Parent=2Chr.g1.t1
2Chr    AUGUSTUS        stop_codon      11864   11866   .       -       0       Parent=2Chr.g1.t1
2Chr    AUGUSTUS        CDS     11864   12940   1       -       0       ID=2Chr.g1.t1.cds;Parent=2Chr.g1.t1
2Chr    AUGUSTUS        start_codon     12938   12940   .       -       0       Parent=2Chr.g1.t1
2Chr    AUGUSTUS        transcription_start_site        13481   13481   .       -       .       Parent=2Chr.g1.t1
# start gene 2Chr.g2
2Chr    AUGUSTUS        gene    22876   31223   0.04    +       .       ID=2Chr.g2
2Chr    AUGUSTUS        transcript      22876   31223   0.04    +       .       ID=2Chr.g2.t1;Parent=2Chr.g2
2Chr    AUGUSTUS        transcription_start_site        22876   22876   .       +       .       Parent=2Chr.g2.t1
2Chr    AUGUSTUS        exon    22876   23456   .       +       .       Parent=2Chr.g2.t1
2Chr    AUGUSTUS        exon    23515   24451   .       +       .       Parent=2Chr.g2.t1
2Chr    AUGUSTUS        start_codon     23519   23521   .       +       0       Parent=2Chr.g2.t1

我想替换1Chr.g1,1Chr.g2,2Chr.g1和2Chr.g2基因的id，就像从g1开始到id结束的顺序，就像在这种情况下的g4。

预期输出

# start gene g1
1Chr    AUGUSTUS        gene    3636    5916    0.1     +       .       ID=g1
1Chr    AUGUSTUS        transcript      3636    5916    0.1     +       .       ID=g1.t1;Parent=g1
1Chr    AUGUSTUS        transcription_start_site        3636    3636    .       +       .       Parent=g1.t1
1Chr    AUGUSTUS        exon    3636    3913    .       +       .       Parent=g1.t1
1Chr    AUGUSTUS        start_codon     3760    3762    .       +       0       Parent=g1.t1
1Chr    AUGUSTUS        intron  3914    3995    1       +       .
1Chr    AUGUSTUS        CDS     3760    3913    1       +       0       ID=g1.t1.cds;Parent=g1.t1
1Chr    AUGUSTUS        stop_codon      5628    5630    .       +       0       Parent=g1.t1
1Chr    AUGUSTUS        transcription_end_site  5916    5916    .       +       .       Parent=g1.t1
# start gene g2
1Chr    AUGUSTUS        gene    5938    8761    0.17    -       .       ID=g2
1Chr    AUGUSTUS        transcript      5938    8761    0.17    -       .       ID=g2.t1;Parent=g2
1Chr    AUGUSTUS        transcription_end_site  5938    5938    .       -       .       Parent=g2.t1
1Chr    AUGUSTUS        exon    5938    6594    .       -       .       Parent=g2.t1
1Chr    AUGUSTUS        stop_codon      6428    6430    .       -       0       Parent=g2.t1
1Chr    AUGUSTUS        intron  6595    7156    0.8     -       .       Parent=g2.t1
1Chr    AUGUSTUS        CDS     6428    6594    0.89    -       2       ID=g2.t1.cds;Parent=g2.t1
# start gene g3
2Chr    AUGUSTUS        gene    11612   13481   0.09    -       .       ID=g3
2Chr    AUGUSTUS        transcript      11612   13481   0.09    -       .       ID=g3.t1;Parent=g3
2Chr    AUGUSTUS        transcription_end_site  11612   11612   .       -       .       Parent=g3.t1
2Chr    AUGUSTUS        exon    11612   13481   .       -       .       Parent=g3.t1
2Chr    AUGUSTUS        stop_codon      11864   11866   .       -       0       Parent=g3.t1
2Chr    AUGUSTUS        CDS     11864   12940   1       -       0       ID=g3.t1.cds;Parent=g3.t1
2Chr    AUGUSTUS        start_codon     12938   12940   .       -       0       Parent=g3.t1
2Chr    AUGUSTUS        transcription_start_site        13481   13481   .       -       .       Parent=g3.t1
# start gene g4
2Chr    AUGUSTUS        gene    22876   31223   0.04    +       .       ID=g4
2Chr    AUGUSTUS        transcript      22876   31223   0.04    +       .       ID=g4.t1;Parent=g4
2Chr    AUGUSTUS        transcription_start_site        22876   22876   .       +       .       Parent=g4.t1
2Chr    AUGUSTUS        exon    22876   23456   .       +       .       Parent=g4.t1
2Chr    AUGUSTUS        exon    23515   24451   .       +       .       Parent=g4.t1
2Chr    AUGUSTUS        start_codon     23519   23521   .       +       0       Parent=g4.t1

我写了下面的bash脚本，但是它花了太长时间，因为我试图计算它的时间，所以对于一个sed它花了1秒，如果有28000迭代，它将花费大约8小时，这是太多的时间。有什么有效的方法吗?

awk '$3 == "gene"' $1 |cut -f9 |grep -o "=.*" |sed -e 's/=//g' >LIST.txt

COUNTER=0
cat LIST.txt | while read line; do
COUNTER=$(expr $COUNTER + 1)
echo "sed -i 's/$line/g$COUNTER/g' $1" |bash
done
rm LIST.txt

另一件事，生成一个文件sedTG45，这是非常烦人的。

Answer 1

您在同一个文件上运行sed多达28,000次。只要稍加预处理，就不难把它压缩到一次。

这是未经测试的，但至少应该为您指明Awk的大致方向。这是一种很小的语言;你可以在不到一个小时的时间里学会它，并在几个星期内变得相当熟练。

awk -F 't' '$3 == "gene" {
g=$9; sub(/^[^=]*=/, "", g); gsub(/=/, "", g);
a[g] = "g" ++n }
{ for(k in a) gsub(k, a[k]) }1' "$1"

简而言之，n维持计数器，a是我们取出的所有基因的关联数组。

这假设定义在应该被替换的出现之前。如果这不是一个有效的假设，您将需要对文件进行两次迭代，但第一次将是只读的，因此这仍然应该比您的暴力尝试要快得多。

---

附录:如果您执意要坚持使用(已更新的)代码，那么在done之后将管道移到Bash(或者仅仅是sed)将会显著改进它。这是一个简单的重构。用Awk还是会更漂亮。

# Use lower case for private variable
counter=0
# Note quotes around $1
# and use of pipe instead of temp file
awk '$3 == "gene"' "$1" |
cut -f9 |
grep -o "=.*" |
sed -e 's/=//g' |
# note IFS='' and read -r
while IFS='' read -r line; do
# avoid paleolithic expr
((counter++))
echo "s/$line/g$COUNTER/g"
done |
# pipe output to sed -i
sed -i -f - "$1"

不是所有的seds都允许-f -。当sed -i运行时，不可避免地会有一个临时文件，但它会在完成后被删除。

演示，带时序:https://ideone.com/GghqiW

Answer 2

首先生成sed命令可能会有所帮助。输入文件将被读取两次，但这应该不是问题。请检查演出。
你想要的sed命令看起来像

s/1Chr.g1.t1/g1/g
s/1Chr.g2.t1/g2/g
s/2Chr.g1.t1/g3/g
s/2Chr.g2.t1/g4/g

当数字大于9时，这些命令可能是错误的，所以稍微改变一下:

awk '
($0 != last) {
n++;
printf("s/([^0-9])%s([^0-9]|$)/\1g%s\2/gn", $0, n)
}
{
last=$0
}' <(grep -Eo "[0-9]+Chr.g[0-9]+" < file)

返回

s/([^0-9])1Chr.g1([^0-9]|$)/1g12/g
s/([^0-9])1Chr.g2([^0-9]|$)/1g22/g
s/([^0-9])2Chr.g1([^0-9]|$)/1g32/g
s/([^0-9])2Chr.g2([^0-9]|$)/1g42/g

现在在sed中使用这些命令(将文件替换为$1在2个地方，将-i添加到sed，当它看起来ok时):

sed -r -f <(
awk '
($0 != last) {
n++;
printf("s/([^0-9])%s([^0-9]|$)/\1g%s\2/gn", $0, n)
}
{
last=$0
}' <(grep -Eo "[0-9]+Chr.g[0-9]+" < file)
) file

Answer 3

sed 's/([= ])[0-9]+Chr.(g[0-9]+)/12/g' file.gff

如果你敢做sed -i，
我建议你先做一个副本。

还注意:将只与基于数字的染色体工作所以没有X,Y, mt,I,II,III, IV…等。