Tag 57: Endlich mal wieder scripten

Eigentlich hatte ich etwas ganz anderes vor, aber mich haben Emails meiner Kollegen erreicht, die gerade an einem Problem mit unseren Daten arbeiten. Nachdem das entfernt auch etwas mit meinen Analysen zu tun hatte, habe ich also angefangen, ein Script zu schreiben, das die Daten auswertet. 

Ich versuche das mal zu erklären: 

Ich habe DNA-Sequenzen. Das sind ca 100 Basen lange "Wörter", die zum Großteil aus den Basen A, T, C und G bestehen. Es gibt auch Lücken und unklare Stellen, wo mehrere der Basen möglich sind. Von diesen Sequenzen habe ich 44.000 von verschiedenen Stellen im Genom und von bis zu 186 verschiedenen Individuen (Pflanzen, aber von der gleichen Art, Adonis). Wenn ich jetzt eine der 44.000 Stellen im Genom anschaue, sind die sich ähnlich, können sich aber an mehreren Basen unterscheiden.

Beispiel:

Pflanze1: AAATGCCG-T

Pflanze2: AATTGCCGAT

Pflanze3: AAATGCCG-T

Pflanze4: AAATCCCGAT

Pflanze5: AATTGCCGAT

Das sieht erst mal ähnlich aus, ist aber leicht unterschiedlich.

Was ich jetzt wissen will: Wie oft kommt das exakt gleiche "Wort " vor? Im Beispiel gibt es eins der "Wörter" zweimal, falls Du suchen magst.

In der Realität habe ich pro Individuum zwei dieser Wörter pro Stelle im Genom, weil meine Art so wie Menschen diploid ist, also zwei Kopien hat. Diese Kopien können identisch oder verschieden sein.

Mein Script geht jede in den Daten vorhandene Stelle im Genom durch, sucht nach identischen Sequenzen, ordnet die identischen Sequenzen einer Nummer zu und trägt das dann in eine Masterliste ein, wo alle 44.000 Stellen sich in je einer Spalte der Tabelle niederschlagen.

Beispiel von oben:

Pflanze1: 1

Pflanze2: 2

Pflanze3: 1

Pflanze4: 3

Pflanze5: 2

Ich habe 12 h daran gesessen und es sind eh nur 10 Zeilen, aber ich bin sehr stolz auf mich. Das ist immerhin mein erstes bash shell, das ich selbst geschrieben habe. Etwas suboptimal ist, dass es auf meinem lahmen Rechner läuft, aber was soll man machen. Kurz vorm Einschlafen ist mir aufgefallen, dass ich einen Fehler gemacht habe, aber das werde ich dann später korrigieren.

Die Laufzeit ist sicher auch nicht optimal. In R wüsste ich, worauf ich besonders achten müsste. Zum Beispiel sollte man da eher über Reihen als Spalten arbeiten, weil das schneller geht. Und man erstellt besser gleich die komplette Liste und füllt dann Werte ein als ständig was dazu zu kopieren. Das dazu kopieren erstellt nämlich jedes Mal eine neue Tabelle, was bei 44.000 mal 186 mal 2 Einträgen schon auf die Rechenzeit geht. Aber ich bin kein Informatiker, meine Lösung ist für meine Verhältnisse schon sehr gut.

Das Script ist auch nach 12 h noch nicht fertig gelaufen.

Meine Gang-Pflanzen

Es gab Toast Hawaii zu Mittag. Und abends hat Nikolai Eier gekocht. Ich glaube, er wollte mich besänftigen, nachdem er es gewagt hat, meine Schokolade zu essen, die er gekauft hat.