2016-09-15 2 views
0

Ich versuche, die Dinuc Anzahl und Frequenzen aus einer Sequenz in einer Textdatei zu finden, aber mein Code gibt nur einzelne Nukleotid zählt.Dinukleotid Anzahl und Häufigkeit

e = "ecoli.txt" 

ecnt = {} 

with open(e) as seq: 
    for line in seq: 
     for word in line.split(): 
      for i in range(len(seqr)): 
       dinuc = (seqr[i] + seqr[i:i+2]) 
       for dinuc in seqr: 
        if dinuc in ecnt: 
         ecnt[dinuc] += 1 
        else: 
         ecnt[dinuc] = 1 

for x,y in ecnt.items(): 
    print(x, y) 

Probeneingabe: "AAATTTCGTCGTTGCCC"

Beispielausgabe: AA: 2 TT: 3 TC: 2 CG: 2 GT: 2 GC: 1 CC: 2

Im Moment bekomme ich nur einzelne Nukleotide für meine Ausgabe:

C 83550600 A 60342100 T 88.192.300 G 92834000

Für die Nukleotide, die also „AAA“ zu wiederholen, hat die Zählung alle möglichen Kombinationen von aufeinander folgenden ‚AA‘ zurückzukehren, so sollte der Ausgang 2 eher als 1. Es spielt keine Rolle, was die Bestellung Dinukleotide sind aufgeführt, ich brauche nur alle Kombinationen und für den Code die korrekte Anzahl für die wiederholten Nukleotide zurückzugeben. Ich fragte meine TA und sie sagte, dass mein einziges Problem darin bestand, meine "for" -Schleife dazu zu bringen, die Dinukleotide zu meinem Wörterbuch hinzuzufügen, und ich denke, dass meine Reichweite falsch ist oder nicht. Die Datei ist sehr groß, daher ist die Sequenz in Zeilen aufgeteilt.

Vielen Dank im Voraus !!!

+1

zeigen einen kurzen Abschnitt der Probeneingabe und den entsprechenden gewünschten Ausgang. – John1024

+0

was ist 'seqr'? Es ist nicht in dem Snippet definiert, das Sie gepostet haben –

+0

Ihr Code ist auf viele Arten kaputt. Was ist 'Seqr'? Warum teilst du hier die Zeile durch Leerzeichen? 'Für word in line.split():', soll es nicht eine DNA-Sequenz oder etwas sein? Sie entfernen das Newline-Symbol nicht. –

Antwort

0

Ich habe mir Ihren Code angeschaut und einige Dinge gefunden, die Sie sich ansehen könnten.

Für meine Lösung zu testen, da ich nicht ecoli.txt hatte, erzeugte ich einen meiner eigenen mit zufälligen Nukleotiden mit der folgenden Funktion:

import random 
def write_random_sequence(): 
    out_file = open("ecoli.txt", "w") 
    num_nts = 500 
    nts_per_line = 80 
    nts = [] 
    for i in range(num_nts): 
     nt = random.choice(["A", "T", "C", "G"]) 
     nts.append(nt) 
    lines = [nts[i:i+nts_per_line] for i in range(0, len(nts), nts_per_line)] 
    for line in lines: 
     out_file.write("".join(line) + "\n") 
    out_file.close() 
write_random_sequence() 

Beachten Sie, dass diese Datei eine einzelne Sequenz von 500 Nukleotiden getrennt in Linien von jeweils 80 Nukleotiden. Um Dinukleotide zu zählen, bei denen Sie das erste Nukleotid am Ende einer Zeile und das zweite Nukleotid am Anfang der nächsten Zeile haben, müssen wir alle diese separaten Zeilen in eine einzige Zeichenfolge ohne Leerzeichen zusammenführen. Lassen Sie uns das zuerst tun:

seq = "" 
with open("ecoli.txt", "r") as seq_data: 
    for line in seq_data: 
     seq += line.strip() 

Versuchen Sie den Druck aus „seq“ und feststellen, dass es enthält alle Nukleotide eine Riesen Zeichenfolge sein sollte. Als nächstes müssen wir die Dinukleotide in der Sequenzkette finden. Wir können das mit Slicing machen, was ich gesehen habe. Für jede Position in der Kette betrachten wir sowohl das aktuelle als auch das nachfolgende Nukleotid.

for i in range(len(seq)-1):#note the -1 
    dinuc = seq[i:i+2] 

Wir können dann die Zählung der Nukleotide tun und Speicherung von ihnen in einem Wörterbuch „ECNT“ sehr ähnlich wie Sie hatten. Der endgültige Code sieht wie folgt aus:

ecnt = {} 
seq = "" 
with open("ecoli.txt", "r") as seq_data: 
    for line in seq_data: 
     seq += line.strip() 
for i in range(len(seq)-1): 
    dinuc = seq[i:i+2] 
    if dinuc in ecnt: 
     ecnt[dinuc] += 1 
    else: 
     ecnt[dinuc] = 1 
print ecnt 
0

eine perfekte Gelegenheit, ein defaultdict zu verwenden:

from collections import defaultdict 

file_name = "ecoli.txt" 

dinucleotide_counts = defaultdict(int) 

sequence = "" 

with open(file_name) as file: 
    for line in file: 
     sequence += line.strip() 

for i in range(len(sequence) - 1): 
    dinucleotide_counts[sequence[i:i + 2]] += 1 

for key, value in sorted(dinucleotide_counts.items()): 
    print(key, value) 
Verwandte Themen