Manchmal macht die Arbeit einfach so Spaß, dass ich jedem davon erzählen muss: Dem armen Menschen der den Fehler machte in die grobe Richtung meines Posters auf einer Konferenz zu schauen, dem Immunologen, den ich im Fahrstuhl treffe, und den Lesern dieses Blogs.
Nebenstehendes Bild zeigt eine Auswertung während meiner momentanen Arbeit, dem Assembley (“Zusammenpuzzeln”) und der erste Analyse aller in meinem Studien-Organismus exprimierten Gene. Die Sequenzen werden von Roche’s “454 Genome Sequenzer” generiert (dazu später mehr). Diese Analyse ist, da ich vor wenigen Monaten erst begonnen habe Programmieren zu Lernen bisher noch ein Spiel mit den Daten:
Eintauchen und ein paar Züge darin Schwimmen, mit einem selbst gebastelten Boot aus Perl-Code in See stechen, BLAST als Außenborder anbringen und dann mit R als Fernrohr nach neuen Ufern suchen.
Was mich hier poetisch werden lässt möchte ich kurz erklären: Die DNA-Sequenzen wurden wie gesagt mit einem Sequenziergerät der sogenannten zweiten Generation (Tobias schreibt hier über die dritte Generation) generiert. Bei dem Datensatz vor dem ich gerade sitze, handelt sich um 100,000 durchschnittlich jeweils 250 Basen lange DNA-Abschnitte1. Um nur die Protein-codierenden Gene zu sequenzieren, nicht den bedeutend größeren nicht-codierenden Teil des Genoms habe ich zum Sequenzieren sogenannte cDNA hergestellt. Man bedient sich dazu einer reversen Transcriptase, die aus einem RNA-Virus stammt und den Prozess der Transkription eben “rückwärts”, also von RNA zu DNA, katalysieren kann.
Nun sind natürlich die meisten Gene länger als 250 Basenpaare. Außerdem liegen sie unterschiedlich stark exprimiert vor, was sich in der cDNA und damit auch in den Daten wiederspiegelt. So habe ich vom am stärkstem exprimierten Gen, der Cytochrome c Oxidase Untereinheit I (CO1), gleich ~4000 (~4%) Sequenzen. Das Gen ist ungefähr 1400 Basen lang und so überschneiden sich natürlich die Sequenzen. Ziel des Assembleys ist es nun anhand dieser Überschneidungen eine gemeinsame Sequenz zu bilden. Dies geschieht im Moment mit der von Roche zu ihrem Sequenzierer mitgelieferten Software (mehr schlecht als recht).
Nun habe ich also die Sequenzen. Eine Datei mit etwa 4000 sogenannten Contigs (aus mehreren Einzelstückchen zusammen-gepuzzelte Sequenzen), eine Datei mit den 100000 “Rohsequenzen” und eine Datei, die mir sagt ob eine Rohsequenz in einem Contig verwendet wurde. Ich Filtere dann mit einem kleinen Perl-script die Einzelsequenzen heraus, die nicht in ein Contig einflossen, und packe sie mit diesen zusammen in eine gemeinsame Datei.
Diese Datei blaste ich dann gegen 4 lokal auf meinem Rechner installierte Datenbanken, das funktioniert ähnlich wie von Alex beschrieben, ist aber schneller da man gegen Datenbanken von geringerer Größe sucht 2. Ein weiterer Unterschied besteht darin, dass man verschiedene Blast-Programme von der Linux/UNIX “command-line” ausführen und damit auch in einem kleinen Perl-Programm automatisieren kann. Die vier Datenbanken sind in meinem Fall Nematoden rRNA, Fisch rRNA, Nematoden-Protein und Fisch-Protein. Die ersten beiden Datenbanken verwende ich um den Datensatz nach rRNA (nicht translatierter, ribosomaler RNA) zu untersuchen, die eigentlich unerwünscht ist aber leider als Verunreinigung in meinen cDNA Präparationen auftaucht. Die beiden anderen Datenbanken verwende ich aus folgendem Grund: Eigentlich möchte ich einen Nematoden (Rundwurm) sequenzieren, der als Parasit im Aal lebt. Nun sind diese Würmer aber leider mit Blut ihres Wirtes gefüllt.
Es kann also durchaus zu Verunreinigungen der Daten durch Wirts-Gene kommen. Zwar wird RNA sehr schnell abgebaut, was mir in dem Fall der Verdauung des Aalblutes ausnahmsweise recht ist, mit einer Rest-Verunreinigung ist trotzdem zu rechnen.
Also sammle ich die mehrere hunderttausend Zeilen lange Ausgabe der Blast-Suchen gegen die verschiedenen Datenbanken ein und vergleiche jeweils die Ergebnisse. Natürlich mache ich das nicht von Hand sondern das macht mein Programm. Die Ergebnisse gibt es dann in einer Komma-separierten Datei (.csv) aus. Diese Dateien öffne ich dann in R und produziere die beiden hier abgebildeten Plots.
In diesem ersten Kreisdiagramm wird ein Problem deutlich: Eine sehr hohe rRNA Kontamination (die Daten stammen noch aus einem Probelauf mit nur 9000 Sequenzen).
Die Lösung des Problems ist, nachdem wir an der Molekularbiologie (der Präparation der cDNA) erfolglos getüftelt haben einfach: Mehr Sequenzieren. Das zweite, wesentlich interessantere Problem sind die Wirts-Gene. Mit einem Anteil von 4% recht viel Konatmination, die sicher als solche identifiziert werden muss. Falsche Identifikation wäre ärgerlich, da Nematoden Gene auch durchaus “wirtsähnlich” sein können, nämlich dann wenn sie mit dem Immunsystem interagieren. Diese Gene gehören für mein Projekt zu den interessantesten.
Hier kommt der zweite bereits in der Einleitung abgebildete Plot ins Spiel: Eine typische Eigenschaft einer Art ist der GC-Gehalt ihrer DNA, d.h. die relative Häufigkeit der komplementären Basen G-C im Gegensatz zu A-T. Diese Häufigkeit variiert (näherungsweise Normal-verteilt) über das Genom verteilt um einen art-typischen Mittelwert.
Die Grafik zeigt nun diese Verteilung. Auf der X-Achse findet man den GC Gehalt in %, auf der Y-Achse die prozentuale Häufigkeit der Beobachtung getrennt für die Sequenzen, die Treffer zu den verschiedenen Datenbanken hatten. Man erkennt: Die eindeutig als Nematoden Gene identifizierten Sequenzen (knallrot) haben im Mittel ~39% GC, die als Fisch-ähnlich identifizierten Gene haben einen GC Mittelwert von ~51%. Für mich sehr beruhigend ist die Verteilung der Sequenzen ohne Blast-Ähnlichkeit (dunkelbraun), sie entspricht in ihrer GC-Verteilung etwa der Verteilung der Nematoden-Gene. Wäre ein großer Teil Fisch-Gene unidentifiziert, hätte die Verteilung beim Fisch-Mittelwert eine kleine Beule. Wahrscheinlich handelt es sich bei den unidentifizierten Genen um neue Nematoden Gene, mit wenig Ähnlichkeit zu bereits in Datenbanken vorhandenen Sequenzen. Dies ist durchaus nichts ungewöhliches in Nematoden. Vertebraten zeigen dagegen einen wesentlich geringeren “Erfindungsreichtum” was neue Gene angeht und man findet fast immer Orthologe in anderen Arten.
Im Gegensatz dazu findet man in der GC-Verteilung der als Fisch ähnlich identifizierten Gene (hellblau) eine Beule an der Stelle, an der der Mittelwert für Nematoden-Gene liegt. Wahrscheinlich wurden einige an den Fischwirt adaptierte Nematoden Gene als Fisch-Gene fehl identifiziert. Ich muss mir die Sequenzen die Blast-Treffer zu Fisch Genen erbrachten und einen niedrigen GC-Gehalt (um 40%) haben nochmal sehr genau anschauen.
Mein Programm produziert, nachdem es die Blast-Suchen durchgeführt hat, indem es dann ein R-script schreibt und dieses dann im R BATCH Mode ausführt, die beiden Plots (und einigen anderen nützlichen Output) automatisch. Sobald ich das Programm soweit fertig gestellt habe, dass es dem Benutzer das Definieren von eigenen Datenbanken erlaubt und damit portabler wird, werde ich nochmal darüber schreiben und dann einen Link zum Programm und einem Manual bereitstellen.
1 Es handelt sich um eine Viertel Picoliter-Platte und noch um die alte “FLX-Chemie” die nächsten Daten (ich bekomme noch 9 weitere Datensätze) werden für den gleichen Preis mit der “Titanium-Chemie” länger (450 Basen) und umfangreicher (200.000 reads auf einer Viertel- Platte).
2 Alle Blast-Suchen brauchen auf meinem Laptop trotzdem noch sehr lange: Etwa 1 Woche bis mein Script gegen die 4 Datenbanken fertig ist. Daher auch noch die alten Plots aus dem Probelauf. Blast Suchen brauchen Prozessor Power, ich muss das ganze hier auf 8 Prozessoren (einem 2x Quad Core mit 16MB GB RAM) laufen lassen, dann hab ich meine Plots über Nacht.
Kommentare (7)