☰

シーケンスの長さに基づいて fastq ファイルをサブセット化するために awk を使用する

tag-icon • tag-icon awk
シーケンスの長さに基づいて fastq ファイルをサブセット化するために awk を使用する

fastqファイルがあります。それが何であるかを説明します。これは次のようなものです

@SRR1024120.7 DBRHHJN1:259:D0PM7ACXX:1:1101:1386:1189 length=100
GATACAGGATGCCTGGGTCTAGGCTGTGTGACCTTGGGCCAGTTCCTCTC
+SRR1024120.7 DBRHHJN1:259:D0PM7ACXX:1:1101:1386:1189 length=100
DDDFFDDBGFEHEHGIGC9F>HG9EH8?DF4?:DF<?3:D?DHIGGDDFH
@SRR1024120.25 DBRHHJN1:259:D0PM7ACXX:1:1101:1752:1149 length=100
CTGCTGCTCATGCTCAT
+SRR1024120.25 DBRHHJN1:259:D0PM7ACXX:1:1101:1752:1149 length=100
BDDDDD<<CC:C+AFFE
@SRR1024120.42 DBRHHJN1:259:D0PM7ACXX:1:1101:2482:1096 length=100
AGCGTGTGCCACCCTACGCCGGC
+SRR1024120.42 DBRHHJN1:259:D0PM7ACXX:1:1101:2482:1096 length=100
DD>DAA@AA@@?2C8AB)?@:DD
@SRR1024120.1 DBRHHJN1:259:D0PM7ACXX:1:1101:1200:1120 length=100
AGACAGAAGGGGAGTACAGCTCTCTGGAACATGAGAGTGCAAGGGGTTGAGTGTTT
+SRR1024120.1 DBRHHJN1:259:D0PM7ACXX:1:1101:1200:1120 length=100
DDDFFFCFGEHI@CGFADFGCCFFGHFGCFFFHGGDGHIFHDFGGI<BF=DHIHHH

4行が1回の読み取りに対応するので

@SRR1024120.7 DBRHHJN1:259:D0PM7ACXX:1:1101:1386:1189 length=100
GATACAGGATGCCTGGGTCTAGGCTGTGTGACCTTGGGCCAGTTCCTCTC
+SRR1024120.7 DBRHHJN1:259:D0PM7ACXX:1:1101:1386:1189 length=100
DDDFFDDBGFEHEHGIGC9F>HG9EH8?DF4?:DF<?3:D?DHIGGDDFH

1回の読み取りに相当し、ガタカAGGATGCCTGGGTCTAGGCTGTGTGACCTTGGGCCAGTTCCTCTC

上でfastqファイルを示しました。私がやりたいのは、リードシーケンスの長さが25以下のリードだけを抽出したいということです。そのため、出力は次のようになります。

@SRR1024120.25 DBRHHJN1:259:D0PM7ACXX:1:1101:1752:1149 length=100
CTGCTGCTCATGCTCAT
+SRR1024120.25 DBRHHJN1:259:D0PM7ACXX:1:1101:1752:1149 length=100
BDDDDD<<CC:C+AFFE
@SRR1024120.42 DBRHHJN1:259:D0PM7ACXX:1:1101:2482:1096 length=100
AGCGTGTGCCACCCTACGCCGGC
+SRR1024120.42 DBRHHJN1:259:D0PM7ACXX:1:1101:2482:1096 length=100
DD>DAA@AA@@?2C8AB)?@:DD

この目的のために awk を使用したいと思います。

こんなことを試してみた

awk 'NR % 2 == 0 {if(length($1) <= 25) print $0}; NR % 2 == 1' test.fastq

しかし、これは次のように印刷されます

@SRR1024120.7 DBRHHJN1:259:D0PM7ACXX:1:1101:1386:1189 length=100
+SRR1024120.7 DBRHHJN1:259:D0PM7ACXX:1:1101:1386:1189 length=100
@SRR1024120.25 DBRHHJN1:259:D0PM7ACXX:1:1101:1752:1149 length=100
CTGCTGCTCATGCTCAT
+SRR1024120.25 DBRHHJN1:259:D0PM7ACXX:1:1101:1752:1149 length=100
BDDDDD<<CC:C+AFFE
@SRR1024120.42 DBRHHJN1:259:D0PM7ACXX:1:1101:2482:1096 length=100
AGCGTGTGCCACCCTACGCCGGC
+SRR1024120.42 DBRHHJN1:259:D0PM7ACXX:1:1101:2482:1096 length=100
DD>DAA@AA@@?2C8AB)?@:DD
@SRR1024120.1 DBRHHJN1:259:D0PM7ACXX:1:1101:1200:1120 length=100
+SRR1024120.1 DBRHHJN1:259:D0PM7ACXX:1:1101:1200:1120 length=100

明らかに私は望んでいない

@SRR1024120.7 DBRHHJN1:259:D0PM7ACXX:1:1101:1386:1189 length=100
+SRR1024120.7 DBRHHJN1:259:D0PM7ACXX:1:1101:1386:1189 length=100
@SRR1024120.1 DBRHHJN1:259:D0PM7ACXX:1:1101:1200:1120 length=100
+SRR1024120.1 DBRHHJN1:259:D0PM7ACXX:1:1101:1200:1120 length=100

私の出力では。

ご協力いただければ幸いです

ありがとう

答え1

レコードとフィールドの区切り文字を使用すると、awk が改行とスペースを使用しないようにすることができます。代わりに、レコードの区切り文字として「\n@」を使用し、フィールドの区切り文字として「\n」を使用することができます。

$ awk 'BEGIN {RS="\n@";FS="\n"} {if (length($2) <= 25) {print "@"$0} }' fastq
@SRR1024120.25 DBRHHJN1:259:D0PM7ACXX:1:1101:1752:1149 length=100
CTGCTGCTCATGCTCAT
+SRR1024120.25 DBRHHJN1:259:D0PM7ACXX:1:1101:1752:1149 length=100
BDDDDD<<CC:C+AFFE
@SRR1024120.42 DBRHHJN1:259:D0PM7ACXX:1:1101:2482:1096 length=100
AGCGTGTGCCACCCTACGCCGGC
+SRR1024120.42 DBRHHJN1:259:D0PM7ACXX:1:1101:2482:1096 length=100
DD>DAA@AA@@?2C8AB)?@:DD

各行は異なるフィールドになるので、$2 で 2 行目の長さを確認できます。レコード区切り文字によって「@」が消えてしまうため、印刷時に「@」を再度追加する必要がありました。

関連情報