删除具有低质量碱基调用的读数

Question

我有一个 fastq 格式的序列数据文件 https://en.wikipedia.org/wiki/FASTQ_format 其中第一行是序列 ID，第二行是序列 [ACGT]，第三行是“+”，第四行是质量值。

@M01610:118:000000000-D49F3:1:1101:14523:2546 1:N:0:CTTGTA
GTACACCTTCATGAAGAACTCCATCACCTTCATCTCCAGGATGCGGTCCTGGGTGCTGTTCCTGGCGATCTCGATCAGCTCGATGTACTCGTGGGGCACGTACTTCAGCTTGTGCCGCAGCTCGGACTTCTTCTCCTCCAGCTCGCTCTTCACCAGCTGGGATCCCCGCAGGTGTATCTTGGTATGCTTGTTCAGGTTGGAGCGGTGGGCAAATTTCCTCCCACAAATGTCACAGGCAAAAGGCTTCTC
+
CCCCCFFFFFFFGGGGGGGGGGHHHHHHHHHHHHHHGHHHGHHHGGGGGHHHGFGGHHHHHHHFHHGGGGHHHGGHGHHHHGGHHHHHHHGHGGGGGGGGHHHHHHHGHHHHHHHHGGGGGGHGGGGGHHHHHHHHHHHHHHHHGGGHHHHHHHHFHHHGEHFGHHGGGGGGGHGHFHHHHHFFHHGGGGGGGGGGFFF?FGGGGFGGGFFFFFFFFFFFFEFFF?FFFFFFEFFEFFFFBFFFFFBFF
@M01610:118:000000000-D49F3:1:1101:9569:5713 1:N:0:CTTGTA
CAAGGAAGGCACGGGGGAGGGGCAAACAACAGATGGCTGGCAACTAGAAGGCACAGGCTAGCCAGGCGGGGAGGCGGCCCAAAGGGAGATCCGACTCGTCGGAGGCCGAAAGCGAAGACGCGGGAGAGGCCGCAGAACCGGCAGAAGGCCTCGGGAAGGGAGGTCCGCTGGATTGAGAGCCGAAGGGACGTAGCAGAAGGACGTCCCGCGCAGGATCCAGTTGGCAACACAGGCGAGCAGCCAAGG
+
CDCCDFFFDCFFGGGGGGGGGGGGGHHHHHHHHHHHHHGGGHHHHHHHHHHHHHHHHGHHHHHHHGHGGGGGCGGFGGGGDHHHHGHGGHHHHGGGGFHGFGAGGGGGAAGFFDBF-DFFF>DF;DFAFDF=CA>CFBE>FFCFEFBFFF0FDDFAFFFFEDC.BFFFDBF.FFEBFFFEFAAC=FFE?>AEFEBFBFFFFFFDFFFFC>-9>=ABFFFFBFFFFFFFFFEFFFCFFA9BBEAFEF

我想删除第 4 行中任何字符与任一符号都不匹配的所有条目 ?@ABCDEFGHIJK

以上示例的输出将是

@M01610:118:000000000-D49F3:1:1101:14523:2546 1:N:0:CTTGTA
GTACACCTTCATGAAGAACTCCATCACCTTCATCTCCAGGATGCGGTCCTGGGTGCTGTTCCTGGCGATCTCGATCAGCTCGATGTACTCGTGGGGCACGTACTTCAGCTTGTGCCGCAGCTCGGACTTCTTCTCCTCCAGCTCGCTCTTCACCAGCTGGGATCCCCGCAGGTGTATCTTGGTATGCTTGTTCAGGTTGGAGCGGTGGGCAAATTTCCTCCCACAAATGTCACAGGCAAAAGGCTTCTC
+
CCCCCFFFFFFFGGGGGGGGGGHHHHHHHHHHHHHHGHHHGHHHGGGGGHHHGFGGHHHHHHHFHHGGGGHHHGGHGHHHHGGHHHHHHHGHGGGGGGGGHHHHHHHGHHHHHHHHGGGGGGHGGGGGHHHHHHHHHHHHHHHHGGGHHHHHHHHFHHHGEHFGHHGGGGGGGHGHFHHHHHFFHHGGGGGGGGGGFFF?FGGGGFGGGFFFFFFFFFFFFEFFF?FFFFFFEFFEFFFFBFFFFFBFF

此处，seq ID @M01610:118:000000000-D49F3:1:1101:9569:5713 1:N:0:CTTGTA 的第 4 行包含 ?@ABCDEFGHIJK 以外的符号，因此被删除。

seq ID 的第 2 行和第 4 行的长度相同，但对于范围从 (200 到 250) 的不同 seqID 不同。

一个单元由4行组成(seq ID, sequence, + sing, quality for sequence)。我想删除所有单元，其中序列（第 2 行）中每个字符的序列质量（第 4 行）与模式中任一字符以外的任何字符匹配？@ABCDEFGHIJK。我已经尝试过这段代码，并且仍在努力

cat file.fq | awk 'NR%4==0' | xargs -n1 awk '{ for(i=0; ++i <= length([=12=]);) printf "%s\n" }'

如有任何帮助，我将不胜感激。

Answer 1

为缓冲区中的每个单元收集行；在每个新的 header 行处理前一个单元（检查缓冲区中的最后一行，是否打印）并清空缓冲区

use warnings;
use strict;
use feature 'say';

sub process_unit {
    my ($rbuf) = @_;
    if (not $rbuf->[-1] =~ tr/?@ABCDEFGHIJK//dr ) { #/ no extra chars
        say for @$rbuf;
    }
}

my $file = shift // die "Usage: [=10=] filename\n";   #/    
open my $fh, '<', $file or die "Can't open $file: $!";

my @buf;
while (<$fh>) {
    chomp;
    if (/^\@/ and @buf) {
        process_unit(\@buf);
        @buf = (); 
    }   

    push @buf, $_; 
}
process_unit(\@buf);  # the last unit

使用tr检查缓冲区最后一行的说明（记录 in perlop）：

tr/..//dr 从其目标字符串中删除所有列出的字符，返回修改后的字符串，同时保持原始字符串不变（由于 "non-destructive" /r 修饰符）。因此，如果在删除允许的字符后还剩下任何东西，我们将丢弃该单元（不要打印它）。

关于选择tr和效率的说明

可以将正则表达式及其匹配运算符与否定字符一起使用 class,

if (not /[^?\@ABCDEFGHIJK]/) { ... }

而不是音译运算符 tr（不是正则表达式）。

但是，即使在匹配运算符的最佳情况下，我也对 tr 方法进行了基准测试，使其 25% 更快。在所有其他情况下，tr 至少比正则表达式的匹配好 2-4 倍。

匹配运算符的 "best case" 是当字符串中的第一个位置出现不可接受的字符时，它会立即匹配并且不会扫描字符串的其余部分。至少可以说，这是相当不现实的。相反，统计上许多（大多数？）字符串不会包含任何这些字符，并且匹配运算符将出现最坏的情况，即扫描整个字符串。请注意，在所有这些中，正则表达式引擎的启动可能是最昂贵的部分。

Answer 2

$ awk '{unit=unit [=10=] ORS} NR%4==0{if (/^[?@ABCDEFGHIJK]+$/) printf "%s", unit; unit=""}' file
@M01610:118:000000000-D49F3:1:1101:14523:2546 1:N:0:CTTGTA
GTACACCTTCATGAAGAACTCCATCACCTTCATCTCCAGGATGCGGTCCTGGGTGCTGTTCCTGGCGATCTCGATCAGCTCGATGTACTCGTGGGGCACGTACTTCAGCTTGTGCCGCAGCTCGGACTTCTTCTCCTCCAGCTCGCTCTTCACCAGCTGGGATCCCCGCAGGTGTATCTTGGTATGCTTGTTCAGGTTGGAGCGGTGGGCAAATTTCCTCCCACAAATGTCACAGGCAAAAGGCTTCTC
+
CCCCCFFFFFFFGGGGGGGGGGHHHHHHHHHHHHHHGHHHGHHHGGGGGHHHGFGGHHHHHHHFHHGGGGHHHGGHGHHHHGGHHHHHHHGHGGGGGGGGHHHHHHHGHHHHHHHHGGGGGGHGGGGGHHHHHHHHHHHHHHHHGGGHHHHHHHHFHHHGEHFGHHGGGGGGGHGHFHHHHHFFHHGGGGGGGGGGFFF?FGGGGFGGGFFFFFFFFFFFFEFFF?FFFFFFEFFEFFFFBFFFFFBFF

Answer 3

这是一种使用 (1) 修改后的输入记录保存程序和 (2) 使用 tr///c c 恭维开关修改的音译运算符的方法。

（我模拟了脚本顶部的文件）

#!/usr/bin/perl
use strict;
use warnings;

my $file =<<'EOF';
@M01610:118:000000000-D49F3:1:1101:14523:2546 1:N:0:CTTGTA
GTACACCTTCATGAAGAACTCCATCACCTTCATCTCCAGGATGCGGTCCTGGGTGCTGTTCCTGGCGATCTCGATCAGCTCGATGTACTCGTGGGGCACGTACTTCAGCTTGTGCCGCAGCTCGGACTTCTTCTCCTCCAGCTCGCTCTTCACCAGCTGGGATCCCCGCAGGTGTATCTTGGTATGCTTGTTCAGGTTGGAGCGGTGGGCAAATTTCCTCCCACAAATGTCACAGGCAAAAGGCTTCTC
+
CCCCCFFFFFFFGGGGGGGGGGHHHHHHHHHHHHHHGHHHGHHHGGGGGHHHGFGGHHHHHHHFHHGGGGHHHGGHGHHHHGGHHHHHHHGHGGGGGGGGHHHHHHHGHHHHHHHHGGGGGGHGGGGGHHHHHHHHHHHHHHHHGGGHHHHHHHHFHHHGEHFGHHGGGGGGGHGHFHHHHHFFHHGGGGGGGGGGFFF?FGGGGFGGGFFFFFFFFFFFFEFFF?FFFFFFEFFEFFFFBFFFFFBFF
@M01610:118:000000000-D49F3:1:1101:9569:5713 1:N:0:CTTGTA
CAAGGAAGGCACGGGGGAGGGGCAAACAACAGATGGCTGGCAACTAGAAGGCACAGGCTAGCCAGGCGGGGAGGCGGCCCAAAGGGAGATCCGACTCGTCGGAGGCCGAAAGCGAAGACGCGGGAGAGGCCGCAGAACCGGCAGAAGGCCTCGGGAAGGGAGGTCCGCTGGATTGAGAGCCGAAGGGACGTAGCAGAAGGACGTCCCGCGCAGGATCCAGTTGGCAACACAGGCGAGCAGCCAAGG
+
CDCCDFFFDCFFGGGGGGGGGGGGGHHHHHHHHHHHHHGGGHHHHHHHHHHHHHHHHGHHHHHHHGHGGGGGCGGFGGGGDHHHHGHGGHHHHGGGGFHGFGAGGGGGAAGFFDBF-DFFF>DF;DFAFDF=CA>CFBE>FFCFEFBFFF0FDDFAFFFFEDC.BFFFDBF.FFEBFFFEFAAC=FFE?>AEFEBFBFFFFFFDFFFFC>-9>=ABFFFFBFFFFFFFFFEFFFCFFA9BBEAFEF
EOF

{
    local $/ = '@';     # set input record separator in this scope to '@'
    open my $fh, '<', $file;
    <$fh>;              # discard first read (will only contain '@')

    while (<$fh>) {
        chomp;
        my ($test) = /\+\n^(.+)$/m; # grab the fourth line

        # print record (with leading @ prepended back to beginning of record)
        #  unless there  are unwanted characters
        print "\@$_" unless $test =~ tr/?@ABCDEFGHIJK//c;
    }   
}

Answer 4

如果您的数据在 'd' 文件中，可以查看：

 perl -ne 'if (/@M0\w*:\d\d\d:/) {$s=$_;$s.=<> for 1..2;$r=<>; if ($r =~ /[^\s?A-K]/) {next} else {$s.=$r;print $s}}' d

删除具有低质量碱基调用的读数

Remove reads that have low quality base calls

linux

perl

awk

fastq