对 fastq 文件进行排序并保持序列长度为 15-17 bp

Question

我有几个非常大的 fastq 文件，我正在使用 cutadapt trim 关闭转座子末端序列，这应该会导致剩余 15-17 个碱基对的基因组 DNA。使用 cutadapt 后，fastq 文件的很大一部分是 15-17 个碱基对，但有些序列要长很多（表明它们上面没有转座子末端序列，它们是我实验的垃圾读取）。

我的问题：在 Linux 中是否有我可以使用的命令或脚本，以便我对这些 fastq 文件进行排序并输出一个新的 fastq，其中仅包含 15-17 个碱基对长的读取，同时仍然保留通常的 fastq 格式？

作为参考，fastq 格式如下所示：

@D64TDFP1:287:C69APACXX:2:1101:1319:2224 1:N:0:
GTTAGACCGGATCCTAACAGGTTGGATGATAAGTCCCCGGTCTAT
+
DDHHHDHHGIHIIIIE?FFHECGHICHHGH>BD?GHIIIIFHIDG
@D64TDFP1:287:C69APACXX:2:1101:1761:2218 1:N:0:
GTTAGACCGGATCCTAACAGGTTGGATGATAAGTCCCCGGTCTAT
+
FFHHHHHJIJJJJJIIJJJIJHIJJGIJIIIFJ?HHJJJJGHIGI

我发现了一个类似的问题 here，但似乎从未找到正确的解决方案。有没有人有任何解决方案？

Answer 1

一次将四行读入一个数组。当读取长度在您的阈值之间时，打印出这四行。

这里有一个如何用 Perl 实现的例子，但是在 Python 或任何其他脚本语言中原理是相同的：

#!/usr/bin/env perl

use strict;
use warnings;

my $fastq;
my $lineIdx = 0;
while (<>) {
    chomp;
    $fastq->[$lineIdx++] = $_;
    if ($lineIdx == 4) {
        my $readLength = length $fastq->[1];
        if (($readLength >= 15) && ($readLength <= 17)) {
            print "$fastq->[0]\n$fastq->[1]\n$fastq->[2]\n$fastq->[3]\n";
        }
        $lineIdx = 0;
    }
}

使用，例如：

$ ./filter_fastq_reads.pl < reads.fq > filtered_reads.fq

这会按照找到的顺序打印出读数。这只是过滤，应该很快。否则，如果您需要根据某些条件进行排序，请在您的问题中指定排序条件。

在Python中：

#!/usr/bin/env python

import sys

line_idx = 0
record = []
for line in sys.stdin:
    record[line_idx] = line.rstrip()
    line_idx += 1
    if line_idx == 4:
        read_length = len(record[1])
        if read_length >= 15 and read_length <= 17:
            sys.stdout.write('{}\n'.format('\n'.join(record)))
        line_idx = 0

对 fastq 文件进行排序并保持序列长度为 15-17 bp

sort fastq file and keep sequences 15-17 bp in length

bioinformatics

fastq