UNIVERZA V LJUBLJANI FAKULTETA ZA MATEMATIKO IN FIZIKO Finančna matematika – 1. stopnja Anže Mramor Algoritmi za optimizacijo rabe kodonov Delo diplomskega seminarja Mentorica: izr. prof. dr. Arjana Žitnik Ljubljana, 2021

(1)

UNIVERZA V LJUBLJANI

FAKULTETA ZA MATEMATIKO IN FIZIKO Finančna matematika – 1. stopnja

Anže Mramor

Algoritmi za optimizacijo rabe kodonov

Delo diplomskega seminarja

Mentorica: izr. prof. dr. Arjana Žitnik

Ljubljana, 2021

(2)

Kazalo

1. Uvod 5

2. Biokemijsko ozadje problema optimalnega zaporedja kodonov 5

2.1. Proteini 6

2.2. Proces translacije mRNA, oziroma sinteza proteina 6 3. Biokemijski problem – optimalno zaporedje kodonov 8

3.1. Predstavitev podatkov 8

3.2. Matematični model 9

4. Dinamično programiranje 11

4.1. Primer dinamičnega programiranja – segmentirani najmanjši kvadrati 11

4.2. Časovna zahtevnost 13

5. Reševanje problema 14

5.1. Rekurzivna formula 14

5.2. Računalniška implementacija algoritma 17

6. Testiranje algoritma in rezultati 21

6.1. Osnovno preverjanje uspešnosti in začetno testiranje algoritma 21

6.2. Testiranje algoritma na realnih podatkih 22

6.3. Testiranje večje količine zaporedij in modeliranje rezultatov 23

6.4. Vpeljava šuma v podatke 27

7. Zaključek 30

Literatura 31

(3)

Algoritmi za optimizacijo rabe kodonov Povzetek

Delo diplomskega seminarja se ukvarja z matematično rešitvijo biokemijskega problema – želimo sestaviti algoritem, ki bo v realnem času poiskal optimalno zaporedje kodonov proteina za izražanje v celicahE. coli. Natančneje, za eksperimentalno dobljene podatke sinteze proteina želimo določiti takšno zaporedje kodonov, da bi vrednosti časov translacije posamezne aminokisline, ki jih izračunamo po formuli lokalnega povprečja treh zaporednih kodonov, od njih odstopale čim manj.

Problem rešujemo z dinamičnim programiranjem. Najprej izpeljemo rekurzivno formulo, ki nam bo rešila optimizacijski problem. Nato rekurzivno formulo razši- rimo v splošen računalniški algoritem, ki ga tudi implementiramo v programskem jeziku Python. Najprej preverimo delovanje algoritma na podatkih treh obstoje- čih proteinov. Zaporedje želenih časov translacije dobimo iz dejanskih proteinov – na tak način zagotovimo, da se začetni aproksimirani podatki ujemajo z rezultati algoritma, ter tako preverjamo uspešnost algoritma. Algoritem se na takšni vrsti podatkov izkaže za delujoč v vseh treh primerih, zato bazo testiranih podatkov raz- širimo na 900 naključno generiranih proteinov. Tudi tu algoritem za vsak protein s 100 % uspešnostjo najde ustrezno zaporedje kodonov. Poleg tega preverjamo še časovno zahtevnost algoritma, ki se izkaže za linearno.

V zadnjem delu naloge spremenimo podatke s šumom (s šumi velikosti 1 %, 5 % in 10 %), tako da so želeni časi drugačni od podanih translacijskih časov za translacijo kodonov. Ponovno preverimo uspešnost algoritma, ki je z večanjem šuma vedno manjša. Izkaže se, da ima ključno vlogo pri uspešnosti algoritma izbira norme znotraj rekurzivne formule. Preizkušamo prvo, drugo in neskončno normo, pri čemer se druga izkaže za najbolj uspešno, prva in neskončna pa sta v deležu uspešno določenih kodonov identični.

Math. Subj. Class. (2010): 92C40, 68W05

Ključne besede: kodoni, proteini, translacijski čas, dinamično programiranje, norme, časovna zahtevnost

(4)

Algorithms for optimization of codon usage Abstract

This thesis is dealing with mathematical solution of a biochemical problem – we wish to develop an algorithm, which will find optimal sequence of codons in a protein for expression in E. coli in real time. To be more exact, we wish to determine such sequence of codons that the difference between time values of experimentally obtained data of times of synthesis of a protein and time values of translation of each individual aminoacid, which are computed using formula of local average of three consecutive codons, would be minimal.

We are solving the problem with dynamic programming. We first derive recursive formula that solves the optimization problem. We then expand it into a computer algorithm, which we implement in Python programming language. We first check the functionality of algorithm on a small sample of three existing proteins. We get the sequence of target times of translation from the chosen proteins themselves – this way we can be sure, whether the inputed data that is being aproximated matches results of the algorithm and therefore how successful it is. With this kind of sample the accuracy rate of algorithm is 100 %. We then expand our base of test data on 900 randomly generated proteins. Algorithm’s accuracy rate remains 100 %. We are also analyzing the running time of the algorithm, which is linear.

In the last part of our thesis we apply noise ratio (with rates of 1 %, 5 % and 10 %) to our data. This way the target times become different from the values algorithm can compute out of given translation times for translation of codons. We again test the accuracy of algorithm, which is getting lower as we are raising the values of noise ratio. We conclude that the key part in accuracy rate of algorithm is played by the choice of norm inside the recursive formula. We test first, second and infinity norm. We find that algorithm has the highest accuracy rate when second norm is chosen, while first and infinity norm have identical accuracy.

Math. Subj. Class. (2010): 92C40, 68W05

Keywords: codons, proteins, translation time, dynamic programming, norms, running time

(5)

1. Uvod

Delo diplomskega seminarja združuje izbrane teme iz področja biokemije ter ma- tematike. Cilj dela diplomskega seminarja je postaviti matematično podlago in izpeljati algoritme, ki bodo v realnem času rešili optimizacijski problem - vrniti mo- rajo optimalno zaporedje kodonov proteina, za njegovo izražanje v celicah E. coli.

Izbira optimalnih kodonov ne vpliva na kemijsko sestavo proteina, vpliva pa na iz- koristek sinteze, zaradi česar je optimizacija rabe kodonov eden ključnih korakov pri sintezi rekombinantnih proteinov (v praktični rabi tip proteinov, ki jih lahko vza- memo iz enega organizma in vstavimo v drugega). Iz dobljenega zaporedja kodonov lahko sestavimo verigo DNA, ki bo nosila zapis posameznega proteina. Če poznamo optimalno zaporedje kodonov v določeni verigi, lahko te podatke uporabimo za pro- izvodnjo, na primer, človeških proteinov v bakterijah. V praksi namreč že obstajajo postopki, s katerimi je možno iz znanega zaporedja kodonov sestaviti DNA. Tega lahko nato vgradimo oziroma vstavimo v bakterije (kot na primer E. coli), ki so z danim dednim zapisom sposobne proizvajati želen protein. Produkcija poljubnega proteina je v biokemiji izjemno uporabna.

Natančneje se bomo osredotočili na biokemijski proces sinteze proteinov oziroma načrtovanje najboljšega nukleotidnega zaporedja za sintezo želenega proteina. Op- timalno rešitev bomo poiskali z modeliranjem prek matematičnih metod oziroma algoritmov. Ideja za diplomsko delo je nastala kot posledica diplomskega dela [7], v kateri avtor razvija algoritem za načrtovanje optimiziranega zaporedja kodonov proteina iz aminokislinskega zaporedja. V omenjenem delu se je avtor osredotočil predvsem na napovedovanje translacijskih časov proteinov, medtem ko se je reše- vanja problema konstruiranja najboljšega zaporedja kodonov pri danih časih lotil s hevrističnimi metodami. V našem delu diplomskega seminarja bomo poskušali slednji problem rešiti z matematičnimi modeli oziroma algoritmi, saj se je v praksi pokazala potreba po hitrejših in bolj natančnih algoritmih.

Diplomsko delo vsebuje več razdelkov. Najprej si bomo pogledali biokemijsko ozadje problema, prek katerega bomo lažje spoznali, kakšen problem pravzaprav rešujemo. Sledila bo postavitev matematičnega modela biokemijskega problema, kjer bomo problem natančno formulirali ter s tem dobili tudi idejo o njegovi rešitvi.

Problem bomo reševali z dinamičnim programiranjem, pristopom v optimizaciji, ki bo natančno predstavljen v četrtem poglavju. Na koncu si bomo ogledali reševanje problema prek rekurzivne formule in oris uporabe algoritma ter algoritem tudi preizkusili na različnih podatkih. Da lahko matematični model za naš problem pravilno zastavimo, si je najprej potrebno podrobneje ogledati dogajanje v celicah, ki pripelje do njega.

2. Biokemijsko ozadje problema optimalnega zaporedja kodonov Življenje na zemlji predstavljajo organizmi. Njihova osnovna gradbena in funk- cionalna enota je celica. Za lažjo predstavo biokemijskega ozadja problema optimalnega zaporedja kodonov si bomo ogledali molekularno biologijo celice oziroma osnovne biokemijske lastnosti organizmov preko enoceličnega modelnega organizma – E. coli. S primernimi modifikacijami in širšim znanjem biologije in kemije, lahko ugotovitve posplošimo tudi na kompleksnejše organizme. Natančneje, ker nas ne zanima celotno biokemijsko ozadje celice, si bomo podrobneje pogledali predvsem najštevilčnejšo skupino makromolekul – proteine. Ti v celici opravljajo glavno kata- lično, strukturno in regulacijsko funkcijo. Osredotočili se bomo predvsem na sintezo

(6)

proteinov v procesu imenovanemtranslacija. Translacija se nanaša na prevajanje zapisa kodonov v aminokisline, ki so osnovni gradniki proteina. Pri procesu sinteze se proteini zvijajo, način zvitja pa vpliva na njihove biokemijske lastnosti. Teoretično biokemijsko ozadje v tem razdelku je povzeto po [1] in [8].

2.1. Proteini. Proteini so makromolekule, sestavljene iz zaporedno vezanih aminokislin. Aminokisline so med seboj povezane s peptidno vezjo, zato takemu zaporedju vezanih aminokislin rečemo tudi polipeptidna veriga. Veriga lahko vsebuje med nekaj deset in več tisoč aminokislin. V našem delu se bomo omejili na proteine, v katerih nastopa med 50 in 800 aminokislin. Proteini krajši od 50 aminokislin imajo namreč preprosta zvitja, pri čemer hitrost translacije ne vpliva na zvijanje. Pri opa- zovani E. Coli, so proteini, ki bi imeli več kot 800 aminokislin redki, zato je privzetek takšne zgornje meje logičen. Podobno se je omejil tudi avtor v [7], zato se zdita ti robni vrednosti smiselni.

Obstaja 20 različnih aminokislin. Polipeptidna veriga je sestavljena iz zaporedja aminokislin, vzetih iz nabora teh dvajsetih različnih. Zaporedje teh aminokislin predstavlja informacijo o sami strukturi in posledično funkciji proteina. Posame- zne aminokisline so različnih velikosti ter imajo različne fizikalno-kemijske lastnosti (polarnost, hidrofobnost oziroma hidrofilnost, zmožnost tvorbe vodikovih vezi in podobno). S svojimi lastnostmi vplivajo na način zvijanja verige v prostoru oziroma celici. Nativno zvitjeproteina, kakor rečemo zvitju, v katerem je protein prisoten in funkcionalen, določa njegovo funkcijo. To pomeni, da lahko majhne spremembe v aminokislinskem zaporedju vplivajo na funkcijo proteina, s čimer dobimo popolnoma drugačen protein.

Zaporedju aminokislin rečemoprimarno zaporedje. To zaporedje določa oziroma zapisuje za vsa nadaljna stanja. V procesu zvijanja se preko vodikovih vezi med glavno verigo proteina tvorijo sekundarne strukture. Hidrofobne aminokisline se umaknejo od vode in tvorijo hidrofobno jedro. Preostali deli polipeptidne verige se nato umestijo v strukturo s tvorbo termodinamsko ugodnih interakcij. Na koncu se preostali deli proteina prilagodijo v prostoru in tako dobimo terciarno strukturo proteina. Izjemno pomembno je, da se protein zvije pravilno, sicer ne more opravljati svoje funkcije. V skrajnih primerih lahko zaradi napačnega zvitja proteinov, na primer v možganskih celicah, nastanejo toksične strukture, ki povzročijo odmiranje celic. Primer tega sta Alzheimerjeva in Parkinsonova bolezen, kjer zaradi nalaganja nepravilno zvitih proteinov pride do propadanja možganskih celic.

2.2. Proces translacije mRNA, oziroma sinteza proteina. V nadaljevanju si bomo pogledali, kako proteini pravzaprav nastanejo, se sintetizirajo. Razumevanje postopka sinteze proteina je izjemno pomembno za razumevanje biokemijskega problema, saj je ključno vprašanje diplomskega dela, kako izvesti translacijo, da bo imel želeni protein ravno dovolj časa, da se pravilno zvije. Oglejmo si, kako translacija mRNA, pri kateri gre za prepis dednega materiala, sploh poteka.

Dedni material celic je shranjen v deoksiribonukleinski kislini (v nadaljevanju DNA). Dedni material nosi zapis za komponente celice, regulacijo izražanja teh komponent in regulacijo delovanja celice kot celote. DNA je dvojna vijačnica, sestavljena iz zaporedno vezanih nukleinskih baz – adenin (A), timin (T), gvanin (G) in citozin (C). Baze prve verige DNA tvorijo komplementarne pare z drugo verigo DNA, ki poteka v nasprotni smeri glede na prvo. V bazni par sta lahko povezana adenin in timin (A-T) ali gvanin in citozin (G-C). Možni so tudi drugačni bazni

(7)

pari, vendar je njihovo parjenje pogojeno z drugačno strukturo DNA, ki ni predmet diplomskega dela.

Vsaka aminokislina je zapisana kot trojček (triplet) nukleotidnih baz (kot recimo ATG). Triplete nukleinskih baz imenujemo kodoni. Iz štirih različnih baz lahko se- stavimo 64 različnih zaporedij dolžine 3, oziroma 64 kodonov. Vsi razen treh izmed tripletov zapisujejo za točno določeno aminokislino. Ostali trije kodoni signalizirajo konec sinteze proteina (pravimo jim tudi stop kodoni). Število kodonov presega število različnih aminokislin, torej lahko za posamezno aminokislino zapisuje več kodonov. Te se po večini razlikujejo le v zadnji nukleinski bazi. Slednje lahko razberemo tudi iz spodnje tabele, v kateri je prikazanih vseh 64 kodonov s pripa- dajočimi aminokislinami za katere zapisujejo. Vsaka aminokislina je predstavljena z angleško tročrkovno kratico, katere najpogosteje uporabljamo tudi v praksi.

U C A G

U

UUU Phe UCU

Ser

UAU Tyr UGU

Cys U

UUC UCC UAC UGC C

UUA Leu UCA UAA Stop UGA Stop A

UUG UCG UAG Stop UGG Trp G

C

CUU

Leu

CCU

Pro

CAU His CGU

Arg

U

CUC CCC CAC CGC C

CUA CCA CAA

Gln CGA A

CUG CCG CAG CGG G

A

AUU Ile

ACU

Thr

AAU Asn AGU

Ser U

AUC ACC AAC AGC C

AUA ACA AAA

Lys AGA

Arg A

AUG Met ACG AAG AGG G

G

GUU

Val

GCU

Ala

GAU Asp GGU

Gly

U

GUC GCC GAC GGC C

GUA GCA GAA

Glu GGA A

GUG GCG GAG GGG G

Tabela 1. Tabela kodonov z aminokislinami, za katere zapisujejo.

Deli DNA, imenovani geni, se s pomočjo encima RNA polimeraze prepišejo v molekule mRNA (messenger RNA). RNA je enovijačna ribonukleinska kislina, kjer so vsi timini zamenjani z uracili (U). Takšni zapisi se lahko prepisujejo v proteine.

Molekule mRNA prepoznajo ribosomi, ki so odgovorni za sintezo proteinov. Ri- bosomi so markomolekule, ki se vežejo na mRNA, in začnejo sintezo proteina pri aminokislini imenovani metionin. Kodon metionina je zapisan s tripletom AUG in je najpogostejši začetni kodon. Ribosom se nato premika po en kodon naprej in vsakega prevede v ustrezno aminokislino. Temu procesu rečemo sinteza proteina oziroma translacija mRNA. Novo nastajajoča polipeptidna veriga se podaljša v notranjosti ribosoma in prihaja na njegovo površje preko tunela. Nemudoma po pri- hodu na površje se začnejo posamezni deli proteina zvijati. To je ključno za končno zvitje proteina, saj je zvitje njegovih posameznih delov lahko oteženo ob prisotnosti ostale polipeptidne verige.

Želene čase translacije posameznih kodonov lahko podamo kot pozitivna realna števila. Ti časi so eksperimentalno izračunani. Pri tem časi niso izraženi v stan- dardnih časovnih enotah (sekunde, milisekunde ipd.), temveč gre za relativne enote,

(8)

ki so med seboj v razmerjih. V splošnem so časi odvisni od stopnje elongacije, ki nam pove, koliko kodonov ene vrste naj bi se prevedlo v eni sekundi. Dejanski časi za naše delo so vzeti iz članka [6], v katerem sta avtorja s pomočjo Markovskih verig določila čase translacije kodonov za E. coli. V splošnem je hitrost translacije med 10 in 20 aminokislin na sekundo odvisno od različnih stanj v celici. Zaradi različnih časov translacije pride do različnih hitrosti potovanja ribosoma po mRNA, kar pomeni, da so hitrosti izhajanja polipeptidnih verig lahko za iste aminokisline različni. Hitrost translacije proteina (čas translacije, ki nas zanima) opišemo kot vsoto lokalnih povprečji časov translacije kodonov oziroma kot vsoto povprečnih vrednosti časov translacije treh zaporednih kodonov. Tako dobimo rezultate, ki so najbolj skladni z eksperimentalno dobljenimi podatki.

3. Biokemijski problem – optimalno zaporedje kodonov

Po natančni obravnavi ozadja dogajanja med sintezo lahko predstavimo problem, ki ga želimo rešiti. Za eksperimentalno dobljene podatke sinteze proteina (oziroma translacije mRNA) želimo določiti takšno zaporedje kodonov, da bi po formuli (lokalno povprečje treh zaporednih kodonov) izračunane vrednosti časov translacije od njih odstopale čim manj. Predpostavljamo, da bo hitrost translacije pridobljenega zaporedja kodonov optimalna za pravilno zvitje proteina. Napačno zvitje proteina, torej če na nekaterih mestih v zaporedju izberemo napačne kodone, lahko vodi v nefunkcionalnost proteina.

V nadaljevanju si bomo najprej pogledali kakšne podatke pravzaprav modeliramo, ter formalno opredelili matematični model, ki bo problem reševal. Izkazalo se bo, da je eden ključnih dejavnikov, ki jih moramo upoštevati pri oblikovanju algoritma, izbira norme.

3.1. Predstavitev podatkov. Vsako izmed 20 aminokislin si lahko predstavimo kot množico diskretnih vrednosti (časi translacije kodonov):

A₁, A₂, . . . , A₂₀.

Vsaka izmed množic vsebuje vsaj enega, vendar ne več kot 6 elementov (časov):

A_i ={aⁱ₁, . . . , aⁱ_n_i}; n_i ∈N, n_i ≤6, i= 1, . . . ,20.

Ker točno vemo, na koliko različnih načinov (iz koliko različnih kodonov z različnimi translacijskimi časi) se lahko zapiše posamezna aminokislina, lahko podamo tabelo vrednosti, v kateri aminokisline razdelimo po skupinah, glede na število možnih časov.

Število aminokislin Število časov translacije

2 1

9 2

1 3

5 4

3 6

Tabela 2. Tabela števil aminokislin s posameznim številom časov translacij.

(9)

Podano imamo aminokislinsko zaporedje, ki določa čas translacije mRNA za že- leni protein. V zaporedju je med 50 in 800 aminokislin, običajno okoli 300 (aminokisline se lahko ponavljajo). Translacijski čas, ki bo nastopal za posamezno aminokislino na določenem mestu, lahko izbiramo iz množice njenih možnih translacijskih časov (izberemo lahko poljuben kodon, ki bo zapisoval za aminokislino na njenem mestu v zaporedju). Na primer, za aminokislino za katero zapisuje natanko en kodon, lahko izberemo zgolj en translacijski čas. Po drugi strani imamo za aminokislino, za katero zapisuje šest kodonov, na voljo kateregakoli izmed šestih translacijskih časov. Na dveh različnih mestih v zaporedju lahko za isto aminokislino izberemo poljuben (drug) kodon, ki bo zanjo zapisoval ter s tem drug translacijski čas iz množice njenih časov. Z določanjem translacijskih časov vsem aminokislinam v zaporedju dobimo translacijski čas proteina. Izbran čas za aminokislino na i-tem mestu označimo z a_i.

Podan je tudi dejanski čas translacije mRNA za želeni protein, kot zaporedje pozitivnih realnih vrednosti:

X₁, X₂, . . . , X_m; m∈N, 50≤m≤800.

Število m ustreza številu aminokislin v podanem zaporedju. Vrednosti X_i so pri- dobljene eksperimentalno. Sami lahko sedaj sestavimo zaporedje, prav tako dolžine m, časov translacij posameznih aminokislin:

Y₁, Y₂, . . . , Y_m; m∈N, 50≤m ≤800.

VrednostiYiso izračunane kot lokalno povprečje translacijskih časov treh zaporednih kodonov, torej po formuli:

Y_i = ai−1+a_i +a_i+1

3 ; i= 2, . . . , m−1.

Problem, ki ga rešujemo pa je naslednji: Za podano znano zaporedje aminokislin (dolžine m) in znan čas translacije mRNA za želeni protein (zaporedje Xi), bi radi določili zaporedje kodonov (oziroma izbrali tako zaporedje časov translacije aminokislin), da bo odstopanje vrednosti Y_i (izračunanih po formuli in iz izbranih podatkov) od vrednosti X_i (podanih vrednosti) čim manjše za vsak i.

3.2. Matematični model. Sedaj lahko iz znanih podatkov in splošno opisanega teoretičnega problema sestavimo matematični model.

Podatki:

• m – dolžina zaporedja,

• j₁, j₂. . . , j_m; j_i ∈ {1,2, . . . ,20} – zaporedje m aminokislin, iz nabora vseh dvajsetih možnih,

• X₁, X₂, . . . , X_m; X_i ∈R – zaporedje eksperimentalnih podatkov.

Iščemo a₁, a₂, . . . , am; ai ∈Aji, da bo ∥X−Y∥ najmanjša, kjer je Yi = a_i−1+a_i +a_i+1

3 ; i= 2, . . . , m−1.

Za robne vrednosti je treba problem premisliti za vsak primer posebej.

Za norme bi lahko vzeli:

• prvo normo: ∥X−Y∥₁ =^∑︁^m−1_i=2 |X_i−Y_i|; i= 2, . . . , m−1,

• drugo normo: ∥X−Y∥2 =^√︂^∑︁^m−1_i=2 (X_i−Yi)²; i= 2, . . . , m−1,

(10)

• neskončno normo: ∥X−Y∥∞ = max{|X_i−Y_i|; i= 2, . . . , m−1}.

Najpomembnejši povezavi med normami sta podani s spodnjim izrekom [5].

Izrek 3.1. Za vsak vektor v∈R^m veljata naslednji neenakosti:

||v||∞ ≤ ||v||₂ ≤ ||v||₁

||v||₁ ≤√

m||v||₂ ≤m||v||∞.

Seveda se postavlja vprašanje, katero izmed naštetih norm je najbolj smiselno uporabiti. Poglejmo si nekaj prednosti in slabosti vsake izmed njih [5]. V prak- tičnih aplikacijah teorije aproksimacije se najpogosteje uporablja druga norma, saj ima veliko ugodnih lastnosti. Druga norma je odvedljiva in nam, v primeru linearne aproksimacijske funkcije, problem, ki ga rešujemo, prevede na problem najmanjših kvadratov. To se zdi zelo ugodno tudi za naš problem, saj bi drugo normo prav zaradi te lastnosti lahko enostavno vključili, ter nato iskali najmanjše skupno odstopanje po metodi najmanjših kvadratov. Poleg tega vemo, da je druga norma strogo konveksna, kar nam zagotavlja, da najboljša aproksimacija obstaja in je hkrati eno- lično določena. Odvisna je samo od funkcije, ki jo aproksimiramo. Vemo tudi, da druga norma najbolj kaznuje tako imenovane osamelce (ang. outliers). To je zelo pomembno pri reševanju našega problema, saj red velikosti posameznih oziroma naj- večjih napak, ne sme vplivati na optimalno rešitev – ključno je le skupno, najmanjše odstopanje. Statistično gledano prav tako vemo, da je druga norma najboljša izbira pri glajenju podatkov, ko so aditivne napake podatkov porazdeljene z normalno porazdelitvijo N(0, σ²).

Druga možna izbira je prva norma, ki nam preračuna celotno skupno odstopanje.

Če minimizramo prvo normo, bo algoritem na koncu vsakega koraka izbral kodon s takšnim translacijskim časom, da bo skupno odstopanje od želenih vrednosti, do tega koraka minimalno. Prva norma je ugodna, saj prav tako kaznuje osamelce.

Hkrati pa prva norma ni strogo konveksna, torej najboljša aproksimacija z uporabo prve norme ni nujno enolična. To nam predstavlja težavo, saj bi se lahko zgodilo, da bi prva norma rešila problem, vendar bi predlagala eno izmed več različnih poti za enako končno skupno odstopanje. To bi lahko vodilo v napačno izbiro zadnjega kodona v zaporedju in posledično, rekurzivno, v napačno, morda ne optimalno, sestavljeno zaporedje.

Zadnja očitna izbira bi bila neskončna norma. Neskončna norma se v praksi običajno uporablja, če so napake podatkov zelo majhne v primerjavi z aproksima- cijskimi napakami [5]. Neskončna norma nam minimizira maksimalno odstopanje med vsemi mesti v zaporedju. To se zdi manj ugodno, kot druga norma, saj je optimizacija algoritma precej bolj omejena. Če se zgodi, da bo na nekem mestu v zaporedju odstopanje veliko, se norma sploh ne bo ukvarjala z drugimi, in bodo tako lahko vsa odstopanja precej velika. Natančneje, ko določimo člen z največjim odstopanjem, imamo največ 5 kodonov s katerimi lahko kodon na danem mestu za- menjamo. Če nobena druga izbira ne ponudi manjšega odstopanja na danem mestu v zaporedju, je neskončna norma že našla svojo vrednost. To ni najbolj primerno, saj lahko obstaja mestoiv zaporedju, kjer absolutna razlika želenih in poračunanih časov ni največja v zaporedju, hkrati pa imamo na voljo drug kodon, ki bi zmanjšal odstopanje, tako i-tega mesta, kot celotne verige.

(11)

4. Dinamično programiranje

Dinamično programiranje je pristop k reševanju problemov. Problem zapore- doma razdelimo na manjše podprobleme. Ti podproblemi so lahko med seboj odvisni. Med postopkom iskanja rešitve z dinamičnim programiranjem rešimo manjše podprobleme in iz njihovih rešitev sestavimo rešitev celotnega problema [2, 4].

Dinamično programiranje temelji na pravilu optimalnosti, ki pravi, da mora biti vsak del optimalne rešitve tudi zase optimalen. Na osnovi pravila optimalnosti sestavimo rekurzivno enačbo, ki nam bo določila optimalno rešitev. Tej enačbi pravimo Bellmanova enačba. Vsakemu problemu, ki ga rešimo s pomočjo metode dinamičnega programiranja, pripada lastna Bellmanova enačba. Kako deluje metoda dinamičnega programiranja si bomo pogledali na primeru, ki je nekoliko podoben našemu, vendar enostavnejši.

4.1. Primer dinamičnega programiranja – segmentirani najmanjši kva- drati. Problem segmentiranih najmanjših kvadratov je znan problem v statistiki in podatkovnem rudarjenju. Problem je natančneje opisan v knjigi [3]. Pri toč- kah, (x₁, y₁),(x₂, y₂), . . . ,(x_n, y_n), narisanih v ravnini, običajno želimo najti premico najboljšega ujemanja. V tem primeru definiramo napako, ki jo pri aproksi- maciji naredi premica y = ax+b, kot vsoto kvadriranih razdalj točk od premice:

e=^∑︁ⁿ_i=1(y_i −ax_i−b)².

Pogosto se nam zgodi, da točk ne moremo povezati z eno samo premico, saj bi bila napaka pri tem očitno prevelika. Kot primer si oglejmo spodnjo sliko. Na sliki je precej očitno, da bi bila najboljša aproksimacija dosežena z dvema premicama.

Slika 1. Primer točk v ravnini, ki bi jih lahko aproksimirali z dvema premicama.

V našem zgledu bomo dovolili aproksimacijo z več premicami. Če imamo podano zaporedje točk, želimo torej odkriti tiste točke zaporedja, pri katerih se zgodi tako imenovana sprememba. Zaradi te spremembe se problemu reče tudi problem

(12)

zaznavanja sprememb. V našem primeru se bomo osredotočili na aproksimacijo točk z nekaj linearnimi premicami, zato je sprememba, ki jo iščemo, sprememba iz ene linearne aproksimacije v drugo.

Sedaj formalno opišimo problem. Podano imamo množico oziroma zaporedje točk P ={(x₁, y₁),(x₂, y₂), . . . ,(x_n, y_n)}, kjer je x₁ < x₂ <· · ·< x_n. S p_i označimo točko (x_i, y_i). Najprej bi radi razdelili množico P v neko manjše število segmentov oziroma odsekov. Vsak odsek mora biti podmnožica točk iz množice P oblike {pi, pi+1, . . . , pj−1, pj} za neka 1 ≤ i ≤ j ≤ n. Nato bi radi za vsak segment S iz naše razdelitve množice P na segmente izračunali premico, ki minimizira napako glede na točke v S. Pri tem naj bo:

• y=ax+b – aproksimacijska premica za posamezen segment,

• a= ⁿ

∑︁

ixiyi−(∑︁

ixi)(∑︁

iyi) n∑︁

ix²_i−(∑︁

ixi)² – smerni koeficient premice,

• b=

∑︁

iyi−a∑︁

ixi

n – odsek na ordinatni osi.

Hkrati definiramo kazen particije kot

kazen =nC +^∑︂

i

e_i,

kjer je n število segmentov, v katere razdelimo P, C dan fiksni koeficient in e_i vrednost napake optimalne premice skozi i-ti segment. Kazen definiramo, ker bi drugače lahko za vsak segment izbrali par zaporednih točk in potegnili premico med njima, ter tako povezali vse točke. S tem bi popolnoma izničili napako aproksimacijske premice, vendar tega ne želimo. Kazen nam tako omogoča, da s spreminjanjem vrednosti C kaznujemo večje število segmentov particije P. Naš cilj je minimizirati kazen. Poglejmo si, kako lahko ta minimum poiščemo z uporabo dinamičnega programiranja.

Kot vemo iz definicije dinamičnega programiranja, moramo problem razdeliti na več podproblemov, poiskati optimalno rešitev za vsakega izmed njih in iz njih z rekurzijo sestaviti optimalno rešitev originalnega problema. Pri iskanju rešitev z dinamičnim programiranjem nam je lahko v veliko pomoč, če opazimo kakšno zani- mivo lastnost problema. V danem primeru lahko hitro ugotovimo naslednje: zadnja točka p_n v optimalni particiji pripada točno enemu segmentu, ki se začne z neko prejšnjo točko p_i. Tako ugotovimo, da če poznamo zadnji segment točk p_i, . . . , p_n, lahko odstranimo te točke in rešujemo problem na preostalih točkah p₁, . . . , pi−1.

Z OPT(i) označimo optimalno rešitev za točkep₁, . . . , p_i, sCdan fiksni koeficinet in z e_i,j velikost napake za optimalno premico skozi točke p_i, . . . , p_j. Potem lahko po zgornjem premisleku zapišemo enačbo za robni primer i=n:

OPT(n) = e_i,n+C+ OPT(i−1).

Za robni primer i = 0 pa nastavimo OPT(0) = 0. Če ponovno uporabimo gornji premislek, lahko sedaj zapišemo splošno rekurzivno enačbo za naš problem:

OPT(j) = min

1≤i≤j{e_i,j+C+ OPT(i−1)}.

Zapisana rekurzivna enačba nam omogoča, da poiščemo optimalno rešitev. Rešitve podproblemov OPT(i) pridobivamo postopoma z večanjem i. Na vsakem koraku za vsak par i ≤ j izračunamo vse velikosti napake za optimalno premico e_i,j za segmentpi, . . . , pj. Na koncu vsakega koraka z izpeljano rekurzivno enačbo poiščemo

(13)

optimalno delitev do tega koraka. Seveda moramo biti pri sestavljanju algoritma posebej pozorni na robne primere.

Celotnega algoritma za rešitev problema ne bomo navedli, ker ni ključen za razumevanje reševanja problemov z dinamičnim programiranjem. Kar si moramo o dinamičnem programiranju zapomniti (in kar nam zelo lepo ilustrira gornji primer), je naslednje:

• pri reševanju problemov z dinamičnim programiranjem je ključna dobra de- finicija in podrobno razumevanje problema samega,

• iskanje vzorcev pred konstruiranjem algoritma nam lahko zelo olajša reševa- nje problema,

• pomembna je natančna specifikacija in premislek o rešitvi robnih primerov,

• želimo najti splošno rekurzivno formulo, ki bo rešila problem, in jo nato lahko implementiramo v računalniški algoritem – ko jo enkrat najdemo, je problem tako rekoč rešen.

4.2. Časovna zahtevnost. Časovna zahtevnost je podatek o tem, koliko časa se bo program oziroma algoritem izvajal pri danih vhodnih podatkih, preden bo vrnil rešitev. Časovno zahtevnost podamo kot funkcijo velikosti vhodnih podatkov (kot so na primer velikosti tabel). Poznamo tri vrste časovne zahtevnosti:

(1) f_B – najboljša možnost oziroma spodnja meja zahtevnosti, (2) f_W – najslabša možnost oziroma zgornja meja zahtevnosti, (3) f_E – pričakovana zahtevnost pri povprečnih podatkih.

Običajno se gleda najslabšo možnost, torej iščemo zgornjo mejo zahtevnosti. Oce- njujemo, koliko časa najdlje bo potreboval algoritem, da bo prišel do rešitve danega problema. Logično je, da nas to najbolj zanima, saj bi želeli zgornjo mejo kar se da minimizirati, da algoritem ne bi tekel predolgo časa, ne glede na konkretne podatke. Običajno je natančno število operacij zelo težko, oziroma celo nemogoče določiti. Zaradi tega uporabljamo tako imenovano O-notacijo, ki označuje red rasti problema. Definicijo O-notacije smo povzeli po [10].

Izrek 4.1. Naj bosta f in g funkciji, ki slikata iz naravnih števil v pozitivna realna števila. Pravimo, da je f(x) = O(g(x)), če obstajata α, n0 ∈ N taka, da za vsak n ≥n0 velja:

f(n)≤αg(n).

V spodnji tabeli so z omenjeno notacijo podane zahtevnosti od najugodnejše, do najneugodnejše. Oznaka n v tabeli nam predstavlja vrednost, ki oriše velikost podatkov (na primer število bitov potrebnih za predstavitev podatkov).

Notacija Zahtevnost O(1) konstantna O(logn) logaritemska

O(n) linearna O(nlogn) vmesna O(n²) kvadratna O(n^c), c >1 polinomska

O(cⁿ) eksponentna

Tabela 3. Tabela časovnih zahtevnosti.

(14)

5. Reševanje problema

Problem rešujemo z dinamičnim programiranjem, saj se zdi ta način najbolj smi- selen, da ga bo rešil pravilno ter z majhno časovno zahtevnostjo. Kot smo že videli v primeru uporabe dinamičnega programiranja, je zelo pomembno dobro razumevanje problema, kar smo uspešno dosegli prek podrobne analize biokemijskega ozadja.

Sledilo je iskanje zanimivih lastnosti problema, ki bi nam olajšale reševanje. To se je pri našem problemu izkazalo za precej težjo nalogo, kot pri zgledu dinamičnega programiranja. Na vsakem koraku mora namreč algoritem preveriti vse kombinacije za vse možne izbire kodona in si zapomniti najboljše ter na koncu izbrati tiste, pri katerih so odstopanja minimalna.

V splošnem je podrobnejša ideja za algoritem sledeča: Če je dolžina zaporedja m, je potrebno narediti m korakov. Na vsakem koraku je potrebno za vsak par števil iz nabora vrednosti za mestii−1 in ipoiskati tako število iz nabora vrednosti za mesto i − 2, da bo odstopanje časa translacije, izračunanega po formuli, od želene vrednosti čim manjše. Nato si je treba vse dobljene možne vrednosti in posamezna odstopanja za določeno mesto zapomniti ter enak postopek ponoviti na naslednjem mestu (torej za vsako kombinacijo vrednosti iz mest i in i+ 1, določiti vse najboljše vrednosti za mesto i−1). Pri tem je treba omeniti, da je podatke najbolje oštevilčiti tako, da imamo v vsakem koraku vse potrebne vrednosti (v i- tem koraku imamo znane vse vrednosti do a_i ∈ A_j_i, torej lahko s tem izračunamo največYi−1). Vi-tem koraku torej računamo odstopanje na mestui−1. Ko pridemo do konca verige, izberemo tisti a_m ∈A_j_m, za katerega je odstopanje najmanjše, kar bo implicitno določalo optimalni am−1 ∈ A_j_m−1. Tako lahko rekurzivno zgradimo celotno optimalno zaporedje izbranih časov translacije za celotno verigo.

5.1. Rekurzivna formula. Zgornjo idejo bomo uporabili za izpeljavo rekurzivne formule, ki jo bo mogoče implementirati v računalniškem programu. Kot vhodne podatke algoritma imamo za zaporedje dolžine m podano:

• zaporedje aminokislin j₁, j₂, . . . , j_m,

• zaporedje želenih časovX₁, X₂, . . . , X_m (eksperimentalno določeni).

Z o^b,c_i bomo v nadaljevanju označevali odstopanje zaporedja doi-tega člena v zaporedju, če na i-tem koraku izberemo kombinacijo kodonov (b, c), z a^b,c_i pa optimalen a, ki ga pri kombinaciji kodonov (b, c) izbere algoritem na i-tem koraku. Ideja za rekurzivno formulo je, da na vsakem koraku gledamo minimalno odstopanje zaporedja, o^b,c_i−1, torej do člena i−1 v zaporedju. Pri tem se odstopanje računa po prvi normi in znani formuli za izračun časa translacije kodona na mestu i−1:

o^b,c_i−1 = min{|Xi−1− a+b+c

3 |+o^a,b_i−2; a ∈A_j_i−2, b ∈A_j_i−1, c∈A_j_i, i= 2, . . . , m}

Razlaga simbolov:

• X_i ∈R – eksperimentalno določeni želeni čas,

• Y_i = ^a+b+c₃ – izračunan translacijski čas (vendar bomo v splošnem uporabljali kar formulo na desni strani enačaja, saj so tako enačbe bolj intuitivne),

• A_j_i – množica vseh kodonov oziroma njihovih časov, ki zapisujejo za i-to aminokislino v zaporedju iz množice aminokislin j_i ∈ {1,2, . . . ,20},

• a – translacijski čas iz množice A_j_i−2 (in enako b ∈A_j_i−1 in c∈A_j_i).

Opomnimo, da bi lahko v rekurzivni formuli uporabili katerokoli izmed omenjenih norm, vendar se je pri prvotnem reševanju problema prva zdela najenostavnejša za

(15)

vključitev v algoritem in je hkrati tudi najbolj pregledna izmed vseh. V tem razdelku bomo eksplicitno vedno uporabljali prvo normo, vendar se zavedamo, da enake formule, s primernimi prilagoditvami, držijo tudi za drugo in neskončno normo. Sle- dnji bomo v algoritem ponovno eksplicitno vključili v šestem poglavju, pri testiranju algoritma.

Algoritem bomo razdelili na dva dela – prvi je algoritem, ki na vsakem koraku poišče vsa odstopanja, o^b,c_i , in vse potencialno optimalne vrednosti, a^b,c_i in ga poimenujemo Algoritem za optimizacijo rabe kodonov. Drugi je algoritem, ki na koncu sestavi zaporedje in ga poimenujemo Sestavi zaporedje. Preden navedemo prvi algoritem, omenimo še, da so množice A_j₁, . . . , A_j_m znotraj vsakega koraka konstantne.

Najprej navedimo Algoritem za optimizacijo rabe kodonov v osnovni verziji, z najbolj očitnimi rekurzivnimi formulami.

Algoritem 1 Algoritem za optimizacijo rabe kodonov - osnovni.

1: procedure Dinamični kodon(X, j₁, . . . , j_m)

2: for i∈(3, m)do

3: for c∈A_j_i do

4: for b ∈A_j_i−1 do

5: a =a₁ ∈A_j_i−2

6: o^b,c_i−1 =|Xi−1−^a+b+c₃ |+o^a,b_i−2

7: for a∈A_j_i−2;a̸=a₁ do

8: if |Xi−1− ^a+b+c₃ |+o^a,b_i−2 < o^b,c_i−1 then

9: o^b,c_i−1 =|Xi−1− ^a+b+c₃ |+o^a,b_i−2

10: a^b,c_i =a

11: end if

12: end for

13: Potencialni_i[id_b][id_c] =a^b,c_i

14: end for

15: end for

16: end for

17: end procedure

Med zadnjo zanko f or se skonstruira matrika Potencialni_i. Da lahko smislno vpeljemo to matriko, moramo najprej kodone znotraj vsake izmed množic A_j_i ošte- vilčiti z zaporednimi številkami. Oznaka id označuje zaporedno številko kodona, znotraj posamezne množice, ki ji pripada. Za dano kombinacijo kodonov lahko tako na mesto [id_b][id_c] v matriki zabeležimoa^b,c_i , ki nam pri tej kombinaciji da najmanjše odstopanje. Za lažjo predstavo si poglejmo, kako bi matrika izgledala na spodnjem primeru. Predpostavimo b ∈A_j_i−1 in c∈A_j_i, pri čemer je |A_j_i−1|= 4 in |A_j_i|= 3.

Potencialni_i =

⎡

⎢

⎣

a^(b_i ¹^,c¹⁾ a^(b_i ¹^,c²⁾ a^(b_i ¹^,c³⁾ a^(b_i ²^,c¹⁾ a^(b_i ²^,c²⁾ a^(b_i ²^,c³⁾ a^(b_i ³^,c¹⁾ a^(b_i ³^,c²⁾ a^(b_i ³^,c³⁾ a^(b_i ⁴^,c¹⁾ a^(b_i ⁴^,c²⁾ a^(b_i ⁴^,c³⁾

⎤

⎥

⎦

.

(16)

V splošnem bi torej matriki Potencialni_i in Potencialnii−1 izgledali kot sledi:

Potencialni_i =

⎡

⎢

⎣

c1 ... cj ... c6

b1 ...

... ...

bi . . . . . . . . . a^(b_i ⁱ^,c^j⁾ ...

b6

⎤

⎥

⎦

,

Potencialnii−1 =

⎡

⎢

⎣

b1 ... bi ... b6

a1 ...

... ...

a⁽_i^bi,cj⁾ . . . . . . a^(a

(bi,cj) i ,bi)

.. i−1

.

a6

⎤

⎥

⎦

.

Opombe k algoritmu za optimizacijo rabe kodonov. Osnova algoritma je, kot že omenjeno, sestavljena iz zunanje zanke, ki teče po zaporedju aminokislin dolžine m. Ker moramo na vsakem koraku poznati vsaj 3 zaporedne vrednosti, mora algori- tem teči od tretjega člena naprej. Rezultat za robne vrednosti se premisli na podoben način, kot za splošni algoritem, vendar za vsako vrednost posebej. Nato sledijo tri notranje zanke, vsaka teče po vseh različnih razpoložljivih vrednostih translacijskih časov za dano aminokislino na tem mestu. V i-tem koraku poračunamo odstopanja in potencialne optimalne vrednosti za mesto i−1 v zaporedju. Znotraj tretje zanke izračunamo o^(b,c)_i−1, ki nam predstavlja odstopanje translacijskega časa aminokisline od želene vrednosti, če na treh zaporednih mestih izberemo par kodonov (b, c) v kombinaciji s kodonom a. Na vsakem koraku moramo izračunati vse možne kombi- nacije. Ker je število različnih kodonov, ki zapisujejo za neko aminokislino, navzgor omejeno s 6, to v najslabšem primeru pomeni 216 kombinacij, kar v smislu dina- mičnega programiranje ni veliko. V algoritmu smo vrednosti preštevilčili na način, da imamo do i-tega koraka znano vse, kar potrebujemo za izračun odstopanja na mestu i−1 v zaporedju, saj je taka formulacija iz programerskega vidika in zaradi aplikacije algoritma precej bolj smiselna in intuitivna.

Na koncu vsake iteracije algoritma shranimo dve spremenljivki:

• v o^(b,c)_i shranimo najmanjše odstopanje Y_i za dan X_i pri nekem paru (b, c),

• v a^(b,c)_i shranimo tak a za katerega je odstopanje med Y_i in X_i pri danem paru (b, c) najmanjše.

Matriko Potencialni_i definiramo zato, da lahko vse dobljene potencialno optimalne a zai-to iteracijo shranimo na isto mesto, saj bomo to pri rekurzivnem klicu potre- bovali.

Oglejmo si še časovno zahtevnost algoritma. Na prvi pogled opazimo, da ima algoritem več vgnezdenih zank, zaradi katerih bi lahko predvidevali, da bo časovna zahtevnost algoritma velika. Vendar ob natančnejši analizi algoritma opazimo, da

(17)

je število korakov navzgor omejeno z neko konstanto. V vsakem koraku algoritma, oziroma v vsaki zanki f or, imamo namreč na voljo največ 6 različnih časov – vsaka aminokislina v zaporedju, ki jo predstavimo z množico možnih časov translacije A_j_i, A_j_i−1 in A_j_i−2, ima največ 6 različnih kodonov, ki zapisujejo zanjo. Torej so vse tri notranje zanke navzgor omejene z 216 različnimi kombinacijami. Tako lahko časovno zahtevnost algoritma ocenimo kot linearno, 216×O(m).

Za konec ponovno omenimo, da bi Algoritem 1 enostavno lahko prevedli na drugo ali neskončno normo. Pri tem bi se spremenil le način računanja odstopanja na vsakem koraku, vse ostalo pa bi ostalo enako. To pomeni, da bo imel algoritem linearno časovno zahtevnost ne glede na izbiro norme.

Poglejmo si sedaj še drugi algoritem Sestavi zaporedje. To je rekurzivni algoritem, ki nam na koncu, ko imamo poračunane vse vrednosti odstopanjao^b,c_i in potencialne optimalne vrednosti a^b,c_i , vrne optimalno zaporedje.

Algoritem 2 Sestavi zaporedje.

1: i=m

2: procedure zaporedje(b, c, i)

3: if i=m then

4: najboljši =a^b_m

5: else

6: najboljši = Potencialni_i[id_b][id_c]

7: end if

8: i−= 1

9: optimalno_zaporedje += najboljši

10: return ZAPOREDJE(najboljši, b, i)

Opombe k algoritmu Sestavi zaporedje. Algoritem se začne na koncu verige, kjer upoštevamo le zadnji, m-ti člen. Najprej določimo najboljši a tako, da pogle- damo vse dobljene a na zadnjem koraku in izberemo tistega, pri katerem je skupno odstopanje od želene vrednosti najmanjše. Ta bo zaradi Algoritma za optimizacijo rabe kodonov shranjen v spremenljivki a^b_m. Algoritem nato deluje rekurzivno – na vsakem koraku za kodon b vzame najboljši aiz prejšnje iteracije, za kodon cpab iz prejšnjega koraka (torej najboljšiaiz predzadnjega koraka). To najlažje naredi prek številčenja kodonov in matrike Potencialni_i, kar smo že vpeljali. Zaporedje, dobljeno s tem algoritmom, je optimalna izbira kodonov za dano zaporedje aminokislin, da bodo čim manj odstopale od želenih časov.

5.2. Računalniška implementacija algoritma. Splošno idejo v obliki rekurziv- nih formul želimo preoblikovati v algoritem tako, da se ga bo dalo uporabiti tudi v praksi oziroma najlažje zapisati v kakšnem računalniškem programu. S tem name- nom smo se odločili napisati kodo za računalniško implementacijo algoritma. Pri tem je prišlo do nekaj manjših, vendar ključnih sprememb – prilagoditev vhodnih podatkov in sprememba shranjevanja in modeliranja z izbranimi vrednosti časov translacije kodonov ter izračunanimi vrednostmi odstopanj. Poglejmo si najprej prilagoditev vhodnih podatkov.

Program ima vključen seznam, ki ga poimenujemo Seznam aminokislin, vseh aminokislin, razporejenih po oznakah. Poleg tega seznam vključuje še kombinacije kodonov, ki zapisujejo za posamezno aminokislino in pripadajoče čase translacije za

(18)

posamezen triplet kodonov. Za jasnejšo predstavo je tu vključena celotna tabela, v katero smo dodali še polna imena aminokislin. Tabela je razvrščena po angleškem abecednem vrstnem redu glede na oznake aminokislin, saj se za njih tudi pri nas v praksi uporabljajo nekatere angleške črke.

Seznam aminokislin Ime

aminokisline

Oznaka aminokisline

Kodoni, ki zapisujejo za aminokislino

Pripadajoč čas translacije

Alanin A GCA 0.0232

Alanin A GCC 0.0962

Alanin A GCG 0.0236

Alanin A GCT 0.0239

Cistein C TGC 0.0379

Cistein C TGT 0.0383

Aspartat D GAC 0.0256

Aspartat D GAT 0.0270

Glutamat E GAA 0.0157

Glutamat E GAG 0.0158

Fenilalanin F TTC 0.1235

Fenilalanin F TTT 0.1316

Glicin G GGA 0.0433

Glicin G GGC 0.0179

Glicin G GGG 0.0272

Glicin G GGT 0.0186

Histidin H CAC 0.0840

Histidin H CAT 0.0847

Izolevcin I ATA 0.1818

Izolevcin I ATC 0.0168

Izolevcin I ATT 0.0174

Lizin K AAA 0.0813

Lizin K AAG 0.0820

Levcin L CTA 0.0833

Levcin L CTC 0.0532

Levcin L CTG 0.0166

Levcin L CTT 0.0552

Levcin L TTA 0.0645

Levcin L TTG 0.0226

Metionin M ATG 0.0917

Asparagin N AAC 0.0602

Asparagin N AAT 0.0610

Prolin P CCA 0.2564

Prolin P CCC 0.0690

Prolin P CCG 0.1031

Prolin P CCT 0.0595

Glutamin Q CAA 0.0625

Glutamin Q CAG 0.0613

Arginin R AGA 0.0629

(19)

Arginin R AGG 0.0943

Arginin R CGA 0.0159

Arginin R CGC 0.01655

Arginin R CGG 0.1075

Arginin R CGT 0.01658

Serin S AGC 0.0599

Serin S AGT 0.0602

Serin S TCA 0.0382

Serin S TCC 0.0813

Serin S TCG 0.0333

Serin S TCT 0.0307

Trenonin T ACA 0.0441

Trenonin T ACC 0.0617

Trenonin T ACG 0.0304

Trenonin T ACT 0.0311

Valin V GTA 0.0207

Valin V GTC 0.0433

Valin V GTG 0.022

Valin V GTT 0.0177

Triptofan W TGG 0.0461

Tirozin Y TAC 0.0295

Tirozin Y TAT 0.0299

Tabela 4. Tabela aminokislin s pripadajočimi kodoni in časi translacije.

V tabeli so navedeni tripleti s timinom, kar se seveda zelo enostavno zamenja z ura- cilom, če je to potrebno (odvisno v kakšni obliki dobimo podatke za aminokisline).

Iz biokemijskih razlogov bi bilo pomembno, da bi v zaporedju kodonov nastopal uracil, saj je pomemben prehod na enovijačno mRNA, ki jo zna prebrati ribosom.

Vendar je za potrebe formuliranja računalniškega algoritma vseeno ali v podatkih nastopa T ali U, saj dogajanje v računalniškem programu ni odvisno od biokemijskih lastnosti, važno je le, da v podatkih nastopa bodisi T, bodisi U. Poleg tega, če podrobneje preučimo tabelo, lahko opazimo, da je v njej navedenih le 61 tripletov.

To je zato, ker trije izmed kodonov zapisujejo za aminokislino, ki predstavlja konec proteina in se pri računanju ne upošteva, torej za naš problem niso pomembne. Ob takšni tabeli so torej vhodni podatki za algoritem naslednji:

• AK_zaporedje– zaporedje aminokislin dolžine m, podano z oznakami ami- nokislin,

• X – zaporedje želenih časov X₁, X₂, . . . , X_m (eksperimentalno določeni).

Zgornji Seznam aminokislin pa nam služi kot parameter, ki je pri nas vedno enak.

Drugo spremembo, shranjevanja in modeliranja vrednosti, si bomo ogledali nepo- sredno v algoritmu, oziroma jo komentirali v nadaljevanju. Na tem mestu ponovno omenimo, da lahko algoritem deluje le za i∈(3, m), ker potrebujemo vsaj 3 zaporedne kodone, da lahko izračunamo vse vrednosti. Robne primere premislimo posebej, kar bo prav tako komentirano v nadaljevanju. Omenimo še, da bo v algoritmu ponovno implementirana le prva norma, saj se je njena predstavitev znotraj algoritma zdela najbolj smiselna. Vključitev druge in neskončne norme je precej podobna, zato ju ne bomo posebej navajali.

(20)

Algoritem 3 Algoritem za optimizacijo rabe kodonov – podroben.

1: procedure Dinamični kodon(AK_zaporedje, X, Seznam aminokislin)

2: for i∈(3, m)do

3: for a∈A do

4: for b ∈B do

5: for c∈C do

6: 3D_tabela[ida][idb][idc] = 2D_tabela[ida][idb] +|Xi−1−^a+b+c₃ |

7: end for

8: end for

9: end for

10: for b ∈B do

11: for c∈C do

12: 2D_tabela[id_b][id_c] = min(3D_tabela[:,id_b,id_c])

13: AK_tabela[id_b][id_c] = argmin(3D_tabela[:,id_b,id_c])

14: end for

15: end for

16: Shrani 2D_tabelo in AK_tabelo za naslednjo iteracijo

17: end for

Opombe k algoritmu za optimizacijo rabe kodonov. Zunanja zanka ostaja enaka kot prej – algoritem teče po zaporedju aminokislin. Na vsakem koraku nato skonstruira tabelo 3D_tabela, velikosti A×B ×C, pri čemer so A, B in C moči množic vseh različnih časov kodonov, ki zapisujejo za neko aminokislino. Z oznakami id_a,id_b in id_c označimo indekse vrednosti a, b in c v posamezni tabeli. To nam omogoča, da je tabela 3D_tabela pravih dimenzij in da vanjo zapišemo primerne vrednosti na pravo mesto. Program namreč skonstruirano tabelo 3D_tabela napolni z odstopanji od želenih časov, če za vsako kombinacijo (a, b) ∈ A×B, izberemo vsako od vrednostic∈C. Na tej točki si pomagamo s tabelo 2D_tabela v kateri so shranjena odstopanja iz prejšnje iteracije. Pri tem opomnimo, da je Yi−1 = ^a+b+c₃ .

Nadalje skonstruiramo tabeli 2D_tabela in AK_tabela, enakih dimenzij (B × C). Nato za vsako kombinacijo (b, c) v tabeli 3D_tabela poiščemo takšen a, pri katerem je odstopanje za fiksno kombinacijo (b, c) najmanjše. Na mesto (b, c) v tabeli 2D_tabela shranimo vrednost a, v tabeli AK_tabela pa id_a. Na koncu vsake iteracije si moramo na poseben seznam (v prejšnjem algoritmu Potencialnii−1) shraniti vse najboljše izbire za a in njihove indekse – torej tabeli 2D_tabela in AK_tabela.

Robne primere moramo premisliti posebej. Po krajšem premisleku se izkaže, da lahko uporabimo popolnoma enak pristop, kot pri splošnem i, le da bodo dimenzije tabel primerno manjše – če imamo na voljo le 2 kodona, izračunamo recimo Y₁ =

a1+a2

2 ter podobno Y_m = ^a^m−1₂^+a^m. Problematičen je tudi prvi korak v iteraciji, saj odstopanja za i = 2 še niso znana, torej ne bi mogli izračunati koraka i = 3 po formulaciji algoritma, zato izračunamo posebej še Y₂ = ^a¹^+a₃²^+a³, kar pa je točno enako kot v splošnem algoritmu. V primerui= 1 ini=mbomo imeli torej namesto tabele 2D_tabela pravzaprav vektor (tabelo z enim stolpcem), kar je logično in tudi izjemno ugodno, saj lahko na ta način v zadnjem koraku dejansko izberemo kodon, ki nam bo dal najmanjše odstopanje. To pa naredimo z že znanim algoritmom Sestavi zaporedje, ki je bil opisan v prejšnjem razdelku in ga tu ne bomo ponavljali.

(21)

6. Testiranje algoritma in rezultati

Končno imamo natančno opisan celoten problem in razčlenjen algoritem, torej ga lahko sedaj implementiramo v programskem jeziku. Program smo preizkusili v programskem jeziku Python. Testiranje algoritma se je zdelo smiselno razdeliti na več delov. V prvem delu si bomo na nekaj naključnih, izmišljenih proteinih pogledali pravilnost delovanja algoritma. Nato bomo enako naredili še na pridobljenih podatkih treh dejanskih proteinov, ki se v praksi pogosto uporabljajo. Sledilo bo testiranje algoritma na veliki količini naključno generiranih verig. Želeni časi bodo pri teh verigah pridobljeni iz konkretnega zaporedja kodonov za vsako verigo. Na koncu si bomo pogledali še, kaj se zgodi, če v omenjene želene čase (podatke) vpeljemošum ter algoritem ponovno preizkusili. Pri vseh testih nas je zanimal predvsem čas trajanja algoritma (čas, ki ga algoritem potrebuje, da najde rešitev za določeno verigo), uspešnost algoritma (kolikšen delež kodonov v verigi je pravilno določil), odzivanje algoritma na s šumom zamaknjene podatke, ter vpliv različnih norm na algoritem.

Pripomnimo še, da smo drugo normo vključili kot vsoto kvadratov odstopanj (brez korenjenja), neskončno normo pa standardno, po definiciji.

6.1. Osnovno preverjanje uspešnosti in začetno testiranje algoritma. Naj- prej smo generirali podatke za tri različne izmišljene proteine in algoritem testno preizkusil na njih. Da program ne bi tekel predolgo časa, smo si za začetek izbrali nekoliko krajše proteine, z dolžino verige med 60 in 80 kodoni. Da bi dejansko preverili pravilnost delovanja programa, smo podatke izbrali kot sledi:

• zaporedje aminokislinA_j_i – zaporedje nukelinskih kislin za izbrane proteine,

• časi translacije kodonov – eksprimentalno določeni, pobrani iz tabele Seznam aminokislin,

• želeni časi X_i – preračunani kot lokalno povprečje treh zaporednih časov translacije kodonov.

Za vrednosti X_i smo torej namesto naključnih oziroma poljubnih vrednosti za potrebe preverjanja pravilnosti delovanja algoritma uporabili funkcijo, ki je kot ar- gument prejela opazovano zaporedje kodonov in zanj preračunala lokalna povprečja za vsako mesto. Tako je funkcija pravzaprav določila kakšne čase bo moral izra- čunati algoritem, da bo lahko izbral točno tisto zaporedje kodonov, ki ga iščemo.

Naloga algoritma je bila nato s postopki dinamičnega programiranja pregledati vse opcije, izračunati vse potrebne čase in odstopanja ter nato določiti, kateri kodon izbrati na katerem mestu, da bo odstopanje čim manjše.

Pravilnost algoritma smo preverjali tako, da smo primerjali želene čase X_i z izbranimi časi iz tabele Seznam aminokislin in iz izračunanimi vrednostmi Yi (po tem, ko je algoritem že sestavil zaporedje). Če sta se časa nai-tem mestu popolnoma ujemala, torej sta bila identična, kar je pomenilo da je algoritem res izbral pravi kodon (po predpostavki nobena dva kodona nimata enakega časa translacije), smo to beležili kot uspeh, če ne, pa kot neuspeh, torej napačno izbrani kodon. Na koncu smo preračunali delež uspešno izbranih kodonov.

Pri osnovnem testiranju funkcionalnosti algoritma (s tremi različnimi izmišljenimi proteini) je prišlo do dveh manjših težav, kar se je takoj opazilo, saj je imel program približno 95 % natančnost. Ob natančnejšem pregledu kode, smo ugotovili, da se je prva težava pojavila pri izbiri prvega, začetnega kodona, torej pri robnem primeru. Napaka na sam algoritem ni imela vpliva, ker se bolj pojavljala kot napaka

(22)

v programu samem – zaradi načina zapisa kode je program občasno izbral napačen začetni kodon.

Druga napaka, ki jo lahko razdelimo na dva dela in je bolj zanimiva, se je pojavila zaradi časov translacije pri nekaterih kodonih. V obeh primerih je šlo za problem pri zaokroževanju, enkrat na preveč, drugič na premalo decimalnih mest. V primeru, ko se je v danem zaporedju pojavil kodon GTG, ki ima dan čas translacije 0.022, je pri načinu računanja želenih translacijskih časov zanj prišlo do odstopanj, saj je program preračunal čas kot 0.02200002, kar je nato pri preverjanju pravilnosti izbire vodilo v napako. Ko smo popravili zaokrožanje na manjše število decimalnih mest, se je pojavil drug problem podobnega tipa. Kodona CGC in CGT, ki oba zapisujeta za isto aminokislino, arginin, imata skoraj identičen čas translacije, prvi 0.01655, drugi 0.01658. Ker tega sprva nismo opazili in smo zaokroževali na 4 decimalna mesta, je nato algoritem, ko je v zaporedju naletel na enega izmed njiju, vedno izbral prvega, CGC. Zaradi podatkov je bilo torej zelo pomembno zaokroževanje na primerno število mest. Po nastavitvi zaokroževanja translacijskih časov na 5 decimalnih mest, je bil algoritem v vseh treh testnih primerih 100 % uspešen.

6.2. Testiranje algoritma na realnih podatkih. Preden smo se lotili množič- nega testiranja algoritma, smo ga želeli preizkusiti tudi na resničnih proteinih. Po- datke o dejanskih proteinih smo pridobili iz spletne strani UniProt [9]. Izbrali smo si tri proteine, ki se v realnosti pogosto uporabljajo – laktozni transporter (ang.

lactose permease), ki je transmembranski protein, ki omogoča prehajanje laktoze skozi celično membrano, represor LacI (ang. DNA-binding transcriptional repressor LacI), ki je DNA vezavni protein, ki inhibira transkripcijo genov za uporabo laktoze, kot vira hranila ter fruktokinazo (ang. fructokinase), ki katalizira prenos fosfatne skupine iz adenozin trifosfata (ATP) na fruktozo.

Tudi pri vseh treh resničnih proteinih se je algoritem izkazal za 100 % uspešnega – točno je zadel vse kodone pri vseh treh proteinih. Poleg tega smo primerjali čase trajanja algoritma vseh treh znanih proteinov v odvisnosti od uporabljene norme.

Graf je prikazan na sliki 2.

Slika 2. Graf časov trajanja algoritma pri znanih proteinih v odvisnosti od norme.