按长度有效分割 fastq 文件

Question

我正在尝试找到一种更省时的方法来按序列长度拆分 fastq 文件，即将一个大的 fastq 文件拆分为多个仅包含相同长度序列的文件。输入是一个普通的 fastq 文件（每个序列 4 行，实际序列在每个四重奏的第二行），具有不同的序列长度：

@HISEQ:28:H8P69ADXX:1:1101:1462:2036 1:N:0:CTTGTA
NCCATAAAGTAGAAAGCACT
+
#00<FFFFFFFFFIIFIIFF
@HISEQ:28:H8P69ADXX:1:1101:1419:2156 1:N:0:CTTGTA
TGGAGAGAAAGGCAGTTCCTGA
+
BBBFFFFFFFFFFIIIIIIIII
@HISEQ:28:H8P69ADXX:1:1101:1378:2223 1:N:0:CTTGTA
TCCTGTACTGAGCTGCCCCGA
+
BBBFFFFFFFFFFIIIIIIII
@HISEQ:28:H8P69ADXX:1:1101:1585:2081 1:N:0:CTTGTA
AAACCGTTACCATTACTGAGT
+
BBBFFFFFFFFFFIIIIFIII

现在我正在使用 awk 过滤掉特定长度或特定范围内的序列：

awk 'BEGIN {OFS = "\n"} {header = [=11=] ; getline seq ; getline qheader ; getline qseq ; if (length(seq) == 22) {print header, seq, qheader, qseq}}'

如果我希望每个序列长度都有一个输出文件，我可以使用 for 循环来管理：

for i in {16..33};
awk -v var=$i 'BEGIN {OFS = "\n"} {header = [=12=] ; getline seq ; getline qheader ; getline qseq ; if (length(seq) == var) {print header, seq, qheader, qseq}}'
done

问题是，虽然它工作正常，但相当耗时，因为我猜我要分别检查整个文件的每个长度。另外我需要事先检查最长和最短的序列。

谁能帮我找到比我的循环更有效的解决方案？如果可能的话，我不必指定一个范围的解决方案，而是一个检查最小和最大长度并自动拆分它们的解决方案。我想用 awk 来做，但我对一切都持开放态度。谢谢本尼迪克特

Answer 1

是这样的吗？

$ awk        '{rec=rec sep [=10=]; sep=ORS} 
       !(NR%4){print rec > fn; rec=sep=""} 
       NR%4==2{fn = length([=10=])".seq"}' file

将生成内容为

的这3个文件

==> 20.seq <==
@HISEQ:28:H8P69ADXX:1:1101:1462:2036 1:N:0:CTTGTA
NCCATAAAGTAGAAAGCACT
+
#00<FFFFFFFFFIIFIIFF

==> 21.seq <==
@HISEQ:28:H8P69ADXX:1:1101:1378:2223 1:N:0:CTTGTA
TCCTGTACTGAGCTGCCCCGA
+
BBBFFFFFFFFFFIIIIIIII
@HISEQ:28:H8P69ADXX:1:1101:1585:2081 1:N:0:CTTGTA
AAACCGTTACCATTACTGAGT
+
BBBFFFFFFFFFFIIIIFIII

==> 22.seq <==
@HISEQ:28:H8P69ADXX:1:1101:1419:2156 1:N:0:CTTGTA
TGGAGAGAAAGGCAGTTCCTGA
+
BBBFFFFFFFFFFIIIIIIIII

由于会有少量这些输出文件，因此无需明确关闭它们。

说明

{rec=rec sep [=12=]; sep=ORS} build the record line by line with ORS in between lines, with lazy initialization of the separator we can eliminate the dangling first separator.

!(NR%4) if the line number is a multiple of 4

{print rec > fn; rec=sep=""} print the record to the file and reset record and separator

NR%4==2 if the line number is a 2 of 4.

{fn = length([=16=])".seq"} set the filename

按长度有效分割 fastq 文件

efficient splitting of fastq files by length

awk

fastq