perl - multifasta ファイルの特定の行の末尾に文字を追加する

Question

互いに隣接するものが仲間になるように、ペアエンド読み取りを含むmultifastaファイルがあります（それらは同じ読み取り名を持っています）。ファイル全体で、最初と 2 番目の読み取りにそれぞれ「/1」と「/2」を追加したいと考えています。ファイル内の読み取り回数がわかりません。ファイルは次のようになります (読みやすくするために読み取りの間に空白行を追加しています)。

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

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 TTTGGGAGATTTTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTC CTATGCTCCCACCCAGGGTCCCTTGAGCTGCTTCCCATTCC

そして、これは私がそれをどのように表示したいかです：

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

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 TTTGGGAGATTTTTAATCAGGGCAAAACCTGTACTAGTAACCACATGTCCAGACTCCTC CTATGCTCCCACCCAGGGGTCCCTTGAGCTGCTTCCCATTCC

次にそれを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

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ます。これを行うには、exchange コマンドを使用して、パターン スペース (ワーキング バッファー) の内容をホールド スペースの内容とx交換します。次に、空の行を に置き換えてから\1、スペースを再度交換します。

中括弧の 2 番目のペアでグループ化されたコマンドは、 で始まるすべての行に対して実行されます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使用して現在の行をそれぞれのファイルに書き込み、次にその行を使用して次の行をパターン空間に読み取り、最後に新しく読み取った行が で始まるかどうかを確認します。そうでない場合は、分岐コマンドを実行してラベルに戻り、ループの別の反復を開始できるようにします。wnHWIbforward

行が で始まる場合はHWI、もう一方のループに移動する必要があります。ただし、その前に、行に を追加する必要があります\2。ループは前のループに似ていますが、別のHWI行が見つかったときにループを終了するとき、前のループに戻るには、コマンドでstart_forwardラベルに戻る必要があります。b

これが役立つことを願っています=）

Answer 2

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 3

これにより、読み取りに/0およびが追加されます。/1

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

Answer 4

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 5

別の解決策：

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 

Answer 6

シャッフル マルチファスタと呼ばれるものがあります。を使用GNU awkしてシャッフルを解除し、2 つのファイルを作成できます。後処理を使用grepまたは実行する必要がないことに注意してください。sedこのコードにより、次の 2 つのファイルが作成されます。

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