Tetrapeptydy

W tym przykładzie sprawdzimy, jakie cztero-aminokwasowe fragmenty najczęściej pojawiają się w białkach.

Białka to liniowe łańcuchy składające się z jednostek zwanych aminokwasami. Ponieważ w przyrodzie występuje jedynie 20 “standardowych” aminokwasów, każdy można oznaczyć inną literą alfabetu (tzw. kod jednoliterowy). Zatem sekwencja aminokwasowa każdego białka może zostać przedstawiona jako długi ciąg liter. Typowa sekwencja aminokwasowa wygląda jak poniżej:

>1dtjC
MKELVEMAVPENLVGAILGKGGKTLVEYQELTGARIQISTRNRRVTITGS
PAATQAAQYLISQRVT

Sekwencję zapisano w formacie FASTA: linijka zaczynająca się od > zawiera identyfikator sekwencji, następnie sama sekwencja może zająć więcej niż jedną linijkę. Koniec sekwencji rozpoznajemy po (i) pustej linii lub (ii) początku nowej sekwencji, czyli kolejnym nagłówku.

Jak można zauważyć, niektóre litery występują częściej niż inne a ich układ jest nieprzypadkowy. Celem tego ćwiczenia jest sprawdzenie, jakie czteroliterowe fragmenty pojawiają się najczęściej?

Słownik o czteroliterowych kluczach

Aby odpowiedzieć na to pytanie, wykorzystamy słownik. Poniższy program:

  • wczytuje plik w formacie FASTA zawierający sekwencje białek (linie 5,6); oprócz linii zawierających sekwencję w pliku fasta można spotkać linie puste oraz linie będące nagłówkiem (zaczynają się od znaku >), te dwie ostatnie możliwości wykluczamy w linii 8

  • każdą sekwencję dzielimy na cztero-literowe fragmenty (linia 10,11)

  • w słowniku T zapamiętujemy, ile razy trafił się każdy z tetrapeptydów

  • na koniec przebiegamy przez cały słownik w poszukiwaniu klucza odpowiadającego największej wartości

 1import math
 2
 3T = {}
 4cnt = 0
 5with open('../../_static/chains_from_db-uniq1000.fasta') as file:
 6    for linia in file:
 7        linia = linia.strip()
 8        if len(linia) > 0 and linia[0] != ">":
 9            # --- tu tniemy sekwencje na kawalki po 4 aminokwasy (4 litery)
10            for i in range(len(linia)-3):
11                pept = linia[i:i + 4]   # --- 4 literowy fragment
12                if not (pept in T):     # --- czy mamy taki w słowniku?
13                    T[pept] = 0         # --- nie, wstawiamy z licznikiem 0
14                T[pept] = T[pept] + 1   # --- zwiększamy zliczenia o 1
15                cnt += 1                # --- liczba wszystkich czwórek
16
17# --- Teraz przejdziemy słownik raz jeszcze i poszukamy najpopulaniejszych tetrapeptydów
18max_v = 0               # --- największa wartość w słwniku
19max_k = ""              # --- klucz odpowiadający największej wartości
20for k, v in T.items():  # --- iterujemy po parach klucz-wartość
21    if v > max_v:
22        max_v = v
23        max_k = k
24print(max_k, max_v/cnt)
25# print(T)
26
27

Uruchom powyższy program dla plików z przykładowymi danymi wejściowymi: mały, duży. Kiedy już zrozumiesz, jak on działa, wprowadź następujące zmiany:

  • W danych wejściowych każda sekwencja zapisana jest w jednej linii, co istotnie upraszcza wczytywanie danych. Zmodyfikuj program tak, aby każda sekwencja mogła być w kilku linijkach. Dopuszczalna jest tylko jedna linia nagłówka dla jednej sekwencji

  • Oblicz prawdopodobieństwo wystąpienia każdego aminokwasu, tzn ile razy pojawia się litera W czy L. Następnie oblicz oczekiwane prawdopodobieństwo każdego tetrapeptydu jako iloczyn odpowiednich prawdopodobieństw i wybierz tetrapeptydy pojawiające się niespodziewanie, czyli istotnie częściej niż wartość oczekiwana.