将字符附加到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 ACAGCTATACACGAAGAATCTCAGCCCTTGTACTTTT GCATAGTCTCATACACGTATCAG AAGCCTCCACCTGGCTAACAGGAATTTGGGGCTTTGGGAGA

HWI-ST1018：1：1101：10016：6488＃0 TTTGGGAGATTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTC CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

这是我多么希望它出现：

HWI-ST1018： 1：1101：10007：34134＃0/1 ACTAGTAACCACATGTCCAGACTCCTCCTATGCTCCCACCCAGGGTCCCTTGAGCTGCTT CCCATTCCCCTAGGGCTGAGACCCA ATATCCTCTATCCCTG

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

然后，我会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做，但我还没有想出如何。请帮忙。提前致谢。

Answer 1

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 
> 

Answer 2

这增加/0和/1到读取：

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

Answer 3

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 

Answer 4

你有什么所谓改组multifasta。您可以使用GNU awk取消清除并创建两个文件。请注意，不需要使用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 

Answer 5

使用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。最后，我们还原空间的内容。忽略所有行，直到我们找到

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

下面我们开始：

你也可以写一个脚本做的一切，将它们分成文件：

#!/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 

在一个较短的形式一条以HWI开头的行。然后我们必须循环，一个用于写入正向数据，另一个用于反向数据。在循环之间有命令用于在将行写入相应文件之前追加相应的\1或\2。循环是类似的，它们简单地将行写入到它们各自的文件中，从输入中加载新行并检查它是否是以HWI开始的行，指示它应该进入下一个循环。

更透彻的解释：

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

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

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

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

希望这有助于=）

Answer 6

另一种解决方案：

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