1
votes

Supprimer les lectures qui ont des appels de base de mauvaise qualité

J'ai un fichier de données de séquence au format fastq https://en.wikipedia.org/wiki/ FASTQ_format où la première ligne est l'ID de séquence, la deuxième ligne est la séquence [ACGT], la troisième ligne est «+» et la 4ème ligne est les valeurs de qualité.

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

Je veux supprimer toutes les entrées où un caractère de la 4ème ligne ne correspond à aucun des symboles ? @ABCDEFGHIJK

La sortie de l'exemple ci-dessus sera

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

linux awk perl fastq

4 commentaires

Une unité est composée de 4 lignes (SEQ ID, sequence, + sing, quality for sequence). Je souhaite supprimer toutes les unités où la qualité de la séquence (4ème ligne) pour chaque caractère de la séquence (2ème ligne) correspond à autre chose que l'un ou l'autre des caractères du motif? @ABCDEFGHIJK. J'ai essayé ce code et j'y travaille toujours {cat file.fq | awk 'NR% 4 == 0' | xargs -n1 awk '{for (i = 0; ++ i <= length ($ 0);) printf "% s \ n"}'}

Avez-vous examiné les différents modules FASTA qui sont disponibles?

@Tanktalus J'ai essayé plusieurs outils utilisés pour couper les fichiers fastq, notamment Trim à gogo et Sickle. Ces outils suppriment les lectures de mauvaise qualité, mais il y a encore des lectures qui ont des bases de mauvaise qualité que je souhaite supprimer.

Attention, s'il est vrai que les fichiers FASTQ sont généralement formatés de telle sorte que l'enregistrement couvre exactement 4 lignes, ce n'est pas sa définition. Il existe des fichiers FASTQ valides dans lesquels des enregistrements individuels s'étendent sur plus de 4 lignes, ce qui peut être nécessaire lors de leur analyse. Quoi qu'il en soit, des outils comme Trim Galore peuvent gérer ce cas avec une configuration appropriée.

4 Réponses :

3
votes

Collectez des lignes pour chaque unité dans un tampon; à chaque nouvelle ligne d'en-tête traiter l'unité précédente (vérifier la dernière ligne du tampon, imprimer ou non) et vider le tampon

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

Explication du contrôle de la dernière ligne du tampon avec tr (documenté dans perlop ):

Le tr /..// dr supprime tous les caractères répertoriés de sa chaîne cible, retournant la chaîne modifiée tout en gardant l'original inchangé (en raison du / r "non destructif" modificateur). Donc, s'il reste quelque chose après la suppression des caractères autorisés, nous supprimons l'unité (ne l'imprimez pas).

Une note sur le choix de tr et efficacité

Il est possible d'utiliser regex et son opérateur de correspondance avec une classe de caractères annulée,

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

au lieu de l'opérateur de translittération tr (qui n'est pas regex).

Cependant, même dans le meilleur des cas pour l'opérateur de correspondance, je compare l'approche tr pour être 25% plus rapide. Dans tous les autres cas, tr surpasse la correspondance de l'expression régulière d'un facteur de 2 à 4, au moins.

Le "meilleur cas" de l'opérateur de correspondance est lorsqu'un caractère inacceptable se trouve en première position dans la chaîne, afin qu'elle corresponde immédiatement et ne scanne pas le reste de la chaîne. C'est pour le moins irréaliste. Au contraire, statistiquement, la plupart des chaînes (la plupart?) N'auront aucun de ces caractères et l'opérateur de correspondance aura son pire cas, en parcourant toute la chaîne. Notez que dans tout cela, le démarrage du moteur regex est probablement la partie la plus chère.

0 commentaires

3
votes

$ 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

0 commentaires

2
votes

Voici une approche utilisant (1) un économiseur d'enregistrement d'entrée modifié et (2) l'opérateur de translittération modifié avec le complément tr /// c c switch.

(J'ai simulé le fichier en haut du script)

#!/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 commentaires

0
votes

si vos données dans le fichier 'd' peuvent être extraites:

 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

0 commentaires