0

Universitatea Politehnica din București

Facultatea de Automatică și Calculatoare

Departamentul de Automatică şi Ingineria Sistemelor

LUCRARE DE DIPLOMĂ

Aplicații ale modelelor Markov ascunse în analiza secvențelor biologice

Absolvent Conducător științific

Istrate Alina-Maria Prof. univ. dr. ing. Cătălin Buiu

București 2013

1

CUPRINS

1. Introducere 2

2. Modele Markov Ascunse (MMA) 4

2.1. Lanț Markov de ordin k 4 2.2. Ce este un Model Markov Ascuns 4

2.3. Arhitectura Modelelor Markov Ascunse 7

2.4. Cauze și efecte 7 2.5. Problema evaluării și algoritmul Forward 9 2.6. Problema decodării și algoritmul Viterbi 9 2.7. Alinierea de secvențe folosind Modele Markov Ascunse 10 2.8.Găsirea de gene cu ajutorul modelelor Markov Ascunse 11 2.9.Modele Markov Ascunse Circulare 12 2.10. Modelul Markov Ascuns Generalizat (MMAG) 13 2.11. Familii de proteine 14 2.12. PFAM 15

3. Studiu de caz. Aplicații ale Modelelor Markov Ascunse în studiul 16

receptorilor de miros 3.1 Receptorii de miros 16 3.2 MMA pentru segmentarea receptorilor de miros 19

3.3 Modele Markov ascunse de profil 32 3.4 Arbore filogenetic 50

4.Concluzii și dezvoltări posibile 52

5.Bibliografie 59

2

Introducere

Bioinformatica este un domeniu interdisciplinar, ce necesită competențe în statistică, informatică, inginerie, chimie analitică, biologie și fizică. Ea este folosită în biologia moleculară, genomică, proteomică, farmacologie, medicină (Pavel Ana Brândușa, Vasile Cristian, Buiu Cătălin. 2011).

Genomica se ocupă cu studiul genomurilor organismelor vii. Genomul

reprezintă totalitatea informației ereditare a unui organism. El este reprezentat fie de molecule de ADN, fie de molecule de ARN (Pavel Ana Brândușa, Vasile Cristian, Buiu Cătălin. 2011).

Proteomica se ocupă cu studiul proteinelor, a structurii și a funcției acestora.

Proteinele sunt molecule ce intră în alcătuirea organismelor vii, fiind componentele principale ale celulei .

Molecula de ADN este numită și molecula vieții deoarece ea conține informația

genetică a organismelor vii. Aceasta are 2 roluri esențiale : codifică proteinele și transmite informația genetică de la o generație la alta. Molecula de ADN are o structură de dublu helix. O nucleotidă din ADN este formată dintr-o pentoză, o bază azotată și o grupare fosfat. Bazele azotate ce intră în alcătuirea moleculei de ADN sunt : Adenina, Guanina,Timina și Citozina (Pavel Ana Brândușa, Vasile Cristian, Buiu Cătălin. 2011).

Proteinele și acizii nucleici sunt două molecule extrem de importante în

bioinformatică. Proteinele sunt formate din aminoacizi. Aminoacizii care intră în componența proteinelor sunt cei 20 de aminoacizi uzuali.

Bioinformatica include studierea structurilor și funcțiilor, evoluția genelor, a

proteinelor și a întregului genom. De o importanță fundamentală în bioinformatică sunt modele matematice concepute pentru a modela secvențe biologice.

Secvențele de ADN sunt adesea modelate ca fenomene probabilistice. O

secvență ADN este modelată ca un lanț Markov de variabile aleatoare care iau valorile (A, C, T, G) (Kal Renganathan Sharma. 2009).

Modelele Markov Ascunse sunt modele probabilistice și sunt construite cu

ajutorul unor concepte precum probabilitatea condiționată. Ele sunt folosite în foarte multe aplicații din bioinformatică.

Sistemul olfactiv include o mare familie de proteine numită receptori de miros. Ei fac posibilă recunoașterea a peste 10 000 de mirosuri diferite. Acești receptori de miros sunt atașați la suprafața celulelor în pasajul nostru nazal. Receptorii de miros sunt cea mai mare familie de gene din genomul uman, cu aproximativ 1000 de membrii. Pentru ca receptorii de miros (OR) să simtă moleculele din afara celulei și să semnalizeze interiorul celulei, trebuie să traverseze membrana celulelor (Richard Axel, Linda Buck, Fred Hutchinson. 2004). Repertoriul uman este cel mai mic dintre cele 4 repertorii complet cunoscute ale mamiferelor (om, câine, șobolan, șoarece). Repertoriul câinelui este mai mare decât cel al șobolanului. Numărul mare de subfamilii din repertoriul câinelui (300) indică un nivel de diversificare mare (Quignon,

3

Pascale, Mathieu Giraud, Maud Rimbault, Patricia Lavigne, Sandrine Tacher, Emmanuelle Morin, Elodie Retout . 2005.).

Modelele Markov Ascunse sunt folosite pentru segmentarea proteinei receptorilor de miros. Modele Markov Ascunse de profil codifică informații despre frecvența anumitor aminoacizi, precum și frecvența inserărilor și ștergerilor din aliniere. Ele sunt folosite pentru a decide dacă o secvență a receptorului de miros aparține unei familii sau nu.

În capitolul 2 Modelele Markov Ascunse sunt explicate în detaliu. Se definește pentru început un Lanț Markov, pentru ca mai apoi să se facă trecerea la Modelele Markov Ascunse. Sunt prezentate exemple cu construcția Modelelor Markov de ordin 0,1,2 și 3 precum și arhitecturile acestor modele. În cadrul acestui capitol sunt discutate atât problemele de evaluare și decodare cât și algoritmii Viterbi și Forward. Aplicații ale Modele Markov Ascunse, precum indentificarea de gene, Modele Markov Ascunse circulare folosite pentru a prezice periodicitatea în ADN și MMA generalizat, sunt și ele discutate. Toate aceste subiecte pot fi găsite în cărți de specialitate precum (Richard Axel , Linda Buck , Fred Hutchinson. 2004), (Selzer Paul, Marhöfer Richard, Rohwer Andreas. 2008), (Sperschneider Volker. 2008), (Stephen A. Krawetz, David D. Womble. 2003), (Kal Renganathan Sharma. 2009), (Warren Ewens, Grant Gregory. 2005).

Capitolul 3 cuprinde un studiu de caz : Aplicații ale Modelelor Markov Ascunse

în studiul receptorilor de miros. Sunt folosite secvențele de ADN de la câine, șobolan și om. Pentru început sunt definiți receptorii de miros apoi se folosește un Model Markov Ascuns de două stări pentru segmentarea proteinei receptorilor de miros. Segmentarea este o modalitate de a descoperii regiunile hidrofobice și hidrofilice caracteristice proteinei 7-TM. În cele din urmă se folosesc Modele Markov Ascunse de profil pentru a testa dacă o secvență face parte sau nu dintr-o anumită familie. Sunt prezentate rezultatele modelării acestor trei secvențe. Acest subiect poate fi găsit în cărți de specialitate precum (Richard Axel , Linda Buck , Fred Hutchinson. 2004 ), (Pascale Quignon, Giraud Mathieu, Rimbault Maud, Lavigne Patricia, Tacher Sandrine, Morin Emmanuelle, Retout Elodie. 2005.), (Zozulya Sergey, Echeverri Fernando, Nguyen Trieu. 2001).

Capitolul 4 este destinat concluziilor și dezvoltărilor posibile. O primă concluzie este aceea că MMA este o unealtă puternică în analiza secvențelor biologice. Sunt evidențiate diferențele dintre repertoriile de gene olfactive de la om, câine și șobolan. S-a observat că numărul genelor codificatoare este semnificativ mai mare la câine și șobolan decât la om. A fost testată și existența unor gene identice în cele trei repertorii. S-a constatat că nu există gene identice la cele trei specii. De asemenea, a fost creat și dotplot pentru perechi din secvențele analizate. O dezvoltare posibilă este repetarea analizei prezentată în lucrare și pentru repertoriul de gene olfactive de la maimuță.

4

2 Modele Markov Ascunse

2.1.Lanț Markov de ordin k

Un lanț Markov este o secvență de variabile aleatoare ale căror probabilități depind de valoarea numărului la momentul anterior. Parametrul de control într-un lanț Markov este probabilitatea de tranziție. Aceasta este o probabilitate condiționată pentru ca, având în vedere starea actuală, sistemul să treacă într-o stare nouă. Într-

un lanț Markov de ordin k, distribuția lui depinde de valorile k precedente (Kal Renganathan Sharma. 2009.). Probabilitatea de tranziție a lui (2.1)

Probabilitățile de tranziție pentru , într-un Model Markov de ordin 1, vor depinde doar de o valoare anterioară . Probabilitățile de tranziție pentru , într-un model Markov de ordin 0, nu vor depinde de valorile anterioare. Similar, într-un

model Markov de ordin 2, probabilitățile de tranziție pentru vor depinde de două valori anterioare si . Un lanț Markov de ordin k este staționar pentru toți t și u : = (2.2)

Aceasta este, pentru lanțul Markov staționar, distribuția lui , independentă de valoarea lui t dar dependentă de variabilele k anterioare.

Pentru a reprezenta structura secvenței primare de ADN cu perechile de bază

de început și de sfârșit, probabilitățile de tranziție pentru un model Markov de ordin 1 sunt reprezentate sub forma unei diagrame. Această diagramă este un graf orientat cu conexiuni nenule și poate fi numită arhitectura unui lanț Markov. Săgețile

indică următoarea stare a perechii de bază (Kal Renganathan Sharma. 2009.).

2.2. Ce este un Model Markov Ascuns

La începutul anului 1990, Krogh, la Universitatea California din Santa Cruz a descris rezultatele preliminarii ale modelării alinierilor multiple de secvențe de proteine. Modelările au fost făcute folosind Modele Markov Ascunse probabilistice. Modelele Markov Ascunse sunt încă privite ca și cutii negre în loc de modele naturale ale problemelor de aliniere de secvență. Multe dintre documentele cheie în care Modelele Markov Ascunse sunt descrise sunt în domeniul de recunoaștere a vorbirii. Aceste modele pot fi aplicate la o mulțime de probleme, cum ar fi : modelarea structurii proteinelor, găsirea de gene, analiza filogeneticii, modelarea seriilor de timp, recunoașterea vorbirii, modelarea regiunilor codificatoare și necodificatoare ale ADN-ului, subfamilii de proteine etc. (Kal Renganathan Sharma. 2009).

5

O secvență de ADN poate fi reprezentă folosind un Model Markov Ascuns. Un astfel de exemplu este prezentat în figura 2.1(Kal Renganathan Sharma. 2009) .

Figura 2.1 : Model Markov pentru o secvență de ADN cu o stare de început și o stare de sfârșit (Kal Renganathan Sharma. 2009)

Un Model Markov Ascuns (MMA) este similar cu un lanț Markov, dar este mai

general, și prin urmare mai flexibil, permițându-ne să modelăm fenomene pe care nu le puteam modela suficient de bine cu un lanț Markov normal. Un MMA este un Model Markov de timp discret, cu unele caracteristici suplimentare. Principalul plus este că atunci când o stare este vizitată de către lanțul Markov, ea "emite" o literă ce aparține unui alfabet fix și independent de timp. Literele sunt emise prin intermediul unei distribuții de probabilitate, independentă de timp, dar de obicei dependentă de stare ( Warren Ewens, Grant Gregory. 2005.).

Există, în primul rând, o secvență de stări vizitate, notate cu , și în al doilea rand, o secvență de simboluri emise, notate cu ,.. . Generarea lor poate fi vizualizată ca un proces cu două etape :

Notăm întreaga secvență de cu Q și întreaga secvență de cu O și scriem

"secvența observată ,. . . . " și" secvența de stări ......" .

De multe ori știm secvența O, dar nu știm secvența Q. Într-un astfel de caz, secvența Q se numește "ascunsă". O caracteristică importantă a Modelelor Markov Ascunse este că putem răspunde la câteva întrebări despre O și Q. Una dintre aceste întrebări se referă la estimarea secvenței stării ascunse cu cea mai mare probabilitate. Am ilustrat acest lucru printr-un exemplu simplu. Se consideră un lanț Markov cu două stări S1 și S2, cu o distribuție inițială uniformă și matricea de tranziții

(2.3)

6

Fie A un alfabet format doar din numerele 1 și 2. Starea emite 1 sau 2 cu o

probabilitate de ½, starea emite 1 cu o probabilitate de ¼ și 2 cu o probabilitate de ¾ . Presupunem că secvența observată este O=2 , 2 , 2 . Cât este ? (2.4)

Există opt posibilități pentru Q. Acestea pot fi scrise și probabilitatea lor calculată. Se

constată că răspunsul la întrebarea de mai sus este Q = . Secvența Q conține mai mulți , chiar dacă este mai probabil să producă un 2 atunci când e vizitat (probabilitate 3/4) decât (probabilitate 1/2). Motivul este că este mult mai

probabil să fie vizitat decât ( = 0.9 și = 0.8). De asemenea, se poate calcula Prob(O)= (2.5)

Acest calcul este folositor pentru a distinge care dintre modele este mai probabil să fi produs O.

În exemplul de mai sus, toate aceste calcule se pot face manual. Cu toate acestea, modelele care apar în practică au multe stari, uneori sute, și un alfabet cu multe simboluri. În aceste cazuri, calculele din exemplul anterior devin imposibile chiar și pentru cele mai rapide calculatoare. Din fericire, există abordări de programare dinamică ce rezolvă această problemă ( Warren Ewens, Grant Gregory. 2005.). Este necesar să se introducă unele notații specifice. Un MMA va consta în cinci componente:

(1) Un set de N stări

(2) Un alfabet de M simboluri de observare distincte A={

(3) Matricea probabilității de tranziție P=( ) , unde

(2.6)

(4) Probabilitățile de emisie: pentru fiecare stare și a A, (2.7)

Probabilitățile formează elementele unei matrici de dimensiune NxM notată

B=( ) (2.8)

(5) Un vector al distribuției inițiale , (2.9) Componentele 1 și 2 descriu structura modelului, iar componentele 3-5 descriu parametrii. reprezintă setul complet de parametrii. (2.10)

7

2.3. Arhitectura Modelelor Markov Ascunse

Modelele Markov Ascunse, sunt grafuri orientate, ale căror noduri reprezintă stări ascunse, iar arcele dintre noduri sunt etichetate cu probabilități de tranziție nenule (Sperschneider Volker. 2008.). Notăm probabilitatea de tranziție de la starea p la starea q cu , iar probabilitatea de emisie pentru caracterul x, în starea q, cu

ecuația (x).

Aspectele liniare ale unei secvențe pot fi surprinse de arhitecturi stânga-

dreapta. O arhitectură este stânga-dreapta, dacă se previne revenirea la o stare odată ce tranziția de la acea stare la orice altă stare a avut loc. În arhitectura standard a modelului Markov ascuns, în plus față de stările de început și de sfârșit, există și alte clase de stări: stări principale, stări de ștergere și stări de inserare (Kal Renganathan Sharma. 2009.).

Mai există de asemenea și arhitecturile MMA buclă și MMA circular (Kal

Renganathan Sharma. 2009.).

2.4.Cauze si efecte

Scopul unui Model Markov Ascuns este de a genera probabilistic șiruri de caractere emise, începând cu nodul q (0) și parcurgând apoi nodurile (stările ascunse) conform probabilităților de tranziție. Separarea între stările ascunse și caracterele observate este principala sursă de flexibilitate pe care modelele Markov Ascunse le ofertă unui utilizator. Opțiunea de modelare a interdependențelor locale dintre pozițiile adiacente, într-un șir emis, este a doua sursă de flexibilitate a modelului. Trebuie reținut faptul că opțiunea de a modela interdependențe pe distanțe lungi nu este disponibilă în Modelele Markov Ascunse. Există cinci distribuții de probabilitate standard, care joacă un rol esențial în utilizarea modelului Markov Ascuns (Sperschneider Volker. 2008). Numim secvențele de stări ascunse" cauze" și secvențele de caractere emise "efecte". Distribuțiile de probabilitate au scopul de a măsura cu ce probabilitate:

• o cauză apare împreună cu un efect (probabilitate comună)

• o cauză dată produce efecte (probabilitate condiționată de efecte, fiind dată o cauză)

• o cauză apare ( probabilitate marginală )

• un efect apare ( probabilitate marginală )

• o cauză apare, fiind dat un efect ( probabilitate condiționată de cauze, fiind dat un efect )

8

Notațiile și definițiile formale sunt: pentru orice secvență de stări Q = q (0) q(1)........ ...q (n), cu starea inițială q (0) și orice secvență de caractere S= x1.....xn (presupunem, pentru simplitate, că toate valorile care apar la numitori sunt mai mari decât zero):

P(Q)= – distribuția marginală a cauzelor (2.11)

P(S|Q)=

( ) – distribuția condiționată a efectelor, fiind dată cauza (2.12)

P(S,Q)=P(Q)P(S|Q) – probabilitatea comună (2.13) P(S)= - distribuția marginală a efectelor (2.14)

P(Q|S)=

– distribuția condiționată a cauzelor, fiind dat efectul (2.15)

Distribuția condiționată a efectelor, cu cauze date, descrie modelul. Distribuția marginală descrie apriori cunoștințele legate de producerea cauzelor. Distribuția condiționată a cauzelor, cu efecte date, descrie, după ce s-a făcut o observație, cunoștințele despre prezența unei anumite cauze. De exemplu, în observarea șirurilor de caractere, o secvență probabilă de stare ar putea defini o aliniere multiplă de șiruri de caractere (Sperschneider Volker. 2008). Probabilitățile condiționate sunt legate de regula lui Bayes:

P(Q|S)=

(2.16)

În teoria Modelelor Markov Ascunse se formulează următoarele ipoteze (Kal

Renganathan Sharma. 2009.):

1. Presupunerile Markov: Starea următoare depinde doar de starea curentă. Modelul rezultat este un model Markov de ordin 1. Se pot folosi și modele Markov ascunse cu ordin mai mare și complexitate mai ridicată.

2. Presupunerea staționării: Probabilitățile de tranziție ale stărilor sunt

independente de timpul real la care tranziția are loc.

3. Presupunerea independenței ieșirii: observatorul curent este independent de

cel anterior. Pentru un model Markov ascuns, descrie secvența O=

(2.17)

9

2.5. Problema evaluării și algoritmul Forward

Fiind dată o secvență O și Modelul Markov Ascuns, pentru λ, ce este utilizat pentru a reprezenta secvența O, găsiți P (O / λ). Să se calculeze numărul de operații

necesare. Se definește o variabilă auxiliară numită variabila forward . Variabila forward este definită ca fiind probabilitatea secvenței de observare parțială O, atunci când secvența se termină la starea i (Kal Renganathan Sharma. 2009.): (2.18) Se observă că următoarea relație recursivă este validă:

(2.19)

unde

și (2.20)

unde N este numărul stărilor MMA și T este lungimea secvenței

Timpul necesar pentru a completa task-ul este . 2.6. Problema decodării și algoritmul Viterbi

Problema constă în găsirea celei mai probabile secvențe de stări pentru o secvență de observații O și un model Markov ascuns λ. Soluția depinde de definiția " celei mai probabile secvențe de stări ". O abordare este aceea de a găsi starea cea

mai probabilă la t = t și concatenarea tuturor stărilor . Uneori, soluția acestei metode nu este semnificativă din punct de vedere fizic. O altă metodă este algoritmul Viterbi (Kal Renganathan Sharma. 2009.). Se găsește întreaga secvență de stări cu probabilitate maximă. Se definește o variabilă auxiliară: (2.21)

Această variabilă auxiliară reprezintă cea mai mare probabilitate pe care secvența de observare parțială și secvența de stări până la t o poate avea atunci când starea actuală este i. Se poate observa că următoarea relație recursivă va fi: (2.22)

Unde

Starea j * se găsește acolo unde este arg max [ (j)]. Pornind de la această stare, secvența de stări revine pe ruta indicată de pointer în fiecare stare . Acest lucru dă setul necesar de stări . Algoritmul Viterbi este ca un graf de căutare, ale cărui noduri sunt formate din stările modelului Markov ascuns în fiecare moment de timp t în intervalul (1, T).

10

2.7.Alinierea de secvențe folosind Modele Markov Ascunse

Similitudinile dintre două secvențe pot fi marcate utilizând modele probabilistice. Procesul de aliniere a spațiilor poate fi convertit în Model Markov Ascuns. Pentru a reprezenta alinierea de perechi cu penalități afine de interval se poate folosi un automat cu trei stări finite. Potrivirea este dată de starea M iar inserția este dată de stările X și Y (Figura 2.2) (Kal Renganathan Sharma. 2009.) .

MMA este derivat din diagrama masină din figura 2.3. Pentru emiterea unei

perechi aliniate a: b, starea M are distribuția de probabilitate .

Figura 2.2 Diagrama masină de stări finite pentru alinierea intervalelor afine (Kal Renganathan Sharma. 2009)

Figura 2.3 Model probabilistic pentru alinierea intervalelor afine (Kal Renganathan Sharma. 2009)

11

Stările X și Y vor avea distribuțiile pentru emiterea simbolului a. reprezintă

starea X. Starea X emite simbolul din secvența x. Parametrii modelului sunt prezentați în figura 2.3. Tranziția de la M la o stare de inserție este dată de δ, iar probabilitatea de a rămâne într-o stare de insertie este dată de ε. La figura 2.3 pot fi adăugate și stările de început și de sfârșit. Adăugarea unei stări de sfârșit poate

introduce modelului un alt parametru . Astfel, Modelul Markov Ascuns emite o aliniere de perechi. Discuțiile despre Modelul Markov Ascuns pot fi extinse la perechi de Modele Markov Ascunse. Ele au nevoie de un spațiu de căutare mai mare pentru a stoca secvența suplimentar emisă. O pereche de Modele Markov Ascunse poate fi folosită pentru a genera o pereche de secvențe aliniate (Kal Renganathan Sharma. 2009.) .

2.8. Găsirea de gene cu ajutorul Modelelor Markov Ascunse O utilizare importantă a Modelelor Markov Ascunse în analiza secvenței ADN este găsirea genelor. Tehnici, cum ar fi găsirea ORF-urilor nu sunt suficient de flexibile pentru a fi utilizate în identificarea genelor . Metodele similare sunt limitate de o serie de probleme, precum: incapacitatea de a detecta gene foarte scurte, care sunt imposibil de distins de zgomotul de fond, incapacitatea de a detecta gene netranslatate. Metodele de identificare de gene pot fi îmbunătățite prin integrarea unor semnale de secvență diferite. MMA face acest lucru în mod natural (Richard Axel, Linda Buck, Fred Hutchinson. 2004). Intronii și exonii

Genele eucariote pot fi foarte mari, dar proteinele pe care le codifică necesită

doar jumătate din numărul de nucleotide pentru a specifica toți aminoacizii. Motivul este faptul că aproape toate genele eucariote sunt împărțite în introni și exoni. În timp ce o porțiune foarte lungă de ADN trece printr-un proces de transcripție pentru fiecare genă, ARN-ul rezultat va fi editat înainte ca ARN mesager să fie trimis spre translație. Părțile genei care sunt eliminate se numesc introni iar celelalte porțiuni se numesc exoni (Richard Axel, Linda Buck, Fred Hutchinson. 2004).

Una dintre cele mai importante probleme în găsirea de gene eucariote, este

găsirea ORF-urilor. Deși există, de multe ori, situsuri de îmbinare pentru a marca poziția intronilor, acestea nu sunt secvențe rigid definite. Anumite semnale din introni pot ajuta în identificarea situsurilor de îmbinare și pot identifica, de asemenea, situsurile de ramificare necesare pentru îmbinare.

12

2.9. Modelele Markov Ascunse Circulare

Modelele Markov Ascunse Circulare sunt folosite pentru a modela periodicitatea în ADN. Modelele periodice din exoni și introni pot fi identificate prin utilizarea unor noi arhitecturi ale MMA, cum ar fi MMA buclă și MMA Circular. Modelele Markov Circulare au o capacitate mai bună de a descoperii tipare periodice în prezența zgomotului. Arhitectura stânga-dreapta nu este ideală pentru reprezentarea exonilor din cauza variației de lungime. A fost proiectat un model atât pentru secvențele de exoni cât și pentru secvențele de introni. Arhitectura MMA avea o formă circulară și un anumit număr de stări principale. Punctul de intrare al unor secvențe adăugate ulterior poate fi determinat utilizând programări dinamice. Periodicitatea poate fi descoperită prin modele cu număr diferit de stări (Kal Renganathan Sharma. 2009.).

Arhitectura modelului circular este afișată în figura 2.4. Săgețile din afară

reprezintă probabilitatea de pornire de la starea respectivă.

Figura 2.4 MMA circular cu 8 stări folosit pentru modelarea periodicității ADN-ului (Kal Renganathan Sharma. 2009)

Formând modele circulare de lungimi diferite, s-a constatat că modelele de lungime 10 au fost cele mai potrivite. Acest lucru este confirmat de faptul că probabilitățile de a sări peste stări în aceste modele sunt mici. Repetarea unei stări într-un model circular de 9 stări este neechivalentă cu săritul peste o stare într-un model circular de 10 stări. Aceste modele circulare nu conțin stări insert independente (precum arhitectura stânga-dreapta) (Kal Renganathan Sharma. 2009).

13

2.10. Modelul Markov Ascuns Generalizat (MMAG)

Un program folosit pentru găsirea genelor este Genie. Acest program este bazat pe Modelul Markov Ascuns Generalizat. Una dintre cele mai dificile probleme în găsirea genelor este determinarea structurii complete a genei. Pentru a rezolva această problemă se utilizează două rețele neurale noi bazate pe frecvențele dinucleotidelor. S-au obținut îmbunătățiri semnificative în sensibilitatea și specificitatea identificării structurii unei gene. Rezultatele experimentale folosind un set standard de gene adnotate arată că Genie a identificat 82% din nucleotidele codificatoare, cu o specificitate de 81% față de 74 și 81% în vechiul sistem (Kal Renganathan Sharma. 2009).

Genie este o implementare a MMAG. Stările Modelului Markov Ascuns

Generalizat sunt submodele arbitrare care emit secvențe de lungime variabilă. Un MMAG este definit în figura 2.5 având ca exemplu o sintaxă simplă a structurii unei gene.

Un MMAG descrie un model mai general, în care fiecare stare poate emite unul

sau mai multe simboluri în conformitate cu o distribuție arbitrară. Fiecare stare reprezintă un submodel independent, care poate fi el însuși un MMA. Un MMAG simplu care modelează structura genei eucariote este prezentat în figura 2.6. Arcele reprezintă stări ce emit șiruri de baze iar nodurile reprezintă tranziții între stări (Kal Renganathan Sharma. 2009).

MMAG este reprezentat ca un graf. Nodurile din graf reprezintă tranziții între

stări. Fiecare stare corespunde unui submodel al caracteristicii unei gene abstracte cum ar fi un exon (E) sau un intron (I). Pentru orice secvență de baze x și stare q, submodelul asociat cu starea q definește o probabilitate pentru secvența x. Această probabilitate este notată cu P (x / q). Când MMAG este privit ca un model statistic generativ, aceasta este probabilitatea secvenței emise când procesul Markov este în starea q. Aceste funcții probabilistice, una pentru fiecare stare, fac parte din definiția MMAG (Kal Renganathan Sharma. 2009).

.

Figura 2.5 Un MMAG simplu pentru o secvență cu multiple gene exon (Kal Renganathan Sharma. 2009)

14

Figura 2.6 Un MMAG ce include constrângeri (Kal Renganathan Sharma. 2009)

MMAG din figura 2.6 reprezintă doar ordonarea de bază a caracteristicilor genei.

2.11.Familii de proteine

Modelele Markov Ascunse au fost aplicate cu succes la multe familii de proteine, cum ar fi globine, immunoglobuline, kinaze și receptori cuplați la proteina G (GPCRs). GPCRs sunt o familie de proteine transmembranare capabile să transporte o varietate de semnale extracelulare purtate de hormoni, neurotransmițători, de miros și de lumină. Modele Markov Ascunse au fost folosite pentru a modela structurile secundare de proteine. FORESST, baza de date ce conține structurile secundare ale familiilor de proteine și Pfam, baza de date ce conține familii de proteine, au fost disponibile în 1997. În 1997, Pfam conținea 527 de familii verificate manual, constând în 39113 alinieri de secvențe și 6.8 milioane de reziduuri. Bazele de date de familii de proteine de obicei se bazează pe mai multe alinieri de secvențe ale membrilor cunoscuți ai familiei (Kal Renganathan Sharma. 2009).

Alinierea de referință este o aliniere multiplă de seturi reprezentative de secvențe, verificată manual.

15

2.12. Pfam

Baza de date Pfam clasifică familii de proteine în funcție de profile. Un profil este un model ce evaluează probabilitatea de apariție a unui aminoacid dat, a unei inserții sau a unei ștergeri pe fiecare poziție într-o secvență de proteine (Selzer Paul, Marhöfer Richard, Rohwer Andreas. 2008).

Pozițiile conservate sunt ponderate mai puternic decât pozițiile mai puțin

conservate. Pfam se bazează pe alinieri de secvențe. Alinierile verificate manual, având o calitate foarte bună, servesc ca puncte de plecare pentru construirea Modelelor Markov ascunse (MMA). Secvențele ulterioare sunt adăugate automat la alinierile individuale ale bazei de date Swissprot.

Alinierile rezultate ar trebui să reprezinte structuri funcționale. Datorită

construcției parțial automată a alinierilor este posibil ca alinierile de secvență să nu aibă nici o relație evolutivă una cu cealaltă. Prin urmare, rezultatele unei căutări în baza de date Pfam ar trebui să fie analizate cu atenție.

16

3. Studiu de caz. Aplicații ale modelelor Markov ascunse în studiul receptorilor de miros

În 2004, Premiul Nobel pentru Fiziologie sau Medicină, a fost câștigat de Richard Axel de la Universitatea Columbia și de Linda Buck de la Fred Hutchinson Cancer Research Center pentru elucidarea sistemului olfactiv. Sistemul olfactiv include o mare familie de proteine numită receptori de miros. Aceștia, în combinație, fac posibilă recunoașterea a peste 10 000 de mirosuri diferite. Acești receptori de miros sunt atașati la suprafața celulelor în pasajul nostru nazal, detectând molecule de miros când acestea sunt inhalate (Richard Axel, Linda Buck, Fred Hutchinson. 2004) .

Receptorii de miros sunt cea mai mare familie de gene din genomul uman, cu

aproximativ 1000 de membri. Cu toate acestea, numai 40% din gene sunt funcționale. Peste 600 de receptori de miros umani sunt pseudogene sau descendenți inactivi ai genelor funcționale. Când genele nu mai sunt necesare pentru miros, de exemplu atunci când culoarea devine simțul principal, selecția naturală nu mai menține funcția (Richard Axel, Linda Buck, Fred Hutchinson. 2004).

Pentru ca receptorii de miros (OR) să simtă moleculele din afara celulei și să

semnalizeze interiorul celulei, trebuie să traverseze membrana celulelor. Pentru a face acest lucru, receptorii de miros conțin șapte domenii transmembranare.

Pseudogenele conțin codoni de stop sau indeluri care rezultă din producția

proteinelor nefunctionale. Analiza receptorilor de miros necesită instrumente mai avansate. Similitudinile

între diferite gene OR sunt prea mici pentru a fi depistate printr-o aliniere simplă de perechi, în schimb dacă le luăm pe toate într-o aliniere multiplă se va genera un semnal foarte puternic.

3.1. Receptorii de miros

Mirosul este o funcție neurosenzorială foarte importantă. În interiorul nasului se alfă milioane de receptori olfactivi ce detectează mirosurile pe care le vor trimite la creier. La rândul său, creierul va spune dacă mirosul este bun sau rău. Aceasta este o modalitate simplă de informare a unui animal de ceea ce ar putea fi dăunător sau benefic pentru el. Acest mecanism a evoluat de-a lungul timpului pentru a oferi șanse de supraviețuire speciilor vii. În ceea ce privește oamenii, capacitățile lor olfactive sunt mult mai puțin dezvoltate decât la câine sau la un animal sălbatic, dar ele rămân totuși esențiale. Ani de zile, cercetătorii s-au întrebat dacă percepția reală a unui miros are o bază genetică. S-a scris mult pe această temă și răspunsul pare a fi "da". Recent, oamenii de știință au descoperit că percepția mirosului pentru iarba

17

tăiată are o bază genetică și este sub influența unui receptor olfactiv numit "receptorul de miros 2J3 " sau OR2J3 (http://web.expasy.org/spotlight/).

Primul pas în indentificarea mirosurilor este interacțiunea unei molecule cu

receptorii de miros aflați pe suprafața neuronilor olfactivi chemosenzoriali. Receptorii de miros ai mamiferelor aparțin familiei de receptori cuplați la proteina G (GPCR). Această familie include și receptori găsiți în retină pentru a simți lumina precum și receptori pentru hormoni și neurotransmițători, cum ar fi serotonina și dopamina. Fiecare neuron olfactiv exprimă un singur tip de receptor de miros ceea ce implică un mecanism sofisticat de alegere a genei receptorului de miros. O altă caracteristică interesantă a simțului olfactiv este recunoașterea combinatorială de mirosuri. Fiecare receptor recunoaște multiple mirosuri și fiecare miros se leagă la mai mulți receptori (Zozulya Sergey, Echeverri Fernando, Nguyen Trieu. 2001. ).

Genele ce codifică receptorii de miros nu conțin introni în regiunile

codificatoare. La mamifere genele receptorilor de miros sunt organizate în grupe de 10 sau mai mulți membrii și sunt situate pe mai mulți cromozomi. La om repertoriul receptorilor de miros conține un număr mare de pseudogene. Studii recente arată că 70% din genele receptorilor de miros pot fi pseudogene.

Identificarea, clonarea, clasificarea bazată pe secvență și analiza receptorilor

de miros umani sunt cerințe esentiale în studierea acestei familii de receptori. Trebuie identificate toate genele receptorilor de miros ce codifică receptori de lungime completă (Zozulya Sergey, Echeverri Fernando, Nguyen Trieu. 2001. ).

Strategia pentru a găsi aceste gene se bazează pe absența intronilor din

secvențele codificatoare ale receptorilor de miros ai mamiferelor, pe similaritatea secvențelor precum și pe prezența secvențelor puternic conservate la toți receptorii de miros cunoscuți ai mamiferelor.

Primul pas este identificarea tututor secvențelor receptorilor de miros la oameni

prin cuvinte cheie extinse și căutări omologice în baze de date de proteine și de ADN. Secvențele rezultate se compară între ele cu BLAST și prin alinieri multiple de secvențe. Se elimină apoi pseudogenele ce conțin ștergeri și alte defecte incompatibile cu funcția receptorului (Zozulya Sergey, Echeverri Fernando, Nguyen Trieu. 2001. ).

Următorul pas este identificarea de membri suplimentari ai familiei. Secvențele

genomice ce conțin zone cu omologie ridicată sunt supuse căutărilor de ORF-uri. Odată identificate ORF-urile, ele sunt translatate și comparate cu receptorii de miros cunoscuți. Un criteriu de recunoaștere a genei receptorului de miros ce codifică un receptor funcțional și de lungime completă este prezența unui ORF neîntrerupt ce începe cu un codon ATG și o unitate șapte- transmembranara completă.

Pentru a valida secvențele obținute prin extragerea de date din genom, toate

secvențele codificatoare ale genelor receptorilor de miros se clonează prin folosirea mai multor secvențe de ADN. De exemplu dacă se folosesc secvențe ADN de la opt persoane se obține o medie de 8 clone independente pentru fiecare receptor. Aceste secvențe sunt izolate și secvențiate în întregime. Secvențierea datelor confirmă identificarea corectă a ORF-urilor codificatoare ale receptorilor de miros de lungime completă. Recent au fost publicate studii despre faptul că secvențele genelor

18

receptorilor de miros și secvențele pseudogenelor acoperă aproape 150 de gene ale receptorilor de miros de lungime completă.

Câinii și șobolanii au o capacitate extrem de dezvoltată în a detecta și identifica moleculele de miros. Simțul mirosului a fost dezvolat pe parcursul evoluției lor pentru a comunica cu lumea exterioara, făcând posibilă identificarea prăzii și pentru a putea evita pericolele. Genele receptorilor de miros au fost identificate pentru prima dată de Axel și Buck . Extragerea datelor din secvența genomului a dus la indentificarea și caracterizarea a 650-900 de gene la om și 1200-1500 de gene la șoareci. Repertoriul olfactiv la șobolani conține între 1700-2000 de gene, în timp ce câinele are aproximativ 1300 de gene (Quignon, Pascale, Mathieu Giraud, Maud Rimbault, Patricia Lavigne, Sandrine Tacher, Emmanuelle Morin, Elodie Retout . 2005.).

Translația secvențelor de gene la câini și șobolani a facut posibilă identificarea

pseudogenelor și identificarea numărului de mutații de închidere a ORF-ului. S-a constatat că 20.3% din genele receptorilor de miros la câine și 19.5 % la șobolan sunt pseudogene. Un cadru simplu de mutație de închidere a detectat la 78 din cele 222 de pseudogene ale câinelui. S-au găsit 43 de pseudogene ce conțineau două mutații și 101 ce conțineau 3 mutații. La șobolan, 153 de pseudogene conțineau o mutație, 48 conțineau 2 mutații și 91 aveau 3 sau mai multe mutații. Pseudogenele cu mai mult de o mutație de închidere a ORF-ului sunt cu siguranță pseudogene reale. În schimb, nu toate pseudogenele ce conțin o singură mutație sunt pseudogene reale .

Pentru a determina gradul de variație a poziției fiecărui aminoacid se aliniază

toate secvențele aminoacizilor pentru receptorii de miros la câini și șobolani.

Tabelul 3.1 prezintă distribuția genelor corespunzătoare repertoriilor de miros la cele 4 mamifere la care ele sunt complet cunoscute:

Numărul de clase

Numărul de familii

Numărul de Subfamilii

Om 2 17 300

Șoarece 2 Nd 241

Câine 2 23 300

Șobolan 2 21 282

Tabel 3.1 distribuiția genelor în familii și subfamilii

Repetoriul câinelui este mai mare decât cel al șobolanului și conține mai multe

poziții puternic conservate (23 de poziții în proteinele receptorilor de miros ale câinelui și 21 de poziții la șobolan ). Nivelul de conservare și numărul mare de subfamilii identificate la câine indică faptul că acesta are un repertoriu mai divers decât șobolanul. Domeniile transmembranare IV și V conțin cea mai mare proporție de aminoacizi extrem de variabili (Quignon, Pascale, Mathieu Giraud, Maud Rimbault, Patricia Lavigne, Sandrine Tacher, Emmanuelle Morin, Elodie Retout . 2005.).

Repertoriul uman este cel mai mic din cele patru repertorii ale mamiferelor și

are cel mai mai mic număr de familii, 17. La fel ca și repertoriul câinelui, el poate fi împărțit în 300 de subfamilii. Repertoriul șobolanului conține doar 282 de subfamilii

19

deși este cel mai mare din cele patru repertorii. Numărul mare de subfamilii, la om, reflectă mai degrabă numărul mare de pseudogene decât diversificarea repertoriului. Aproximativ 126 de subfamilii sunt formate în întregime din pseudogene. În schimb numărul mare de subfamilii din repertoriul câinelui indică un nivel de diversificare mare.

Există pseudogene atât în clasele I și II, cât și în familii și subfamilii, dar nu

sunt distribuite egal. Clasa I include mai puține pseudogene (17% pentru câine și 13 % pentru șobolan ) decât clasa II ( 23 % pentru câine și 20 % pentru șobolan). O variabilitate mai mare există la familii și subfamilii. De exemplu în familia 6 ce aparține clasei 2, 34 % la câini și 10% la șobolani sunt pseudogene iar familia 10 are 13% pentru câini și 20% pentru șobolani (Quignon, Pascale, Mathieu Giraud, Maud Rimbault, Patricia Lavigne, Sandrine Tacher, Emmanuelle Morin, Elodie Retout . 2005.).

Doar 22 de subfamilii ale câinelui și 11 subfamilii ale șobolanului au fost găsite

pe mai mult de un cromozom. Din modul în care genele șobolanului sunt numite reiese că ordinul genelor din genom tinde să respecte ordinea filogenetică.

Vom folosi un MMA de 2 stări pentru segmentarea proteinei receptorilor de

miros, ca o modalitate de a descoperi regiunile hidrofobice și hidrofilice caracteristice proteinei 7-TM. Apoi vom folosi un Model Markov Ascuns de profil ce reprezintă această familie importantă de proteine , obținută din baza de date on-line Pfam pentru a decide dacă o anumită secvență a receptorului de miros aparține familiei sau nu. Ca o comparație, vom încerca să luăm aceeași decizie bazându-ne pe alinierea globală între aceeași secvență și un receptor tipic (număr GenBank NP_002368) , ce a fost folosit pentru a crea modelul Markov Ascuns de profil . Vom ilustra puterea de aliniere a MMA de profil . Vom vedea , de asemenea , modul în care MMA de profil poate fi utilizat pentru a obține o aliniere multiplă . (Richard Axel, Linda Buck, Fred Hutchinson. 2004) .

3.2. Modele Markov ascunse pentru segmentarea receptorilor de miros

Secvențele de gene și proteine pot conține regiuni distincte ale căror proprietăți chimice diferă mult. MMA ne poate ajuta în definirea limitelor exacte ale acestor regiuni. Segmentarea este folosită pentru a defini întinderi mult mai mari de utilizare a nucleotidei eterogene în secvențe genomice . Ea mai poate fi utilizată pentru a identifica caracteristicile biologice responsabile de această eterogenitate.

Receptorii de miros sunt proteine 7-transmembranare , ceea ce înseamnă că

fiecare trebuie să traverseze membrana celulei de șapte ori. Proteina trebuie să aibă șapte segmente hidrofobice (care nu reacționează cu membrana celulelor grase) și șapte segmente hidrofilice (care nu reacționează cu citoplasma apoasă și medii extra-celulare). Fiecare aminoacid are propiul său nivel caracteristic de hidrofobicitate, de la extrem de hidrofob la extrem de hidrofil. Dar pentru că nu toți aminoacizii într-o regiune hidrofobă a proteinei vor fi extrem de hidrofobi (și invers pentru regiunile hidrofile), trebuie să utilizăm Modele Markov Ascunse pentru a

20

segmenta receptorii de miros în segmente hidrofobe și segmente hidrofile (Richard Axel, Linda Buck, Fred Hutchinson. 2004).

Segmentarea receptorilor olfactivi de la câine

Vom construi un MMA pentru a prezice regiunile hidrofobice ale unei proteine codificate de una dintre genele corespunzătoare portofoliului de receptori olfactivi care conține: 1121 de gene, dintre care 222 pseudogene. Citim în Matlab secvențele de aminoacizi corespunzătoare. Folosim funcția fastaread care citește un format fasta și returnează datele din fișier ca o structură.

orseqs = fastaread('prot-dog.fasta'); orseqs

orseqs =

848x1 struct array with fields:

Header

Sequence

Iată secvența de AA corespunzătoare primei intrări din repertoriu:

>>or1= orseqs(1).Sequence

ans =

MELMDGNYTLVTEFILLGFPTRPELQIVLFLVFLTLYGIILTGNIGLMMLIRT....

și lungimea ei:

>> length(or1)

ans = 314

Vom folosi această primă secvență și funcția Matlab proteinplot pentru a afișa diverse proprietăți corespunzătoare secvenței.

proteinplot(or1)

21

Figura 3.1 Diverse proprietăți ilustrate cu funcția proteinplot

22

Din lista de proprietăți selectăm hydrophobicity (Kyte & Doolittle) și obținem:

Figura 3.2 Figura anterioară după selectarea proprietății de hidrofobicitate (Kyle & Doolittle)

Exportăm figura, o salvăm cu numele 'hydro.fig' și apoi o deschidem.

>> open('hydro.fig')

>> hold on

23

Figura 3.3 Graficul corespunzător hidrofobicității primei secvențe din fișierul fasta

Considerăm această secvență de aminoacizi ca fiind generată de un proces stochastic care alternează între două stări: 1. în afara membranei 2. în membrană. Vom putea folosi astfel un MMA pe care îl antrenăm cu primele 30 de secvențe din repertoriu: >> for i=1:20

intseqs(i) = {aa2int(orseqs(i).Sequence)};

end

-aa2int convertește un șir de aminoacizi într-un vector de numere întregi pe baza unui tabel de mapare.

>> intseqs{1}

ans =

Columns 1 through 15

13 7 11 13 4 8 3 19 17 11 20 17 7 14 10

24

Columns 15 through 30

11 11 8 14 15 17 2 15 7 11 6 10 20 11 14 ...

Matricea de tranziție a MMA va avea dimensiunea 2*2 (avem 2 stări ascunse). O generăm astfel:

>> T = [0.95 0.05;

0.05 0.95]

T = 0.9500 0.0500

0.0500 0.9500

Matricea de emisie corespunzătoare MMA va avea dimensiunea 2*20 (sunt 20 de AA). Iată o estimare initială a ei:

>> E = [0.018 0.067 0.067 0.067 0.018 0.067 0.067 0.067 0.067 0.01 0.01 0.067 0.018 0.018 0.067 0.067 0.067 0.067 0.067 0.01;

0.114 0.007 0.007 0.007 0.114 0.007 0.007 0.025 0.007 0.114 0.114 0.007 0.114 0.114 0.025 0.025 0.025 0.025 0.025 0.114];

Acum vom determina matricile exacte folosind funcția hmmtrain care estimează probabilitățile de tranziție și emisie pentru un model Markov ascuns. Rezultatul poate fi un vector ce conține o singură secvență sau o matrice ce conține câte o secvență pe fiecare rând.

[estT , estE] =hmmtrain(intseqs , T ,E)

estT =

0.7257 0.2743

0.1420 0.8580

estE =

Columns 1 through 12

0.0470 0.0834 0.0686 0.0492 0.0060 0.0516 0.0402 0.0258 0.0152 0.0507 0.1086 0.0956

0.0643 0.0042 0.0159 0.0143 0.0475 0.0126 0.0131 0.0521 0.0296 0.0975 0.1618 0.0000

25

Columns 13 through 20

0.0858 0.0000 0.0816 0.0586 0.0424 0.0110 0.0318 0.0470

0.0205 0.1148 0.0171 0.1085 0.0748 0.0031 0.0494 0.0989

Apoi folosim algoritmul Viterbi implementat în Matlab pentru a determina traiectoria stărilor și vom afișa stările estimate peste figura generată anterior (hidro.fig).

>> estimatedStates = hmmviterbi(aa2int(or1),estT,estE);

plot(estimatedStates)

hold off

Iată și graficul obținut:

Figura 3.4 Hidrofobicitatea după netezire și segmentare

26

Segmentarea receptorilor olfactivi de la șobolan

Vom construi un MMA pentru a prezice regiunile hidrofobice ale unei proteine codificate de una dintre genele corespunzătoare portofoliului de receptori olfactivi care conține: 1493 de gene, dintre care 293 pseudogene. Folosim funcția fastaread pentru a citi în Matlab secvențele corespunzătoare.

>>orseqs = fastaread('prot-rat.fasta'); >>orseqs orseqs = 1201x1 struct array with fields: Header Sequence >>or1=orseqs(2).Sequence or1 = MEAPVPIPTANRTSVLEFLLLGVTEKRGLQLLLFGILLITYTLTLLGNLFI.... >>length(or1) ans = 324 proteinplot(or1)

Figura 3.5 Diverse proprietăți ilustrate cu funcția proteinplot

27

Selectăm hydrophobicity (Kyle& Doolittle) și obținem:

Figura 3.6 Figura anterioară după selectarea proprietății de hidrofobicitate (Kyle & Doolittle)

Se salvează figura cu numele ’ hydro2.fig’.

>> open('hydro2.fig')

>> hold on

28

Figura 3.7 Hidrofobicitatea primei secvențe a receptorilor de miros la șobolan

Vom putea folosi astfel, ca și la câine, un MMA pe care îl antrenăm cu primele 30 de secvențe din repertoriu:

>> for i=1:30

intseqs(i) = {aa2int(orseqs(i).Sequence)};

end

Matricea de tranziție a MMA va avea dimensiunea 2*2 (avem 2 stari ascunse). O generăm astfel:

>> T = [0.95 0.05;

0.05 0.95]

T= 0.9500 0.0500

0.0500 0.9500

29

O estimarea a matricii de emisie corespunzătoare MMA :

>> E = [0.018 0.067 0.067 0.067 0.018 0.067 0.067 0.067 0.067 0.01 0.01 0.067 0.018 0.018 0.067 0.067 0.067 0.067 0.067 0.01;

0.114 0.007 0.007 0.007 0.114 0.007 0.007 0.025 0.007 0.114 0.114 0.007 0.114 0.114 0.025 0.025 0.025 0.025 0.025 0.114];

[estT , estE] =hmmtrain(intseqs , T ,E)

estT =

0.9999 0.0001

0.0000 1.0000

estE =

Columns 1 through 12

0.0613 0.0367 0.0341 0.0227 0.0402 0.0227 0.0208 0.0392 0.0205 0.0828 0.1440 0.0301

0.1628 0.0000 0.0051 0.0000 0.3417 0.0000 0.0000 0.2483 0.0000 0.0000 0.0000 0.0000

Columns 13 through 20

0.0328 0.0696 0.0432 0.0907 0.0697 0.0057 0.0433 0.0900

0.0000 0.0000 0.0000 0.0000 0.2421 0.0000 0.0000 0.0000

Folosim algoritmul Viterbi implementat în Matlab pentru a determina traiectoria stărilor și vom afișa stările estimate peste figura generată anterior (hidro2.fig).

>> estimatedStates = hmmviterbi(aa2int(or1),estT,estE);

plot(estimatedStates)

hold off

30

Graficul obținut:

Figura 3.8 Hidrofobicitatea după netezire și segmentare

31

Segmentarea receptorilor olfactivi de la om

Vom construi un MMA pentru a prezice regiunile hidrofobice ale unei proteine codificate de una dintre genele corespunzătoare portofoliului de receptori olfactivi care conține 347 de gene. Se parcurg aceeași pași ca la câine și șobolan.

Prima secvență din fișierul fasta pentru receptorii de miros la om este :

Or1= MPNSTTVMEFLLMRFSDVWTLQILHSASFFMLYLVTLMGNILIVTVTTCD...

Folosim această secvență pentru a afișa diverse proprietăți ale receptorilor de miros

Figura 3.9 Hidrofobicitatea

Putem netezi graficul utilizând o fereastră alunecatoare de 20 de aminoacizi și de asemenea prin raportarea la nivelul mediu al hidrofobicității fiecărei regiuni.

Această netezire ar trebui să facă mai ușoară vizualizarea segmentelor hidrofobice.

Figura 3.10 Hidrofobicitatea după netezire și segmentare

32

3.3. Modele Markov Ascunse de profil

O problemă importantă este caracterizarea seturilor de proteine omoloage (familii de gene) în funcție de modelele comune din secvența lor . Acest lucru ne permite, de exemplu, să determinăm dacă o nouă proteină aparține unei familii anume sau nu. O abordare clasică pentru a caracteriza familiile este de a porni de la o aliniere multiplă a tuturor elementelor din familie. Rezumatul alinierii poate fi caracterizat prin frecvența cu care unele simboluri apar în anumite poziții. Uneori, aceste frecvențe pot fi reprezentate printr-o matrice de punctare a poziției, precizând probabilitatea de a observa un simbol specific la o anumită poziție. Această metodă, totuși, nu funcționează bine pentru cazurile care includ spații de lungime variabilă. Modele Markov Ascunse pot fi utilizate în această situație pentru a oferi o caracterizare mai flexibilă a modelelor de secvență . Acest lucru se face prin definirea unui MMA de profil (MMA de profil). Un mod de a gândi MMA de profil este ca fiind o descriere abstractă a unei familii de proteine, un alt mod este ca fiind situații statistice ale alinierii multiple de secvențe (Richard Axel, Linda Buck, Fred Hutchinson. 2004).

Mai formal, MMA de profil codifică informații despre frecvența anumitor

aminoacizi, precum și frecvența inserărilor și ștergerilor din aliniere. Acestea sunt realizate din mai multe alinieri de secvențe omoloage. Ele conțin stări de potrivire, care descriu distribuția de aminoacizi pe fiecare poziție, precum și stări de inserare și ștergere care permit adăugarea sau îndepărtarea reziduurilor. Există o stare de potrivire, o stare de inserție și o stare de ștergere pentru fiecare coloană a alinierii multiple. Cu alte cuvinte, pentru o aliniere de 250 de aminoacizi, în MMA de profil există 250 de potriviri, 250 de inserții și 250 stări de ștergere. Pentru fiecare stare de potrivire și inserție există o probabilitate de emisie. Nici un aminoacid nu este emis de stările de ștergere.

Fiecare poziție a profilului are diferite probabilități de emisie asociate cu fiecare

simbol și fiecare margine a grafului are propriile probabilități de tranziție. Orice cale de la starea de început la starea de final poate corespunde cu o instanță valabilă a modelului.

Alinierea unei secvențe cu MMA de profil este echivalentă cu alinierea acesteia la sutele de secvențe utilizate pentru a produce modelul.

33

MMA de profil la câine

Din alinierea multiplă a membrilor unei familii de proteine se poate obține MMA de profil. Vom folosi pentru aceasta baza de date PFAM.

>> web('http://www.sanger.ac.uk/Software/Pfam')

Figura 3.11 Baza de date Pfam

Vom genera o versiune aleatoare de AA pe baza primei secvențe din repertoriu:

>> randor1 = randseq(length(or1), 'fromstructure', aacount(or1))

unde randseq generează aleator o secvență iar aacount numără aminoacizii dintr-o secvență.

randor1 =

SSSNPLENTCDPYKSDVKPATSNKIKLLCGAYPVDYTMLGSERSPSFLLTIT.....

Putem căuta în întreaga bază de date PFAM folosind secvența de AA corespunzătoare primei secvențe din repertoriul nostru (vezi mai jos):

34

Figura 3.12 Secțiunea din Pfam destinată căutării familiei corespunzătoare unei anumite secvențe

or1 =

MELMDGNYTLVTEFILLGFPTRPELQIVLFL...

Vom copia secvența de AA și o vom introduce în căsuța de căutare :

Figura 3.13 Aceeași secțiune de mai sus în care a fost introdusă prima secvență de aminoacizi .

Iată rezultatele obținute :

(http://pfam.sanger.ac.uk/search/sequence/results/8A467D64-E222-11E2-9A03-

BDE4780A6372), cea mai buna potrivire cu familia 13853:

35

Figura 3.14 rezultatele returnate în urma căutării în baza de date

Figura 3.15 Familia secvenței ce a fost introdusă în secțiunea de căutare

Alternativ, putem realiza și în Matlab această căutare pentru primele 60 de familii din Pfam. Se folosește funcția hmmprofalign care întoarce scorul alinierii optime a secvenței de aminoacizi cu MMA de profil. De asemenea, se folosește și funcția gethmmprof pentru a returna modelul.

>> seqs = {or1,randor1};

for i = 1:60

for j = 1:2

[score(i,j)] = hmmprofalign(gethmmprof(i), seqs(j));

end

end

score

score = 111.6303 -53.6475 -52.1349 -53.7212 -36.9297 -44.4756

36

-23.2696 -28.5125 -38.3675 -31.0312 -49.1357 -52.5800 -14.5213 -3.3388 -8.3739 -14.2729 -41.3303 -28.8950 -11.5414 -9.7119 -20.1695 -8.8635 -413.4147 -411.4147 -17.3759 -12.2411 -6.6870 -12.0361 -57.0009 -39.3764 -87.1938 -77.4234 -15.4252 -16.2869 -21.6247 -15.1793 -14.1947 -15.6088 -14.1431 -15.5939 -2.9992 -2.5814 -154.2015 -163.1153 -4.0557 -0.3789 -17.1447 -14.3371 -38.7431 -33.6060 -100.0165 -99.3197 -8.7972 -19.0439 -14.3674 -11.5374 -18.2473 -9.2198 -9.4353 -17.6704 -15.0932 -23.1484 -15.3614 -13.3109 -27.1863 -36.9958 -17.0496 -13.2726 -9.5009 -16.1573 -10.5556 -6.6097 -13.6748 -11.1062 -101.0705 -89.9091 -8.9421 -12.4306 -7.4489 -7.7473 -16.1694 -15.8224 -20.2487 -14.1285 -3.4784 -17.1826 -33.7205 -28.8425 -11.0987 -9.4504 -18.0708 -12.0879 -11.0073 -9.0649 -13.8493 -15.7857 -21.1720 -16.1620 -12.9487 -3.5971 -12.1487 -10.4700 -24.4452 -24.6346 -6.9352 -3.0578 -28.8861 -23.6426 -68.8823 -54.8723

37

-25.5320 -26.4901 -11.8119 -9.5748 -7.4897 -14.3176 -5.9573 -16.5193 -24.4975 -21.7178 Se alege potrivirea semnificativă. Vom folosi în continuare familia 13853 detectată mai sus.

>> hmm7tm = gethmmprof(13853);

Funcția gethmmprof returnează modelul corespunzător numărului de accesare.

Regăsim alinierile multiple corespunzătoare:

>> seqs = gethmmalignment(13853, 'type', 'seed');

>> seqs

seqs =

118x1 struct array with fields:

Header

Sequence

Gethmmalignment întoarce doar alinierile folosite pentru a genera MMA când tipul e ‚seed’ și întoarce toate alinierile când tipul e ‚full’. Valoarea default e ’ full ’.

Afișăm aceste rezultate (parțial):

>> disp([char(seqs.Header) char(seqs.Sequence)])

Figura 3.16 O parte din alinierile folosite pentru a genera MMA

38

Vom alinia profilul corespunzător cu primele 50 de secvențe din repertoriul nostru:

>> for i=1:50

[Score(i), Seqs(i).Aligned] = hmmprofalign(hmm7tm, orseqs(i).Sequence);

end

hmmprofmerge(Seqs,Score)

Funcția hmmprofmerge(seqs,score) - sortează secvențele folosind scorurile

Figura 3.17 Secvențele aliniate sortate cu ajutorul scorurilor

39

Trebuie evidențiat de ce e mai bine să folosim MMA în loc de alinierea de perechi . De fiecare dată cand se aliniază o secvență cu un MMA, este ca și cum secvența ar fi aliniată la sutele de secvențe ce au fost folosite pentru a crea MMA. Pentru a confirma acest lucru se consideră următoarele: familia la care aliniem receptorul de miros este familia rhodopsinei. Se poate încerca alinierea cu una din secvențele folosite pentru dezvoltarea MMA. Pentru început trebuie returnată o secvență pentru rhodopsină și una pentru receptorii de miros.

>>rhod = getgenpept('NP_002368','sequenceonly',true) >>or = orseqs(300).Sequence rhod =

MDGSNVTSFVVEEPTNISTGRNASVGNAHRQIPIVHWVIMSISPVGF .......

or =

MDAANDSVVSEFVLIGLSNSWDMHLFLFLFFSVFYVG...

Acum se poate face alinierea globală de perechi a celor două secvențe. BLOSSUM 30 este matricea de punctare. Penalizările pentru deschiderea și extinderea unui spațiu în aliniere sunt setate la 5. >>[Score, Alignment] = nwalign(or, rhod, 'scoringmatrix', 'blosum30', 'gapopen', 5, 'extendgap',5) >>perm = randperm(length(or)); >>randor = or(perm); >>[Score, Alignment] = nwalign(randor, rhod, 'scoringmatrix', 'blosum30', 'gapopen', 5, 'extendgap',5) Funcția nwalign face alinierea globală Needleman-Wunsch pentru 2 secvențe [Score, Alignment] = nwalign(randor, rhod,Gapopen ,penalty)- definește penalizarea pentru deschiderea unui spatiu în aliniere. [Score, Alignment] = nwalign (randor,rhod,'scoringmatrix','blosum30') – definește matricea de punctare ce va fi folosită pentru aliniere [Score, Alignment] = nwalign(randor, rhod,’extendgap’ ,penalty-definește penalizarea pentru extinderea unui spațiu în aliniere. Dacă extendgap nu e specificat, extinderile spațiilor sunt punctate cu aceeași valoare ca gapopen. randor = LANHIVIPGINFVGSRFYMLIAFIHKSSAEFIVTDHIFL

40

Figura 3.18 Rezultatele (parțiale) alinierii rhodopsinei cu ultima secvență a receptorilor de miros folosind algoritmul Needleman și Wunsch

Figura 3.19 Rezultatele (parțiale) alinierii rhodopsinei cu o secvență aleatoare a receptorilor de miros folosind algoritmul Needleman și Wunsch

41

Daca se aliniază ambele secvențe cu MMA, semnificația alinierii este mai evidentă.

[score_or, align_or] = hmmprofalign(hmm7tm,or) [score_rand, align_rand] = hmmprofalign(hmm7tm,randor)

Figura 3.20 rezultatele (parțiale) alinierii ambelor secvențe cu MMA

.

42

MMA de profil la șobolan

Generăm o versiune aleatoare de AA pe baza primei secvențe din repertoriu:

>> randor1 = randseq(length(or1), 'fromstructure', aacount(or1))

randor1 = FLLIHALQLPRTVTLFLLAVALLLVYASLLTWVCWGLTILAFPGIPVPA.......

Or2 = MPPSWAAPPWANQSRTHELEFVLLGFAHVPSLRPMLAALFLAA

La fel ca și la câine căutăm familia de proteine în Pfam cu ajutorul secvenței or2

Figura 3.21 Secțiunea din Pfam destinată căutării în care a fost introdusă secvența or2

Rezultatele obținute sunt :

Figura 3.22 Rezultatele returnate în urma căutării în baza de date

43

Figura 3.23 Familia secvenței ce a fost introdusă în secțiunea de căutare

Această căutare poate fi realizată și în Matlab, pentru primele 59 de familii din Pfam:

>> seqs = {or2,randor1};

for i = 2:60

for j = 2:2

[score(i,j)] = hmmprofalign(gethmmprof(i), seqs(j));

end

end

score

score = 119.8519 -71.5610 -39.5903 -56.9205 -46.7771 -58.2262 -20.1624 -22.5645 -25.9889 -27.5829 -41.4862 -51.7628 -11.3996 -7.8657 -4.7461 -16.1892 -51.8786 -41.7697 -16.3782 -12.3931 -12.6144 -16.4628 -403.3688 -414.6230

44

-11.3620 -17.1977 -6.9425 -11.1805 -41.9125 -52.7274 -92.9339 -93.4134 -15.0501 -17.2962 -16.1952 -18.1741 -15.4882 -15.2801 -4.4304 -11.0762 -8.8403 -9.1098 -166.9274 -170.5167 -8.1228 -5.8912 -10.7062 -15.5480 -38.3652 -41.6987 -96.6683 -95.9923 -5.1546 -15.1763 -2.7689 -5.8790 -19.1551 -25.8321 -14.7972 -21.3430 -18.5159 -12.7044 -18.8163 -12.4073 -32.2345 -41.5606 -22.1871 -17.9247 -11.2730 -11.3524 -7.5160 -11.2305 -16.4537 -17.7670 -106.6456 -91.8016 -7.0615 -13.7007 -6.3464 -13.4060 -13.2034 -14.0388 -25.0833 -21.0919 -20.8497 -22.0326 -32.0251 -34.1810 -12.8920 -3.2061 -13.7950 -17.8878 -15.6184 -15.3751 -13.2075 -6.2590 -11.8900 -5.8779 -9.5014 -5.5183 -9.0971 -13.5692 -19.9636 -28.4693 -4.5630 3.8388 -18.9590 -31.6370 -72.3591 -53.0879 -24.4002 -28.3283 -9.9745 -7.1494 -11.4843 -17.5292 -15.5195 -13.4378 -19.7734 -26.9954 Din scorul de mai sus se alege cea mai semnificativă potrivire.

45

Vom folosi în continuare familia 13853 detectată mai sus.

>> hmm7tm = gethmmprof(13853);

Alinierile multiple corespunzătoare:

>> seqs = gethmmalignment(13853, 'type', 'seed');

>>seqs

seqs =

118x1 struct array with fields:

Header

Sequence

Afișăm aceste rezultate (parțial):

>> disp([char(seqs.Header) char(seqs.Sequence)])

Figura 3.24 O parte din alinierile folosite pentru a genera MMA

46

Vom alinia profilul corespunzător cu primele 50 de secvențe din repertoriul nostru:

>> for i=2:50

[Score(i), Seqs(i).Aligned] = hmmprofalign(hmm7tm, orseqs(i).Sequence);

end

hmmprofmerge(Seqs,Score)

Figura 3.25 Secvențele aliniate sortate cu ajutorul scorurilor

47

Se face aceeași aliniere, ca și la câine, a secvențelor or și randor cu rhodopsina. Se face de asemenea și alinierea celor două secvențe or și randor cu MMA de profil. Pentru început trebuie returnată o secvență pentru rhodopsină și una pentru receptorii de miros.

>>rhod = getgenpept('NP_002368','sequenceonly',true) >>or = orseqs(282).Sequence rhod =

MDGSNVTSFVVEEPTNISTGRNASVGNAHRQIPIVHWVIMSISPVGFVENGIL...

or =

MAGKNYTFVTQFFLTAFTEHPEWGLPLFLLFLSFYLSTLLGNT

Acum se poate face alinierea globală de perechi a celor două secvențe. BLOSSUM 30 este matricea de punctare. Penalizările pentru deschiderea și extinderea unui spațiu în aliniere sunt setate la 5. [Score, Alignment] = nwalign(or, rhod2, 'scoringmatrix', 'blosum30', 'gapopen', 5, 'extendgap',5) perm = randperm(length(or)); randor = or(perm); [Score,Alignment]= nwalign(randor,rhod2, 'scoringmatrix', 'blosum30', 'gapopen', 5, 'extendgap', 5) randor = IFLEIPCNCFFTIVPGADFYESPGIREYKQIAYLSDKSGNGVP

Figura 3.18 Rezultatele (parțiale) alinierii rhodopsinei cu ultima secvență a receptorilor de miros folosind algoritmul Needleman și Wunsch

48

Figura 3.27 Rezultatele (parțiale) alinierii rhodopsinei cu o secvență aleatoare a receptorilor de miros folosind algoritmul Needleman și Wunsch

Dacă se aliniază ambele secvențe cu MMA, semnificația alinierii este mai evidentă.

[score_or, align_or] = hmmprofalign(hmm7tm,or) [score_rand, align_rand] = hmmprofalign(hmm7tm,randor)

Figura 3.28 Rezultatele (parțiale) alinierii ambelor secvențe cu MMA

49

MMA de profil la om

Se urmează aceeași pași prezentați mai sus. Rezultatul alinierii secvenței receptorului de miros cu rhodopsina, precum și rezultatul alinierii unei secvențe aleatoare cu rhodopsina sunt afișate în figura de mai jos:

Figura 3.29 Rezultatele (parțiale) alinierii folosind algoritmul Needleman și Wunsch

50

3.4 Arborele filogenetic

Arborii filogenetici la câine și șobolani

Vom creea un arbore filogenetic din membrii unei familii de proteine. Receptorii de miros fac parte, de fapt, dintr-o familie mult mai mare de proteine cunoscută sub numele de Receptori Cuplați la Proteina G. GRPC are cinci grupuri principale: Adhesion, Secretin, Glutamate, Frizzled/TAS2, Rhodopsin. Pentru a creea arborele sunt folosite câteva din aceste grupuri. Secvențele pot fi returnate din baza de date GenBank folosind funcția getgenbank. >>data = {'Adhesion 1' 'NP_001775'; 'Adhesion 2' 'NP_001965'; 'Glutamate 1' 'NP_000830'; 'Glutamate 2' 'NP_000836'; 'Rhod-Alpha 1' 'NP_001051'; 'Rhod-Alpha 2' 'NP_000946'; 'Rhod-Delta 1' 'NP_002368'; 'Rhod-Delta 2' 'NP_473372'}; >>for prot = 1:8 seqs(prot).Header = data{prot,1}; seqs(prot).Sequence = getgenpept(data{prot,2},'sequenceonly','true'); end Funcția Getgenpept-întoarce informații despre secvența din baza de date NCBI GenPept Getgenpept(data{prot,2},' sequenceonly ',true)-întoarce doar secvența intrării GetPept ca un vector de caractere

Pentru a construi arborele se calculează distanțele UPGMA folosind corecția Jukes-Cantor. >>distances = seqpdist(seqs,'Method','Jukes-Cantor'); >>tree = seqlinkage(distances,'UPGMA',seqs) Arborele este : >>h = plot(tree,'orient','bottom'); >>ylabel('Evolutionary distance') Funcția seqpdist selectează metoda ce va fi folosită pentru a calcula distanțele dintre fiecare pereche de secvențe iar funcția Seqlinkage construiește un arbore filogenetic din distanțele perechilor de bază. Tree=seqlinkage(distances,'UPGMA',seqs) trece o listă de nume pentru a eticheta nodurile frunză ale arborelui filogenetic

51

Figura 3.30 Arborele filogenetic Dacă se adaugă două secvențe ale receptorilor de miros trebuie refăcut arborele. data2 = {'Olfactory 1';'Olfactory 2'}; for prot = 1:2 seqs(prot+8).Header = data2{prot,1}; seqs(prot+8).Sequence = orseqs(prot).Sequence; end distances = seqpdist(seqs,'Method','Jukes-Cantor'); tree = seqlinkage(distances,'UPGMA',seqs) h = plot(tree,'orient','bottom'); ylabel('Evolutionary distance');

Figura 3.31 Arborele filogenetic după adăugarea celor două secvențe

52

4. Concluzii și dezvoltări posibile

1. Utilizarea MMA se dovedește a fi o unealtă puternică în analiza secvențelor biologice.

2. Repertoriul de gene olfactive (potențial codificatoare de receptori olfactivi) la câine și la șobolani este mult mai vast decât la om (347 la om, 1121 la câine, 1493 la șobolan).

>> orseqsH=fastaread('347OR.fasta')

orseqsH =

347x1 struct array with fields:

Header

Sequence

>> orseqsD=fastaread('Seq-dog-fasta.txt')

orseqsD =

1121x1 struct array with fields:

Header

Sequence

>> orseqsR=fastaread('Seq-rat-fasta.txt')

orseqsR =

1493x1 struct array with fields:

Header

Sequence

Și numărul genelor codificatoare este semnificativ mai mare la câine și șobolan:

>> orseqsH=fastaread('347OR.fasta')

orseqsH =

347x1 struct array with fields:

Header

Sequence

53

>> orseqsD=fastaread('Prot-dog.tfa')

orseqsD =

848x1 struct array with fields:

Header

Sequence

>> orseqsR=fastaread('Prot-rat.txt')

orseqsR =

1201x1 struct array with fields:

Header

Sequence

3. Cu codul Matlab de mai jos am testat existența unor gene identice în repertoriile celor 3 specii:

for i=1:347

for j=1:848

if strcmp(orseqsH(i).Sequence, orseqsD(i).Sequence)

aHD=aHD+1;

end

end

end

Și am găsit 0 identități între genele de la om și cele de la câine:

>> aHD

aHD =

0

54

Între om și șobolan:

>> aHR=0

for i=1:347

for j=1:1201

if strcmp(orseqsH(i).Sequence,orseqsR(i).Sequence)

aHR=aHR+1;

end

end

end

aHR =

0

Iar între câine și șobolan:

>> aDR=0;

for i=1:848

for j=1:1201

if strcmp(orseqsD(i).Sequence,orseqsR(i).Sequence)

aDR=aDR+1;

end

end

end

aDR =

0

55

4. Am creat și un dotplot pentru perechi din secvențele analizate.

>> seqdotplot(orseqsH(1),orseqsD(2))

Figura 4.1 Dotplot pentru o secvență de la om și o secvență de la câine

56

Sau

>> seqdotplot(orseqsH(1),orseqsR(2))

Figura 4.2 Dotplot pentru o secvență de la om și o secvență de la șobolan

57

Sau cu o fereastră alunecătoare de lungime 11 și stringență 7:

>> seqdotplot(orseqsH(30),orseqsR(10),11,7)

Figura 4.3 Dotplot pentru o secvență de la om și o secvență de la șobolan cu o fereastră alunecătoare de lungime 11 și strigență 7

58

Cu o stringență mai redusă:

>> seqdotplot(orseqsH(30),orseqsR(10),11,3)

Figura 4.4 Dotplot pentru o secvență de la om și o secvență de la șobolan cu o fereastă alunecătoare de lungime 11 și strigență 3

5. Se poate repeta aceeași analiză prezentată în lucrare și pentru repertoriul de gene olfactive de la maimuță. Acesta trebuie identificat din baze de date biologice.

59

Bibliografie

Cristianini Nello, Hahn Matthew. 2006. “ Introduction to Computational Genomics.” , http://www.computational-genomics.net/case_studies/olfactoryreceptors_demo.html,

accesat la data de 15.05.2013

Gerritsen Baillie Vivienne. 2013. Protein Spotlight http://web.expasy.org/spotlight/, accesat la data de 20.06.2013.

Pascale Quignon, Giraud Mathieu, Rimbault Maud, Lavigne Patricia, Tacher Sandrine, Morin Emmanuelle, Retout Elodie 2005 „Genome Biology”. “Dog Genetics and Whole Genome Mapping” ,

http://dogs.genouest.org/ORrepertoire.html, accesat la data de 3.05.2013. Pascale Quignon, Giraud Mathieu, Rimbault Maud, Lavigne Patricia, Tacher Sandrine, Morin Emmanuelle, Retout Elodie. 2005. ”Genome Biology” “The Dog and Rat Olfactory Receptor Repertoires.”, http://genomebiology.com/2005/6/10/R83 , accesat la data de 02.05.2013

Pavel Ana Brândușa, Vasile Cristian, Buiu Cătalin. 2011. Biomatematică și Bioinformatică. Concepte și Aplicații. p 113-115

Richard Axel, Linda Buck, Fred Hutchinson. 2004 . The Boulevard of Broken Genes Hidden Markov Models . p.61-77

Selzer Paul, Marhöfer Richard, Rohwer Andreas. 2008. Applied Bioinformatics an Introduction. p.5,61

Sperschneider Volker. 2008. Bioinformatics Problem Solving Paradigms . p.265-267

Stephen A. Krawetz, David D. Womble. 2003. Introduction to Bionformatics. p.357

Kal Renganathan Sharma. 2009. Bioinformatics Sequence Alingment and Markov Models . p.133-165

Warren Ewens, Grant Gregory. 2005. Statisical Methods in Bioinformatics : An Introduction, second edition. p.409-411

Zozulya Sergey, Echeverri Fernando, Nguyen Trieu. 2001. ”Genome Biology ” “The Human Olfactory Receptor Repertoire.” Vol. 2 , http://genomebiology.com/2001/2/6/research/0018