Эффективное разбиение файлов 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 из 4.

{fn = length($0)".seq"} установить имя файла