0

我有几个非常大的 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

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

4

1 回答 1

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
于 2015-07-31T22:39:24.133 回答