Universitaire biologie — jaar 3 · Bachelor Year 3

5Bio-informatica en sequentieanalyse

Een bioloog die zojuist een gen van een diepzeeworm heeft gesequenced plakt zijn $300$ aminozuren in een webformulier en verneemt drie seconden later dat het eiwit een verre neef van een menselijk kinase is, met $31\,\%$ identiteit over $280$ residuen en een kans $10^{-40}$ dat de gelijkenis toeval is. Achter die drie seconden zitten een algoritme voor dynamisch programmeren uit 1970, een statistische theorie van toevallige uitlijningen, substitutiematrices die uit duizenden eiwitfamilies zijn gedestilleerd, en een databank van zo’n honderd miljard residuen. Dit hoofdstuk gaat over de redenering in die zwarte doos: hoe twee sequenties zo worden uitgelijnd dat de uitlijning aantoonbaar de beste is, hoe de score iets gaat betekenen, hoe een overeenkomst van een toeval wordt onderscheiden, en hoe patronen worden gevonden in een genoom waar nog niemand naar heeft gekeken. De wiskunde is elementair — een recursie, een logaritme, een poissonverdeling — en zij is de moeite waard, omdat elke conclusie uit een sequentievergelijking erop rust.

5.1 Twee sequenties uitlijnen

Definitie 5.1 (Uitlijning en score)

Een uitlijning van twee sequenties schrijft ze boven elkaar, met hiaten (–) ertussen zodat de kolommen een residu met een residu of een residu met een hiaat paren, en geen enkele kolom twee hiaten paart. Haar score is de som over de kolommen van een substitutiescore $s(a,b)$ voor elk paar residuen en van een hiaatstraf voor elk hiaat: een lineaire straf $-d$ per hiaatpositie, of, realistischer, een affiene straf $-d - (k-1)e$ voor een reeks van $k$ hiaten, waarbij de openingskosten $d$ groter zijn dan de verlengingskosten $e$ , omdat één insertie van verscheidene residuen één enkele evolutionaire gebeurtenis is. Een globale uitlijning beslaat beide sequenties van begin tot eind; een lokale uitlijning zoekt het hoogst scorende paar deelreeksen en negeert de rest, en dat is wat men wil wanneer een gedeeld domein in twee verder ongerelateerde eiwitten zit.

Stelling 5.2 (Needleman–Wunsch)

Zij $x = x_{1}\dots x_{m}$ en $y = y_{1}\dots y_{n}$ , met lineaire hiaatstraf $d$ . Definieer $F(i,j)$ als de beste score van een globale uitlijning van de prefixen $x_{1}\dots x_{i}$ en $y_{1}\dots y_{j}$ . Dan geldt $F(i,0) = -id$ , $F(0,j) = -jd$ , en voor $i,j \ge 1$

F(i,j) = \max\bigl\{\,F(i-1,j-1) + s(x_{i},y_{j}),\; F(i-1,j) - d,\; F(i,j-1) - d\,\bigr\}.

$F(m,n)$ is de optimale globale score, een optimale uitlijning wordt teruggevonden door vanaf $(m,n)$ de keuzen terug te volgen die elk maximum opleverden, en de hele berekening kost $mn$ stappen. De variant van Smith–Waterman voor lokale uitlijning voegt $0$ als vierde mogelijkheid aan het maximum toe, zet de randen op $0$ en leest het antwoord af bij de grootste waarde in de tabel.

Bewijs. Beschouw de laatste kolom van een willekeurige uitlijning van de twee prefixen. Zij is een van drie dingen: $x_{i}$ boven $y_{j}$ , $x_{i}$ boven een hiaat, of een hiaat boven $y_{j}$ . Haar weghalen laat een uitlijning over van $(x_{1}\dots x_{i-1}, y_{1}\dots y_{j-1})$ , van $(x_{1}\dots x_{i-1}, y_{1}\dots y_{j})$ , of van $(x_{1}\dots x_{i}, y_{1}\dots y_{j-1})$ , waarvan de score hoogstens $F$ van dat paar is; en omgekeerd kan elk van die optimale uitlijningen met de bijbehorende laatste kolom worden verlengd. De beste score die op elk soort kolom eindigt is dus $F$ van het kortere paar plus de score van de kolom, en het optimum is de grootste van de drie. De randen liggen vast (tegenover een leeg prefix zijn alleen hiaten mogelijk). Inductie naar $i + j$ vult de tabel; het aantal cellen is $(m+1)(n+1)$ . Bij lokale uitlijning betekent de extra mogelijkheid $0$ “begin hier een nieuwe uitlijning”, wat van $F(i,j)$ de beste score maakt van een uitlijning die eindigt op $(i,j)$ , en de beste lokale uitlijning eindigt ergens. ∎

Voorbeeld 5.3 (Een tabel van vier bij drie)

Lijn GAT uit tegen GCAT, met $+1$ voor een overeenkomst, $-1$ voor een verschil en $d = 1$ . De randen zijn $0, -1, -2, -3, -4$ langs de bovenkant en $0, -1, -2, -3$ langs de zijkant. Rij voor rij invullen geeft $F(\text{G},\text{G}) = 1$ , $F(\text{G},\text{C}) = 0$ , $F(\text{G},\text{A}) = -1$ , $F(\text{G},\text{T}) = -2$ ; $F(\text{A},\text{G}) = 0$ , $F(\text{A},\text{C}) = 0$ , $F(\text{A},\text{A}) = 1$ , $F(\text{A},\text{T}) = 0$ ; $F(\text{T},\text{G}) = -1$ , $F(\text{T},\text{C}) = -1$ , $F(\text{T},\text{A}) = 0$ , $F(\text{T},\text{T}) = 2$ . Het optimum is $2$ , en terugvolgen — diagonaal vanaf (T,T), diagonaal vanaf (A,A), dan naar links van (G,C) naar (G,G), dan diagonaal — geeft

\begin{array}{c} \texttt{G-AT}\\ \texttt{GCAT} \end{array}

drie overeenkomsten en één hiaat: $3 - 1 = 2$ .

De tabel van Needleman–Wunsch voor GAT tegen GCAT (overeenkomst +1, verschil -1, hiaat -1). Elke cel is de beste score voor de twee prefixen die daar eindigen; het rode pad dat vanuit de hoek is teruggevolgd is de optimale uitlijning. — De tabel van Needleman–Wunsch voor GAT tegen GCAT (overeenkomst $+1$ , verschil $-1$ , hiaat $-1$ ). Elke cel is de beste score voor de twee prefixen die daar eindigen; het rode pad dat vanuit de hoek is teruggevolgd is de optimale uitlijning.

Methode 5.4 (Twee sequenties uitlijnen)

(1) Kies de scoring: een substitutiematrix die past bij de verwachte divergentie (BLOSUM62 voor eiwitten op onbekende afstand, overeenkomst/verschil voor DNA) en affiene hiaatstraffen (met BLOSUM62 meestal openen $-11$ , verlengen $-1$ ). (2) Beslis globaal of lokaal: globaal voor twee sequenties die over hun hele lengte homoloog worden geacht, lokaal in de overige gevallen. (3) Vul de tabel met de recursie, en bewaar bij elke cel een verwijzing naar de keuze die haar maximum gaf. (4) Volg terug vanaf $(m,n)$ (globaal) of vanaf de maximale cel tot een nul (lokaal), en schrijf de uitlijning van rechts naar links. (5) Beoordeel het resultaat niet naar zijn ruwe score maar naar zijn statistische significantie (hieronder), en kijk ernaar: lange hiaten, reeksen met weinig complexiteit en uitlijningen die tot een herhaling beperkt blijven zijn waarschuwingen.

5.2 Scoren: wat een overeenkomst waard is

Definitie 5.5 (Substitutiematrices)

Een substitutiematrix geeft $s(a,b)$ voor elk paar aminozuren als een log-oddsscore:

s(a,b) = \frac{1}{\lambda}\,\log\frac{q_{ab}}{p_{a}\,p_{b}},

waarin $q_{ab}$ de frequentie is waarmee $a$ en $b$ uitgelijnd worden aangetroffen in betrouwbare uitlijningen van verwante eiwitten, $p_{a} p_{b}$ de frequentie waarmee zij door toeval zouden worden gepaard, en $\lambda$ een schaal die zo is gekozen dat de waarden handige gehele getallen worden. Een positieve score betekent dat het paar vaker in homologen voorkomt dan door toeval; de identiteitsscores zijn het grootst voor zeldzame aminozuren (tryptofaan $+11$ , cysteïne $+9$ in BLOSUM62) en het kleinst voor veelvoorkomende (leucine $+4$ , alanine $+4$ ), en behoudende substituties (isoleucine–valine $+3$ ) scoren positief terwijl ingrijpende (tryptofaan–glycine $-2$ ) negatief scoren. De PAM-matrices (Dayhoff, 1978) werden afgeleid uit nauw verwante eiwitten en door matrixvermenigvuldiging naar grotere afstanden geëxtrapoleerd; de BLOSUM-matrices (Henikoff en Henikoff, 1992) werden rechtstreeks geteld in blokken uitgelijnde sequenties die op een gegeven identiteit waren geclusterd — BLOSUM62 uit blokken bij $62\,\%$ — en zij zijn de standaard omdat zij gemeten zijn en niet geëxtrapoleerd, op de afstand waar zij worden gebruikt.

Propositie 5.6 (Waarom log-odds)

Wil een scoringsschema bij lokale uitlijning bruikbaar zijn, dan moet de verwachte score van een willekeurig gepaarde kolom, $\sum_{a,b} p_{a} p_{b}\, s(a,b)$ , negatief zijn en moeten sommige scores positief zijn; anders zouden toevallige uitlijningen onbegrensd groeien en zou het hoogst scorende segment de hele sequentie zijn. Onder die voorwaarde is elk zulk schema gelijkwaardig aan een log-oddsschema voor zekere doelfrequenties $q_{ab}$ — de uitlijningen die het als optimaal zal vinden zijn die waarvan de residuparen als $q_{ab}$ zijn verdeeld. De matrix kiezen is dus de divergentie kiezen die men verwacht op te sporen: een matrix voor nauwe verwanten (BLOSUM80, PAM30) heeft scherpere positieve en hardere negatieve waarden, een matrix voor verre verwanten (BLOSUM45, PAM250) is vlakker.

Bewijs. Op dit niveau zonder bewijs aangenomen. ∎

Voorbeeld 5.7 (Identiteit, gelijkenis en de schemerzone)

Twee willekeurige eiwitsequenties die optimaal met hiaten worden uitgelijnd halen door toeval ongeveer $15\text{ tot }20\,\%$ identiteit. Boven $35\,\%$ identiteit over honderd residuen zijn twee eiwitten vrijwel zeker homoloog; tussen $20\,\%$ en $35\,\%$ ligt de schemerzone, waar identiteit alleen niets kan beslissen en de statistiek hieronder dat moet doen. Homologen kunnen ver onder de zone vallen: de subeenheden van hemoglobine en myoglobine delen $25\,\%$ identiteit, lysozym en $\alpha$ -lactalbumine $40\,\%$ , en veel paren eiwitten met dezelfde vouwing delen minder dan $15\,\%$ , alleen op te sporen door profielen of structuren te vergelijken.

5.3 Een databank doorzoeken

Definitie 5.8 (BLAST)

Een zoekvraag van $300$ residuen tegen een databank van $10^{11}$ uitlijnen met volledig dynamisch programmeren zou $3\times 10^{13}$ celbewerkingen per zoekactie kosten. BLAST (Altschul en collega’s, 1990) ruilt een beetje gevoeligheid in voor duizendvoudige snelheid met drie stappen: (1) maak een lijst van de woorden van de zoekvraag (drie residuen voor eiwitten, elf basen voor DNA) en van hun hoogscorende buren; (2) doorzoek de databank op exacte woordtreffers — de seeds; (3) verleng elke seed in beide richtingen zonder hiaten tot de score een vastgestelde hoeveelheid onder haar beste waarde zakt, en houd de hoogscorende segmentparen (HSP’s) over, en verbind daarna nabije HSP’s met gehiaat dynamisch programmeren in een smalle band. Een echte homoloog bevat vrijwel altijd minstens één exact woord van drie residuen; een toevallige gelijkenis zelden, en die wordt nooit verlengd.

De heuristiek van BLAST. Korte exacte woorden die de zoekvraag en een databankvermelding delen (rood) zijn seeds; elk wordt langs zijn diagonaal verlengd zolang de score blijft stijgen, en alleen verlengingen die hoog blijven worden hoogscorende segmentparen.

Stelling 5.9 (De statistiek van een toevallige treffer)

Voor een zoekvraag van lengte $m$ die wordt doorzocht tegen een databank van totale lengte $n$ , met een scoringsschema waarvan de verwachte score negatief is, is het aantal ongehiaatte lokale uitlijningen met een score van minstens $S$ dat door toeval ontstaat poissonverdeeld met gemiddelde

E = K\,m\,n\,e^{-\lambda S},

waarin $\lambda$ en $K$ alleen van het scoringsschema en de residufrequenties afhangen ( $\lambda$ is de schaal van de log-oddsmatrix). $E$ is de verwachtingswaarde van de score $S$ . Schrijft men de score in bits, $S' = (\lambda S - \ln K)/\ln 2$ , dan wordt de formule $E = m n\, 2^{-S'}$ , en de kans dat minstens één toevallige uitlijning $S$ haalt is $P = 1 - e^{-E}$ , wat gelijk is aan $E$ wanneer $E$ klein is.

Gedeeltelijk bewijs. De exponentiële staart is de stelling van Karlin–Altschul en wordt aangenomen: de maximale segmentscore van een toevalswandeling met negatieve drift heeft een verdeling waarvan de staart als $e^{-\lambda S}$ afneemt, met $\lambda$ de positieve wortel van $\sum_{a,b} p_{a} p_{b} e^{\lambda s(a,b)} = 1$ — en dat is precies de vergelijking die de log-oddsmatrix consistent maakt. Gegeven die staart is de rest tellen. Hoogscorende segmenten kunnen bij elk van de $mn$ paren posities beginnen, zij zijn zeldzaam en zij zijn vrijwel onafhankelijk; het aantal dat $S$ overtreft is dus poissonverdeeld met een gemiddelde dat evenredig is met $mn$ en met de staartkans, $E = Kmn\,e^{-\lambda S}$ . De kans op geen enkele is $e^{-E}$ . De vervanging door de bitscore is algebra: $e^{-\lambda S} K = 2^{-(\lambda S - \ln K)/\ln 2}$ . Voor gehiaatte uitlijningen geldt dezelfde vorm met $\lambda$ en $K$ uit simulatie geschat. ∎

Voorbeeld 5.10 (Een E-waarde lezen)

Een zoekvraag van $250$ residuen tegen een databank van $5\times 10^{10}$ residuen heeft $mn = 1.25\times 10^{13} \approx 2^{43.5}$ . Een treffer met een bitscore van $60$ heeft $E = 2^{43.5 - 60} = 2^{-16.5} \approx 10^{-5}$ : vrijwel zeker een homoloog. Een treffer met $S' = 40$ heeft $E = 2^{3.5} \approx 11$ : elf zulke scores zijn door toeval te verwachten, en de treffer betekent niets. Dezelfde uitlijning, met dezelfde bitscore, tegen een tien keer grotere databank doorzocht, heeft een tien keer grotere $E$ — significantie is een eigenschap van de zoekactie, niet van het paar. De gangbare drempel is $E < 10^{-3}$ voor een betrouwbare homoloog; $E \approx 0.01$ – $1$ verdient een tweede blik met een profielmethode.

E = mn\,2-S': elke extra bit halveert het verwachte aantal toevallige treffers, en een tien keer grotere databank kost 3.3 bits significantie voor dezelfde uitlijning. — $E = mn\,2^{-S'}$ : elke extra bit halveert het verwachte aantal toevallige treffers, en een tien keer grotere databank kost $3.3$ bits significantie voor dezelfde uitlijning.

5.4 Profielen, verborgen toestanden en motieven

Definitie 5.11 (Meervoudige uitlijning en profielen)

Een meervoudige sequentie-uitlijning zet een familie van sequenties in kolommen van homologe residuen. Exact dynamisch programmeren over $k$ sequenties kost $n^{k}$ en is voorbij drie onmogelijk; praktische programma’s lijnen stapsgewijs uit, eerst het naaste paar volgens een gidsboom, daarna sequenties en groepen bij de groeiende uitlijning, met ronden van verfijning. Een voltooide uitlijning wordt samengevat als een profiel: per kolom de frequentie van elk residu en van hiaten. Een profiel-HMM, een verborgen markovmodel, formaliseert dit als een keten van overeenkomsttoestanden, één per geconserveerde kolom, die elk residuen met hun eigen kansen uitzenden, met invoeg- en verwijdertoestanden die op elke positie extra of ontbrekende residuen toelaten; het model van een familie (een Pfam-vermelding) scoort een nieuwe sequentie met de kans van het beste pad door de toestanden, en vindt homologen ver onder de schemerzone van de paarsgewijze vergelijking, omdat een kolom die alleen hydrofobe residuen duldt dat ook zegt, terwijl één enkele sequentie dat niet kan.

Een profiel-HMM van een familie met vier kolommen. Elke overeenkomsttoestand M zendt een residu uit met de eigen frequenties van de kolom; invoegtoestanden I (met zelflussen) laten extra residuen toe, verwijdertoestanden D slaan een kolom over. Een sequentie scoren is haar meest waarschijnlijke pad zoeken.

Definitie 5.12 (Motieven en informatie-inhoud)

Een motief is een kort patroon — een bindingsplaats voor een transcriptiefactor, een splitssignaal, een fosforyleringsplaats — weergegeven door een positiegewichtsmatrix van de frequentie $f_{i}(b)$ van elke base of elk residu $b$ op elke positie $i$ . De informatie-inhoud van positie $i$ is $R_{i} = 2 - H_{i}$ bits voor DNA, waarbij $H_{i} = -\sum_{b} f_{i}(b)\log_{2} f_{i}(b)$ haar entropie is: $2$ bits voor een onveranderlijke base, $0$ voor een positie waar alle vier even waarschijnlijk zijn. Het totaal $R = \sum_{i} R_{i}$ wordt getekend als een sequentielogo, waarbij elke positie een stapel letters is waarvan de totale hoogte $R_{i}$ is en waarvan de letters op frequentie zijn geschaald.

Propositie 5.13 (Hoeveel informatie een plaats nodig heeft)

Een plaats die $\gamma$ keer in een genoom van $G$ posities moet worden gevonden, en nergens anders, heeft ongeveer $R_{\text{nodig}} = \log_{2}(G/\gamma)$ bits informatie-inhoud nodig: het motief moet de $G$ kandidaatposities tot de $\gamma$ echte terugbrengen, en elke bit halveert het aantal kandidaten. De waargenomen motieven van goed onderzochte bacteriële regulatoren voldoen aan die voorspelling — plaatsen in E. coli voor een repressor die enkele tientallen plekken in een genoom van $4.6\,\mathrm{Mb}$ bindt dragen $16\text{ tot }18$ bits; motieven van eukaryote transcriptiefactoren, met $8\text{ tot }12$ bits in een genoom van $3\times 10^{9}$ , kunnen hun doelwitten niet alleen aanwijzen, en daarom werken zij in combinaties en in het open chromatine van Hoofdstuk 1.

Bewijs. Een willekeurige positie past op een motief met informatie-inhoud $R$ met kans ongeveer $2^{-R}$ (elke bit specificiteit halveert de kans), dus is het verwachte aantal toevallige treffers in $G$ posities $G\,2^{-R}$ . Willen de echte plaatsen eruit springen, dan moet dit van de orde $\gamma$ of kleiner zijn: $G\,2^{-R} \le \gamma$ , dat wil zeggen $R \ge \log_{2}(G/\gamma)$ . ∎

Voorbeeld 5.14 (Verwachte toevallige treffers)

Een restrictieplaats van zes vaste basen heeft $R = 12$ bits en past op een willekeurige positie met kans $4^{-6} = 2^{-12}$ : ongeveer $1100$ keer in een genoom van E. coli van $4.6\,\mathrm{Mb}$ op beide strengen gelezen (de plaats is palindroom, dus één keer per positie), en $7\times 10^{5}$ keer in het menselijk genoom. Een eukaryote factor waarvan het motief $10$ bits draagt past op $3\times 10^{9}\times 2^{-10} \approx 3$ miljoen posities in het menselijk genoom, enkele duizenden keren meer dan de genen die zij reguleert. Een motief alleen is een zwakke voorspeller in een groot genoom; de toestand van het chromatine, de naburige motieven en de conservering van de plaats tussen soorten maken pas een voorspelling.

Een sequentielogo van een promotormotief van het TATA-boxtype. De hoogte van elke stapel is de informatie-inhoud van die positie, 2 - H_i bits; de eerste vier posities zijn vrijwel onveranderlijk en dragen het grootste deel van de ongeveer 12 bits van het motief. — Een sequentielogo van een promotormotief van het TATA-boxtype. De hoogte van elke stapel is de informatie-inhoud van die positie, $2 - H_{i}$ bits; de eerste vier posities zijn vrijwel onveranderlijk en dragen het grootste deel van de ongeveer $12$ bits van het motief.

5.5 Van sequentie naar functie

Methode 5.15 (Een onbekend eiwit annoteren)

Gegeven een nieuwe coderende sequentie: (1) vertaal haar in het juiste leesraam en kijk of er een signaalpeptide, transmembraansegmenten en gebieden met weinig complexiteit zijn; (2) doorzoek de eiwitdatabanken met BLAST en lees de treffers met $E < 10^{-3}$ , en let erop of de uitlijning het hele eiwit beslaat (een echte ortholoog) of een segment (een gedeeld domein); (3) doorzoek de domeindatabanken met profiel-HMM’s, die families vinden die BLAST mist en die het eiwit in domeinen verdelen; (4) leid orthologie af, niet enkel gelijkenis, door na te gaan of de beste treffer in het andere genoom de zoekvraag als zijn beste treffer heeft (wederzijds beste treffers) of door het eiwit in een genboom te plaatsen (Hoofdstuk 25); (5) draag de functie van orthologen behoedzaam over — een geconserveerd katalytisch residu pleit voor dezelfde chemie, een ontbrekend ertegen — en voorspel de structuur; (6) behandel elke voorspelling als een hypothese voor de werkbank.

Propositie 5.16 (Structuur uit sequentie)

De vouwing van een eiwit wordt bepaald door zijn sequentie (Hoofdstuk 7), en haar uit de sequentie berekenen was vijftig jaar lang het centrale onopgeloste probleem van het vak. Drie benaderingen slaagden na elkaar. Homologiemodellering bouwt de structuur van een eiwit op die van een opgeloste homoloog, en dat gaat betrouwbaar boven $30\,\%$ identiteit. Co-evolutieanalyse maakt gebruik van het feit dat twee residuen die in de vouwing contact maken over een diepe meervoudige uitlijning geneigd zijn samen te muteren, zodat statistisch gekoppelde paren kolommen voorspelde contacten zijn, en genoeg contacten leggen een vouwing vast. Methoden uit het diepe leren, getraind op de honderdduizend opgeloste structuren en op zulke uitlijningen, voorspellen de meeste bolvormige eiwitstructuren nu met bijna experimentele nauwkeurigheid (de CASP-beoordelingen van 2020), en databanken bevatten een voorspelde structuur voor vrijwel elke bekende eiwitsequentie. Wat zij minder goed voorspellen is wat één enkele structuur niet vat: wanordelijke gebieden, andere conformaties, het effect van een puntmutatie en complexen.

Links: een voorspelde eiwitstructuur, gekleurd naar het vertrouwen van het model, van hoog (blauw) tot laag (oranje) in een wanordelijke lus. Rechts: een bio-informaticakantoor — genoombrowsers en bomen op de schermen, en geen natte werkbank te bekennen.

Opmerking 5.17 (De grenzen van de gevolgtrekking)

De meeste functionele annotaties in de databanken zijn nooit getoetst; zij zijn overgenomen van een homoloog, die op haar beurt door overname was geannoteerd. Fouten planten zich voort en vermenigvuldigen zich, en een verkeerde annotatie op een goed verbonden eiwit kan een hele familie besmetten. De remedies zijn de bovenstaande: onderscheid orthologie van homologie, lees de uitlijning, zoek de katalytische residuen, en onthoud dat “hypothetisch eiwit” een eerlijk etiket is dat een derde van de genen in de meeste genomen nog verdient.

5.6 Opgaven

Oefening 5.1 ★

Definieer globale en lokale uitlijning en geef voor elk één biologische situatie die erom vraagt.

Oplossing

Oplossing van Oefening 5.1.

Globaal: beide sequenties van begin tot eind uitgelijnd, elk residu in een kolom — voor twee eiwitten die over hun hele lengte homoloog worden geacht, zoals orthologen van een huishoudenzym. Lokaal: het hoogst scorende paar deelreeksen, de rest genegeerd — om een gedeeld domein te vinden (een SH2-domein in twee verder ongerelateerde signaaleiwitten), of een gen in een lange genoomsequentie.

Oefening 5.2 ★

Vul de tabel van Needleman–Wunsch voor AGC tegen AAC in met overeenkomst $+1$ , verschil $-1$ , hiaat $-1$ , en geef de optimale uitlijning en de score.

Oplossing

Oplossing van Oefening 5.2.

Randen $0,-1,-2,-3$ beide kanten op. Rij A: $1, 0, -1$ . Rij G: $0, 0, -1$ . Rij C: $-1, -1, 1$ . Optimum $F(3,3) = 1$ : AGC boven AAC zonder hiaten (overeenkomst, verschil, overeenkomst: $1 - 1 + 1 = 1$ ).

Oefening 5.3 ★

In BLOSUM62 scoort tryptofaan–tryptofaan $+11$ en leucine–leucine $+4$ . Leg met de log-oddsformule uit waarom de identiteit van het zeldzamere residu meer waard is.

Oplossing

Oplossing van Oefening 5.3.

$s(a,a) = \lambda^{-1}\log\bigl(q_{aa}/p_{a}^{2}\bigr)$ . Tryptofaan is zeldzaam ( $p_{W} \approx 0.013$ ), dus de kans dat twee tryptofanen toevallig worden uitgelijnd, $p_{W}^{2}$ , is uiterst klein, en een geconserveerd tryptofaanpaar is een veel sterker teken van homologie dan een geconserveerd leucinepaar ( $p_{L} \approx 0.1$ ); de log-oddsverhouding is navenant groter.

Oefening 5.4 ★

Wat is een E-waarde? Een zoekactie levert een treffer met $E = 3$ op. Wat betekent dat getal, en is de treffer een homoloog?

Oplossing

Oplossing van Oefening 5.4.

De E-waarde is het aantal uitlijningen met een score van minstens die hoogte dat door toeval te verwachten is bij een zoekactie van deze zoekvraag tegen een databank van deze omvang. $E = 3$ betekent dat drie zulke scores door toeval worden verwacht: de treffer is geen bewijs van homologie (zij kan het nog steeds zijn, maar de zoekactie kan het niet zeggen).

Oefening 5.5 ★★

Een zoekvraag van $400$ residuen wordt doorzocht tegen $2\times 10^{11}$ residuen. Bereken de E-waarde van treffers met bitscores $45$ , $55$ en $65$ . Welke bitscore geeft $E = 10^{-3}$ ? Hoe verandert het antwoord als de zoekvraag $40$ residuen lang is?

Oplossing

Oplossing van Oefening 5.5.

$mn = 400\times 2\times 10^{11} = 8\times 10^{13} = 2^{46.2}$ . $E(45) = 2^{1.2} \approx 2.3$ ; $E(55) = 2^{-8.8} \approx 2\times 10^{-3}$ ; $E(65) = 2^{-18.8} \approx 2\times 10^{-6}$ . Voor $E = 10^{-3}$ is $S' = 46.2 + 10.0 = 56$ bits nodig. Een zoekvraag van $40$ residuen heeft een tien keer kleinere $mn$ , $2^{42.9}$ : $53$ bits volstaan — maar een korte zoekvraag haalt zelfs dat zelden.

Oefening 5.6 ★★

Bereken de informatie-inhoud van een motief waarvan de vier posities basenfrequenties (A, C, G, T) van $(1,0,0,0)$ , $(0.5,0,0.5,0)$ , $(0.25,0.25,0.25,0.25)$ en $(0.7,0.1,0.1,0.1)$ hebben. Hoeveel toevallige treffers heeft het in een genoom van $4.6\,\mathrm{Mb}$ ?

Oplossing

Oplossing van Oefening 5.6.

Informatie-inhouden: $2$ , $1$ , $0$ , en $2 - H$ met $H = -(0.7\log_{2} 0.7 + 3\times 0.1\log_{2} 0.1) = 0.36 + 1.00 = 1.36$ , dus $0.64$ . Totaal $R = 3.64$ bits. Toevallige treffers: $9.2\times 10^{6}$ posities op twee strengen $\times 2^{-3.64} \approx 7\times 10^{5}$ — het motief is alleen vrijwel nutteloos.

Oefening 5.7 ★★

Leg uit waarom affiene hiaatstraffen realistischer zijn dan lineaire, en waarom zowel een zeer hoge als een zeer lage straf op het openen van een hiaat slechte uitlijningen geeft.

Oplossing

Oplossing van Oefening 5.7.

Een insertie van verscheidene residuen is één mutatiegebeurtenis, dus zouden haar kosten niet lineair met haar lengte moeten groeien: een openingskost plus een kleine verlengingskost bootst dat na. Een te hoge openingsstraf dwingt verschillen af waar een hiaat hoort en lijnt alles na een echte insertie verkeerd uit; een te lage straf strooit overal hiaten, laat residuen toevallig overeenkomen en blaast de identiteit op.

Oefening 5.8 ★★

Een BLAST-zoekactie van een menselijk eiwit tegen een vliegendatabank geeft een beste treffer met $E = 10^{-30}$ die de residuen 50–180 van de zoekvraag van $600$ residuen beslaat. Is het vliegeneiwit de ortholoog van het menselijke? Welke aanvullende toets ligt voor de hand?

Oplossing

Oplossing van Oefening 5.8.

Niet noodzakelijk: de uitlijning beslaat een segment van $130$ residuen, en dat is de handtekening van een gedeeld domein en niet van een ortholoog die over zijn hele lengte is uitgelijnd. Toets: doorzoek met het vliegeneiwit het menselijk eiwitbestand terug (is de zoekvraag zijn beste treffer, over de hele lengte?), bepaal het domein met een profiel-HMM, en bouw een genboom van de familie in verscheidene soorten.

Oefening 5.9 ★★

Waarom sporen profielmethoden homologen op die de paarsgewijze uitlijning mist? Geef een voorbeeld van een kolompatroon dat een profiel vat en één enkele sequentie niet.

Oplossing

Oplossing van Oefening 5.9.

Een profiel legt kolom voor kolom vast wat de familie duldt: een positie die altijd hydrofoob is maar nooit hetzelfde residu, een onveranderlijk katalytisch residu, een positie die bij de helft van de familie altijd een hiaat is. Een paarsgewijze uitlijning scoort elk residu tegen één ander residu en kan niet weten dat een valine op positie 40 “even goed is als” het isoleucine dat de zoekvraag daar heeft. Het profiel weegt bovendien de geconserveerde kolommen, zodat zwakke gelijkenis die zich juist daar samenbalt significant wordt.

Oefening 5.10 ★★★

Laat zien dat onder een scoringsschema met positieve verwachte score de lokale uitlijning volgens Smith–Waterman van twee lange willekeurige sequenties een score heeft die lineair met hun lengte groeit, en leg uit waarom de theorie van de E-waarde daardoor faalt. Wat betekent dit voor het uitlijnen van DNA met overeenkomst $+1$ en verschil $-1$ bij een GC-gehalte van $60\,\%$ ?

Oplossing

Oplossing van Oefening 5.10.

Bij een positieve verwachte score $\mu > 0$ per kolom is de cumulatieve score langs de diagonaal van twee willekeurige sequenties een toevalswandeling met positieve drift: na $n$ kolommen is zij ongeveer $\mu n$ , zodat de beste lokale uitlijning in wezen het geheel is en haar score als $\mu n$ groeit in plaats van als $\log n$ . De theorie van Karlin–Altschul, die een negatieve drift vergt zodat hoge scores zeldzame uitstapjes zijn, geldt dan niet en er bestaat geen $\lambda$ . Voor DNA bij $60\,\%$ GC is de kans op een overeenkomst $2(0.3^{2}) + 2(0.2^{2}) = 0.26$ , zodat de verwachte score $0.26 - 0.74 = -0.48$ is: nog steeds negatief, en de statistiek houdt stand; maar een schema met overeenkomst $+1$ en verschil $-0.3$ zou verwachting $+0.04$ hebben en zou het hele genoom als één uitlijning melden.

Oefening 5.11 ★★★

De tabel van Needleman–Wunsch vergt $mn$ geheugencellen; voor twee chromosomen van $100\,\mathrm{Mb}$ is dat $10^{16}$ . Beschrijf twee ideeën waarmee genoomuitlijners dat vermijden (seeds en ketenen; banden), en wat elk daarvan opgeeft.

Oplossing

Oplossing van Oefening 5.11.

Seeds en ketenen: zoek met een hashtabel exacte of bijna exacte treffers van $k$ -meren tussen de twee sequenties, houd die welke op samenhangende diagonalen liggen, keten ze aaneen, en voer dynamisch programmeren alleen uit in de gaten tussen geketende seeds; dat geeft uitlijningen op in gebieden zonder seed (sterk uiteenlopende stukken). Banden: zijn de twee sequenties bekend als vrijwel collineair, bereken dan alleen de cellen binnen een band van breedte $w$ rond de diagonaal, tegen kosten $wn$ in plaats van $mn$ ; dat geeft elke uitlijning op met een insertie die groter is dan de band.

Oefening 5.12 ★★★

Een HMM voor het vinden van genen bij bacteriën heeft toestanden voor de drie codonposities en voor niet-coderend DNA. Leg uit hoe het model coderende van niet-coderende sequentie kan onderscheiden zonder enige informatie over stopcodons (denk aan het codongebruik), en waarom dezelfde aanpak in een menselijk genoom veel moeilijker is.

Oplossing

Oplossing van Oefening 5.12.

Coderende sequentie heeft een periode van drie: de drie codonposities hebben verschillende basensamenstellingen (de derde is het meest scheef), en het codongebruik is per soort ongelijk. Een model met drie coderende toestanden achter elkaar, die elk basen uitzenden met de samenstelling van die codonpositie, geeft coderend DNA over een venster van enkele tientallen codons een hogere kans dan de niet-coderende toestand doet, zelfs zonder stopcodons. In een menselijk genoom zijn de exons kort ( $150\,\mathrm{bp}$ ) en gescheiden door introns van kilobasen, zodat het coderende signaal kort en onderbroken is; het model moet bovendien splitsplaatsen herkennen, die zwakke signalen zijn, en de enorme hoeveelheid niet-coderende sequentie levert vele valse coderende segmenten op.

5.7 Vraagstuk: een sequentie uit de diepzee

Probleem 5.1

Weekendvraagstuk — een onbekend eiwit met de hand uitgelijnd, tegen de databanken doorzocht met berekening van zijn significantie, zijn regulerende motief in bits gewogen, en zijn gen getoetst aan de statistiek van toevallige open leesramen, met als besluit de E-waarde van de beste treffer, het aantal bits dat een plaats nodig heeft en de lengte die een leesraam moet hebben om geloofd te worden

Gegevens: een eiwit van $300$ residuen uit een diepzeeworm. Eiwitdatabank: $1.2\times 10^{11}$ residuen. Genoom van de worm: $1.6\,\mathrm{Gb}$ , $38\,\%$ GC. Scoring voor uitlijningen met de hand: overeenkomst $+1$ , verschil $-1$ , hiaat $-1$ . Bitscore van de beste BLAST-treffer: $92$ ; van de tiende treffer: $38$ .

Deel I — Met de hand.

Lijn de peptiden KQT en KAQT uit met de recursie van Needleman–Wunsch: schrijf de tabel op en geef de optimale uitlijning en de score.
Doe hetzelfde met Smith–Waterman (lokaal) voor GATCAT tegen ACAT: zoek de beste lokale uitlijning en haar score.
Hoeveel celbewerkingen kost een globale uitlijning van het eiwit van $300$ residuen tegen een eiwit van $450$ residuen? En tegen de hele databank?
Als een computer $10^{9}$ bewerkingen per seconde uitvoert, hoe lang duurt de volledige databankuitlijning van vraag 3 dan? Waarom wordt in plaats daarvan BLAST gebruikt?
Een identiteitsscore in BLOSUM62 is $+4$ voor alanine ( $p_{A} = 0.074$ ) en $+11$ voor tryptofaan ( $p_{W} = 0.013$ ). Bereken met $\lambda = 0.347$ (eenheden van een halve bit) de doelfrequenties $q_{AA}$ en $q_{WW}$ , en de verhouding $q/p^{2}$ voor elk. Geef de betekenis.
Twee eiwitten delen $24\,\%$ identiteit over $250$ residuen. Zeg waarom identiteit alleen hier de homologie niet kan beslechten en wat dat wel zou doen.

Deel II — De zoekactie.

Bereken $mn$ voor de zoekvraag tegen de databank, en $\log_{2}(mn)$ .
Bereken de E-waarde van de beste treffer ( $S' = 92$ ) en van de tiende treffer ( $S' = 38$ ).
Welke bitscore komt bij deze zoekactie overeen met $E = 10^{-3}$ ? En met $E = 1$ ?
Dezelfde beste treffer wordt gevonden wanneer de databank tot $1.2\times 10^{12}$ residuen is gegroeid. Wat is haar E-waarde?
De tiende treffer lijnt de residuen 200–260 van de zoekvraag uit met $40\,\%$ identiteit over $60$ residuen. Zeg met de E-waarde of dit bewijs van homologie is, en wat een profielzoekactie zou kunnen toevoegen.
De beste treffer is een menselijk kinase, uitgelijnd over de residuen 10–290. Zijn beste treffer in het eiwitbestand van de worm is de zoekvraag. Wat stelt deze wederzijdse toets vast, en wat niet?

Deel III — Een motief.

Stroomopwaarts van het gen ligt een kandidaat-bindingsplaats voor een transcriptiefactor van acht posities met informatie-inhouden $2, 2, 1.6, 2, 0.8, 1.2, 0.4, 0.3$ bits. Wat is $R$ in totaal?
Hoeveel toevallige treffers heeft het motief in het genoom van $1.6\,\mathrm{Gb}$ (beide strengen, $3.2\times 10^{9}$ posities)?
De factor reguleert ongeveer $200$ genen. Hoeveel bits zou een motief nodig hebben om in dit genoom $200$ plaatsen alleen aan te wijzen?
Hoeveel van het tekort zou een tweede, aangrenzend motief van $8$ bits kunnen leveren, als de twee binnen een vaste tussenafstand samen moeten voorkomen?
Een positie met frequenties $(0.5, 0.5, 0, 0)$ voor (A, C, G, T): bereken haar entropie en haar informatie-inhoud.
Leg met het informatieargument uit waarom bacteriële transcriptiefactoren doorgaans langere en beter geconserveerde plaatsen hebben dan eukaryote.

Deel IV — Het gen zelf.

Wat is in willekeurig DNA met gelijke basensamenstelling de kans dat een codon een stopcodon is? Wat is het verwachte aantal codons voordat er een stopcodon verschijnt (een geometrische verdeling)?
Het genoom van de worm is voor $38\,\%$ GC. Herbereken met de werkelijke basenfrequenties de kans dat een willekeurig codon een stopcodon is (TAA, TAG, TGA), en de verwachte lengte van het leesraam. Welke kant op duwt een laag GC-gehalte het vinden van genen?
Wat is de kans dat een willekeurig open leesraam minstens $100$ codons lang is? En minstens $300$ ?
In het genoom van $1.6\,\mathrm{Gb}$ geven zes leesramen op twee strengen ongeveer $3.2\times 10^{9}$ codonbeginpunten. Hoeveel willekeurige open leesramen van minstens $100$ codons zijn te verwachten? En van minstens $300$ ?
Leg uit waarom “open leesraam langer dan $100$ codons” bij een bacterie een bruikbare genzoeker is en in dit genoom niet, en wat een eukaryote genzoeker in plaats daarvan gebruikt.
Het gen van de worm heeft zes exons van gemiddeld $150\,\mathrm{bp}$ . Leg uit hoe reads uit RNA-sequencing de exonstructuur vaststellen die de genoomsequentie alleen onduidelijk laat.
Vat samen: de E-waarde van de beste treffer (vraag 8), het aantal bits dat nodig is om $200$ plaatsen aan te wijzen (vraag 15), en het verwachte aantal willekeurige leesramen van $300$ codons in het genoom (vraag 22).

Oplossing

Oplossing van Probleem 5.1.

1. Rijen K, Q, T; kolommen K, A, Q, T; randen $0,-1,-2,-3,-4$ en $0,-1,-2,-3$ . Rij K: $1, 0, -1, -2$ ; rij Q: $0, 0, 1, 0$ ; rij T: $-1, -1, 0, 2$ . Optimum $2$ : K-QT boven KAQT. 2. Beste lokale score $3$ : CAT tegen CAT (residuen 4–6 van GATCAT met 2–4 van ACAT); ATCAT tegen A-CAT scoort eveneens $4 - 1 = 3$ . 3. $300\times 450 = 1.35\times 10^{5}$ bewerkingen; tegen de databank $300\times 1.2\times 10^{11} = 3.6\times 10^{13}$ . 4. $3.6\times 10^{4}$ s, tien uur per zoekvraag; de seeds van BLAST slaan vrijwel de hele tabel over en antwoorden in seconden. 5. $q_{ab} = p_{a}p_{b}e^{\lambda s}$ . Alanine: $e^{1.39} = 4.0$ , $q_{AA} = 0.074^{2}\times 4.0 = 0.022$ , verhouding $4$ . Tryptofaan: $e^{3.82} = 45$ , $q_{WW} = 0.013^{2}\times 45 = 0.0077$ , verhouding $45$ . Een uitgelijnd tryptofaanpaar komt $45$ keer vaker voor in homologen dan door toeval, een alaninepaar maar vier keer; alanineparen zijn in absolute zin niettemin talrijker omdat alanine veel voorkomt. 6. $24\,\%$ ligt in de schemerzone, waar toevallige uitlijningen $15\text{ tot }20\,\%$ halen; de E-waarde van de uitlijning, geconserveerde motieven op de juiste posities, een treffer van een profiel-HMM op een bekende familie, of een gedeelde vouwing zouden het beslechten. 7. $mn = 300\times 1.2\times 10^{11} = 3.6\times 10^{13}$ ; $\log_{2}(mn) = 45.0$ . 8. $E(92) = 2^{45 - 92} = 2^{-47} \approx 7\times 10^{-15}$ ; $E(38) = 2^{7} = 128$ . 9. $E = 10^{-3}$ bij $S' = 45 + 10 = 55$ bits; $E = 1$ bij $45$ bits. 10. Tien keer $mn$ : $E \approx 7\times 10^{-14}$ , nog altijd overweldigend. 11. Met $E = 128$ is de tiende treffer wat het toeval oplevert; $40\,\%$ identiteit over $60$ residuen is geen bewijs. Een profielzoekactie van de residuen 200–260 tegen de domeindatabank zou kunnen laten zien of dat segment een bekend domein is, met een statistiek die de paarsgewijze vergelijking niet heeft. 12. Wederzijds beste treffers over de volle lengte zijn verenigbaar met een-op-een-orthologie; zij bewijzen haar niet — een duplicatie in één lijn na de splitsing geeft twee co-orthologen, en het verlies van de echte ortholoog kan een paraloog als beste treffer achterlaten. Een genboom met verscheidene soorten is de toets. 13. $R = 2 + 2 + 1.6 + 2 + 0.8 + 1.2 + 0.4 + 0.3 = 10.3$ bits. 14. $3.2\times 10^{9}\times 2^{-10.3} \approx 2.5\times 10^{6}$ toevallige treffers. 15. $\log_{2}(3.2\times 10^{9}/200) = \log_{2}(1.6\times 10^{7}) \approx 24$ bits. 16. Samen voorkomen op een vaste tussenafstand telt de bits op: $10.3 + 8 = 18.3$ , wat $8$ van de $13.7$ ontbrekende levert; ongeveer $5.7$ bits (een factor $50$ in toevallige treffers) moet van elders komen — toegankelijkheid van het chromatine, verdere partners. 17. $H = -(0.5\log_{2}0.5 + 0.5\log_{2}0.5) = 1$ bit; $R = 2 - 1 = 1$ bit. 18. Een bacteriële factor moet haar enkele plaatsen vinden in een genoom van $4.6\,\mathrm{Mb}$ zonder hulp van het chromatine: zij heeft zo’n $19$ bits nodig, en haar plaatsen zijn lang en geconserveerd. Een eukaryoot genoom is duizend keer groter en vraagt tien bits meer, en toch hebben de factoren korte plaatsen; zij bereiken specificiteit door combinatie en door de beperking tot toegankelijk chromatine, wat de regulatie ook beter evolueerbaar maakt, omdat een korte plaats gemakkelijk wordt gewonnen of verloren. 19. $3/64 = 0.047$ ; het verwachte aantal codons vóór een stopcodon is $64/3 \approx 21$ . 20. $p_{A} = p_{T} = 0.31$ , $p_{G} = p_{C} = 0.19$ : $P(\text{TAA}) = 0.31^{3} = 0.030$ , $P(\text{TAG}) = P(\text{TGA}) = 0.31^{2}\times 0.19 = 0.018$ ; samen $0.066$ , verwachte leesraamlengte $15$ codons. AT-rijk DNA zit vol stopcodons, dus zijn willekeurige open leesramen korter en springen lange er sterker uit. 21. $(61/64)^{100} = e^{-4.80} = 0.008$ ; $(61/64)^{300} = e^{-14.4} = 5.6\times 10^{-7}$ . 22. Elk maximaal open leesraam eindigt op een stopcodon, en $3.2\times 10^{9}$ codonbeginpunten bevatten $3.2\times 10^{9}\times 3/64 = 1.5\times 10^{8}$ stopcodons: ongeveer $1.5\times 10^{8}\times 0.008 = 1.2\times 10^{6}$ willekeurige leesramen van minstens $100$ codons, en $1.5\times 10^{8}\times 5.6\times 10^{-7} \approx 80$ van minstens $300$ . 23. Een bacterie van $4.6\,\mathrm{Mb}$ heeft zo’n $4\times 10^{5}$ stopcodons en dus ongeveer $3500$ toevallige leesramen van $100$ codons maar vrijwel geen van $300$ ; haar genen tellen gemiddeld $300$ codons en $88\,\%$ van het DNA is coderend, zodat een lang open leesraam vrijwel altijd een gen is. Bij de worm codeert $1.5\,\%$ van het DNA, tellen de exons gemiddeld $50$ codons — korter dan de toevalsdrempel — en overspoelen een miljoen toevallige leesramen van $100$ codons hen. Eukaryote genzoekers gebruiken signalen van splitsplaatsen, codonvoorkeur in een verborgen markovmodel, homologie met bekende eiwitten en, vooral, gesequencede transcripten. 24. Een read uit een gesplitste boodschapper lijnt in twee stukken op het genoom uit, gescheiden door een intron: de splitsing markeert beide splitsplaatsen tot op de base; de readdekking bakent de exons af en gepaarde reads verbinden opeenvolgende exons tot één transcript, wat uitwijst welke van verscheidene kandidaat-splitsplaatsen wordt gebruikt. 25. $E \approx 7\times 10^{-15}$ voor de beste treffer; ongeveer $24$ bits om $200$ plaatsen in het genoom aan te wijzen; zo’n $80$ toevallige leesramen van $300$ codons in het hele genoom.

Begrippen gedefinieerd in dit hoofdstuk

Bekijk alle 479 begrippen in de begrippenlijst