perl - 将字符附加到 multifasta 文件中特定行的末尾

Question

我有一个带有双端读取的 multifasta 文件，因此彼此相邻的那些是配对的（它们具有相同的读取名称）。我想在整个文件中分别将“/1”和“/2”附加到第一次和第二次读取。我不知道文件中有多少读取。这是文件的样子（为清楚起见，在读取之间添加空行）：

HWI-ST1018:1:1101:10007:34134#0 ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTT CCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG

HWI-ST1018:1:1101:10007:34134#0 GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAG GAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG

HWI-ST1018:1:1101:10016:6488#0 ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAG AAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA

HWI-ST1018:1:1101:10016:6488#0 TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTTAGTAACCACATGTCCAGACTCCTC CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

这就是我希望它出现的方式：

HWI-ST1018:1:1101:10007:34134#0/1 ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTT CCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG

HWI-ST1018:1:1101:10007:34134#0/2 GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAG GAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG

HWI-ST1018:1:1101:10016:6488#0/1 ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAG AAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA

HWI-ST1018:1:1101:10016:6488#0/2 TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTTAGTAACCACATGTCCAGACTCCTC CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

然后，我将对其进行 grep，删除“--”分隔符并将正向读取（即以“/1”结尾的那些）和反向读取（即以“/2”结尾的那些）保存在不同的文件中，如下所示：

grep -A 2 "/1" filename.fa | sed '/--/d' > reads_1.fa
grep -A 2 "/2" filename.fa | sed '/--/d' > reads_2.fa

我认为这可以用 sed 和 awk 来完成，但我还没有弄清楚怎么做。请帮忙。提前致谢。

score 2 · Accepted Answer

使用 sed 生成中间文件：

#!/bin/sed -f

1 {
    x
    s/^$/\\1/
    x
}
/^HWI/ {
    G
    s/\n//
    x
    y/12/21/
    x
}

在一行中：

sed -e '1{x;s/^$/\\1/;x};/^HWI/{G;s/\n//;x;y/12/21/;x}'

命令非常简单。第一对大括号中的命令在第一行执行，它们初始化保持空间（辅助缓冲区）\1。为此，我们使用x交换命令将模式空间（工作缓冲区）的内容与保持空间的内容交换。然后我们用替换一个空行\1，然后再次交换空格。

对以 .开头的每一行执行分组在第二对大括号中的命令HWI。首先，我们将保持空间的内容附加到模式空间中。因为它是以换行符开头附加的，所以下一个命令将其删除。现在我们必须将数字从 1 交换为 2，从 2 交换为 1。首先我们再次交换空格的内容，然后使用y命令交换字符。它定义当找到 1 或 2 时，它们必须分别替换为 2 或 1。最后，我们恢复空间的内容。

您还可以编写一个脚本来完成所有操作，将它们分成文件：

#!/bin/sed -f

/^HWI/! d

:start_forward
s/$/\\1/

:forward
w reads_1.fa
n
/^HWI/! b forward

s/$/\\2/

:reverse
w reads_2.fa
n
/^HWI/! b reverse

b start_forward

以较短的形式：

sed -e '/^HWI/!d;:s;s/$/\\1/;:f;wreads_1.fa
    n;/^HWI/!bf;s/$/\\2/;:r;wreads_2.fa
    n;/^HWI/!br;bs'

在这里，我们首先忽略所有行，直到找到以开头的行HWI。然后我们必须循环，一个用于写入前向数据，另一个用于后向数据。在循环之间有用于附加相应\1或\2在将行写入相应文件之前的命令。循环是类似的，它们简单地将行写入各自的文件，从输入中加载新行并检查它是否是以开头的行HWI，表明它应该进入下一个循环。

更彻底的解释：

第一个命令在该行不以开头时执行HWI（我们通过在其后添加一个来否定匹配!）。该命令是d删除一行，同时强制sed加载下一行并重新启动脚本。实际上，我们循环直到找到以 . 开头的字符串HWI。

现在，我们使用:命令来定义一个名为start_forward. 标签只不过是脚本中我们可以跳转到的位置的名称。如果我们一直在标签之间跳转并且永远不会到达脚本的末尾，我们最终将永远不会重新启动脚本，因此在HWI找到以开头的第一行之后，第一个命令永远不会被执行。我们要做的第一件事是将附加\1到行尾。

现在我们定义一个新的标签，叫做forward它，当我们遍历这些行时，它会被用来跳回来。reads_1.fa循环非常简单，首先我们使用命令将当前行写入相应的文件w，然后我们使用该n行将下一行读入模式空间，最后我们检查新读取的行是否以HWI. 如果没有，我们执行b分支命令跳回forward标签，允许我们开始循环的另一个迭代。

如果这些行确实以开头HWI，我们现在必须转到另一个循环。不过，在此之前，我们必须在该行后面加上\2. 该循环类似于前一个循环，除了当我们在HWI找到另一行时退出循环时，我们必须使用命令分支回到start_forward标签b才能切换回前一个循环。

希望这会有所帮助=)

score 1 · Accepted Answer

1

这增加/0了/1以下内容：

perl -pe 'if (/#0/) { $x = 1 - $x; s:#0:#0/$x: }'

于 2012-10-09T09:32:56.503 回答

score 1 · Accepted Answer

awk 单线：

 awk -F'#' 'NF==2{a[$1]=($1 in a)?++a[$1]:1;$0=$0"/"a[$1];}1' file

测试

kent$  cat t.txt
HWI-ST1018:1:1101:10007:34134#0
ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTT 
CCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG

HWI-ST1018:1:1101:10007:34134#0
GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAG 
GAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG

HWI-ST1018:1:1101:10016:6488#0
ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAG 
AAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA

HWI-ST1018:1:1101:10016:6488#0
TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTC 
CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

kent$  awk -F'#' 'NF==2{a[$1]=($1 in a)?++a[$1]:1;$0=$0"/"a[$1];}1' t.txt
HWI-ST1018:1:1101:10007:34134#0/1
ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTT 
CCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG

HWI-ST1018:1:1101:10007:34134#0/2
GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAG 
GAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG

HWI-ST1018:1:1101:10016:6488#0/1
ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAG 
AAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA

HWI-ST1018:1:1101:10016:6488#0/2
TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTC 
CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

score 1 · Accepted Answer

awk 'BEGIN{i=1}{if($0~/#0/){print $0"/"i;if(i==1)i=2;else i=1;}else {print}}' your_file

测试如下：

> cat temp
>HWI-ST1018:1:1101:10007:34134#0
ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTT
CCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG
>HWI-ST1018:1:1101:10007:34134#0
GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAG
GAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG
>HWI-ST1018:1:1101:10016:6488#0
ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAG
AAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA
>HWI-ST1018:1:1101:10016:6488#0
TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTC
CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

执行：

> awk 'BEGIN{i=1}{if($0~/#0/){print $0"/"i;if(i==1)i=2;else i=1;}else {print}}' temp
>HWI-ST1018:1:1101:10007:34134#0/1
ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTT
CCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG
>HWI-ST1018:1:1101:10007:34134#0/2
GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAG
GAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG
>HWI-ST1018:1:1101:10016:6488#0/1
ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAG
AAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA
>HWI-ST1018:1:1101:10016:6488#0/2
TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTC
CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC
>

score 1 · Accepted Answer

你有所谓的 shuffled multifasta。您可以使用GNU awkunshuffle 并创建两个文件。请注意，无需使用grep或sed执行任何后处理。此代码将为您创建两个文件：

awk 'NR%4==1 { getline one; printf "%s/1\n%s\n", $0, one > "reads_1.fa" } NR%4==3 { getline two; printf "%s/2\n%s\n", $0, two > "reads_2.fa" }' file.txt

输入：

HWI-ST1018:1:1101:10007:34134#0
ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG
HWI-ST1018:1:1101:10007:34134#0
GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAGGAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG
HWI-ST1018:1:1101:10016:6488#0
ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAGAAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA
HWI-ST1018:1:1101:10016:6488#0
TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

结果：

内容reads_1.fa：

HWI-ST1018:1:1101:10007:34134#0/1
ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG
HWI-ST1018:1:1101:10016:6488#0/1
ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAGAAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA

内容reads_2.fa：

HWI-ST1018:1:1101:10007:34134#0/2
GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAGGAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG
HWI-ST1018:1:1101:10016:6488#0/2
TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

score 1 · Accepted Answer

另一种解决方案：

awk 'BEGIN{RS=""}{if(NR<3){sub(/#0/,"#0/"NR);print $0,"\n"}else{NR=1;sub(/#0/,"#0/"NR);print $0,"\n"}}' file

结果：

awk 'BEGIN{RS=""}{if(NR<3){sub(/#0/,"#0/"NR); print $0,"\n"}else{NR=1;sub(/#0/,"#0/"NR);print $0, "\n"}}' file 
HWI-ST1018:1:1101:10007:34134#0/1 
ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTT 
CCCATTCCCCTAGGGCTGAGACCCAATATCCTCTATCCCTG 

HWI-ST1018:1:1101:10007:34134#0/2 
GTGCAGGCATGTTGGGGCGTGTCTCAGAGCCTGAACTTCCCTTCCAGTCAGTGCTGGAAG 
GAGGTGGGCAGGGGAATGATAGAAAGGAAGGAGTGGATTGG 

HWI-ST1018:1:1101:10016:6488#0/1 
ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTTGCATAGTCTCATACACGTATCAG 
AAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA 

HWI-ST1018:1:1101:10016:6488#0/2 
TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTC 
CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

perl - 将字符附加到 multifasta 文件中特定行的末尾

6 回答 6

Related

Reference