按长度有效分割 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 = $0 ; 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 = $0 ; getline seq ; getline qheader ; getline qseq ; if (length(seq) == var) {print header, seq, qheader, qseq}}'
done

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

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

Answer 1

这样的？

$ awk        '{rec=rec sep $0; sep=ORS} 
       !(NR%4){print rec > fn; rec=sep=""} 
       NR%4==2{fn = length($0)".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 $0; sep=ORS} 逐行构建记录，行间使用 ORS，通过延迟初始化分隔符我们可以消除悬空的第一个分隔符。

!(NR%4)如果行号是4的倍数

{print rec > fn; rec=sep=""}将记录打印到文件并重置记录和分隔符

NR%4==2 如果行号是 2 of 4。

{fn = length($0)".seq"}设置文件名