Comparaison entre séquences et numéro d'accessions swissprot

**zaboug** · 09/07/2009, 15h37

Bonjour à tous,

comme le titre du post l'indique, je me demandais comment gérer les temps d'attente avec mechanize.
Je m'explique, je veux interroger automatiquement cette page.
J'ai donc un script qui remplie le champs en haut de la page, qui clique sur le bouton "submit", puis qui suis le lien menant vers la page des resultats :

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
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
#!/usr/bin/perl -w
 
 
use strict;
use diagnostics;
use WWW::Mechanize;
use Data::Dumper;
use vars qw($opt_h);
use Getopt::Std;
 
##recuperation des arguments - initialisation
getopts("h") or die("erreur(s) de saisie dans les options. Faire -h pour obtenir l'aide\n");
 
if (defined($opt_h))
{
	&help;
	exit;
}
 
 
my $peptide = "NTDQASMPDNTAAQK";
my $AC = "P14105";
 
my $url = 'http://services.uniprot.org/blast/';
 
#creation du robot qui naviguera sur internet
my $bot = WWW::Mechanize->new();
 
#on va sur la page qui nous interesse
$bot->get($url);
$bot->form_number(1);
 
$bot->set_fields(query => $peptide);
$bot->submit();
 
die $bot->res->status_line unless $bot->success;
 
 
$bot->follow_link(text => 'sequence similarity search');
 
open(TEST, '>', 'test.html') or die("Impossible d'ouvrir test.html");
 
print(TEST Dumper $bot->content());
 
close(TEST);
 
sub help
{
	print qq{
analyse_identification_TMT.pl va sur le site de swissprot, lui entre passe une sequence peptidique et verifie que le fichier TMT contient la bonne identification.
 
Parametres (a venir):
	-p : peptide
	-a : AC
};
}

La première fois que j'ai lancé le script, ça a très bien marché.
J'ai retesté juste après avoir ajouter les variables $peptides et $AC, ce qui ne devrait rien changer pour le script, et ça n'a pas marché.
Enfin, si en quelque sorte : mon script a bien copier une page html, mais pas celle de résultats. En effet, le site met parfois un peu de temps (15-30sec) à répondre à la requête, et dans ce cas, il y a une page intermédiaire disant que la requete est en train de tourner depuis X secondes... C'est cette page que mon script copie.

D'où ma question : comment gérer ce cas de figure avec mechanize? comment le faire patienter jusqu'à ce que la bonne page de résultats arrive?
Je n'ai pas trouvé de fonction prédéfinie pour ça, mais peut être existe-il un moyen détourné?

Bref, la moindre idée sera étudiée

**djibril** · 09/07/2009, 16h30

Il existe des modules plus en adéquation avec ce que tu souhaites faire : Bio::DB::SwissProt. Ensuite, tu devras lire la documentation du module Bio::Seq pour voir les méthodes à utiliser.

Voici un exemple de script que je t'ai fais pour que tu puisse voir comment ça fonctionne.

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
#!/usr/bin/perl
use warnings;
use strict;
use Carp;
 
use Bio::DB::SwissProt;
 
my @ListeAccession = qw / TTK_HUMAN BRCA2_HUMAN /;
foreach my $accession ( @ListeAccession ) {
  my $sp = Bio::DB::SwissProt->new();
  my $ObjetBioseq = $sp->get_Seq_by_acc($accession); 
  print "Accession numbre : ",$ObjetBioseq->accession_number,"\n";
  print "Id : ",$ObjetBioseq->display_id,"\n";
  print "Description : ",$ObjetBioseq->description,"\n";
  print "Sequence : ",$ObjetBioseq->seq,"\n";
  print "========\n\n";
}

**zaboug** · 09/07/2009, 16h37

en fait, si j'ai bien compris, le script que tu me proposes fait l'inverse de ce que je veux :
ton script permet de récupérer la séquence à partir de l'accession number, c'est ça?

Mon but est de vérifier que l'identification d'une liste de peptides a été bien faite. Donc de récupérer l'Accession Number à partir d'une sequence peptidique.

Cela dit, ton exemple est une bonne base de travail. Juste une question : le module Bio:

B::SwissProt utilise la même base que celle qu'il y a sur internet? Les données sont bien à jour?

**djibril** · 09/07/2009, 16h39

La plus part des modules bioperl utilisent les scripts mis à disposition du NCBI, EMBL ou autre et les recherches sont faites directement dans les bases en ligne.

**zaboug** · 09/07/2009, 16h56

je viens de lire la doc du module Swissprot que tu m'as proposé, et je n'ai pas l'impression qu'il existe de fonction qui fait ce que je veux...

J'ai une liste de peptides et la liste d'accession number correspondant. Je veux vérifier que l'identification a été bien faite. Donc sur le site internet, je fournissais le peptide, et vérifiais que l'accession number que j'ai apparaissait bien dans la page de résultats. Bon, ça c'est que je voulais faire automatiquement et qui bloque pour l'instant à cause du temps de chargement de la page.

Je n'ai pas trouvé de solution dans la partie bioinformatique du forum Perl, et je ne vois pas non plus de fonction le permettant dans le module que tu m'a proposé.

du coup, j'aurais tendance à vouloir revenir à la solution internet, mais je me trompe peut être?

**djibril** · 09/07/2009, 17h45

Si tu utilises mon script en fournissant les accessions. Tu peux ensuite récupérer les séquences et les comparer.

Peux tu me donner une accession et sa séquence adéquat que tu souhaites vérifier ?

**zaboug** · 09/07/2009, 17h54

par exemple, le peptide NTDQASMPDNTAAQK a été identifié par l'AC P14105 chez le poulet, ce qui correspond à la myosine-9. Mais sur swissprot, on se rend compte que la myosine-11 est identifiée avec la même p-value que la myosine-9. Enfin bon, dans un premier temps, le but est de vérifier sur swissprot qu'on retrouve bien les bons AC pour nos peptides.

Pour commencer, j'ai donc deux fichiers d'environ 300 lignes à vérifier.

ce sont des résultats de spectrométrie de masse, donc les peptides sont pas très grands...

**djibril** · 09/07/2009, 18h04

En faite, quand tu soumets ton peptide sur uniprot, il fait un blast.
Donc c'est peut être ce que tu devrais faire via un module adéquat.

Fais un blast sur le NCBI avec ton peptide

>seq
NTDQASMPDNTAAQK

sur la base swissprot et tu verras ton résultat. Tu peux le faire via Perl

**zaboug** · 10/07/2009, 09h45

Pas bête du tout, je vais tenter ça

Merci Djibril, je mettrais mes avancées ici, mais du coup, ce post n'est peut être pas au bon endroit dans le forum

**djibril** · 10/07/2009, 10h05

Je le déplace dans le forum bioinformatique et je change son intitulé.

**zaboug** · 10/07/2009, 10h17

C'est super, merci

et ça tombe bien, car je viens de faire une recherche sur CPAN, et que je ne trouve pas de module permettant de faire ce que je veux

Enfin, je fais pas souvent de recherche sur CPAN, donc je suis peut être passé à côté de quelque chose...
Le seul module que j'ai trouvé qui ressemble à ce que je veux, c'est WWW:

DB car il y a une fonction blast() à laquelle on peut passer une sequence peptidique. Mais il faut aussi lui passer un cutoff, une matrice, et là, ça me dépasse complètement

Du coup, ça doit pas être le bon module.
J'ai vu aussi dans ce forum bioinfo qu'il est souvent question de blast, mais toujours à propos de séquences nucléiques. Est-ce la même chose pour les séquences peptidique?
Pourriez vous me donner des liens vers de la doc à ce sujet, car je rame un peu là...

**djibril** · 10/07/2009, 10h39

Une petite recherche de blast sur le CPAN te donne ceci

Bio::Tools::Run::RemoteBlast

**zaboug** · 10/07/2009, 10h40

bon, je viens de trouver sur ce forum une discussion à propos du module RemoteBlast.

Je vais tester ce module et essayer de me dépatouiller avec

**zaboug** · 10/07/2009, 10h42

Envoyé par djibril

Une petite recherche de blast sur le CPAN te donne ceci

oui, c'est la première recherche CPAN que j'ai faite (je suis un boulet, mais pas à ce point là

). Le truc c'est que la plupart de ces modules permettent d'analyse une sortie de blast, pas de faire un blast (ou alors, j'ai pas compris ce que j'ai lu, et je peux être un boulet à ce point là

)

EDIT : je suis un boulet, j'avais pas vu le module remoteBlast lors de ma recherche (enfin si, mais pas compris ce que signifiais remote...

)
Mais le principal c'est de finalement le tester, non?

**djibril** · 10/07/2009, 10h44

C'est bien, tu commence à saisir.
Fais un blast avec Perl, ensuite, via un autre module boulder::blast, tu parses le fichier résultat et ensuite tu pourras faire tes vérifications.

Voilou

**zaboug** · 10/07/2009, 11h08

Bon, j'ai commencé à tester le module remoteBlast à partir du script qu'ils proposent sur la page de présentation du module.
pour rappel, voici le code:

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
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
#!/usr/bin/perl -w
 
#Remote-blast "factory object" creation and blast-parameter initialization
 
use Bio::Tools::Run::RemoteBlast;
use strict;
my $prog = 'blastp';
my $db   = 'swissprot';
my $e_val= '1e-10';
 
my @params = ( '-prog' => $prog,
	 '-data' => $db,
	 '-expect' => $e_val,
	 '-readmethod' => 'SearchIO' );
 
my $factory = Bio::Tools::Run::RemoteBlast->new(@params);
 
#change a paramter
$Bio::Tools::Run::RemoteBlast::HEADER{'ENTREZ_QUERY'} = 'Homo sapiens [ORGN]';
 
#remove a parameter
delete $Bio::Tools::Run::RemoteBlast::HEADER{'FILTER'};
 
my $v = 1;
#$v is just to turn on and off the messages
 
my $str = Bio::SeqIO->new(-file=>'./amino.fa' , '-format' => 'fasta' );
 
while (my $input = $str->next_seq()){
#Blast a sequence against a database:
 
#Alternatively, you could  pass in a file with many
#sequences rather than loop through sequence one at a time
#Remove the loop starting 'while (my $input = $str->next_seq())'
#and swap the two lines below for an example of that.
my $r = $factory->submit_blast($input);
#my $r = $factory->submit_blast('amino.fa');
 
print STDERR "waiting..." if( $v > 0 );
while ( my @rids = $factory->each_rid ) 
{
  foreach my $rid ( @rids ) 
  {
	my $rc = $factory->retrieve_blast($rid);
	if( !ref($rc) ) 
	{
	  if( $rc < 0 ) 
	  {
		$factory->remove_rid($rid);
	  }
	  print STDERR "." if ( $v > 0 );
	  sleep 5;
	} 
	else 
	{
	  my $result = $rc->next_result();
	  #save the output
	  my $filename = $result->query_name()."\.out";
	  $factory->save_output($filename);
	  $factory->remove_rid($rid);
	  print "\nQuery Name: ", $result->query_name(), "\n";
	  while ( my $hit = $result->next_hit ) 
	  {
		next unless ( $v > 0);
		print "\thit name is ", $hit->name, "\n";
		while( my $hsp = $hit->next_hsp ) 
		{
		  print "\t\tscore is ", $hsp->score, "\n";
		}
	  }
	}
  }
}
}
 
# This example shows how to change a CGI parameter:
$Bio::Tools::Run::RemoteBlast::HEADER{'MATRIX_NAME'} = 'BLOSUM25';
 
# And this is how to delete a CGI parameter:
delete $Bio::Tools::Run::RemoteBlast::HEADER{'FILTER'};

J'ai créé le fichier amino.fa avec trois de mes sequences peptidiques, ce qui donne :

Code :

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

1
2
3
4
5
6
>seq1
ALLEVVQSGGK
>seq2
DMPYQPFSK
>seq3
EELSNVLAAMR

puis j'ai lancé le script.
le script fonctionne bien, mais il ne trouve pas de similarité :
voici une partie du fichier de sortie seq1.out :

Code :

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

<b>No significant similarity found.</b> For reasons why, <A HREF = "Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=FAQ"><b>click here</A>.</b><br><br>

J'ai essayé aussi en supprimant le '-expect' des params :

Code :

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

1
2
3
my @params = ( '-prog' => $prog,
	 '-data' => $db,
	 '-readmethod' => 'SearchIO' );

en me disant que mes sequences matchaient avec une p-value plus faible. Mais je n'obtiens toujours pas de résultats.

Y a-t-il des paramètres particulier à régler pour qu'il soit moins stringent, vu que mes peptides sont petits?
Je sais qu'il y a des identifications possibles puisque je les avais d'abord testé sur internet via le site uniprot...

**djibril** · 10/07/2009, 11h31

je viens de tester, j'ai aussi le même souci. Je vais chercher.

**zaboug** · 10/07/2009, 11h35

En comparant mes fichiers de résultats et les résultats obtenus par internet, et je me suis rendue compte que la base de données interrogées n'est pas la même :

- sur internet : UniprotKB
- avec le module : swissprot

j'ai essayé de mettre uniprotkb à la place de swissprot dans mon script, mais j'obtiens le message d'erreur suivant :

--------------------- WARNING ---------------------
MSG: <hr><p id="blastErr"><font color="red">An error has occurred on the server, The server is unable to format right now, please try again in a few minutes. If the problem (Informational Message: No alias or index file found for protein database [UniProtKB] in search path [/export/home/splitd/blastdb/blast2:/blast/db/disk.blast/blast2:]) persists - Contact Blast-help@ncbi.nlm.nih.gov and include your RID: 5DZ1P1S7016</font></p><hr>

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

Est ce que swissprot et uniprotkb sont deux base de données différentes?
Si oui, où peut-on trouver une liste de toutes les bases de données supportée par le module?

**djibril** · 10/07/2009, 11h54

Le module effectue le blast sur le site du NCBI : http://blast.ncbi.nlm.nih.gov/Blast.cgi

Les bases de données disponibles sont celles du site, notamment dans ton cas swissprot.

Le problème est que via le site, on a bien un résultat et vi le module non.

Il faut donc chercher pourquoi.

**zaboug** · 10/07/2009, 11h56

bon, en cherchant sur google, j'ai trouvé un pdf expliquant que le module RemoteBlast permettait d'automatiser les blasts en utilisant la base de données swissprot UniprotKB) :

Automating Blast
►BioPerl
Provides API’s for RemoteBlast and LocalBlast.
Provides API’s for parsing Blast output.
►Bio::Search::HSP::GenericHSP
►Bio::Search::Hit::BlastHit
►Database used –swissprot(UniprotKB)

J'en déduis qu'uniprotKB et swissprot correspondent à la même base de données.
Mais il y a quand même quelque chose de bizarre :
dans le fichier de resultats obtenu avec le module, il est précisé :

Code :

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

1
2
3
4
5
6
7
 
BLASTP 2.2.21+
[...]
Database: Non-redundant SwissProt sequences
    Posted date:  Jul 9, 2009  10:12 PM
  Number of letters in database: 11,275,576
  Number of sequences in database:  20,452

alors que sur le site internet, les informations données dans "job informations" sont :

Code :

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

1
2
3
4
5
Query  	 entry Sequence ;1 AA
Date run: 	2009-07-09 13:34:14 UTC+0100 on services.uniprot.org
Program 	NCBI BLASTP 2.2.20 [Feb-08-2009]
Database 	UniProtKB release 15.4 of 16-Jun-2009
  	9092914 sequences; 2958719688 amino acids.

donc sur internet, il y a beaucoup plus de sequences disponibles qu'avec le module!!!
Je vois pas bien comment contourner le problème du coup

EDIT : le site internet dont je parle c'est http://services.uniprot.org/blast, mais je viens de tester le site NCBI du lien que tu m'as donné et les résultats obtenu sur les deux sites sont équivalent (rassurant n'est ce pas

) et cette fois, la base de données utilisée est la même que le module :

Database Name
swissprot
Description
Non-redundant SwissProt sequences See details
Program
BLASTP 2.2.21+ Citation