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

【问题标题】：sort fastq file and keep sequences 15-17 bp in length对 fastq 文件进行排序并保持序列长度为 15-17 bp
【发布时间】：2020-11-20 04:44:15
【问题描述】：

我有几个非常大的 fastq 文件，我正在使用 cutadapt 从其中修剪转座子末端序列，这应该会导致剩余 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，但似乎从未找到正确的解决方案。有没有人有任何解决方案？

【问题讨论】：

如果bioawk 可以选择以fastq 格式输出，则您链接的问题的答案看起来应该可以解决问题。 OP 从未发布示例输入/输出，或将其标记为已接受。 ://
链接的答案没有回答问题，因为结果不是 FASTQ 格式。我建议避免使用框架，除非它严格遵循规范。

标签： bioinformatics fastq

【解决方案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

【讨论】：