Comment fusionner deux fichiers FASTA (un fichier avec saut de ligne) en Perl ?

Question

J'ai deux fichiers Fasta suivants :

fichier1.fasta

>0
GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT
>1
GTTAAGTTATATCAAACTAAATATACATACTATAAA
>2
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC

fichier2.qual

>0
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 40 15 40 40
>1
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40
40 40 40 40 40 40 40 40 40 40 40
>2
40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4
40 8 3 29 10 19 18 40 19 15 5

Notez le saut de ligne dans le fichier "qual" pour chaque en-tête fasta - marqué par ">".Le numéro d'en-tête du fichier (« > ») est le même pour les deux fichiers.Nombre de qualités numériques = longueur de la séquence.

Ce que je veux faire, c'est ajouter ces deux fichiers donnant :

GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT  40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 15 40 40
GTTAAGTTATATCAAACTAAATATACATACTATAAA  40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40 40 40 40 40 40 40 40 40 40 40 40
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC  40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4 40 8 3 29 10 19 18 40 19 15 5

Mais d’une manière ou d’une autre, mon code ci-dessous ne parvient pas à le faire correctement ?En particulier, la deuxième ligne de chaque entrée du fichier « qual » n'est pas imprimée.

use strict;
use Data::Dumper;        
use Carp;
use File::Basename;      

my $fastafile = $ARGV[0] || "reads/2039F.2.fasta"; 
my $base      = basename( $fastafile, ".fasta" );
my $qualfile  = "reads/" . $base . ".qual";
print "$qualfile\n";

open SEQ, '<', $fastafile or die $!; #Seq
open PRB, '<', $qualfile or die $!; #quality


while (my $seq = <SEQ>) {
     my $qual = <PRB>;
     chomp($seq);
     chomp($qual);

     if ($seq =~ /^>/ || $qual =~ /^>/) {
         next;
     }
     else {
         print "$seq\t$qual\n";      
     }

}

Quelle est la bonne façon de procéder ?

Answer 1

Il vous manque la deuxième ligne (et toutes les lignes suivantes) des scores de qualité et vous manqueriez également des lignes de séquence supplémentaires.À cette fin et à des fins de réutilisation du code, la manière de gérer les séquences FASTA est sous forme d'entrées/enregistrements entiers :

local $/ = "\n>";
while (my $seq = <SEQ>) {
     my $qual = <PRB>;
     chomp($seq);  $seq =~ s/^>*.+\n//;  $seq =~ s/\n//g;
     chomp($qual);  $qual =~ s/^>*.+\n//;  $qual =~ s/\n/ /g;

     print "$seq\t$qual\n";      

}

Vous pouvez également facilement capturer l'en-tête FASTA lors du premier remplacement.

Answer 2

Le problème est que vous avancez dans le fichier en parallèle, donc lorsque la ligne est ">" dans un fichier, il se peut qu'elle ne soit pas ">" dans le suivant.

La façon dont vous lisez les données se fait par paires, comme ceci :

1: >0 
2: >0
1: GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT
2: 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
1: >1
2: 40 40 40 40 40 40 40 40 15 40 40
1: GTTAAGTTATATCAAACTAAATATACATACTATAAA
2: >1
1: >2
2: 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40
1: GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC
2: 40 40 40 40 40 40 40 40 40 40 40
1: EOF
2: >2
1: EOF
2: 40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4
1: EOF
2: 40 8 3 29 10 19 18 40 19 15 5

Le même ensemble de données appliqué à vos règles de boucle ferait ceci :

1: GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT
2: 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
1: GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC
2: 40 40 40 40 40 40 40 40 40 40 40

Vous devez donc soit séparer la logique de boucle, soit trouver un moyen de faire correspondre les fichiers.

Voici une tentative de séparation de la recherche, mais je ne l'ai pas testée.

fileIO: {
  while( 1 ){ 
   my $seq; 
   my $qual  = q{};
   while( 1 ){ 
     $seq = <SEQ>; 
     last fileIO if not $seq;  # stop at end of file
     last if $seq !~ /^>/; 
  }
  while( 1 ){ 
     my $qual_in = <PRB>;
     last fileIO if not $qual_in; # stop at end of file 
     last if $qual_in =~ /^>/ and $qual ne q{}; 
     next if $qual_in =~ /^>/ and $qual eq q{}; 
     $qual .= $qual_in;
  }
  print "$seq \n $qual \n";

 }
}

Mise à jour

J'ai refactorisé le code ci-dessus en une seule fonction qui lira un morceau d'un descripteur de fichier arbitraire selon les besoins, cela semble fonctionner selon les besoins.Notez bien sûr que j'ai expérimenté ici un peu une astuce que je voulais utiliser pour quelque chose de pratique.

use strict;
use warnings;

# 
#  readUntilNext( $fileHandle, \$scalar_ref ); 
#
#  returns 0 when nothing could be read from the fileHandle. 
#  otherwise returns 1; 
#

sub readUntilNext {
    my ($fh)            = shift;
    my ($output)        = shift;
    my ($output_buffer) = '';
    while (1) {
        my $line = <$fh>;
        if ( !$line ) { # No more data
            # No data to flush to user, return false.
            return 0 if $output_buffer eq q{};
            last;  # data to  flush to user, loop exit. 
        }
        if ( $line =~ /^>/ ) {
            # Didn't get anything, keep looking. 
            next if $output_buffer eq q{};
            # Got something, flush data to user. 
            last;
        }
        chomp($line);
        $output_buffer .= $line;
    }
    # Data to flush to user 
    # Write to the scalar-reference 
    $$output .= $output_buffer;
    return 1;
}

open my $m, '<', 'a.txt';
open my $n , '<', 'b.txt';
# Creates 2 scalar references every loop, and only loops as long 
# as both files have data. 
while ( readUntilNext( $m, \my $seq ) && readUntilNext( $n, \my $qual ) ) {
    print "$seq\t$qual\n";
}

Et le code ci-dessus, testé, fait exactement ce que vous voulez faire.

Notez à ce sujet mes affaires

while( readUntilNext( $m, \my $seq ) ) { 
}

est fondamentalement la même chose que

my $seq; 
while( readUntilNext( $m, \$seq ) ) { 
}

Sauf que le premier crée à chaque fois un nouveau scalaire, garantissant que la même valeur ne sera pas visible dans une boucle successive ;

donc ça devient plutôt :

while( 1 ){ 
 my $seq; 
 last if not readUntilNext($m, \$seq);
 do { 
    # loop body here
 }
}

Answer 3

Voici une solution n'utilisant pas Perl, mais des commandes Shell simples :

prompt>grep -v '^>[0-9]' file1.fasta > tmp1
prompt>(tr '\012' ' ' < file2.qual; echo) | sed 's/>[0-9]* /\n/g' | sed 1d > tmp2
prompt>paste tmp1 tmp2
GAATAGATGTTTCAAATGTACCAATTTCTTTCGATT    40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 15 40 40
GTTAAGTTATATCAAACTAAATATACATACTATAAA    40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 20 40 40 40 40 40 40 40 40 40 40 40 40 40 40
GGGGCTGTGGATAAAGATAATTCCGGGTTCGAATAC    40 40 40 40 7 40 40 5 40 40 40 40 40 40 40 40 37 13 31 20 15 40 10 11 4 40 8 3 29 10 19 18 40 19 15 5
prompt>

J'ai cherché pendant de nombreuses années la commande coller (sachant "c'est une opération super basique, quelqu'un doit déjà mis en œuvre quelque chose pour résoudre ce problème").

La deuxième ligne de commande traduit d'abord toutes les nouvelles lignes en espaces, et la commande echo est ajoutée pour ajouter une nouvelle ligne finale à l'entrée (car sed ignorera les lignes dépourvues d'EOL), joignant ainsi toutes les lignes d'entrée en une seule ligne qui sera ensuite exécutée par la commande sed. se divise à nouveau (note de portabilité :tous les programmes sed ne fonctionneront pas avec des longueurs de ligne arbitraires, mais GNU sed le fait).