2012-10-28 22 views
12

Używam Pythona i wyrażenie regularne, aby znaleźć ORF (otwarta ramka odczytu).Jak znaleźć otwartą ramkę odczytu w Pythonie

Znajdź podnapis ciąg, który składa się wyłącznie z liter ATGC (bez spacji i nowych linii), które:

rozpoczyna się ATG, kończy TAG lub TAA lub TGA i należy rozważyć kolejności od pierwszy znak, potem drugi, a następnie trzeci:

Seq= "CCTCAGCGAGGACAGCAAGGGACTAGCCAGGAGGGAGAACAGAAACTCCAGAACATCTTGGAAATAGCTCCCAGAAAAGC 
AAGCAGCCAACCAGGCAGGTTCTGTCCCTTTCACTCACTGGCCCAAGGCGCCACATCTCCCTCCAGAAAAGACACCATGA 
GCACAGAAAGCATGATCCGCGACGTGGAACTGGCAGAAGAGGCACTCCCCCAAAAGATGGGGGGCTTCCAGAACTCCAGG 
CGGTGCCTATGTCTCAGCCTCTTCTCATTCCTGCTTGTGGCAGGGGCCACCACGCTCTTCTGTCTACTGAACTTCGGGGT 
GATCGGTCCCCAAAGGGATGAGAAGTTCCCAAATGGCCTCCCTCTCATCAGTTCTATGGCCCAGACCCTCACACTCAGAT 
CATCTTCTCAAAATTCGAGTGACAAGCCTGTAGCCCACGTCGTAGCAAACCACCAAGTGGAGGAGCAGCTGGAGTGGCTG 
AGCCAGCGCGCCAACGCCCTCCTGGCCAACGGCATGGATCTCAAAGACAACCAACTAGTGGTGCCAGCCGATGGGTTGTA 
CCTTGTCTACTCCCAGGTTCTCTTCAAGGGACAAGGCTGCCCCGACTACGTGCTCCTCACCCACACCGTCAGCCGATTTG 
CTATCTCATACCAGGAGAAAGTCAACCTCCTCTCTGCCGTCAAGAGCCCCTGCCCCAAGGACACCCCTGAGGGGGCTGAG 
CTCAAACCCTGGTATGAGCCCATATACCTGGGAGGAGTCTTCCAGCTGGAGAAGGGGGACCAACTCAGCGCTGAGGTCAA 
TCTGCCCAAGTACTTAGACTTTGCGGAGTCCGGGCAGGTCTACTTTGGAGTCATTGCTCTGTGAAGGGAATGGGTGTTCA 
TCCATTCTCTACCCAGCCCCCACTCTGACCCCTTTACTCTGACCCCTTTATTGTCTACTCCTCAGAGCCCCCAGTCTGTA 
TCCTTCTAACTTAGAAAGGGGATTATGGCTCAGGGTCCAACTCTGTGCTCAGAGCTTTCAACAACTACTCAGAAACACAA 
GATGCTGGGACAGTGACCTGGACTGTGGGCCTCTCATGCACCACCATCAAGGACTCAAATGGGCTTTCCGAATTCACTGG 
AGCCTCGAATGTCCATTCCTGAGTTCTGCAAAGGGAGAGTGGTCAGGTTGCCTCTGTCTCAGAATGAGGCTGGATAAGAT 
CTCAGGCCTTCCTACCTTCAGACCTTTCCAGATTCTTCCCTGAGGTGCAATGCACAGCCTTCCTCACAGAGCCAGCCCCC 
CTCTATTTATATTTGCACTTATTATTTATTATTTATTTATTATTTATTTATTTGCTTATGAATGTATTTATTTGGAAGGC 
CGGGGTGTCCTGGAGGACCCAGTGTGGGAAGCTGTCTTCAGACAGACATGTTTTCTGTGAAAACGGAGCTGAGCTGTCCC 
CACCTGGCCTCTCTACCTTGTTGCCTCCTCTTTTGCTTATGTTTAAAACAAAATATTTATCTAACCCAATTGTCTTAATA 
ACGCTGATTTGGTGACCAGGCTGTCGCTACATCACTGAACCTCTGCTCCCCACGGGAGCCGTGACTGTAATCGCCCTACG 
GGTCATTGAGAGAAATAA" 

Co próbowałem:

# finding the stop codon here 

def stop_codon(seq_0): 

     for i in range(0,len(seq_0),3): 
      if (seq_0[i:i+3]== "TAA" and i%3==0) or (seq_0[i:i+3]== "TAG" and i%3==0) or (seq_0[i:i+3]== "TGA" and i%3==0) : 
       a =i+3 

       break 

      else: 
       a = None 

# finding the start codon here 

startcodon_find =[m.start() for m in re.finditer('ATG', seq_0)] 

Jak znaleźć sposób sprawdzenia kodonu start, a następnie znaleźć kodon pierwszego stopu. Następnie znajdź następny kodon start i kodon stop.

Chciałbym uruchomić to dla trzech klatek. Jak wspomniano wcześniej, trzy ramki będą traktować pierwsze, drugie i trzecie znaki sekwencji jako początek.

także sekwencja musi zostać podzielona na mniejsze części 3. Nie powinno być dla niektórych rzeczą tak:

ATG TTT AAA ACA AAA TAT TTA TCT AAC CCA ATT GTC TTA ATA ACG CTG ATT TGA 

Każda pomoc będzie mile widziane.

Moja ostatnia odpowiedź:

def orf_find(st0): 

    seq_0="" 
    for i in range(0,len(st0),3): 
     if len(st0[i:i+3])==3: 
      seq_0 = seq_0 + st0[i:i+3]+ " " 

    ms_1 =[m.start() for m in re.finditer('ATG', seq_0)] 
    ms_2 =[m.start() for m in re.finditer('(TAA)|(TAG)|(TGA)', seq_0)] 

    def get_next(arr,value): 
     for a in arr: 
      if a > value: 
       return a 
     return -1 




    codons = [] 
    start_codon=ms_1[0] 
    while (True): 
     stop_codon = get_next(ms_2,start_codon) 
     if stop_codon == -1: 
      break 
     codons.append((start_codon,stop_codon)) 
     start_codon = get_next(ms_1,stop_codon) 
     if start_codon==-1: 
      break 

    max_val = 0 
    selected_tupple =() 
    for i in codons: 
     k=i[1]-i[0] 
     if k > max_val: 
      max_val = k 
      selected_tupple = i 

    print "selected tupple is ", selected_tupple 

    final_seq=seq_0[selected_tupple[0]:selected_tupple[1]+3] 

    print final_seq 
    print "The longest orf length is " + str(max_val) 



output_file = open('Longorf.txt','w') 
output_file.write(str(orf_find(st0))) 

output_file.close() 

Powyższa funkcja zapisu nie może mi pomóc na piśmie treść do pliku tekstowego. Wszystko, co dostaję, jest BRAK .. Dlaczego ten błąd .. Czy ktoś może pomóc?

+0

Trzeba powrócić wyjście z funkcji. Właśnie zredagowałem swoją odpowiedź. – MoRe

Odpowiedz

8

Jak oznaczyłeś to Biopython Przypuszczam, że znasz Biopython. Czy sprawdziłeś już dokument? http://biopython.org/DIST/docs/tutorial/Tutorial.html#htoc231 może pomóc.

I dostosowane kod z powyższego linku nieco do pracy w kolejności:

from Bio.Seq import Seq 

seq = Seq("CCTCAGCGAGGACAGCAAGGGACTAGCCAGGAGGGAGAACAGAAACTCCAGAACATCTTGGAAATAGCTCCCAGAAAAGCAAGCAGCCAACCAGGCAGGTTCTGTCCCTTTCACTCACTGGCCCAAGGCGCCACATCTCCCTCCAGAAAAGACACCATGAGCACAGAAAGCATGATCCGCGACGTGGAACTGGCAGAAGAGGCACTCCCCCAAAAGATGGGGGGCTTCCAGAACTCCAGGCGGTGCCTATGTCTCAGCCTCTTCTCATTCCTGCTTGTGGCAGGGGCCACCACGCTCTTCTGTCTACTGAACTTCGGGGTGATCGGTCCCCAAAGGGATGAGAAGTTCCCAAATGGCCTCCCTCTCATCAGTTCTATGGCCCAGACCCTCACACTCAGATCATCTTCTCAAAATTCGAGTGACAAGCCTGTAGCCCACGTCGTAGCAAACCACCAAGTGGAGGAGCAGCTGGAGTGGCTGAGCCAGCGCGCCAACGCCCTCCTGGCCAACGGCATGGATCTCAAAGACAACCAACTAGTGGTGCCAGCCGATGGGTTGTACCTTGTCTACTCCCAGGTTCTCTTCAAGGGACAAGGCTGCCCCGACTACGTGCTCCTCACCCACACCGTCAGCCGATTTGCTATCTCATACCAGGAGAAAGTCAACCTCCTCTCTGCCGTCAAGAGCCCCTGCCCCAAGGACACCCCTGAGGGGGCTGAGCTCAAACCCTGGTATGAGCCCATATACCTGGGAGGAGTCTTCCAGCTGGAGAAGGGGGACCAACTCAGCGCTGAGGTCAATCTGCCCAAGTACTTAGACTTTGCGGAGTCCGGGCAGGTCTACTTTGGAGTCATTGCTCTGTGAAGGGAATGGGTGTTCATCCATTCTCTACCCAGCCCCCACTCTGACCCCTTTACTCTGACCCCTTTATTGTCTACTCCTCAGAGCCCCCAGTCTGTATCCTTCTAACTTAGAAAGGGGATTATGGCTCAGGGTCCAACTCTGTGCTCAGAGCTTTCAACAACTACTCAGAAACACAAGATGCTGGGACAGTGACCTGGACTGTGGGCCTCTCATGCACCACCATCAAGGACTCAAATGGGCTTTCCGAATTCACTGGAGCCTCGAATGTCCATTCCTGAGTTCTGCAAAGGGAGAGTGGTCAGGTTGCCTCTGTCTCAGAATGAGGCTGGATAAGATCTCAGGCCTTCCTACCTTCAGACCTTTCCAGATTCTTCCCTGAGGTGCAATGCACAGCCTTCCTCACAGAGCCAGCCCCCCTCTATTTATATTTGCACTTATTATTTATTATTTATTTATTATTTATTTATTTGCTTATGAATGTATTTATTTGGAAGGCCGGGGTGTCCTGGAGGACCCAGTGTGGGAAGCTGTCTTCAGACAGACATGTTTTCTGTGAAAACGGAGCTGAGCTGTCCCCACCTGGCCTCTCTACCTTGTTGCCTCCTCTTTTGCTTATGTTTAAAACAAAATATTTATCTAACCCAATTGTCTTAATAACGCTGATTTGGTGACCAGGCTGTCGCTACATCACTGAACCTCTGCTCCCCACGGGAGCCGTGACTGTAATCGCCCTACGGGTCATTGAGAGAAATAA") 


table = 1 
min_pro_len = 100 

for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]: 
    for frame in range(3): 
     for pro in nuc[frame:].translate(table).split("*"): 
      if len(pro) >= min_pro_len: 
       print "%s...%s - length %i, strand %i, frame %i" % (pro[:30], pro[-3:], len(pro), strand, frame) 

ORF jest również tłumaczone. Możesz wybrać inną tabelę tłumaczeń. Sprawdź http://biopython.org/DIST/docs/tutorial/Tutorial.html#sec:translation

EDIT: wyjaśnienie kodu:

na samym szczycie utworzyć obiekt sekwencji z Twojego ciąg. Zwróć uwagę na seq = Seq("ACGT"). Dwie pętle for tworzą 6 różnych ramek. Wewnętrzna pętla for tłumaczy każdą ramkę zgodnie z wybraną tablicą translacji i zwraca łańcuch aminokwasów, gdzie każdy kodon stop jest kodowany jako *. Funkcja split dzieli ten ciąg, usuwając te symbole zastępcze, co daje listę możliwych sekwencji białek. Poprzez ustawienie min_pro_len można zdefiniować minimalną długość łańcucha aminokwasowego dla wykrywanego białka. 1 to standardowa tabela. Check out http://www.ncbi.nlm.nih.gov/Taxonomy/Utils/wprintgc.cgi#SG1 Tutaj widzisz, że kodon inicjacyjny to AUG (jest równy ATG), a kodony końcowe (* w sekwencji nukleotydów) to TAA, TAG i TGA, dokładnie tak jak chciałeś. Możesz także użyć innej tabeli tłumaczeń.

Po dodaniu

print nuc[frame:].translate(table) 

prawo wewnątrz sekundę dla pętli można dostać coś takiego:

PQRGQQGTSQEGEQKLQNILEIAPRKASSQPGRFCPFHSLAQGATSPSRKDTMSTESMIRDVELAEEALPQKMGGFQNSRRCLCLSLFSFLLVAGATTLFCLLNFGVIGPQRDEKFPNGLPLISSMAQTLTLRSSSQNSSDKPVAHVVANHQVEEQLEWLSQRANALLANGMDLKDNQLVVPADGLYLVYSQVLFKGQGCPDYVLLTHTVSRFAISYQEKVNLLSAVKSPCPKDTPEGAELKPWYEPIYLGGVFQLEKGDQLSAEVNLPKYLDFAESGQVYFGVIAL*REWVFIHSLPSPHSDPFTLTPLLSTPQSPQSVSF*LRKGIMAQGPTLCSELSTTTQKHKMLGQ*PGLWASHAPPSRTQMGFPNSLEPRMSIPEFCKGRVVRLPLSQNEAG*DLRPSYLQTFPDSSLRCNAQPSSQSQPPSIYICTYYLLFIYYLFICL*MYLFGRPGCPGGPSVGSCLQTDMFSVKTELSCPHLASLPCCLLFCLCLKQNIYLTQLS**R*FGDQAVATSLNLCSPREP*L*SPYGSLREI 

(zauważyć gwiazdki są w pozycjach kodonu stop)

EDIT : Odpowiedz na drugie pytanie:

Musisz zwrócić ciąg, który chcesz zapisać w pliku. Załóż łańcuch wyjściowy i odesłać go na końcu funkcji:

output = "selected tupple is " + str(selected_tupple) + "\n" 
output += final_seq + "\n" 
output += "The longest orf length is " + str(max_val) + "\n" 
return output 
+0

daje mi błąd "STR" nie ma atrybutu "reverse_complement" – Nodnin

+0

zauważyłeś pierwsze 3 wiersze kodu, który opublikowałem? Istnieje konwersja z twojego ciągu znaków do obiektu sekwencji. – MoRe

+0

Ale nie widzę, jak rozpocznie się sekwencja z "ATG" i kończy się z kodonem stop .. jak mogę znaleźć wiele podciągów w jednej ramce odczytu – Nodnin

9

Jeśli chcesz ręcznie kod go:

import re 
from string import maketrans 

pattern = re.compile(r'(?=(ATG(?:...)*?)(?=TAG|TGA|TAA))') 

def revcomp(dna_seq): 
    return dna_seq[::-1].translate(maketrans("ATGC","TACG")) 

def orfs(dna): 
    return set(pattern.findall(dna) + pattern.findall(revcomp(dna))) 

print orfs(Seq)