Удалить чтения, которые имеют низкое качество базовых вызовов - PullRequest
1 голос
/ 01 апреля 2019

У меня есть файл данных последовательности в формате fastq https://en.wikipedia.org/wiki/FASTQ_format, где первая строка - это идентификатор последовательности, вторая строка - это последовательность [ACGT], 3-я строка - это «+», а 4-я строка - это значения качества.

@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

Здесь 4-я строка seq ID @M01610:118:000000000-D49F3:1:1101:9569:5713 1:N:0:CTTGTA содержит символы, отличные от ?@ABCDEFGHIJK, и поэтому удаляется.

Длина 2-й и 4-й строки для идентификатора seq одинакова, но отличается для разных seqID в диапазоне (от 200 до 250).

Единица состоит из 4 строк (seq ID, sequence, + sing, качество для последовательности). Я хочу удалить все единицы, в которых качество последовательности (4-я строка) для каждого символа в последовательности (2-я строка) соответствует чему-либо, кроме символов шаблона? @ABCDEFGHIJK. Я пробовал этот код и все еще работаю над ним

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

Буду признателен за любую помощь.

Ответы [ 4 ]

3 голосов
/ 01 апреля 2019
$ awk '{unit=unit $0 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
3 голосов
/ 01 апреля 2019

Соберите строки для каждой единицы в буфере; при каждой новой строке заголовка обрабатывайте предыдущий блок (проверьте последнюю строку в буфере, печатать или нет) и очистите буфер

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: $0 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 (задокументировано в perlop ):

tr/..//dr удаляет все перечисленные символы из своей целевой строки, возвращая измененную строку, сохраняя оригинал без изменений (из-за «неразрушающего» модификатора /r). Поэтому, если после удаления разрешенных символов остается что-либо, мы отбрасываем модуль (не печатайте его).


Примечание по выбору tr и эффективности

Можно использовать регулярное выражение и его оператор сопоставления с классом отрицанных символов,

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

вместо оператора транслитерации tr (который не регулярное выражение).

Однако, даже в самом лучшем случае для оператора матча я оцениваю подход tr, чтобы он был 25% быстрее. Во всех остальных случаях tr опережает соответствие регулярному выражению, по крайней мере, в 2-4 раза.

Оператор соответствия «лучший случай» - это когда недопустимый символ находится на первой позиции в строке, чтобы он сразу совпадал и не сканировал остальную часть строки. Это довольно нереально, если не сказать больше. Напротив, статистически во многих (большинстве?) Строках не будет ни одного из этих символов, и у оператора совпадения будет наихудший случай сканирования всей строки. (Но большая часть стоимости, вероятно, в любом случае, это запуск двигателя регулярного выражения.)

2 голосов
/ 01 апреля 2019

Вот подход, который использует (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;
    }   
}
0 голосов
/ 16 апреля 2019

если ваши данные в файле '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
...