Dieser Beitrag Blue Collar Bioinformatics; Automatisches Abrufen von Expressionsdaten mit Python und R Es ist eine Übersetzung. Es ist ein bisschen alt, aber ich kann auch für mich selbst lernen. Wir freuen uns über Tippfehler oder Fehlübersetzungen.
Übersetzt unten
Diesen Februar gehe ich nach Japan, um am Bio-Hackason 2010 teilzunehmen, der sich auf die Integration biologischer Daten konzentriert. Also dachte ich über einen Anwendungsfall als Vorbereitung nach.
Welche zusätzlichen Daten würden beim Ringen mit einem Datensatz zur Lösung des Problems beitragen? Ein typisches Beispiel sind Informationen zu Daten auf Genexpressionsniveau. Wenn Sie eine Liste von Genen von Interesse haben und deren Expressionsniveaus wissen möchten.
In einer Welt, in der der Traum vom Semantic Web wahr geworden ist, werden Sie wahrscheinlich eine Reihe von RDF-Tripeln wie diese abfragen:
Dann können Sie das Expressionsniveau des interessierenden Zelltumors (proB) unter normalen Bedingungen erhalten, und dann können Sie das Expressionsniveau des vorliegenden Gens erhalten. Idealerweise möchte ich, dass die Ausdrucksebene Metadaten enthält. Dann können Sie im Kontext des Experiments sehen, was die Zahl 7.65 bedeutet.
Sie können diese Abfrage nicht aus einem Kanal im Semantic Web werfen, sondern ähnliche Probleme automatisch lösen. Es gibt viele Werkzeuge.
Der Gene Expression Omnibus (GEO) von NCBI hostet Expressionsdaten online. Sie können den Datensatz mit der API Eutils von Biopython abfragen. Wenn Sie den Datensatz identifiziert haben, an dem Sie interessiert sind, rufen Sie die Daten unter Bioconductor's GEOQuery ab und entsprechen Sie der UCSC RefGene-ID. Sie können. All dies kann durch Ausführen von R aus Python in Rpy2 verknüpft werden.
Übersetzung: Jetzt ist
pipeRvielleicht besser
Setzen Sie sich zunächst Ziele auf höchster Ebene. Hier suchen wir nach Wildtyp-Expressionsdaten von Maus-proB-Zellen.
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):
"""Nach dem Titel zu urteilen, ob es sich bei der Stichprobe um einen Wildtyp handelt"""
return result.samples[0][0].startswith("WT")
Anschließend wird eine Abfrage an GEO gesendet, um die verfügbaren Daten basierend auf der Art und dem Zelltumor abzurufen. In der Praxis können Sie Ihre private Funktion "_is_wild_type" durch eine Funktion ersetzen, die interaktiv wirkt und Ihren Biowissenschaftler nach Ihrer experimentellen Methode auswählen lässt. Beachten Sie auch, dass das Endergebnis pickle ist und sofort gespeichert wird. Durch nur einmaliges Abfragen auf diese Weise wird die Belastung des externen Servers verringert.
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)
Die Herausforderung, eine Abfrage an das GEO zu senden und die Ergebnisse zu erhalten, wird von der Entrez-Oberfläche von Biopython erledigt. Nach dem Aufpumpen der Abfrage werden die Ergebnisse in ein einfaches Objekt mit einer Beschreibung analysiert, die die Ausdrucksebene zusammen mit dem Titel und der ID für jede Probe enthält.
def cell_type_gsms(organism, cell_types, email):
"""Erhalten Sie GEO für bestimmte Arten und Zelltumoren mit 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
Wenn Sie die Art des Experiments auswählen, durchbrechen Sie damit die erste Barriere. Als nächstes erhalten Sie den Wert der Ausdrucksebenendaten. Dazu wird ein Wörterbuch erstellt, das RefGene-IDs Ausdrucksebenen zuordnet. Mit anderen Worten, das Ergebnis ist:
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])]
Der eigentliche Aufwand besteht darin, das Expressionsniveau zu ermitteln und es der Gen-ID zuzuordnen, die die GEOquery-Bibliothek von Bioconductor für Sie erledigt. Die Genexpressionskorrespondenztabelle und Metadaten höherer Ordnung werden in einer lokalen Datei gespeichert. Ersteres ist eine tabulatorgetrennte Tabelle im R-Stil, und letzteres befindet sich in JSON. Speichern Sie sowohl lokal als auch in einem einfach zu lesenden Format.
def write_gsm_map(self, gsm_id, meta_file, table_file):
"""Daten zum GEO-Expressionsniveau werden unter Verwendung von Bioconductor erhalten und in einer Tabelle gespeichert.
"""
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)
''')
Diese Funktion erstellt eine Korrespondenztabelle zwischen Sondennamen und Ausdrucksebenen. Um es bequemer zu machen, sollte die Sonde des Expressionsmikroarrays an die ID des biologischen Gens angepasst werden. Dies erfolgt durch erneutes Aufrufen der GEO-Bibliothek. Die Daten werden erneut in Form eines R-Datenrahmens abgerufen, und nur die interessierenden Elemente werden in einer durch Tabulatoren getrennten Datei gespeichert.
def _write_gpl_map(self, gpl_id, gpl_file):
"""Verwenden Sie R, um GEO-Plattformdaten abzurufen und in Tabellenform zu speichern"""
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)
''')
Abschließend werde ich die bisherige Verarbeitung zusammenfassen. Laden Sie jede Korrespondenztabellendatei herunter, analysieren Sie sie und verbinden Sie die Korrespondenz mit der Ausdrucksebene-> Sonden-ID-> Transkript-ID. Das endgültige Wörterbuch enthält die ID des Gens, das dem Expressionsniveau zugeordnet ist, und wird daher zurückgegeben.
def get_gsm_tx_values(self, gsm_id, save_dir):
""""Produkt transkribieren"=> 「GEO,Erstellen Sie eine Korrespondenztabelle mit "Ausdrucksebene aus GSM-Datei erhalten"
"""
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
Das vollständige Skript (https://github.com/chapmanb/bcbb/blob/022c50b3132bc3d009ed1d533a3954bd564a9ed3/rest_apis/find_geo_data.py) ist eine Integration der alten Snippets. Wir haben mit der Abfrage begonnen und schließlich die erforderlichen Daten auf Ausdrucksebene erhalten.
Es ist jedoch ein Umweg im Vergleich zum idealen mit RDF geworden. Der Wunsch, den Zugriff auf fragegesteuerte Daten so einfach wie möglich zu gestalten, ist daher eine wiederkehrende Herausforderung für die Entwickler von Bioinformatik-Tools.
Die Übersetzung endet hier
Das Skript selbst lautet https://github.com/chapmanb/bcbb/blob/022c50b3132bc3d009ed1d533a3954bd564a9ed3/rest_apis/find_geo_data.py Es ist in.
Das Repository des Autors bcbb enthält viele andere nützliche Programme für die Bioinformatik.
Recommended Posts