bash 脚本 - 使用 sed 中的模式列表删除子字符串

bash script - Use patterns list in sed to remove substrings

我有这个文件 (adapters.txt),其中包含模式列表:

cactctttccctacacgacgctcttccg
cactctttccctacacgacgctcttccgaatcta
cactctttccctacacgacgctcttccgaatctaatt
cactctttccctacacgacgctcttccgaatctaatta
cactctttccctacacgacgctcttccgaatctag
cactctttccctacacgacgctcttccgaatctagc
cactctttccctacacgacgctcttccgacctcattcc
cactctttccctacacgacgctcttccgacctcattcccaccctcttccg
cactctttccctacacgacgctcttccgatc
cactctttccctacacgacgctcttccgatccaatt
cactctttccctacacgacgctcttccgatttagc
cactctttccctacacgacgctcttccgatttagct
cactctttccctacacgacgctcttccgatttcattc
cactctttccctacacgacgctcttccgatttcattcttcccc
cactctttccctacacgacgctcttccgattttatttc
cactctttccctacacgacgctcttccggatcta
cactctttccctacacgacgctcttccggatctaatt
cactctttccctacacgacgctcttccggatctaattc
cactctttccctacacgacgctcttccggatctaattca
cactctttccctacacgacgctcttccggatctagctt
cactctttccctacacgacgctcttccggttcta
cactctttccctacacgacgctttccgatcta
cactctttccctacacgacgctttccgatctaattc
cactctttccctacacgacgtcttccgatctaattctggaccatagtgcaatgt
cactctttccctacacgcgctcttccgatcta
cactctttccctacacgcgctcttccgatctaattcg
cactctttccctacacgcgctcttccgatctaattcgg
cactctttccctacacgcgctcttccgatctaattcggcgg
cactctttccctacacgcgctcttccgatctagct
cactctttccctaccgacgctcttccgatcta
cactctttccctacacgacg

我需要从“sequences.fasta”文件中找到并删除这些模式:

>seq01
cactctttccctacacgacgctcttccgWANTEDSEQUENCE
>seq01
cactctttccctacacgacgctcttccgaatctaWANTEDSEQUENCE
>seq03
cactctttccctacacgacgctcttccgaatctaattWANTEDSEQUENCE
>seq04
cactctttccctacacgacgctcttccgaatctaattaWANTEDSEQUENCE
>seq05
cactctttccctacacgcgctcttccgatctaattcggWANTEDSEQUENCE
>seq06
cactctttccctacacgcgctcttccgatctaattcggcggWANTEDSEQUENCE
>seq07
cactctttccctacacgcgctcttccgatctagctWANTEDSEQUENCE
>seq08
cactctttccctaccgacgctcttccgatctaWANTEDSEQUENCE

所以想要的输出应该是:

>seq01
WANTEDSEQUENCE
>seq02
WANTEDSEQUENCE
>seq03
WANTEDSEQUENCE
>seq04
WANTEDSEQUENCE
>seq05
WANTEDSEQUENCE
>seq06
WANTEDSEQUENCE
>seq07
WANTEDSEQUENCE
>seq08
WANTEDSEQUENCE

(只是为了示例,我使用了“WANTEDSEQUENCE”而不是真正的序列)

我尝试了以下方法(以及一些变体。我也尝试了 while read):

ADAPS=($(cat adapters.txt))
FASTA="sequences.fasta"


for ADAP in "${ADAPS[@]}";
do
    sed "s/${ADAP}//g" "${FASTA}" > output.fasta
done

但是我得到了这个:

>seq01
ctcttccgWANTEDSEQUENCE
>seq01
ctcttccgaatctaWANTEDSEQUENCE
>seq03
ctcttccgaatctaattWANTEDSEQUENCE
>seq04
ctcttccgaatctaattaWANTEDSEQUENCE
>seq05
cactctttccctacacgcgctcttccgatctaattcggWANTEDSEQUENCE
>seq06
cactctttccctacacgcgctcttccgatctaattcggcggWANTEDSEQUENCE
>seq07
cactctttccctacacgcgctcttccgatctagctWANTEDSEQUENCE
>seq08
cactctttccctaccgacgctcttccgatctaWANTEDSEQUENCE

我该如何解决这个问题?

Sort adapters.txt in reverse order by its line length,从其输出创建一个 sed 脚本,并将其与 bash 的命令替换 <(...) 一起使用,并应用第二个 sed它到 sequences.fasta:

sed -f <(awk '{ print length, [=10=] }' adapters.txt | sort -rn | cut -d" " -f2- | sed -E 's/(.*)/s|&||/') sequences.fasta

输出:

>seq01
WANTEDSEQUENCE
>seq01
WANTEDSEQUENCE
>seq03
WANTEDSEQUENCE
>seq04
WANTEDSEQUENCE
>seq05
WANTEDSEQUENCE
>seq06
WANTEDSEQUENCE
>seq07
WANTEDSEQUENCE
>seq08
WANTEDSEQUENCE

adapters.txt 的排序是必要的,因为它包含同一文件中其他字符串的子字符串。

多行和文件中的相同代码:

awk '{ print length, [=11=] }' adapters.txt | sort -rn | cut -d" " -f2- > adapters_sorted.txt
sed -E 's/(.*)/s|&||/' adapters_sorted.txt > sed.script
sed -f sed.script sequences.fasta

这是一个 POSIX awk 解决方案:

$ awk 'NR==FNR{seq[FNR]=[=10=]; x=FNR; next}
      {for(i=1; i<=x; i++) if ([=10=] ~ "^" seq[i]) {sub(seq[i],""); print [=10=]; next}
      print}
      ' <(awk '{ print length()"\t"[=10=]}' adapters.txt | sort -nr | cut -f2)  sequences.fasta
>seq01
WANTEDSEQUENCE
>seq01
WANTEDSEQUENCE
>seq03
WANTEDSEQUENCE
>seq04
WANTEDSEQUENCE
>seq05
WANTEDSEQUENCE
>seq06
WANTEDSEQUENCE
>seq07
WANTEDSEQUENCE
>seq08
WANTEDSEQUENCE

gawk您可以在其中对序列从长到短进行内部排序:

$ gawk 'BEGIN{PROCINFO["sorted_in"] = "@val_num_desc"}
      NR==FNR { seq[[=11=]] = length([=11=]); next }
      {for (e in seq)  if([=11=]~"^" e) {sub(e,""); print [=11=]; next}
      print}
      ' adapters.txt sequences.fasta
  # same output

使用 GNU awk sorted_in:

$ cat tst.awk
NR==FNR {
    adapters2lengths[] = length()
    next
}
!/^>/ {
    PROCINFO["sorted_in"] = "@val_num_desc"
    for (adapter in adapters2lengths) {
        if ( index([=10=],adapter) == 1 ) {
            [=10=] = substr([=10=],adapters2lengths[adapter]+1)
            break
        }
    }
}
{ print }

.

$ awk -f tst.awk adapters.txt sequences.fasta
>seq01
WANTEDSEQUENCE
>seq01
WANTEDSEQUENCE
>seq03
WANTEDSEQUENCE
>seq04
WANTEDSEQUENCE
>seq05
WANTEDSEQUENCE
>seq06
WANTEDSEQUENCE
>seq07
WANTEDSEQUENCE
>seq08
WANTEDSEQUENCE

这个和@dawg 的 gawk 解决方案在功能上的区别在于,这个解决方案只进行字符串比较,而他们的解决方案进行正则表达式比较——只有当你的“adapters.txt”文件包含正则表达式元字符时才重要,所有否则等于我只是更喜欢使用字符串,除非我需要正则表达式。

使用 GNU edbash

#!/usr/bin/env bash

ed -s sequences.fasta < <(
  printf '%s\n' '1,$-1s/$/\|/' '1,$j' 's/^/,s\//' 's/$/\/\//' '$a' ,p Q . ,p Q |
  ed -s adapters.txt
)