Cet article Blue Collar Bioinformatics; Récupération automatisée des données d'expression avec Python et R C'est une traduction. C'est un peu vieux, mais je peux aussi étudier pour moi-même. Nous nous félicitons de toute erreur typographique ou erreur de traduction.
Traduit ci-dessous
En février, je vais au Japon pour participer au Bio-Hackason 2010, qui se concentre sur l'intégration des données de biologie. Donc, je pensais à un cas d'utilisation en guise de préparation.
Lors de la lutte avec un ensemble de données, quelles données supplémentaires aideraient à résoudre le problème? Un exemple typique est l'information sur les données de niveau d'expression génique. Si vous avez une liste de gènes qui vous intéressent et que vous souhaitez connaître leurs niveaux d'expression.
Dans un monde où le rêve du web sémantique est devenu réalité, je pense que vous interrogerez un ensemble de triplets RDF comme celui-ci:
Ensuite, le niveau d'expression de la tumeur cellulaire (proB) d'intérêt dans des conditions normales peut être obtenu, puis le niveau d'expression du gène à portée de main peut être obtenu. Idéalement, je voudrais que le niveau d'expression défini ait des métadonnées. Ensuite, vous pouvez voir ce que signifie le nombre 7,65 dans le contexte de l'expérience.
Vous ne pouvez pas lancer cette requête à partir d'un canal sur le Web sémantique, mais pour résoudre automatiquement des problèmes similaires. Il existe de nombreux outils.
Gene Expression Omnibus (GEO) de NCBI héberge des données de niveau d'expression en ligne. Vous pouvez interroger l'ensemble de données avec l'API Eutils de Biopython. Une fois que vous avez identifié l'ensemble de données qui vous intéresse, récupérez les données sur [GEOQuery] de Bioconductor (http://www.bioconductor.org/packages/bioc/html/GEOquery.html) et correspondent à l'identifiant UCSC RefGene. Vous pouvez. Tous ces éléments peuvent être liés en exécutant R à partir de python dans Rpy2.
Traduction: Maintenant peut-être que
pipeRest mieux
Tout d'abord, fixez-vous des objectifs de haut niveau. Ici, nous recherchons des données de niveau d'expression de type sauvage des cellules proB de souris.
def main():
organism = "Mus musculus"
cell_types = ["proB", "ProB", "pro-B"]
email = "[email protected]"
save_dir = os.getcwd()
exp_data = get_geo_data(organism, cell_types, email, save_dir)
def _is_wild_type(result):
"""Juger d'après le titre si l'échantillon est de type sauvage"""
return result.samples[0][0].startswith("WT")
Il envoie ensuite une requête à GEO pour obtenir les données disponibles en fonction de l'espèce et de la tumeur cellulaire. Dans le vrai travail, vous pourriez remplacer votre fonction privée _is_wild_type par une autre qui agit de manière interactive, laissant votre bioscientifique choisir en fonction de votre méthode expérimentale.
Notez également que le résultat final est pickle et sauvegardé sous la main. En n'interrogeant qu'une seule fois de cette manière, la charge sur le serveur externe est réduite.
def get_geo_data(organism, cell_types, email, save_dir, is_desired_result):
save_file = os.path.join(save_dir, "%s-results.pkl" % cell_types[0])
if not os.path.exists(save_file):
results = cell_type_gsms(organism, cell_types, email)
for result in results:
if is_desired_result(result):
with open(save_file, "w") as out_handle:
cPickle.dump(result, out_handle)
break
with open(save_file) as save_handle:
result = cPickle.load(save_handle)
Le défi d'envoyer une requête au GEO et d'obtenir les résultats est fait par l'interface Entrez de Biopython. Après avoir pompé la requête, le résultat est analysé dans un objet simple avec une description, contenant le niveau d'expression associé au titre et à l'ID de chaque échantillon.
def cell_type_gsms(organism, cell_types, email):
"""Obtenez GEO pour des espèces spécifiques et des tumeurs cellulaires en utilisant Entrez.
"""
Entrez.email = email
search_term = "%s[ORGN] %s" % (organism, " OR ".join(cell_types))
print "Searching GEO and retrieving results: %s" % search_term
hits = []
handle = Entrez.esearch(db="gds", term=search_term)
results = Entrez.read(handle)
for geo_id in results['IdList']:
handle = Entrez.esummary(db="gds", id=geo_id)
summary = Entrez.read(handle)
samples = []
for sample in summary[0]['Samples']:
for cell_type in cell_types:
if sample['Title'].find(cell_type) >= 0:
samples.append((sample['Title'], sample['Accession']))
break
if len(samples) > 0:
hits.append(GEOResult(summary[0]['summary'], samples))
return hits
Avec cela, si vous choisissez le type d'expérience, vous franchirez la première barrière. Ensuite, obtenez la valeur des données au niveau de l'expression. Cela se fait en créant un dictionnaire qui mappe les ID RefGene aux niveaux d'expression. En d'autres termes, le résultat est:
WT ProB, biological rep1
[('NM_177327', [7.7430266269999999, 6.4795213670000003, 8.8766985500000004]),
('NM_001008785', [7.4671954649999996, 5.4761453329999998]),
('NM_177325', [7.3312364040000002, 11.831475960000001]),
('NM_177323', [6.9779868059999997, 6.3926399939999996]),
('NM_177322', [5.0833683379999997])]
Le vrai problème est d'obtenir le niveau d'expression et de le mapper à l'ID du gène, ce que la bibliothèque GEOquery de Bioconductor fait pour vous. La table de correspondance d'expression génique et les métadonnées d'ordre supérieur sont stockées dans un fichier local. Le premier est un tableau délimité par des tabulations de style R, et le second est conservé dans JSON. Enregistrez à la fois localement et dans un format facile à lire.
def write_gsm_map(self, gsm_id, meta_file, table_file):
"""Les données de niveau d'expression GEO sont obtenues à l'aide de Bioconductor et stockées dans une table.
"""
robjects.r.assign("gsm.id", gsm_id)
robjects.r.assign("table.file", table_file)
robjects.r.assign("meta.file", meta_file)
robjects.r('''
library(GEOquery)
library(rjson)
gsm <- getGEO(gsm.id)
write.table(Table(gsm), file = table.file, sep = "\t", row.names = FALSE,
col.names = TRUE)
cat(toJSON(Meta(gsm)), file = meta.file)
''')
Cette fonction crée une table de correspondance entre le nom de la sonde et le niveau d'expression. Pour le rendre plus pratique, la sonde du microréseau d'expression doit être appariée à l'ID du gène biologique. Ceci est fait en appelant à nouveau la bibliothèque GEO. Les données sont à nouveau récupérées sous la forme d'un bloc de données R et seuls les éléments d'intérêt sont enregistrés dans un fichier délimité par des tabulations.
def _write_gpl_map(self, gpl_id, gpl_file):
"""Utilisez R pour obtenir les données de la plateforme GEO et les enregistrer au format tabulaire"""
robjects.r.assign("gpl.id", gpl_id)
robjects.r.assign("gpl.file", gpl_file)
robjects.r('''
library(GEOquery)
gpl <- getGEO(gpl.id)
gpl.map <- subset(Table(gpl), select=c("ID", "RefSeq.Transcript.ID"))
write.table(gpl.map, file = gpl.file, sep = "\t", row.names = FALSE,
col.names = TRUE)
''')
Enfin, je résumerai le traitement jusqu'à présent. Téléchargez et analysez chaque fichier de table de correspondance et connectez la correspondance avec le niveau d'expression-> ID de sonde-> ID de transcription. Le dictionnaire final aura l'ID du gène associé au niveau d'expression, il sera donc renvoyé.
def get_gsm_tx_values(self, gsm_id, save_dir):
""""Produit de transcription"=> 「GEO,Créer une table de correspondance du "niveau d'expression obtenu à partir du fichier GSM"
"""
gsm_meta_file = os.path.join(save_dir, "%s-meta.txt" % gsm_id)
gsm_table_file = os.path.join(save_dir, "%s-table.txt" % gsm_id)
if (not os.path.exists(gsm_meta_file) or \
not os.path.exists(gsm_table_file)):
self._write_gsm_map(gsm_id, gsm_meta_file, gsm_table_file)
with open(gsm_meta_file) as in_handle:
gsm_meta = json.load(in_handle)
id_to_tx = self.get_gpl_map(gsm_meta['platform_id'], save_dir)
tx_to_vals = collections.defaultdict(list)
with open(gsm_table_file) as in_handle:
reader = csv.reader(in_handle, dialect='excel-tab')
reader.next() # header
for probe_id, probe_val in reader:
for tx_id in id_to_tx.get(probe_id, []):
tx_to_vals[tx_id].append(float(probe_val))
return tx_to_vals
Le script complet (https://github.com/chapmanb/bcbb/blob/022c50b3132bc3d009ed1d533a3954bd564a9ed3/rest_apis/find_geo_data.py) est une intégration des anciens extraits. Nous avons commencé avec la requête et avons finalement obtenu les données de niveau d'expression requises.
Cependant, c'est devenu un détour par rapport à l'idéal en RDF. Ainsi, la volonté de simplifier au maximum l'accès aux données questionnées est un défi récurrent pour les créateurs d'outils bioinformatiques.
La traduction se termine ici
Le script lui-même est https://github.com/chapmanb/bcbb/blob/022c50b3132bc3d009ed1d533a3954bd564a9ed3/rest_apis/find_geo_data.py C'est dedans.
Le référentiel de l'auteur bcbb contient de nombreux autres programmes utiles pour la bioinformatique.