[bioinfo] NCBI et BioPerl et OO

Version imprimable

Bonjour à tous,

Je découvre ce site, il parait très intéressant. Je suis dans la recherche biomoléculaire et j'utilise bioperl afin d'analyser des séquences d'ADN. Je suis bloquée dans un programme et j'aimerais un coup de main de votre part.

J'aimerais pouvoir accéder aux valeurs contenues dans un objet composé d'objets. La classe est get_SeqFeatures et j'ai plusieurs instances de cette classe dont j'aimerais extraire des informations à l'aide d'une boucle.
Voici un exemple d'objet obtenu par : print Dumper($Features)."\n";.

Code:

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
39
40
41
42
43
44
45
46
47
48
$VAR1 = bless( {
 
                 '_gsf_tag_hash' => {
                                      'mol_type' => [
                                                      'mRNA'
                                                    ],
                                      'db_xref' => [
                                                     'taxon:185431'
                                                   ],
                                      'isolate' => [
                                                     '927/4 GUTat10.1'
                                                   ],
                                      'chromosome' => [
                                                        '1'
                                                      ],
                                      'organism' => [
                                                      'Trypanosoma brucei TREU927'
                                                    ],
                                      'note' => [
                                                  'single VAT derivative of TREU927/4 (GPAL/KE/70/EATRO 1534)'
                                                ]
                                    },
 
               '_source_tag' => 'EMBL/GenBank/SwissProt',
 
                 '_gsf_seq' => bless( {
                                        'display_id' => 'XM_001219120',
                                        'primary_id' => '115504656',
                                        'accession_number' => 'XM_001219120',
                                        '_seq_length' => '2097',
                                        '_version' => '1',
                                        'alphabet' => 'dna',
                                        'desc' => 'Trypanosoma brucei hypothetical protein, conserved (Tb927.1.4280) mRNA, complete cds.',
                                        '_root_verbose' => 0,
                                        'seq' => 'ATGGAGAAGGTAGCAACAACAGCGTCTAAGGACAGTTCTGCTGTTAGGATTGTTGAACCAGTGATGTGA'
                                      }, 'Bio:: PrimarySeq' ),
 
                 '_location' => bless( {
                                         '_strand' => 1,
                                         '_seqid' => 'XM_001219120',
                                         '_location_type' => 'EXACT',
                                         '_start' => '1',
                                         '_end' => '2097',
                                         '_root_verbose' => 0
                                       }, 'Bio::Location::Simple' ),
 
  '_primary_tag' => 'source'
               }, 'Bio::SeqFeature::Generic' );

J'aimerais récupérer les valeurs contenues dans "'_gsf_tag_hash'", la référence ('taxon:185431') et l'isolat ('927/4 GUTat10.1'). Mais tous mes objets de classe get_seqFeatures ne possèdent pas "_gsf_tag_hash. Il faut donc que je commence par un test (if exists) et ensuite que j'extraie l'information.
Pourriez vous m'expliquer comment accéder aux informations provenant d'objets composés?

Voici ce que j'ai fait

Code:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
      my($Features)= $seq->get_SeqFeatures;
      foreach my $k (keys %$Features)
      {
        #print "$k : $Features->{$k}\n";
        if ($k eq "_gsf_tag_hash")
        {
                my $SousObjet1=($Features->{$k});
                foreach my $key (keys(%$SousObjet1))
                {
                        #print "$key \t $SousObjet1->{$key}\n ";
                        if ($key eq "db_xref")
                        {
                        @taxon=($SousObjet1->{'db_xref'});
                        print $taxon[0];
                        }
 
                }
 
        }

mais $taxon[0] possède une array.

Pourriez- vous m'aider?

Merci beaucoup,

Jasmine,

26/10/2006, 11h32
djibril

Bienvenue à une bioinformaticienne sur le forum :mrgreen: et un de plus. Tout d'abord, reedite ton message en mettant ton code perl entre balise code (tu selectionne ton code et click sur le #), ce sera plus clair, plus lisible pour tout le monde
26/10/2006, 11h50
stoyak

Bienvenue à toi!

Alors, après indenté 1 peu ton code et l'avoir mis entre balise, peux-tu nous préciser certaines choses:

- donne nous 1 petit bout de script pour voir ta façon de coder et voir sous quelle forme tu appelles ta classe (module, procédures, etc ..)
- de quelle manière sont stockées tes informations? hash? tableau? scalaire?
la façon de récupérer tes infos va dépendre de ça!

Merci beaucoup pour cet accueil. Je débute et j'espère que la façon dont je programme n'est pas trop affreuse. Si vous avez des conseils de style, n'hésitez pas, ils sont les bien venus. Je vous donne tout le code ainsi vous aurez toutes les réponses.

Code:

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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
#!/usr/local/bin/perl
 
 
# Ce script interroge NCBI au sujet des informations suivantes:
#       Acc number (ID)
#       Gi (Primary Id)
#       Type de la séquence (ADN, ARN ou PROT)
#       Description de la séquence
#       Séquence
 
 
 
 
use Bio::Perl;
use Bio::Annotation::Collection;
use strict;
use warnings;
use FileHandle;
use Data::Dumper;  # visualisation des OO
 
 
# Args    : -db       database ('protein' or 'nucleotide')
#           -query    query string
#           -mindate  minimum date to retrieve from
#           -maxdate  maximum date to retrieve from
#           -reldate  relative date to retrieve from (days)
#           -datetype date field to use ('edat' or 'mdat')
#           -ids      array ref of gids (overrides query)  exemple :   -ids=>[195052,2981014,11127914]
#           -maxids   the maximum number of IDs you wish to collect (defaults to 100)
 
 
my ($i)=0;
my ($db)="nucleotide";
my ($mindate)="";
my ($maxdate)="";
my ($reldate)="";
my ($datetype)="";
my ($ids)="";
my($maxids)="";
my ($query_string) = 'XM_001219120';
 
 
 
 
 
 
   my $query = Bio::DB::Query::GenBank->new(-db=>$db,
                                            -query=>$query_string,
					    -mindate => $mindate,
					    -maxdate => $maxdate,
                                            -reldate => $reldate,
                                            -datetype => $datetype,
                                            -ids => $ids,
                                            -maxids => $maxids
                                          );
 
        my $count = $query->count;
        print "Compte du nombre d'entrees = $count\n";
 
 
 
Info();
 
 
close;
 
 
 
 
 
 
sub Info
 
{
 
 
   # get a genbank database handle
   my $gb = new Bio::DB::GenBank;
   my $stream = $gb->get_Stream_by_query($query);
   while (my $seq = $stream->next_seq)
   {
 
      my($string);                        #$seq->revcom->seq();
      $string = $seq->seq();
      my($Acc);
      $Acc= $seq->accession_number();
      #my ($Id);
      #$Id=$seq->display_id();       # the human read-able id of the sequence = Acc
      my ($Type);
      $Type=$seq->alphabet();         # one of 'dna','rna','protein'
      my ($PrimId);
      $PrimId=$seq->primary_id();       # a unique id for this sequence irregardless  = Gi
      my ($Desc);
      $Desc = $seq->desc();             # a description of the sequence
      my (@taxon);
 
 
 
 
       my($Features)= $seq->get_SeqFeatures;      # The 'top level' sequence features
       my ($AllSeq)=$seq->get_all_SeqFeatures;  # All sequence features, including sub-features
 
 
 
      foreach my $k (keys %$Features)
      {
        print "$k : $Features->{$k}\n";
        if ($k eq "_gsf_tag_hash")
        {
                my $SousObjet1=($Features->{$k});
                foreach my $key (keys(%$SousObjet1))
                {
                        print "$key \t $SousObjet1->{$key}\n ";
                        if ($key eq "db_xref")
                        {
                        @taxon=($SousObjet1->{'db_xref'});
                        print $taxon[0];
                        }
 
                }
 
        }
 
     }
 
 
 
 
 
 
 
      #print Dumper($Features)."\n";
 
 
 
 
   }
 
 
 
 
      $i++;
}

Mon but est d'ensuite faire de ce script un module. Est-ce une bonne idée? Il renverrai un tableau indexé contenant toutes les infomations.
Dans mon travail, j'utilise MySQL et j'ai déjà écrit un programme créant une table à partir des informations que j'arrive à prélever, il fonctionne très bien mais j'aimerais plus d'informations. Vaut il mieux réécrire un second programme et utiliser un module me donnant les informations NCBI ou au contraire ne faire qu'un seul script interrogeant NCBI et chargeant la DB?

Merci,

Jasmine,

26/10/2006, 13h34
djibril

Ce forum est visité par de nombreuses personnes bioinformaticien, informaticien, etc. Donc si tu veux avoir plus d'informations et d'aides, il faut que tu sois le plus clair possible.
Tout le monde ne connait pas le NCBI et ne sait pas ce que fait ton module et le fichier que tu veux parser. Donc faudrait que tu donnes des bouts d'exemple de fichier à parser et t'explique un peu le module bioperl que tu utilises. De plus, tu parses pas le NCBI :mouarf: , mais des fichiers GENBANK, donc attention à l'abus de langage ;) .

Sinon je pense pas que tu ai vraiment besoin de faire un module, mais bon pourquoi pas, c'est bien pour l'apprentissage.
En ce qui concerne ton souci de $taxon[0], comme tu le sais, un accesion number peut avoir souvent plusieurs taxonomy, donc c'est pour ça que tu peux te retrouver avec un array, ce qui veut dire qu'il t'a renvoyé une référence à un tableau. pour avoir la premiere case de ce tableau par exemple, c'est ainsi :
${$taxon[0]}[0]. J'ai écrit une doc dans la FAQ sur les références, je te conseil donc de la lire.

Sinon quel est le résultat de ton code? t'as fait des print, qu'est ce que ça donne?
26/10/2006, 13h57
Jasmine80

Recoucou,

Oui, tu as raison, je mélange NCBI et Genbank. Je vais faire attention.
J'ai remplacé dans le code print "\n".${$taxon[0]}[0]."\n"; et cela fonctionne très bien, merci.

Voici ce que mon script me donne:

Compte du nombre d'entrees = 1

_gsf_tag_hash : HASH(0x2198e0c)
mol_type ARRAY(0x2198e48)
db_xref ARRAY(0x2198e24)

taxon:10255
strain ARRAY(0x2198e6c)
isolate ARRAY(0x2198e90)
organism ARRAY(0x2198eb4)

-----------------------
_source_tag : EMBL/GenBank/SwissProt

-----------------------
_gsf_seq : Bio::PrimarySeq=HASH(0x219f320)

-----------------------
_location : Bio::Location::Simple=HASH(0x2198c68)

-----------------------
_primary_tag : source

-----------------------

Je découvre les références depuis ce matin. J'ai trouvé un chouette turorial.
=> http://sylvain.lhullier.org/publicat...hapitre11.html
(http://sylvain.lhullier.org/publications/intro_perl/)

j'utilise aussi http://perl.enstimac.fr/DocFr/perl.html

Merci beaucoup,

Jasmine,

Merci beaucoup maintenant je peux récupérer toutes les informations dont j'ai besoin.

Voici mon script:

Code:

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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
sub Info
 
{
 
 
   # get a genbank database handle
   my $gb = new Bio::DB::GenBank;
   my $stream = $gb->get_Stream_by_query($query);
   while (my $seq = $stream->next_seq)
   {
 
 
 
      my(@Acc);
      $Acc[$i]= $seq->accession_number();
 
 
      $Info{$Acc[$i]}{"Sequence"}= $seq->seq();
      $Info{$Acc[$i]}{"Type"}=$seq->alphabet();
      $Info{$Acc[$i]}{"Gi"}=$seq->primary_id();
      $Info{$Acc[$i]}{"description"}=$seq->desc();
 
 
      my @Taxon="";
      my @Souche="";
      my @Isolat="";
 
      my($Features)= $seq->get_SeqFeatures;      # The 'top level' sequence features
 
      foreach my $k (keys %$Features)
      {
        #print "$k : $Features->{$k}\n";
        if ($k eq "_gsf_tag_hash")
        {
                my $SousObjet1=($Features->{$k});
                foreach my $key (keys(%$SousObjet1))
                {
                        #print "$key \t $SousObjet1->{$key}\n ";
                        if ($key eq "db_xref")
                        {                                   #$Table{$i}{"Acc"}=$Acc[$i];
                        @Taxon=($SousObjet1->{'db_xref'});
                        $Info{$Acc[$i]}{"Taxon"}=${$Taxon[0]}[0];
                        }
                        if ($key eq "db_xref")
                        {
                        @Souche=($SousObjet1->{'strain'});
                        $Info{$Acc[$i]}{"Souche"}=${$Souche[0]}[0];
                        }
                        if ($key eq "db_xref")
                        {
                        @Isolat=($SousObjet1->{'isolate'});
                        $Info{$Acc[$i]}{"Isolat"}=${$Isolat[0]}[0];
                        }
 
                }
 
        }

Est ce correct comme style?
J'ai juste un avertissement pour un nombre impair d'éléments dans mon tableau indexé.

Encore merci, mon problème est quasi résolu, il ne me reste plus qu'à en faire un module.

Jasmine,

26/10/2006, 14h29
djibril

c'est bien. A toi d'ecrire comme ce sera lisible et comprehensible si tu reviens dans 1 mois ou plus dessus avec des commentaires. Il y aura peut etre des cas particulier à traiter dans ton script, mais tu verra cela en fonction de tes tests.

voilà, bon courage
26/10/2006, 14h37
Jasmine80

Oui, merci pour ton aide. Il est très chouette ce site. Je vais essayer d'aider les autres en répondant à leurs questions si du moins elles ne sont pas trop difficiles.
Dernière petite chose.
Que puis-je faire pour le message "odd number of elements in hash asignment" ? Dois-je faire une "if !defined" et assigner une valeur vide?

Merci,
26/10/2006, 18h02
stoyak

Tu m'as l'air d'utiliser régulièrement BioPerl. Sais-tu s'il existe un module pour récupérer les annotations fonctionnelles et si certains modules ont été développés pour l'analyse statistiques des puces a ADN?
Si tu as des infos je suis preneuse!
26/10/2006, 18h56
Jasmine80

Bonsoir Stoyak ,

Je ne sais pas si de tels modules existents. Voici le site à consulter :

http://www.cpan.org/

Pour la liste des modules va directement sur ce lien : http://search.cpan.org/~birney/bioperl-1.2.3/

Je ne fais pas de microarray mais certains modules existent:
Bio-Affymetrix-0.5
Bio::Affymetrix::CHP
Regarde sur CPAN et vois ce qui te convient. Pour ce qui est des annotations fonctionnelles, je ne sais pas. Je n'utilise BioPerl que depuis 3 semaines. Je découvre comme toi. Enfin, peut-être que tout ce que je te dis est basic et que tu t'y connais doix fois mieux que moi. Bonne chance dans tes recherches.

Je te conseille le langage de programmation "R" (l'open source de S). Pour les chips et les biostat, je crois que c'est le meilleur de tous, et de loin. Il y a moyen d'utiliser R via Perl, j'ai visité un site à ce sujet mais son adresse est au boulot. Si cela t'intéresse, demain je poste son lien.
http://www.r-project.org/

Avec Perl regarde le package suivant si cela t'intéresse:
Statistics::R => Controls the R (R-project) interpreter through Perl.
http://search.cpan.org/~gmpassos/Sta...tatistics/R.pm

Biiiiz,

Jasmine,

PS: je te conseille de poser ta question dans un nouveau post afin que d'autres la liste.
27/10/2006, 11h22
stoyak

je connais déjà le langage R et le module perl Statistics::R, je l'utilise depuis presque 6 mois maintenant. J'ai fait le tour de tout les soucis possibles je crois! Tu pourras certainement retrouver mes posts à ce sujet!
Ajouter 1 nouveau post ne me semble pas fort utile, tu verras que les réponses à ce sujet ne sont pas légion: seuls les quelques bionfos qui traînent ici te répondront :-)

Merci de tes renseignements!
27/10/2006, 14h30
Jasmine80

Bonjour,

J'aimerais pouvoir t'aider mais je ne connais pas la réponse à ta question. Si je trouve des modules pouvant t'intéresser, je te les MP.

Voici ce que mon programme va rechercher comme informations

AY223491 => Souche => Khateen
AY223491 => ProtGi => GI:37903922
AY223491 => Taxon => taxon:10255
AY223491 => Type => dna
AY223491 => ProtId => AAP48876.1
AY223491 => Gi => 37903921
AY223491 => Isolat => S/KHA
AY223491 => ProtProduit => 14K membrane protein
AY223491 => Description => Variola virus strain Khateen 14K membrane protein gene, complete cds.
AY223491 => ProtNote => similar to Vaccinia virus strain Copenhagen A27L
AY223491 => Sequence => ATGGACGGAACTCTTTTCCCTGGAGATGATGATCTTGCAA

Si ProtProduit te suffit comme annotation fonctionnelle, je peux t'expliquer comment faire.

Bonne journée,

Jasmine,