Discussion :

[pontarose]

bonjour,

J'ai 2 fichiers un fichier fastq et fichier qui contient des identifiants d'intérêt.
J'aimerais parser le fichier fastq de manière à obtenir un tableau de hash avec pour clé l'identifiant et pour valeur les 3 lignes suivantes.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
@HCDPQ1D0501
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT.
  +HCDPQ1D0501
  !''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65.....

Ensuite je souhaiterai parcourir mon tableau de hash et parcourir les identifiants de mon fichier pour récupérer au final un fastq des identifiants d'intérêt


J'étais parti comme avec un fichier fasta, mais en utilisant Bio::Seq::Quality mais la méthode next_seq n'existe pas dans ce package!




####parser le fichier fastq

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
use Bio::SeqIO;
use Bio::Seq::Quality;
 
my $file = 'file.fastq';
my $in  = Bio::Seq::Quality->new(-file => $file , '-format' => 'fastq');
 
my %hash;
my $nb_seq++;
 
while ( my $seq = $in->next_seq() ){
    $nb_seq++;          
    $hash{$seq->seq} = $seq->id ;
}
 
 
####parser le fichier des identifiants
open (DATA, "identifiant_file") || die("Impossible d'ouvrir le fichier: $! "); 
my  $nb=0;
my @tab;
while ( my $ligne = <DATA> ) {
	chomp $ligne;
        my ( $id, $mot) = split /\s+/, $ligne;
        $nb++;
		$tab[$nb][1]=$id;
}
 
##récupérer fastq
open(FILE, '>final.fastq') || die "pbe: $!";
 
foreach my $seq (keys %hash){   
         for(my $i=1; $i<= $nb; $i++){
                  if($hash{$seq} eq $tab[$i][1]){
                         print (FILE "$hash{$seq}\n$seq\n");
                         last;
                  }
         }
}
close (FILE);

[Jasmine80]

Bonjour,


Je ne connais pas du tout ce module. Pourrais-tu poster un fichier complet en exemple?


merci.

[Beniou]

Bonjour,

Pour cela tu devrais utiliser le module Bio::SeqIO::fastq qui hérite des fonctions next_seq du module Bio::SeqIO

[pontarose]

Merci Beniou,
le module Bio::SeqIO::fastq hérite bien des fonctions next_seq du module Bio::SeqIO.

Cependant d'après mon exemple, j'obtiens:
@HCDPQ1D0501 #$hash{$seq}
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT. #$seq

Alors que je voudrais:
clé =@HCDPQ1D0501
valeur = (les 3 lignes suivantes)
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT.
+HCDPQ1D0501
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>>CCCCCCC65.....

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
use Bio::SeqIO::fastq;
 
my $file = 'file.fastq';
my $in  = Bio::Seq::fastq->new(-file => $file , '-format' => 'fastq');
 
my %hash;
my $nb_seq++;
 
while ( my $seq = $in->next_seq() ){
    $nb_seq++;          
    $hash{$seq->seq} = $seq->id ;
}
foreach my $seq (keys %hash){  
	print "$hash{$seq}\n$seq\n"; 
}

[Beniou]

Pour avoir toutes le informations que tu souhaites il faut utiliser les bonnes méthodes du module Bio::Seq::Quality une fois que tu parcours les séquences : qual(), qual_text() etc. suivant ce que tu veux obtenir et de quelle manière (tableau, texte etc.)

Voici avec ton exemple :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
10
11
12
13
use strict;
use warnings;
use Bio::SeqIO::fastq;
 
my $file = 'file.fastq';
my $fastq  = Bio::SeqIO::fastq->new(-file    => $file,
                                    -format  => 'fastq');
 
while ( my $seq = $fastq->next_seq() ){
    print "ID   : ".$seq->id()."\n";
    print "SEQ  : ".$seq->seq()."\n";
    print "QUAL : ".$seq->qual_text()."\n";
}

Cela affiche

ID : HCDPQ1D0501
SEQ : GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAACTCACAGTTT
QUAL : 0 6 6 9 7 7 7 7 9 9 9 10 8 8 4 4 4 10 10 8 7 4 4 4 4 8 13 16 9 9 9 12 10 9 6 6 8 8 9 9 20 20 34 34 37 29 29 29 29 29 29 34 34 34 34 34 34 34 21 20

Une choss : les qualités sont automatiquement converties en PHRED (voire l'option variant de la méthode "new" pour définir tes types de qualités : solexa, illumina ou sanger)

Si tu veux vraiment les qualités telles quelles et bien peut être qu'il y a une méthode du package Quality mais je ne l'ai pas trouvé.