Comment récupérer les données brutes des GEO Profiles d'un gène donné ?

**Norore** · 09/02/2011, 15h12

Bonjour,

je suis actuellement face à un problème de taille, je cherche à récupérer, pour un gène donné (dont je connais le GeneID), les données brutes de ses GEO Profiles afin de pouvoir les analyser à coup de scripts. J'ai bien trouvé via le site du NCBI comment accéder aux profils pour un gène, malheureusement je ne vois pas comment récupérer les informations de façon automatisée et en données brutes de préférence.

Apparemment, d'après les CPAN, le module Bio:: DB::Expression devrait me permettre de trouver les GEO, mais je ne vois pas du tout comment cela fonctionne ni si il s'agit bien de ce que je cherche à réaliser ! A moins que je ne doive passer par la récupération de la banque de données GEO et par des scripts pour extraire les informations dont j'ai besoin...

L'un de vous a-t-il déjà rencontré cette problématique ? Si oui, comme a-t-il résolu l'affaire ?

Merci d'avance à la ou les personnes qui saura me venir en aide.
Cordialement, Norore.

**Beniou** · 09/02/2011, 17h21

Bonjour,

Je pratique depuisun moment le site du NCBI mais les GEO, je ne connaissais pas. J'ai regardé d'un peu plus près à quoi cela ressemblait mais j'avoue que ce n'est pas du tout ma branche donc je ne vois pas très bien ce que tu veux récupérer qudn tu parles de "données brutes" : les profils ? les datasets ?

Pour obtenir ce que tu veux via le site tu mets une query du style en recherche de profils?

mongene [Gene]

Mets un exemple de requête et de fichiers qui t'intéressent pour que l'on comprenne mieux.

Sinon, le module Bio::DB::Expression peut peut être faire l'affaire mais au vue de la documentation du module ca fait un peu peur...
J'ai vu qu'il y avait un autre module Bio::DB::Expression::geo mais là, la doc est cette fois-ci inexistante...
Vu les exemples du premier module, je ne vois même pas comment faire des requêtes : on a l'impression que le module récupère tout ...

Donc là perso, je ne vois pas vraiment comment arriver à tes fins (sans plus de précisions sur les données brutes etc.) à l'aide de ce ou ces modules.

Peut être via une construction de requêtes comme mentionné sur le site ou après téléchargement mais bon je ne connais le volume de données que cela représente...

Désolé de ne pas faire plus pour le moment

**Norore** · 10/02/2011, 09h15

Bonjour.

Mince, je pensais que les GEO Profiles étaient utilisés de façon régulière en bioinfo. Soit, je ne me laisse pas abattre, je vais tâcher de m'expliquer plus clairement et plus abondamment.

Si l'on se réfère à cette page sur le gène 119 (ADD2), en regardant le menu de droite, on peut repérer, au niveau des Links un lien pointant sur les GEO Profiles relatifs à ce gène. Ce qui intéresse le biologiste avec lequel je travaille (et donc moi aussi par extension), ceux sont les profils sous forme de graphique. Ce que j'entends donc par les données brutes, ceux sont les données relatives à ce graphique, qui me permettront ainsi de repérer les données intéressantes à récupérer (down regulation, grandes différences entre les différentes barres du graphique).
Le problème qui se pose donc est le suivant : comment récupérer ces données sachant que je pars du fichier XML du gène, dont j'extrais certaines données d'intérêt, et que je n'ai aucune idée du numéro d'accession pour les GEO Profiles de ce gène ?
J'ai beau lire attentivement le fichier XML, je ne vois aucun numéro d'accession correspondant aux GEO Profiles et j'ai bien peur que ce soit ce qui me bloque le plus à l'heure actuelle...

Merci pour le lien vers le tutoriel de la construction des requêtes sur GEO, je vais tâcher de voir si je ne peux pas passer par une requête via un autre numéro d'accession

! Oui, la méthode est un peu sioux et tirée par les cheveux, mais je n'ai pas toujours le choix

...

Cordialement.

Edit : je précise juste que je récupère le fichier XML en le téléchargeant via les Eutils grâce à un script Perl.

**Norore** · 10/02/2011, 11h30

Bon, à force de chercher et de tâtonner, je crois que j'ai trouvé ce que je souhaite faire. Je vais encore faire grossir mon disque dur mais dans l'immédiat je ne vois pas de meilleure solution (c'est l'administrateur du serveur qui va être content

!). Voici la solution que j'ai trouvé, ça pourra peut-être servir à quelqu'un d'autre, à coup de script Perl et grâce aux Eutils :

Code perl :

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
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
#!/usr/bin/perl -w                                  -*- Perl -*-
#
# $Id$
#
#
####
# Inspired by: http://www.ncbi.nlm.nih.gov/books/NBK25498/#chapter3.Application_3_Retrieving_large
# Tutorial for use Eutils tool from NCBI
 
use strict;
 
use LWP::Simple;
use LWP::UserAgent;
 
# get GEO Profiles by GeneID access number
&geo_profiles;
 
sub geneid_xml{
 
# initialize some values
    my $gene = $ARGV[0];
    my $geneid = $gene."[uid]";
    my $db = "gene";
 
#assemble the esearch URL
    my $base = 'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/';
    my $url = $base . "esearch.fcgi?db=$db&term=$geneid&usehistory=y";
 
#post the esearch URL
    my $output = get($url);
 
# parse WebEnv, QueryKey and Count (# records retrieved)
    my $web = $1 if ($output =~ /<WebEnv>(\S+)<\/WebEnv>/);
    my $key = $1 if ($output =~ /<QueryKey>(\d+)<\/QueryKey>/);
    my $count = $1 if ($output =~ /<Count>(\d+)<\/Count>/);
 
# open output file for writing
    open OUT, "> $gene.xml" or ( warn "Can't create or rewritte file $gene.xml!\n" and die );
 
# retrieve data in batches of 500
    my $retmax = 500;
    for (my $retstart = 0; $retstart < $count; $retstart += $retmax) {
        my $efetch_url = $base ."efetch.fcgi?db=$db&WebEnv=$web";
        $efetch_url .= "&query_key=$key&retstart=$retstart";
        $efetch_url .= "&retmax=$retmax&rettype=xml&retmode=text";
        my $efetch_out = get($efetch_url);
        print OUT "$efetch_out";
    }
    close OUT;
 
}
 
sub geo_profiles{
 
# initialize some values
    my $gene = $ARGV[0];
    my $geneid = $gene."[Gene]";
    my $db = "gds";
    my $dbFrom = "gene";
 
# assemble the esearch URL
    my $base = 'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/';
    my $url = $base . "esearch.fcgi?db=$db&dbFrom=$dbFrom&term=$geneid&usehistory=y";
 
 
# post the esearch URL
    my $output = get($url);
 
 
# parse WebEnv, QueryKey and Count (# records retrieved)
    my $web = $1 if ($output =~ /<WebEnv>(\S+)<\/WebEnv>/);
    my $key = $1 if ($output =~ /<QueryKey>(\d+)<\/QueryKey>/);
    my $count = $1 if ($output =~ /<Count>(\d+)<\/Count>/);
 
 
# open output file for writing
    open OUT, "> $gene.geo" or ( warn "Can't create or rewritte file $gene.geo!\n" and die );
 
# retrieve data in batches of 500
    my $retmax = 500;
    for (my $retstart = 0; $retstart < $count; $retstart += $retmax) {
        my $efetch_url = $base ."efetch.fcgi?db=$db&WebEnv=$web";
        $efetch_url .= "&query_key=$key&retstart=$retstart&report=docsum";
        $efetch_url .= "&retmax=$retmax&retmode=file";
        my $efetch_out = get($efetch_url);
        print OUT "$efetch_out";
    }
 
    close OUT;
 
    open OUT, "< $gene.geo" or ( warn "Can't read file $gene.geo!\n" and die );
 
    # parsign summaries looking for GEO access number
    my @geo_acc;
    while(<OUT>){
	chomp;
	if ($_ =~ /^\d+: (\w+) / ){
	    my $geo = $1;
	    push(@geo_acc, $geo);
	}
    }
 
    close OUT;
 
    # download files for each GEO access number stored and gunzip them
    my $path_geo = "/home/nolwenn/Bureau/GEO/";
    chdir($path_geo) or ( warn "Can't change directory to $path_geo!\n" and die );
    foreach my $acc (@geo_acc){
	my $request = "ftp://ftp.ncbi.nih.gov/pub/geo/DATA/SOFT/GDS_full/".$acc."_full.soft.gz";
	system("wget ".$request);
    }
    system("gunzip *.gz");
 
}

**Beniou** · 10/02/2011, 13h24

C'est vrai que la méthode du "télécharger tout puis analyser" est un bon moyen à court terme. En plus, une fois que les données sont en local on peut faire plein de choses en plus.

Par contre, juste une mise en garde vu que cela m'est déjà arrivé : quelle est la pérennité de tes données ? Quid des mises à jour ? Durée de vie des analyses qui en résultent ?
Je dis cela parce que dans les bases de données bio (gènes ou SNPs par exemple) sont mises à jour périodiquement,ce qui est bien mais le problème vient que certaines entrées sont supprimées (ca ca va), plusieurs entrées sont "mergées" en 1 seule, d'autres changent carrément d'identifiant etc. et j'en passe. Et ce, sans historique des modifications comme ca l'est dans certaines bases. Du coup, c'est assez difficile de mettre à jour ses données locales. On garde alors plusieurs versions. Ou alors peut être utiliser des logiciels comme biomaj ...

Donc voilà, c'était juste pour préciser ce point qui me paraît important quand on utilise des données locales que l'on télécharge.

**Norore** · 10/02/2011, 13h51

Oui, je suis au courant de ces difficultés majeures, le problème principal réside dans le fait que j'utilise plusieurs bases de données, récupérées en local, que je fais correspondre les unes aux autres, du coup certaines bases de données sont mise à jour quotidiennement, quand c'est précisé ce qui est exceptionnel, d'autres sont vieilles de plusieurs jours à mois (quand elles n'ont pas carrément 2 ans).
Bref, c'est encore un casse-tête de plus sans lequel nous serions au chômage

!
Merci de me prévenir pour les pertes de données suite aux mises à jour, je vais tâcher d'y prêter attention

!

**Beniou** · 10/02/2011, 14h28

Envoyé par Norore

Bref, c'est encore un casse-tête de plus sans lequel nous serions au chômage

Cela fait partie de notre travail il est vrai et est indispensable. Heureusement, une fois que ce travail pénible est passé (du moins jusqu'à de nouvelles bases à utiliser ou à mettre à jour), il nous reste toute la partie analyse-recherche qui, pour ma part, est l'essence même de notre travail.

**Norore** · 10/02/2011, 16h15

Bon, j'ai laissé tomber la première solution que j'ai trouvé, en voyant la tête des fichiers, mes yeux ont failli sortir de leurs orbites

!
Je garde le script sous le coude mais j'en ai réécrit un qui lui est plus proche de ce que je cherche à obtenir et qui, je pense, devrait plaire au biologiste. Ce sera un peu plus léger, pas très optimal parce que du coup je vais avoir plein de petits fichiers, mais je devrais m'y retrouver plus facilement

!

Ma fonction a maintenant cette tête là à peu près :

Code perl :

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
sub geo_profiles_dat{
 
# initialize some values
    my $gene = $ARGV[0];
    my $geneid = $gene."[Gene]";
    my $db = "geo";
    my $dbFrom = "gene";
 
# assemble the esearch URL
    my $base = 'http://eutils.ncbi.nlm.nih.gov/entrez/eutils/';
    my $url = $base . "esearch.fcgi?db=$db&dbFrom=$dbFrom&term=$geneid&usehistory=y";
 
 
# post the esearch URL
    my $output = get($url);
 
 
# parse WebEnv, QueryKey and Count (# records retrieved)
    my $web = $1 if ($output =~ /<WebEnv>(\S+)<\/WebEnv>/);
    my $key = $1 if ($output =~ /<QueryKey>(\d+)<\/QueryKey>/);
    my $count = $1 if ($output =~ /<Count>(\d+)<\/Count>/);
 
 
# open output file for writing
    open OUT, "> $gene.gds.dat" or ( warn "Can't create or rewritte file $gene.geo!\n" and die );
 
# retrieve data in batches of 500
    my $retmax = 500;
    for (my $retstart = 0; $retstart < $count; $retstart += $retmax) {
        my $efetch_url = $base ."efetch.fcgi?db=$db&WebEnv=$web";
        $efetch_url .= "&query_key=$key&retstart=$retstart&report=docsum";
        $efetch_url .= "&retmax=$retmax&retmode=file";
        my $efetch_out = get($efetch_url);
        print OUT "$efetch_out";
    }
 
    close OUT;
 
    open OUT, "< $gene.gds.dat" or ( warn "Can't read file $gene.gds.dat!\n" and die );
 
    # parsign summaries looking for GEO access number
    my %geo_acc;
    while(<OUT>){
    	chomp;
    	if ($_ =~ /^\d+: (\w+) record \| \w+ (\w+) / ){
    	    my $geo = $1;
	    my $ref = $2;
    	    $geo_acc{$geo} = $ref;
    	}
    }
 
    close OUT;
 
    # download files for each GEO access number stored and gunzip them
    my $path_geo = "/home/nolwenn/Bureau/GEOSamples/";
    chdir($path_geo) or ( warn "Can't change directory to $path_geo!\n" and die );
    while ( my ($acc, $ref) = each (%geo_acc) ){
	my $file = $gene."_".$acc."_".$ref.".txt";
	if ( !$file ){
	    my $request = "http://www.ncbi.nlm.nih.gov/geo/gds/getDatum.cgi?ID=".$acc.":".$ref."&type=text";
	    system("wget ".$request);
	    system("mv getDatum.cgi?ID=".$acc.":".$ref."&type=text ".$file);
	}
	else{
	    print "File $file already exists in directory $path_geo.";
	}
    }
 
}

Il ne me reste plus qu'à trouver comment comparer les dates des fichiers et ça devrait rouler

!

**Beniou** · 10/02/2011, 16h53

Ok.

Sinon une précision en regardant ton script : lorsque tu récupères à la fin tes fichiers avec le wget, cela utilise les eutils ?

Code :

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

1
2
3
 
my $request = "http://www.ncbi.nlm.nih.gov/geo/gds/getDatum.cgi?ID=".$acc.":".$ref."&type=text";
system("wget ".$request);

Parce qu'il faut savoir que le site du NCBI n'aime pas trop qu'on les "harcèle" de requêtes type wget. Si mes souvenirs sont bons, dans leur "bonnes pratiques", ils définissaient un maximum de 3 requêtes par seconde pour cette méthode.
Je peux me tromper mais vérifie quand même : ce serait dommage de se faire bloquer son IP parce qu'on a abusé un tout petit peu

**Norore** · 10/02/2011, 17h23

Je ne pense pas que cette commande passer par les Eutils, je n'ai trouvé cet URL qu'à force de fouiner (et en regardant un peu leur source HTML qui insère un tableau par iframe) et de trouver comment leur système fonctionne pour afficher ces données.
J'ai eu beau fouiller les informations concernant les Eutils, je n'ai rien trouvé de concluant pour cette méthode, ou alors je suis passée à côté...

Je vais continuer de chercher mais je crains de devoir leur demander par Email s'ils ont une solution à me proposer. A moins que je ne passe par leur module R, puis que je trouve comment convertir un fichier R afin de l'utiliser en Perl...

**Norore** · 11/02/2011, 12h02

Finalement je vais laisser tomber la seconde méthode, ça évitera les râleries (à raison) du NCBI. De plus, j'aurais bien plus d'informations et je pourrais donc parcourir plus aisément les fichiers récupérés.
Seul défaut, découvert ce matin, si le NCBI fait une maintenance sur les Eutils, le script plante

!

En tout cas, merci pour le renseignement Beniou

!

Comment récupérer les données brutes des GEO Profiles d'un gène donné ?

Bioinformatique Perl

Discussions similaires

Partager

Partager