Gensequenz auslesen - Vernünftiger Ansatz? (Prototyp)

[ne0h]

Hi,

nun denn, gehts mal an etwas eher praktisch/technisches.

Folgendes:

Gegeben ist eine Textdatei mit einem Stück Gen-Sequenz wie z.B.:

Code: Alles auswählen

ACCGACGATGACAGCGATAGC

Nun will ich die Sequenz aus der Datei lesen und dazu das Komplement erzeugen, also zu jedem "A" ein "T", zu jedem "C" ein "G" und natürlich umgekehrt. Dazwischen sollen "Pipes" (|) als Verbindungsglied auf der Konsole angezeigt werden (was ich mit der Zeile...

Code: Alles auswählen

print len(sequence) * "|"

...ausführe) .


Mein Ansatz sieht dazu aus wie folgt:

Code: Alles auswählen

#!/usr/bin/env python
# -*- coding: utf-8 -*-

try:
	data_obj = file("sequenz.dna", "r")
	sequence = data_obj.readline()
	sequence.strip()
	if sequence.isupper():
		print sequence
		print len(sequence) * "|"
		compl_sequence = ""
		for i in sequence:
			if i == "A":
				compl_sequence += "T"
			elif i == "C":
				compl_sequence += "G"
			elif i == "T":
				compl_sequence += "A"
			elif i == "G":
				compl_sequence += "C"
		print compl_sequence
	else:
		print "Sequenz ist falsch implementiert."
except:
	dat_obj = None
	print "Datei konnte nicht geoeffnet werden."

Das ganze soll nur als kleines "Skript" dienen und erstmal relativ minimalistisch aufgebaut sein.

Was haltet Ihr von meinem Ansatz?

Wo habe ich schlechten Code geschrieben?


Gruß

ne0h

[CM]

Hoi,

ich will Dich ja nicht enttäuschen, aber das geht viel einfacher:

Code: Alles auswählen

In [1]: from Bio.Seq import Seq

In [2]: seq = Seq("ACCGACGATGACAGCGATAGC")

In [3]: seq.complement()
Out[3]: Seq('TGGCTGCTACTGTCGCTATCG', Alphabet())

Biopython bietet auch diverse Parser für versch. biol. Fileformate.

Bzw. ohne biopython, ein Auszug aus einem alten Projekt von mir:

Code: Alles auswählen

class DNA:
    """Class representing DNA as a string sequence."""
    _basecomplement = {'A':'T','C':'G','T':'A','G':'C','N':'N','M':'K','K':'M',
		       'R':'Y','Y':'R','W':'W','S':'S','H':'D','V':'B','D':'H',
		       'B':'V'}
	
    def __init__(self, seq):
        self.seq = seq.upper()

    def __repr__(self):
        """Returns DNA string"""
        letters = list(self.seq)
        return ''.join(letters)

    def complement(self):
        """Returns the complementary DNA string.""" 
        letters = [self._basecomplement[base] for base in self.seq] 
        return ''.join(letters)
	
    def reversecomplement(self):
        """Returns the reverse complement of the DNA string."""
        letters = self.seq[::-1]
        letters = [self._basecomplement[base] for base in letters] 
        return ''.join(letters)

(Der Code mag nicht mehr optimal sein - meine Anfänge, lange nicht mehr dran gearbeitet.) Der Link unter dem alten Projekt funktioniert nicht mehr, aber den Code habe ich noch, falls Interesse besteht.

Was ich nicht verstehe: Was sollen diese "Pipes"?
Gruß,
Christian

PS Einrückung korrigiert.

[CM]

Zu Deinem Ansatz: unabhängig davon, ob Du externe Module nehmen willst eine Klasse einsetzen möchtest, usw.: Du solltest doch dicts nutzen. Das erspart Dir die vielen ifs/elifs.

statt

Code: Alles auswählen

if sequence.isupper():
...
else:
...

Könntest Du auch hierbei try/except verwendet, bzw. einen solchen winzigen Fehler schlicht ignorieren und eben string.upper() anwenden. Fehler im Parsen, bzw. Öffnen der Datei haben außerdem nichts mit dem Code zu tun. Auch dies kann man trennen - wodurch ein nur etwas größeres Projekt gleich viel stabiler würde.

Gruß,
Christian

[Leonidas]

Ja, die Lösung mit einem Dictionary ist definitiv hübscher als eine if-elif-Kaskade. Hätte ich auch so gemacht, allerdings hatte ich in meiner DNA nur 4 Basen verwendet Aber das weiß CM bestimmt viel besser als ich.

Aber ich bin neugierig: was sind dann das für tolle Basen, wie N, R, Y, S, H, V, B, D?

[CM]

http://de.wikipedia.org/wiki/Nukleotide#Ambiguity_Code

[EyDu]

Code: Alles auswählen

>>> import string
>>> table = string.maketrans("ACTG", "TGAC")
>>> "ACCGACGATGACAGCGATAGC".translate(table)
'TGGCTGCTACTGTCGCTATCG'

[CM]

Interessante Option, danke.

Aber, wenn ich nochmal Werbung machen darf, die Leistungsfähigkeit von biopython können wir hier nicht im Forum zusammenstückeln .

[ne0h]

Hallo,


@CM

Danke für diesen Tip. Da ich eigentlich vor hatte, ein kleines Programm zu schreiben, welches mir die DNA Sequenzen simuliert (bis hin zu graphischen Darstellungen) hat sich da wohl nun erledigt.

Ich hatte nicht damit gerechnet, dass es bereits ein solches Modul gibt für Python.

Aber dazu mal ne Frage (wenn auch etwas blöd klingend): Wo genau finde ich denn solche Module bzw. wo finde ich eine Übersicht über alle verfügbaren Python Module? Ich meine, Googlen ist ja schön und gut, aber irgendwie weiss ich da nie, ob es sich um etwas aktuelles handelt bzw. um ein wirklich "gutes" Modul.

Gibt es nicht etwas dass dem CPAN bei Perl entspricht?


Zum Thema:

Ok, Deine Lösung ist natürlich um einiges schlanker und übersichtlicher sowie effizienter.

Aber was ich nicht ganz verstehe ist das hier:

Könntest Du auch hierbei try/except verwendet, bzw. einen solchen winzigen Fehler schlicht ignorieren und eben string.upper() anwenden. Fehler im Parsen, bzw. Öffnen der Datei haben außerdem nichts mit dem Code zu tun. Auch dies kann man trennen - wodurch ein nur etwas größeres Projekt gleich viel stabiler würde.

Was genau meinst Du mit:

...winzigen Fehler schlicht ignorieren und eben string.upper() anwenden.

Und wie kann ich den letzten Teil verstehen? Wie genau sollte Deiner Meinung nach die Trennung aussehen (wohin gehört dieser Teil?) ?

Danke für die Verbesserungsvorschläge.


Ach ja, ehrlich gesagt, bei genauerem Betrachten empfinde ich gerade diese Lösung hier von EyDu:

Code: Alles auswählen

>>> import string
>>> table = string.maketrans("ACTG", "TGAC")
>>> "ACCGACGATGACAGCGATAGC".translate(table)
'TGGCTGCTACTGTCGCTATCG'

als noch etwas eleganter, wobei ich jetzt einfach mal ganz pauschal nur von den Standardmodulen ausgehe, die Python mit sich bringt.

Ich hatte diese "Übersetzung" anhand einer Tabelle gestern auch im Sinn, dachte aber, dass es zu viel Aufwand würde (na ja, ich habe das wohl auch nicht wirklich gut durchdacht und mich nur auf simple Konstrukte ie if - elif - else versteift)....


Gruß

ne0h

[ne0h]

Was ich nicht verstehe: Was sollen diese "Pipes"?

Die stellen in der Konsole nur die "Verbindung" der einzelnen Basen dar (also das Komplement bzw. die zugehörige Base).

Gruß

ne0h

[ne0h]

Das mit den Modulen hat sich erledigt. Ich denke, ich habe gefunden was ich gesucht habe: http://pypi.python.org/pypi


ne0h

[Leonidas]

Ich denke, ich habe gefunden was ich gesucht habe: http://pypi.python.org/pypi

Ja, die Sachen im Cheeseshop sind in der Regel auch easy_install'able

[CM]

Hoi,

eine ausführliche Liste wissenschaftsnaher Software findest Du außerdem hier. Es ist ein Wiki. Solltest Du weitere Software kennen: Einfach hinzufügen.

Was genau meinst Du mit:

...winzigen Fehler schlicht ignorieren und eben string.upper() anwenden.

Was ich meine ist folgendes: Angenommen Du hast einen String "aTGcTgAAtCCG", dann kannst Du das natürlich als einen Fehler ansehen - wie Du es oben machst oder .upper() anwenden und alles ist ok. Erst, wenn ein Nicht-Nukleotid in der Sequenz ist, würde ich einen Fehler unterstellen. Je nach Quelle der Sequenzen, kommt es nämlich zu Groß- oder/und Kleinschreibung.

Außerdem: Was ich meinte mit der Fehlerbehandlung ist: Fehler können - mal auf Dein Beispiel bezogen - an verschiedenen Stellen auftreten. Mit der Datei, die Du öffnest könnte etwas nicht stimmen. Oder es könnte ein Fehler in der Sequenz auftreten. Aber Du kontrollierst nur: Gibt es überhaupt einen Fehler oder nicht? Python eröffnet Dir aber die Möglichkeit nach Fehlerquellen zu unterscheiden und entsprechend zu reagieren. (siehe http://docs.python.org/tut/node10.html ) Und zusätzlich kannst Du auch sowohl die Daten einlesen und dort auf Fehler checken als auch die Sequenz bearbeiten - und dort wieder auf Fehler checken und entsprechend reagieren.
Keine Sorge: Das Verständnis dafür kommt mit der Zeit.

Gruß,
Christian

[ne0h]

Ja, die Sachen im Cheeseshop sind in der Regel auch easy_install'able

Ach, das ist also der berühmte Cheeseshop?


@CM

Ok, ich denke ich verstehe Deine Punkte.

Sehe ich auch ein, dass es besser ist, einen Fehler zu Klassifizieren und demnach auch zu behandeln.

Ich werde heute Abend mal das Skript etwas erweitern und ausführlicher ausbauen.

Und vielen Dank für die Links!


ne0h