Technische Universiteit Delft
Faculteit Elektrotechniek, Wiskunde en Informatica
Delft Institute of Applied Mathematics
Het rangschikken van genen (Gene Ranking)
Verslag ten behoeve van het
Delft Institute of Applied Mathematics
als onderdeel ter verkrijging
van de graad van
BACHELOR OF SCIENCE
in
TECHNISCHE WISKUNDE
door
Cindy Boon
Delft, Nederland
Juli 2014
c 2014 door Cindy Boon. Alle rechten voorbehouden.
Copyright BSc verslag TECHNISCHE WISKUNDE
“Het rangschikken van genen”
(Gene Ranking)
Cindy Boon
Technische Universiteit Delft
Begeleider
Martin van Gijzen
Overige commissieleden
M. Keijzer
C. Vuik
C. Kraaijkamp
Juli, 2014
Delft
2
SAMENVATTING
E´en van de grootste uitdagingen in de geneeskunde is op het moment het vinden van onderliggende
oorzaken voor genetische ziektes. Een mogelijk hulpmiddel in dit onderzoek is het rangschikken van genen
op mate van belang. Dit rapport beschrijft een onderzoek naar dit rangschikken van genen. Een menselijk
lichaam bevat 20.000-25.000 genen, die doordat ze in verbinding staan met elkaar, zijn te ordenen naar
mate van belang. Het doel van het maken van deze rangschikking is het vinden van de meest bepalende
genen, die een grote rol spelen bij bepaalde genetische ziektes. Het doel van het onderzoek is begrip
krijgen van hoe genen zijn verbonden, en het verbeteren van het al bestaande GeneRank algoritme, een
algoritme over het rangschikken van genen.
Het GeneRank algoritme dat het fundament is van dit onderzoek, neemt twee componenten van data mee
in de berekening van de rangschikking, zodat de rangschikking zo betrouwbaar mogelijk is. Het eerste
component is de experimentele data, die een rangschikking van genen geeft met eventuele fouten door
experimentele onnauwkeurigheid. Het tweede component is de biologische kennis over verbindingen, ofwel
de verwachte verbindingen, tussen genen, waardoor de experimentele onnauwkeurigheid zo veel mogelijk
wordt ingeperkt. De biologische data is opgebouwd uit de moleculaire functies, de biologische processen,
en de cellulaire componenten van de genen. In de beschreven algoritmes is er even veel gewicht gelegd op
de experimentele data als op de biologische data.
Het rangschikken van genen heeft veel overeenkomsten met het rangschikken van pagina’s, zoals het
PageRank algoritme van Google doet. Zowel genen als pagina’s zijn onderling verbonden, en zowel genen
als pagina’s zijn te modelleren in een netwerk, en te rangschikken op mate van belang. Het grootste verschil
tussen deze twee problemen is dat als een willekeurige pagina op het web met een andere willekeurige
pagina verbonden is, dit andersom niet hoeft te gelden, terwijl wanneer een willekeurig gen met een
ander willekeurig gen verbonden is, dit ook andersom geldt, ofwel het gen netwerk is symmetrisch. Door
de symmetrie van het GeneRank probleem is er na omschrijvingen een snel convergerende methode voor
symmetrische matrices toegepast, die een goede oplossing geeft. We hebben het GeneRank algoritme
hiermee verbeterd.
4
1
Voorwoord
In het rapport dat voor u ligt, leest u over het wiskundig rangschikken van genen. Het project dat dit
rapport tot stand heeft gebracht, is mijn bachelorproject. Ik heb dit mogen doen in de sectie Numerieke
Wiskunde aan de Technische Universiteit in Delft. Het project is uitgevoerd ter afsluiting van de bachelor
Technische Wiskunde die ik heb mogen volgen in de drie collegejaren van 2011-2014. In ongeveer drie
maanden heb ik geprobeerd om de opgedane kennis in het dik twee jaar studeren van wiskunde, en het
leren onderzoeken, samen te brengen in een eindrapport.
Het rapport is interessant voor iedereen die ge¨ınteresseerd is in het onderzoek naar het rangschikken van
genen door het uitvoeren van geschikte algoritmes. Voor iedereen die het PageRank probleem van Google
kent, is het rapport mogelijk interessant, omdat het GeneRank probleem veel overeenkomsten heeft met
het PageRank probleem, en ook op enkele punten verschilt.
Mijn begeleider Martin van Gijzen wil ik bedanken voor alle uren die hij in mijn ontwikkeling en de ontwikkeling van het rapport heeft gestoken. Zowel voor de vele interessante uren gedurende onze afspraken,
als voor het nalezen van verschillende versies van het rapport. Ik kan met zekerheid zeggen dat in het
rapport na iedere iteratie, na iedere versie die Martin nakeek, er een kleinere fout overbleef, zoals ook de
algoritmes in het rapport als kenmerk hebben.
Ook een woord van dank aan de commissieleden die dit rapport beoordelen: Meneer Vuik, Mevrouw
Keijzer en Meneer Kraaijkamp, bedankt!
Ik hoop dat u het rapport met plezier zult lezen, en dat u interessante kennis mag opdoen.
Cindy Boon Delft, 25 juli 2014.
6
Inhoudsopgave
1 Voorwoord
6
2 Inleiding
10
3 De verbinding tussen genen
3.1 Het rangschikken van genen . . . . . . . . .
3.2 De verbinding tussen genen, deel 1 . . . . .
3.2.1 Ontologie en annotatie . . . . . . . .
3.2.2 Componenten binnen een annotatie
3.2.3 Van GO annotatie tot GO matrix .
3.3 De verbinding tussen genen, deel 2 . . . . .
3.3.1 Betekenis experimentele vector ex .
3.3.2 Gen microarray experimenten . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
12
12
13
13
14
15
16
16
16
4 De GeneRank matrix
4.1 Transitiematrix T1 . . . . . . . . . . . . .
4.1.1 Opbouw met connectiematrix C .
4.1.2 Van connectie- tot transitiematrix
4.2 Transitiematrix T2 . . . . . . . . . . . . .
4.3 Gehele transitiematrix T . . . . . . . . . .
4.3.1 Weergave in graaf- en matrixvorm
4.4 Dempingsfactor d . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
18
18
18
19
20
20
21
21
5 Algoritmes om de GeneRank vector te berekenen
5.1 De stelselmethode . . . . . . . . . . . . . . . . . . . . .
5.1.1 Beschrijving methode . . . . . . . . . . . . . . .
5.1.2 Implementatie . . . . . . . . . . . . . . . . . . .
5.2 Eigenwaardeprobleem . . . . . . . . . . . . . . . . . . .
5.2.1 Markov matrix . . . . . . . . . . . . . . . . . . .
5.2.2 De stelling van Perron-Frobenius . . . . . . . . .
5.2.3 Eigenwaardes en stationaire rangschikkingsvector
5.3 De machtsmethode . . . . . . . . . . . . . . . . . . . . .
5.3.1 Beschrijving methode . . . . . . . . . . . . . . .
5.3.2 Implementatie . . . . . . . . . . . . . . . . . . .
5.4 De verschoven machtsmethode . . . . . . . . . . . . . .
5.4.1 Beschrijving methode . . . . . . . . . . . . . . .
5.4.2 Implementatie . . . . . . . . . . . . . . . . . . .
5.5 De methode van Jacobi . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
r
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
22
22
22
23
24
24
24
25
26
26
27
28
28
28
29
6 Het gebruik van symmetrie
6.1 Aanpassing systeem tot symmetrisch computatieprobleem . . .
6.2 De CG methode . . . . . . . . . . . . . . . . . . . . . . . . . .
6.2.1 Voorbereidende kennis . . . . . . . . . . . . . . . . . . .
6.2.2 Beschrijving CG methode . . . . . . . . . . . . . . . . .
6.2.3 Implementatie . . . . . . . . . . . . . . . . . . . . . . .
6.3 De PCG methode toegepast op het GeneRank probleem . . . .
6.3.1 Beschrijving PCG toegepast op het GeneRank probleem
6.3.2 Implementatie . . . . . . . . . . . . . . . . . . . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
30
30
30
30
31
31
32
32
32
.
.
.
.
.
.
.
8
6.4
Begrenzing eigenwaardes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
7 Numerieke experimenten
7.1 Dataset met recente data . . . . . . . . . . . . . .
7.2 Convergentie in de machtsmethode . . . . . . . . .
7.2.1 Het residu . . . . . . . . . . . . . . . . . . .
7.2.2 Implementatie van het residu . . . . . . . .
7.2.3 Convergentiefactor machtsmethode . . . . .
7.3 Convergentie in de verschoven machtsmethode . .
7.3.1 Het residu . . . . . . . . . . . . . . . . . . .
7.3.2 Implementatie van het residu . . . . . . . .
7.3.3 Resultaten . . . . . . . . . . . . . . . . . .
7.4 Convergentie in de methode van Jacobi . . . . . .
7.5 Convergentie in de PCG methode toegepast op het
7.5.1 Het residu . . . . . . . . . . . . . . . . . . .
7.5.2 Implementatie van het residu . . . . . . . .
7.5.3 Convergentiefactor PCG methode . . . . .
7.6 Vergelijken resultaten machts- en PCG methode .
7.6.1 Variatie dempingsfactor d . . . . . . . . . .
7.6.2 Tijdsduur en aantal iteraties . . . . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
GeneRank
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
probleem
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
33
35
35
35
35
35
36
37
37
37
38
40
40
40
41
41
42
42
46
8 Conclusie
48
9 Discussie
49
10 Appendix
10.1 Methodes . . . . . . . . . . . . . . . . . . . . . . . . . . . .
10.1.1 Stelselmethode . . . . . . . . . . . . . . . . . . . . .
10.1.2 Machtsmethode . . . . . . . . . . . . . . . . . . . . .
10.1.3 Verschoven machtsmethode . . . . . . . . . . . . . .
10.1.4 PCG methode toegepast op het GeneRank probleem
10.2 Gen Ontologie (GO) databank . . . . . . . . . . . . . . . .
10.3 Constructie matrix met meest recente data . . . . . . . . .
10.3.1 Inlezen . . . . . . . . . . . . . . . . . . . . . . . . . .
10.3.2 “Deel van”-relatie ⇒ “is een”-relatie . . . . . . . . .
10.3.3 Symmetrisch maken van de matrix . . . . . . . . . .
10.4 Inlezen willekeurige datasets van de GO website . . . . . . .
10.4.1 Constructie matrix uit gedownloade informatie . . .
9
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
51
51
51
51
52
53
54
55
55
55
56
56
56
2
Inleiding
Het menselijk lichaam bestaat uit een groot aantal organen. Deze organen zijn opgebouwd uit verschillende
weefsels, die weer bestaan uit miljarden cellen. In iedere cel zit een kopie van erfelijkheidsmateriaal, in
de vorm van chromosomen. En op de chromosomen zitten de genen. Er zijn 46 chromosomen, en 20.00025.000 genen in het menselijk lichaam.
Het rangschikken van deze 20.000-25.000 genen op mate van belang, is een hulpmiddel in de analyse van
gen expressie datasets in de moleculaire biologie. Biomedische onderzoekers probeerden al signalen die
uit de gen microarray experimenten komen, te verbinden met onderliggende genetische factoren van een
ziekte. De GeneRank methode gaat verder in op het mogelijk maken van een goede rangschikking van
genen. Het idee achter de GeneRank methode is om netwerk data, ofwel biologische kennis, te gebruiken
om een rangschikking van genen te geven in een microarray experiment. De rangschikking van genen
is dan minder gevoelig voor variatie, veroorzaakt door experimentele bias, dat een onderzoek slechts
gebaseerd op fout-gevoelige experimenten, wel is. De GeneRank methode is een nogal nieuwe wiskundig
methode, dat in 2005 naar buiten kwam naar het idee van Morrison et al., zie [6]. In essentie is GeneRank
een intu¨ıtieve aanpassing van de PageRank methode van Google, waarbij vele wiskundige eigenschappen
zijn behouden. In dit rapport zal enkele keren de vergelijking worden getrokken met Google PageRank,
ten behoeve van het verhelderen van het GeneRank principe. De eerste vergelijking volgt nu meteen.
Net zoals pagina’s op het web, zijn genen weer te geven in een netwerk, met genen als knopen en verbindingen tussen de genen als paden. De eerste doelstelling van dit rapport is het verhelderen van hoe
genen met elkaar in verbinding staan. Genen staan op verschillende manieren met elkaar in verbinding,
bijvoorbeeld het gezamenlijk maken van bepaalde stoffen, prote¨ınes, om cellen aan te sturen. Een ander
voorbeeld van een verbindende factor, is dezelfde functie hebben, of bepaalde componenten gemeenschappelijk hebben. De tweede doelstelling is het verbeteren van de GeneRank methode door het bekijken van
verschillende methodes toegepast op het GeneRank probleem; het vinden van snel convergerende methodes. Er zal worden gekeken naar de afname van de fout per iteratie, en daarmee de convergentiesnelheid
van de methodes. De laatste doelstelling is te kijken naar wat de voordelen zijn van de symmetrische
connectiematrix, een eigenschap van de GeneRank methode. Deze eigenschap maakt dat we de geconjugeerde gradi¨enten methode kunnen toepassen. Ook zorgt deze eigenschap ervoor dat er gemakkelijk
uitspraken kunnen worden gedaan over de ligging van de eigenwaardes.
De onderzoeksvragen die achtereenvolgens zullen worden behandeld zijn:
1. Wat houdt de verbinding tussen genen in?
2. Welke effici¨ente en snelle methodes zijn er om de rangschikkingsvector voor het GeneRank probleem
te vinden?
3. Zijn er voordelen voor het vinden van de rangschikkingsvector door de symmetrie van de connectiematrix?
Het volgende hoofdstuk begint met het behandelen van de eerste deelvraag, over de verbinding tussen
genen en de data die in dit rapport daarvoor zullen worden gebruikt. Vervolgens zal in hoofdstuk 3
worden uitgelegd hoe de transitiematrix, waarop het GeneRank algoritme is toegepast, en verscheidene
methodes zullen worden toegepast, is opgebouwd. Dan, wanneer het duidelijk is hoe de matrix waarmee zal
worden gewerkt, er uit ziet, kan in hoofdstuk 4 op de tweede deelvraag worden ingegaan. Er zullen enkele
methodes worden uitgelegd. De methodes die worden behandeld zijn: de stelselmethode, ofwel een directe
methode, de (verschoven) machtsmethode, en de methode van Jacobi. In hoofdstuk 5 wordt de derde
deelvraag behandeld, over de symmetrie van de connectiematrix. Hier zal onder andere de laatste methode,
de geconjugeerde gradi¨entenmethode met preconditionering, worden uitgelegd. Het laatste hoofdstuk,
10
hoofdstuk 6, bestaat uit de numerieke experimenten. De data zal worden beschreven, en de resultaten
van de behandelde methodes zijn weergegeven, op het gebied van tijdsduur van het algoritme, het aantal
iteraties, en convergentie, om na te gaan welke van de methode het meest effici¨ent is.
Ten slotte, is er de conclusie, de discussie, waarin suggesties worden gegeven voor mogelijk vervolg onderzoek, de referenties, waarin de geraadpleegde literatuur is opgenomen, en de appendix, waarin de matlab
codes zijn weergegeven die gebruikt zijn om dit rapport tot stand te brengen.
11
3
De verbinding tussen genen
In dit hoofdstuk zal basis informatie worden gegeven over genen en netwerken, zodat het duidelijk is wat
het doel is van het rangschikken van genen. Achtereenvolgens wordt er uitleg gegeven over wat genen
zijn en doen, wordt de verbinding tussen genen toegelicht, en wordt de gen ontologie (GO) database
ge¨ıntroduceerd.
3.1
Het rangschikken van genen
Ons genenbestand bevat erfelijk materiaal opgeslagen in lange DNA strengen. De genen bestaan uit
stukken DNA, zie figuur 1 verkregen uit [1]. Alle genen samen bepalen het functioneren van de cellen
waaruit het organisme is opgebouwd. De genen vertellen door middel van een code aan elk van de cellen
wat en wanneer iets te doen. Elke cel bevat exact dezelfde genen, maar lang niet alle genen zijn in alle
cellen actief. Elk type cel heeft een eigen patroon van genexpressie. In de ene cel komen genen tot expressie
die de cel tot een zenuwcel maken; in een spiercel zijn die genen stil en zijn weer hele andere genen actief.
Welke genen wel, en welke genen niet tot expressie komen, wordt geregeld door een zeer complex systeem.
Figuur 1: Schematische afbeelding van het DNA en de genen
Genen zijn op verschillende manieren met elkaar verbonden. Er is een mate van belang van genen die te
bepalen is door te kijken naar de verbindingen tussen de genen. Het doel is om de meest bepalende genen
te vinden die een rol spelen bij bepaalde genetische ziektes. In de gerangschikte vector staat bovenaan
het gen dat het “meest belangrijk” is, dus het meest in contact staat met andere genen, en onderaan het
gen dat het “minst belangrijk” is, ofwel het minst in contact staat met andere genen.
Vele bestaande methodes zijn empirisch, met een mogelijk grote bias, aangezien deze onderzoeken slechts
op empirische experimenten zijn gebaseerd, en daarmee meer onderhevig aan variatie veroorzaakt door
experimentele afwijkingen. Om een kleinere bias te verkrijgen, wordt er gemodelleerd waarbij gekeken
wordt naar de verwachte verbindingen tussen genen. Er wordt gebruik gemaakt van grafentheorie, waarbij
12
er twee dominante componenten van informatie uiteenvallen bij het algoritme: de experimentele informatie
van de DNA-microarray data, en de kennis over het onderliggende netwerk van genen; de biologie. Het
GeneRank algoritme combineert deze twee componenten.
De genen zijn weer te geven in een ongerichte graaf van een eindige verzameling knopen. Elk gen is een
knoop. Elke verbinding tussen twee genen is een pad; bepaalde vooraf bekende informatie.
In de volgende paragraaf zal nader worden ingegaan op de twee onderscheide componenten binnen de
transitiematrix:
1. Deel 1: de kennis over de biologie van genen; de kans dat een gen met een ander gen in contact
staat, bijvoorbeeld rekening houdend met functie of prote¨ıne constructie;
2. Deel 2: de informatie van de DNA-microarray data.
3.2
De verbinding tussen genen, deel 1
Het eerste deel van de transitiematrix is te vergelijken met de hyperlinks, de verbindingen tussen pagina’s,
bij Google PageRank. De data die de verbinding tussen genen inhoudt, kan verschillend zijn en op
verschillende plekken worden gevonden, zoals:
• Rat Genoom Databank (RGD);
• Biomedische ontologie¨en, van de volwassen hersenen van de muis tot zebravis anatomie en ontwikkeling, bijvoorbeeld te vinden in [2];
• Gen Ontologie (GO) annotatie: gemeenschappelijke factoren op het gebied van biologisch proces,
cellulair component, en moleculaire functie.
In dit rapport is het laatst genoemde, Gene Ontology (GO), te vinden in [10], beschouwd als databank
van verbindingen tussen genen. GO interacties beschrijven op moleculair niveau of een paar van genen
in verbinding staan met elkaar.
3.2.1
Ontologie en annotatie
Een ontologie is een gecontroleerde woordenschat of -lijst, van goed gedefinieerde termen met gespecificeerde relaties tussen hen, en is in staat ge¨ınterpreteerd te worden door zowel mensen als computers.
Annoteren is biologische informatie voorzien van een geschikte context. Specifiek voor DNA-annotatie
wordt ontdekt genetisch materiaal bijvoorbeeld zowel gerelateerd aan de positie van het gen in het genoom,
en aan de intronen en exonen, ofwel de stukjes DNA die zich bevinden in het gen, maar respectievelijk
wel en niet worden gebruikt om prote¨ıne te coderen, met andere woorden de delen van het gen die respectievelijk wel en niet in het uiteindelijke mRNA terechtkomen. Als ook aan de regulerende sequenties,
aan herhalingen, aan de naam en de genproducten, zoals RNA en prote¨ınes.
De annotaties worden opgeslagen in databanken om annoteren eenvormig te maken. GO is een belangrijke
ontologie binnen de bio-informatica als het gaat om het beschrijven van eigenschappen van de onderzochte
genen in een cel. Biologen verliezen vaak veel tijd en moeite aan het zoeken naar beschikbare informatie
over elk gebied van onderzoek. Er is ook een brede variatie aan terminologie door de jaren heen ontstaan,
dat effectief zoeken moeilijk maakt voor zowel computers als wetenschappers. GO vergemakkelijkt dit en
is dus een datastructuur door gezamenlijke inspanning, dat alle relevante entiteiten en hun onderlinge
relaties en regels binnen dat domein bevat. Het vervult de vraag naar consistente beschrijvingen van gen
producten in verschillende databases. Het doel ervan is om de weergave van de eigenschappen van genen
13
en de genproducten te standaardiseren, onafhankelijk van het soort of ras waaruit ze komen. Duidelijk
wordt nog vermeld op de GO website dat GO niet een manier is om biologische databanken te unificeren.
3.2.2
Componenten binnen een annotatie
Gen ontologie annotatie geeft een samenvatting van wat bekend is over een gen of gen product. De GO
databank uit [10] bevat drie componenten van informatie. Het GO project heeft deze drie gestructureerde
gecontroleerde ontologie¨en ontwikkeld, drie GO secties, die gen producten beschrijven. De ontologie¨en zijn
in termen van:
• de moleculaire functies van het gen (de elementaire activiteiten van een gen product op het moleculaire niveau, zoals binding, transportatie, als mutatie). Dit is anders dan biologische processen
die meer dan ´e´en activiteit inhouden);
• het biologisch proces waarin het gen is verwikkeld (operaties of verzamelingen van moleculaire gebeurtenissen met een gedefinieerd begin en einde, toebehorend aan het functioneren van ge¨ıntegreerde
levende eenheden, zoals cellen, weefsels, organen en organismes, bijvoorbeeld cellulaire processen);
• de locatie van het gen binnen de cel, ofwel het cellulaire component.
Samenvattend: GO is een platform om genen en gen producten te beschrijven in termen van hun moleculaire functies, hun biologische processen, en hun subcellulaire localisaties of celullaire componenten.
Binnen het GO project is de informatie over de genen soort-overschrijdend.
Elke ontologie-term kan worden gerepresenteerd door een unieke ontologie specifieke identificatie, een
GO-getal, ofwel een ontologie ID. De aangeleverde data binnen een ontologie, uit [15], ziet er per gen als
volgt uit:
[Term]
id: GO:0000001
name: mitochondrion inheritance
namespace: biological_process
def: "The distribution of mitochondria, including
the mitochondrial genome,
into daughter cells after mitosis or meiosis,
mediated by interactions between mitochondria
and the cytoskeleton." [GOC:mcc, PMID:10873824, PMID:11389764]
exact_synonym: "mitochondrial inheritance" []
is_a: GO:0048308 ! organelle inheritance
is_a: GO:0048311 ! mitochondrion distribution
De relatie “is a” beschrijft dat dit gen, nummer 0000001, verbonden is met gen nummer 0048308 en
0048311.
Enkele componenten die buiten het bereik van de GO databank vallen:
• Processen, functies of componenten die uniek zijn voor mutanten of ziektes;
• prote¨ıne-prote¨ıne interactie;
• omgeving, evolutie en expressie;
14
• anatomische of histologische kenmerken boven het niveau van cellulaire componenten, inclusief cel
types.
In figuur 2 is te zien hoe gen nummer 0065010, de ‘extracellular membrane-bounded organelle’, zich
verhoudt tot voorvaderen en afstammelingenvan het gen.
Figuur 2: Familie van het gen: ‘extracellular membrane-bounded organelle’.
Een grafisch voorbeeld van het netwerk van gen nummer 0005840, ‘ribosome’, is te zien in figuur 3.
Figuur 3: Netwerk van het gen: ‘ribosome’.
3.2.3
Van GO annotatie tot GO matrix
De GO matrix is een gerichte acyclische graaf. Elke term heeft relaties met ´e´en of meerdere andere
termen. Genen zijn verbonden wanneer ze een annotatie delen, gedefinieerd door de GO databank. GO
definieert drie netwerken, ´e´en voor elk van de drie eerder genoemde GO secties: moleculaire functie,
biologische proces en cellulaire component. In het GeneRank algoritme wordt een variant gebruikt van
15
de GO matrix, zoals in [6]. De acyclische gerichte graaf wordt aangepast. In de GO matrix staan 0’en,
1’en en 2’en. Daar waar een 0 staat, betekent het dat twee genen geen door GO benoemde verbinding
hebben. Daar waar een 1 staat, zijn de drie door GO genoemde GO secties hetzelfde voor een paar van
genen. Daar waar een 2 staat, is het ene gen gedeeltelijk gelijk aan het andere gen wat betreft de drie
GO secties. Daar waar er 2’en staan in de GO matrix, worden bladknopen, ofwel knopen zonder kind,
toegewezen als annotaties voor elk gen, zodat de connectiematrix slechts 0’en en 1’en bevat. Verder wordt
ervoor gezorgd dat de matrix symmetrisch is, door te zorgen dat daar waar een 1 staat in de aangepast
GO matrix op plek cij , ook een 1 staat in de getransponeerde aangepaste GO matrix.
We noemen de symmetrische aangepaste GO matrix: connectiematrix C.
3.3
De verbinding tussen genen, deel 2
Het tweede deel van de transitiematrix is opgebouwd uit experimentele empirische resultaten. Om precies
te zijn: gen microarray expressie data. Gebruik makend van gen microarrays, is het mogelijk het gen
expressie profiel uit meer dan 30.000 genen t verkrijgen. Nu eerst informatie over wat de experimentele
vector inhoudt, daarna informatie over hoe de experimentele vector gevonden wordt.
3.3.1
Betekenis experimentele vector ex
De experimentele vector ex = (exi ), met i = 1, 2, · · · , n, bevat de absolute waarde van de expressieverandering van alle genen i, met i = 1, 2, · · · , n. De expressieverandering is een maat van “relevantie” van dat
gen verkregen door microarray experimenten. De experimentele vector is de gen expressie karakterisering.
Deze vector is te vergelijken met de gepersonaliseerde PageRank vector, waar door teleportatie in het
random walk proces wordt gestuurd in de richting van de geliefde pagina’s van een gebruiker. Zoals in
artikel [6], pag. 2, beschreven is, heeft het PageRank algoritme een random walk interpretatie waar de
rangen corresponderen met de invariante maat van een teleporterende random walk op het web. Dit is
equivalent met dat de rang van een webpagina proportioneel is met de bestede tijd op een webpagina al
surfend op het web. Nu de vergelijking naar de GeneRank methode. In de GeneRank methode wordt ook
gestuurd, maar dan volgens expressieniveau.
In het veld van moleculaire biologie, is genexpressiekarakterisering de maat van de activiteit, de expressie
van duizenden genen op hetzelfde moment, om een globaal beeld te cree¨eren van de celullaire functie.
Vele experimenten van deze soort meten en geheel genoom tegelijkertijd, ofwel, ieder gen aanwezig in een
bepaalde cel. Van de paar meest significante genen wordt verwacht dat ze worden ge¨ıdentificeerd door de
experimentele data.
3.3.2
Gen microarray experimenten
De experimentele vector wordt gevonden door experimenten te doen met gen microarray expressies. Een
gen microarray, zoals beschreven in [4] bestaat uit een groot aantal bekende DNA controle rijen die zich
bevinden op verschillende locaties op een dia. De DNA in een cel verandert niet, maar zekere afgeleides
van de DNA, de mRNA en de tRNA, worden geproduceerd als stimulatie plaatsvindt. De messenger RNA
(mRNA), ofwel ‘boodschapper’ RNA, speelt een centrale rol in het tot expressie brengen van genetische
informatie. mRNA is een vorm van RNA, welke als ‘boodschapper’ twee processen met elkaar verbindt.
Een impressie van dit proces waarin mRNA een rol speelt is: begonnen wordt met DNA, dan is er
transcriptie, waarbij een stuk DNA, een gen, overgeschreven wordt tot mRNA, en dan de translatie,
16
waarbij het mRNA wordt vertaald naar een keten van aminozuren, een prote¨ıne, ofwel een eiwit. Ieder
gen is een recept hoe een bepaald prote¨ıne te maken. De prote¨ınes geven aan wat de cel behoort te doen,
en bepalen hiermee hoe het organisme vorm wordt gegeven.
Het proces ziet er globaal gezien uit als in figuur 4 uit [5].
Figuur 4: Gen expressie
Samenvattend, is genexpressie het proces waarbij informatie in een gen ‘tot expressie komt’, doordat het
gen afgelezen wordt en RNA en prote¨ınes worden gemaakt. Deze experimenten leveren geen data op het
absolute niveau van een specifiek gen, ofwel de werkelijke concentratie mRNA, maar zijn te gebruiken
om het expressieniveau tussen condities en genen te vergelijken. Microarrays kunnen gelijktijdig toezien
op de expressieniveaus van duizenden genen.
Het principe achter het aangeven van expressieniveau, zoals beschreven in [3] is dat de hoeveelheid van
een bepaalde stof, fluorescentie, die wordt gemeten op iedere rij-specifieke locatie direct proportioneel
is met de hoeveelheid mRNA aanwezig in de geanalyseerde steekproef. Door deze metingen bevat de
experimentele vector die wordt gebruikt in de GeneRank methode de meest actieve, belangrijke genen
bovenaan, en de minst actieve, belangrijke genen, onderaan. Meer hierover staat geschreven in [3].
17
4
De GeneRank matrix
In dit hoofdstuk zal de matrix waarmee in de hierna volgende methodes zal worden gewerkt, opgebouwd,
en wordt er uitleg gegeven over de dominante eigenwaarde en eigenvector.
De matrix waarop algoritmes zullen worden toegepast, is een transitiematrix. Een matrix die kans bevat
dat een gen met een ander gen verbonden is. Zoals toegelicht in het vorige hoofdstuk, bevat de transitiematrix, die in dit rapport zal worden gebruikt, twee soorten data, die worden gerepresenteerd door twee
transitiematrices T1 en T2 .
1. Verbindingsmatrix T1 ; bevat het logisch verwachte verband tussen genen, in ons geval de annotaties
in de GO matrix.
2. Teleportatiematrix T2 ; bevat het verband tussen genen met betrekking tot de experimentele data.
4.1
Transitiematrix T1
De matrix T1 die nu zal worden ge¨ıntroduceerd modelleert een netwerk met verbindingen tussen genen.
Deze verbindingen kunnen worden weergegeven door een graaf. In de uitleg die volgt wordt gewerkt met
een voorbeeldgraaf en een bijbehorende voorbeeldconnectie- en transitiematrix.
4.1.1
Opbouw met connectiematrix C
De transitiematrix wordt opgebouwd met kennis van de connectiematrix C. C is een symmetrische matrix
met een 1 waar een gen verbonden is met een ander gen door GO-annotatie, en een 0 waar er geen
verbinding is. Een voorbeeld van een connectiematirx met 5 genen is:








C=







0
1
1
1
1
1
0
0
1
0
1
0
0
0
0
1
1
0
0
1
1
0
0
1
0
















In de bijbehorende voorbeeldgraaf in figuur 5 staan de pijlen voor een link, ofwel het getal 1 in de matrix.
Merk op dat de graaf ongericht is, en de matrix daarmee symmetrisch. Dit is het geval, aangezien het
geldt voor i 6= j, dat wanneer gen i met gen j in verbinding staat, het ook geldt dat gen j met gen i in
verbinding staat. Dit in tegenstelling tot hoe de Google matrix is opgebouwd. Daar namelijk, wanneer
pagina i naar pagina j verwijst door middel van een hyperlink, ofwel er is een gerichte verbinding van
pagina i naar pagina j, dan hoeft pagina j niet naar pagina i te verwijzen.
18
Figuur 5: Graaf bij connectiematrix C
4.1.2
Van connectie- tot transitiematrix
Van de connectiematrix C kan een transitiematix T1 worden gemaakt, die nodig is voor de methodes
die zullen worden gebruikt om de rangschikkingsvector te vinden. De stap van connectiematrix C →
transitiematrix T1 gaat als volgt: de transitiematrix T1 is gedefinieerd door T1 = CD−1 , met C = C T en
diagonaalmatrix D gedefinieerd door: D = diag(d1 , . . . , dn ), met di = deg(i) de som van de getallen in
rij i; de graad. In dit geval zijn het enen die worden opgeteld.
De zojuist geconstrueerde transitiematrix T1 bevat de volgende informatie: laat pij (t) de kans aanduiden
dat de keten in toestand Sj is op tijd t, gegeven dat het in toestand Si was op tijd t − 1. Ofwel,
pij (t) = P (X(t) = Sj |X(t − 1) = Si ),
met pij (t) is de transitiekans. De transitiematrix T1 is dan T1 (t) = [pij (t)], een n × n matrix, met n het
aantal toestanden, en elke kolom van T1 (t) telt op tot 1. Het laatste maakt dat T1 een kolom-stochastische
matrix is.
De matrix representatie is nu als volgt:
• Laat n = 5 het totaal aantal genen;
• Laat k de uitgaande verbindingen van gen j met j = 1, . . . , 5;
• Laat T1 de matrix zijn met elementen:
T1ij =
1
k
als gen i in contact staat met gen j;
0 anders.
19
De transitiematrix T1 behorend bij de graaf in figuur 5, en daarmee behorend bij de matrix C, is:

T1 = C ∗ D−1
4.2







=







0
1
1
1
1
1
0
0
1
0
1
0
0
0
0
1
1
0
0
1
1
0
0
1
0
 
 
 
 
 
 
 
 
∗
 
 
 
 
 
 
 
4
0
0
0
0
0
2
0
0
0
0
0
1
0
0
0
0
0
3
0
0
0
0
0
2
−1























=







0
1
2
1
1
3
1
4
0
0
1
3
1
4
0
0
0
1
4
1
2
0
0
1
4
0
0
1
3
1
2



0 



0 



1 
2 


0 
Transitiematrix T2
De transitiematrix T2 is opgebouwd uit de experimentele vector ex, en ziet er als volgt uit: T2 =
met vectoren van lengte n, ofwel:







ex
T2 =
∗ eT = 

||ex||1





ex1
||ex||1
ex1
||ex||1
ex2
||ex||1
···
ex1
||ex||1
ex2
||ex||1
···
ex2
||ex||1
..
.
..
.
..
.
exn
||ex||1
exn
||ex||1
···
ex
T
||ex||1 ∗e ,







.. 

. 



exn

||ex||1
Omdat de transitiematrices verbindingen tussen genen onderling representeren, zouden we verwachten
dat er slechts 0’en op de diagonaal staan. Op de diagonaal in T2 staan echter geen 0’en. Wanneer we
kijken naar een matrix die in dit werk wordt gebruikt als data, zien we dat deze matrix C 4047 genen
representeert. Op de diagonaal van T2 staat daarom een nogal klein getal, namelijk 1/4047 = 0.00024710.
4.3
Gehele transitiematrix T
De twee zo juist verkregen transitiematrices opgeteld, geeft de transitiematrix T waarmee zal worden
gewerkt in het volgende hoofdstuk. De graaf bij transitiematrix T is nu een graaf waarbij alle knopen
met elkaar zijn verbonden, zoals te zien in de graaf in figuur 6.
20
Figuur 6: Schematische afbeelding dna
4.3.1
Weergave in graaf- en matrixvorm
Om T1 en T2 op te kunnen tellen tot T wordt de scalar 0 < d < 1 ge¨ıntroduceerd die zorgt dat het geheel
convex blijft. Een voorbeeld van een gehele transitiematrix T met in totaal 5 genen:

0


 0


 1

2
T = dT1 + (1 − d)T2 = d 

 1

 3

 1

6
0
1
2
1
3
0
1
6
1
2
1
6
0
0
1
2
0
0
1
3
1
3
1
6
1
6





1 

3 





1 

3  + (1 − d) 




1 

6 






0
ex1
||ex||1
ex1
||ex||1
ex2
||ex||1
···
ex1
||ex||1
ex2
||ex||1
···
ex2
||ex||1
..
.
..
.
..
.
exn
||ex||1
exn
||ex||1
···







.. 

. 



exn

||ex||1
Hoe dichter d bij 0 ligt, des te meer benadrukt het model de experimentele data. Hoe dichter d bij 1 ligt,
des te meer benadrukt het model de netwerk data van de GO databank.
4.4
Dempingsfactor d
De keuze van de dempingsfactor d be¨ınvloed de prestatie van gen rangschikking. In artikel[6], pag. 11,
wordt na berekeningen op een waarde van d = 0.5 uitgekomen, anders dan bij de PageRank methode,
waar deze waarde rond de d = 0.85 ligt. De optimale keuze van d is een interessant onderwerp en verdient
nader aandacht te krijgen.
21
5
Algoritmes om de GeneRank vector te berekenen
In dit hoofdstuk over methodes worden verschillende algoritmes bekeken om een rangschikkingsvector te
berekenen. Er worden achtereenvolgens twee soorten methodes bekeken:
• een directe methode;
• iteratieve methodes.
De directe methode is de stelselmethode, de iteratieve methodes zijn: de machtsmethode en de methode
van Jacobi.
Morrison et al. gebruikten in [6] een directe methode om de GeneRank vector te verkrijgen, een methode
die erg ineffici¨ent kan zijn wanneer het probleem heel groot. Voor grote ijle systemen is Gaussische
eliminatie en acherwaartse substitutie niet geschikt. Grote ijle systemen komen vaak op als parti¨ele
differentiaal vergelijkingen of optimalisatieproblemen worden opgelost. Wel geschikt, is een benaderde
oplossing met behulp van iteratieve methodes. Een iteratieve methode is een wiskunde procedure die
een rij van verbeterde benaderingsoplossingen genereert voor een klasse van problemen. Op de directe
methode na, zullen er dus iteratieve methodes worden besproken.
Om een rangschikkingsvector te verkrijgen, moet de methode convergeren. Wanneer er een oplossing van
een systeem van vergelijkingen wordt gezocht, dan gebruikt de iteratieve methode een beginschatting
om een succesvolle benadering voor de oplossing te genereren. Een iteratieve methode heet convergent
als de corresponderende rij voor een gegeven beginbenadering convergeert. Iteratieve methodes zijn vaak
de enige keuze voor niet-lineaire vergelijkingen. Echter, iteratieve methodes zijn vaak nuttig om lineaire
problemen op te lossen die een groot aantal variabelen bevatten, waar directe methodes erg duur zouden
kunnen zijn, soms zelfs onmogelijk, ook met de best beschikbare berekenkracht.
Eerst zal de directe methode worden behandeld, vervolgens zal er algemene informatie volgen over een
eigenwaardeprobleem, dan wordt met (verschoven) machtsmethode behandeld, en als laatste de methode
van Jacobi, die een lineair stelsel oplost. De optimale verschuiving van de verschoven machtsmethode
zal worden berekend, en er zal worden geconcludeerd dat de machtsmethode en de methode van Jacobi
hetzelfde zijn.
5.1
De stelselmethode
5.1.1
Beschrijving methode
Een directe methodes, in dit geval de stelselmethode, probeert het gegeven probleem door een eindige
rij van operaties op te lossen. Wanneer er geen afrondfouten zijn, dan levert een directe methode een
exacte oplossing, zoals wanneer het een lineair systeem van vergelijkingen Ax = b oplost met behulp van
Gaussische eliminatie.
Een directe methode geeft een exacte oplossing, maar is slechts toe te passen op voldoende kleine problemen. De methode kan erg veel geheugenruimte gebruiken, zoals het geval is bij het GeneRank probleem.
Voor grote n, kost de stelselmethode veel tijd. Het oplossen van een stelsel met n vergelijkingen en n
onbekenden kost veel geheugen.
We beginnen met het systeem T r = r met T = dT1 + (1 − d)ex ∗ eT .
22
Enkele omschrijvingen van het systeem geeft:
Tr = r
[dT1 + (1 − d)ex ∗ eT ]r = r
dT1 r + (1 − d)T ex ∗ eT r = r
dT1 r + (1 − d)ex = r
(I − dT1 )r = (1 − d)ex,
waarbij geldt dat eT r = 1, aangezien r een kansvector is.
De vergelijking :
(I − dT1 )r = (1 − d)ex
komt overeen met het systeem:
(I − dCD−1 )r = (1 − d)ex,
(1)
met D de matrix met de som van de rijen van C, met D een diagonaalmatrix, die informatie bevat over
de graad van een knoop, ofwel het aantal paden dat grenst aan de knoop.
Dan volgt dat de rangschikkingsvector r met A = (I − dCD−1 ) en b = (1 − d)ex is:
r = A−1 b.
5.1.2
Implementatie
De beschreven stelselmethode komt uit artikel[6], pag. 13. In het volgende algoritme geldt dat ex is
de experimentele data, onder andere onderdeel van de beginvector, degrees is de som van de rijen van
matrix C, ind is een indicator, D is hier D−1 , ofwel de som van de rijen van C op de diagonaal, en r is
de rangschikkingsvector.
De invoer voor het algoritme:
• de matrix C die het netwerk van genen representeert;
• de experimentele vector ex, ofwel microarray data, ofwel de mate van expressieverandering per gen
in een vector;
• de parameter d ∈ [0, 1], die bepaalt hoeveel nadruk er op de experimentele data ligt en hoe veel op
de transitiematrix.
De uitvoer: rangschikkingsvector r.
Het algoritme is beschreven in algoritme 1.
23
Algorithm 1 Stelsel algoritme
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
procedure
ex = abs(ex)
norm ex = ex/max(ex)
c = sparse(C)
degrees = sum(C)
ind = find(degrees == 0)
degrees(ind) = 1
D1 = sparse(diag(1./degrees))
A = eye(size(c)) − d ∗ (c0 ∗ D1)
b = (1 − d) ∗ norm ex
r = A\b
5.2
Eigenwaardeprobleem
In de volgende paragraaf wordt de machtsmethode behandeld, een methode die een eigenwaardeprobleem
oplost. Een eigenwaardeprobleem is een mogelijkheid waarop het GeneRank probleem te beschrijven is.
Deze paragraaf zal informatie verschaffen over de eigenschappen van de matrix T , die nodig zijn om het
eigenwaardeprobleem op te lossen, de bijbehorende eigenwaardes, en de unieke stationaire eigenvector r.
De stelling van Perron-Frobenius speelt hierbij een grote rol.
5.2.1
Markov matrix
De transitie matrix T van genen heeft alle eigenschappen van een Markov matrix, een eigenschap die zo
meteen van pas komt. Een Markov matrix is een stochastische matrix met de eigenschap dat de matrix
geheugenloos is. Een stochastische matrix is een vierkante matrix T met de eigenschappen:
1. ti,j ≥ 0, voor alle i en j;
PN
2.
i=1 ti,j = 1, voor alle j.
Ofwel, alle elementen van T zijn niet-negatief en alle kolomsommen zijn kansvectoren.
Een Markov proces is een stochastisch proces met een eindige of aftelbare toestandsruimte die voldoet
aan de Markov eigenschap. Een collectie random variabelen X(t)t∈T met toestandsruimte S = S1 , S2 , . . .
voldoet aan de Markov eigenschap als, voor alle integers m en t0 < t1 < · · · < tm < t, geldt:
P (X(t) = Sj |X(tm ) = Sim , . . . , X(t0 ) = Si0 ) = P (X(t) = Sj |X(tm ) = Sim ).
Ofwel, een Markov keten is geheugenloos in de zin van dat de toestand van de keten in de volgende
tijdsperiode alleen afhangt van de huidige tijdsperiode, en niet van de afgelopen toestanden van de keten.
We zijn ge¨ınteresseerd in de discrete-tijd Markov keten met eindige toestandsruimte. Een stochastisch
proces Xn , n ≤ 0 met toestandsruimte S heet een discrete-tijd Markov keten (DTMC) als voor alle i en
j in S geldt:
P (Xn+1 = j|Xn = i, Xn−1 , . . . , X0 ) = P (Xn+1 = j|Xn = i).
In woorden: Gegeven het heden Xn en het verleden X0 , . . . , Xn−1 van het proces hangt de toekomst
Xn+1 alleen af van het heden en niet van het verleden. De conditionele kansen P (Xn+1 = j|Xn = i)
heten 1-staps overgangskansen.
5.2.2
De stelling van Perron-Frobenius
Om de stelling van Perron-Frobenius te kunnen defini¨eren, geven we enkele definities. De re¨ele transitiematrix T is niet-negatief als alle elementen van T niet-negatief zijn. Een niet-negatieve vierkante matrix
24
T heet primitief als er een k bestaat, zodat alle elementen van T k positief zijn. Er geldt dat als T primitief
is, dat alle andere eigenwaardes van T voldoen aan:
|λ| < λmax .
T heet irreducibel als ∀ i, j er een k = k(i, j) bestaat, zodat (T k )ij > 0. Ofwel T irreducibel als T = (tij )
een n × n positieve matrix: tij > 0 voor 1 ≤ i, j ≤ n. In woorden, de matrix T is irreducibel als elk
gen met elk ander gen in verbinding staat, ofwel het is mogelijk van elke willekeurige knoop naar een
willekeurige ander knoop te gaan.
De stelling van Perron-Frobenius, zoals beschreven in [19] heeft als aanname dat T irreducibel is. De
stelling geeft onder andere de volgende drie eigenschappen:
• T heeft een positieve re¨ele eigenwaarde λmax zodat alle andere eigenwaardes van T voldoen aan:
|λ| ≤ λmax ,
dus de spectraal radius ρ(T ) is gelijk aan λmax ;
• λmax heeft algebra¨ısche en geometrische multipliciteit 1, en heeft een eigenvector r = (r1 , r2 , . . . , rn )
van T , zodat alle elementen van r positief zijn: T r = λr, ri > 0 voor 1 ≤ i ≤ n. De bij r behorende
eigenwaarde is λmax ;
• Elke niet-negatieve eigenvector verschillend van r is een veelvoud van r.
Het geldt dat T irreducibel is, ofwel ieder gen dat gerepresenteerd wordt in T is verbonden met elk willekeurig ander gen in T . Dit geldt wegens de constructie van T , zoals in het vorige hoofdstuk. Nu, wegens T
irreducibel, mogen we de stelling van Perron-Frobenius toepassen, en gelden de zojuist opgesomde punten
voor de transitiematrix T van het GeneRank probleem.
5.2.3
Eigenwaardes en stationaire rangschikkingsvector r
Stelling: Voor transitiematrix T is eigenwaarde 1 de grootste eigenwaarde, en deze is uniek.
Bewijs: We bewijzen achtereenvolgens de existentie en uniciteit van de grootste eigenwaarde 1. Existentie:
Wegens
T kolom-stochastisch, geldt dat er een eigenwaarde 1 bestaat, ofwel het geldt dat matrix T met
P
T
=
1 met i = 1, · · · , N , voldoet aan T r = r. Ook zegt de stelling dat er geen eigenwaarde |λ| > 1
j i,j
is. De reden dat er geen |λ| > 1 is, is dat de machten van T k zouden groeien wanneer
|λ| > 1, maar
P
kenmerk van een Markov matrix is dat iedere T k weer een Markov matrix is, met j Ti,j = 1, dus volgt
dat |λ| ≤ 1. Hieruit kunnen we concluderen dat 1 een eigenwaarde, en dat 1 is de grootste eigenwaarde
van T .
Eenduidigheid: De Perron-Frobenius stelling verzekert dat T primitief is, ofwel dat als λmax = 1, dat
deze dan uniek is. Doordat matrix T slechts positieve elementen bevat, geldt wegens Perron-Frobenius dat er een stationaire
kansvector bestaat met eigenwaarde λ = 1. Deze stationaire kansvector is de corresponderende strikt
positieve eigenvector, de rangschikkingsvector. Het i’de element in de rangschikkingsvector, r(k) , is de
kans dat de keten op tijdstip k in toestand Si is, dat is,
(k)
ri
:= P (X(k) = Si ),
Pn
met r = (r1 , r2 , . . . , rn )T de kansverdelingsvector, waarvan de elementen optellen tot 1, ofwel, i=1 ri = 1.
Het i’de element van de unieke stationaire verdeling r is de proportie tijd dat de keten spendeert in
25
toestand Si op de lange termijn. Het Markov proces zal convergeren naar een stationaire vector r. Wanneer
de transitiekansen niet meer vari¨eren in de tijd, dan is de Markov keten stationair.
Het systeem dat wordt opgelost om deze stationaire rangschikkingsvector te verkrijgen is Tr=r, zodat:
rk+1 = T rk .
Verder geldt dat de op ´e´en na grootste eigenwaarde λ2 gelijk is aan de dempingsfactor d. Dit zal worden
bewezen in het hoofdstuk over symmetrie.
5.3
De machtsmethode
5.3.1
Beschrijving methode
De machtsmethode zoals beschreven in [7] begint met een beginschatting r0 en vormt de volgende termen
r1 = T r0 , r2 = T r1 en zo verder, zodat in het algemeen rk+1 = T rk . Iedere stap is een matrix-vector
multiplicatie, en na k stappen geeft de methode: rk = T k r0 .
Veronderstel nu dat T een verzameling van n lineair onafhankelijke eigenvectoren {x1 , x2 , . . . , xn } heeft,
waarvoor geldt: T xi = λi xi , i = 1, . . . , n en waarbij λ1 correspondeert met de grootste eigenwaarde, in
ons geval λ1 = 1, λ2 met de op ´e´en na grootste eigenwaarde, enzovoorts, ofwel:
1 = |λ1 | > |λ2 | ≥ |λ3 | ≥ · · · ≥ |λn |
met de grootste eigenwaarde, λ1 = 1, uniek. Deze eigenwaarde heeft algebra¨ısche multipliciteit gelijk aan
1.
De vereiste beginvector r0 ziet er als volgt uit:
r0 = c1 x1 + c2 x2 + · · · + cn xn
met ||r0 ||1 = 1.
Itereren van r0 met matrix T geeft dan:
T k r0 = c1 T k x1 + c2 T k x2 + · · · + cn T k xn
= c1 λk1 x1 + c2 λk2 x2 + · · · + cn λkn xn
= c1 x1 + c2 λk2 x2 + · · · + cn λkn xn
= rk
Dan zo lang als de veronderstelling r0 een component van de eigenvector x1 bevat, zodat c1 6= 0, en aan
de voorwaarde is voldaan dat de matrix T diagonaliseerbaar is en daarmee de eigenvectoren xi een basis
vormen voor Rn , zal een component geleidelijk aan dominant worden. Er geldt dan dat λrkk = rk met de
1
grootste eigenwaarde λmax = λ1 = λk1 = 1, is:
rk = c1 x1 + c2
λ2
λ1
k
x2 + · · · + cn
= c1 λk1 x1 + c2 λk2 x2 + · · · + cn λkn xn
26
λn
λ1
k
xn
Dan vinden we na k iteraties:
|λi |
= |λi | < 1 → lim |λi |k = 0 → rk ≈ c1 λk1 x1 = c1 x1
k→∞
|λ1 |
Ofwel met andere notatie: rk+1 = ||TTrrkk||1 . Wat gebeurt, is dat met iedere iteratie de vector rk wordt
vermenigvuldigd met de matrix T , en vervolgens genormaliseerd.
5.3.2
Implementatie
Het implementeren in matlab wordt op twee verschillende manieren gedaan:
• De matrix wordt voorafgaand aan het itereren gevormd;
• De matrix wordt impliciet gevormd.
De beginvector waarmee in beide gevallen wordt gerekend is: r0 =
ex
||ex||1 .
De eerste manier, waarbij de matrix voorafgaan aan het itereren wordt gevormd, pakt als transitiematrix
T:
T = d ∗ T1 + (1 − d) ∗ T2 = d ∗ T1 + (1 − d) ∗
ex
∗ ones(1, n)
||ex||1
met T1 en T2 de twee onderdelen van de transitiematrix T . Het algoritme is beschreven in algoritme 2,
met
Algorithm 2 Machtsmethode
procedure
loop:
rk+1 = T ∗ rk
normr = norm(rk+1 , 1)
rk+1 = rk+1 /normr
6: end
1:
2:
3:
4:
5:
De tweede manier, waarbij de matrix T impliciet wordt gevormd, vormt r op de volgende manier:
r = d ∗ (C ∗ (D\r0)) + (1 − d) ∗ norm ex
Het itereren gaat zoals in algoritme 3, met r0 =
ex
||ex||1 .
27
Algorithm 3 Machtsmethode
1:
2:
3:
4:
5:
6:
procedure
loop:
rk+1 = d ∗ (C ∗ (D\rk )) + (1 − d) ∗ norm ex
normr = norm(rk+1 , 1)
rk+1 = rk+1 /normr
end
5.4
De verschoven machtsmethode
5.4.1
Beschrijving methode
De verschoven machtsmethode wordt beschreven, omdat het mogelijk een snellere convergentie van het
residu heeft, ofwel een snellere manier om de rangschikkingsvector te verkrijgen.
Laat (λ, x) een eigenpaar van de matrix T , en σ een constante, dan is (λ − σ, x) een eigenpaar van de
matrix (A − σI). Het getal σ heet een verschuiving, aangezien het een eigenwaarde λ van T verschuift
2|
naar λ − σ van B = A − σI. Kies µ > 0, met µ een eigenwaarde van B. Nu zoeken we de σ zodat |µ
|µ1 |
minimaal, ofwel dat minimaal is:
|(λk −σ|)2
(|λk −σ)1 | .
Nu het berekenen van de optimale verschuiving σ. Bekend is dat λmax = 1. Met het commando ’eig’ in
matlab, wordt gevonden dat in de connectiematrix C van het GeneRank probleem geldt dat −d ≤ λi ≤
0.5, voor i = 2, 3, . . . , n, en wegens Perron-Frobenius hebben we de eigenschap dat T primitief is, ofwel
dat:
1 = |λ1 | > |λ2 | ≥ |λ3 | ≥ · · · ≥ |λn |.
Dat impliceert met d = λ2 en −d = λmin :
|λmin − σ| = |λ2 − σ|
⇒ | − d − σ| = |d − σ|
d+σ =d−σ
⇒ σ = 0.
Ofwel de optimale verschuiving is: σ = 0 voor d = 0.5.
In het algemeen geldt: A − σI convergeert met een optimale ratio, wanneer σ =
in [18].
5.4.2
λ2 +λn
,
2
zoals geschreven
Implementatie
Een kleine aanpassing aan de machtsmethode zonder verschuiving geeft algoritme 4, met r0 =
28
ex
||ex||1 .
Algorithm 4 Verschoven machtsmethode
procedure
T = d ∗ (C ∗ (inv(D))) + (1 − d) ∗ (norm ex ∗ ones(1, n))
eps = 1e − 12
loop
rk+1 = T ∗ rk − s ∗ rk
normr = norm(rk+1 , 1)
7:
rk+1 = rk+1 /normr
8:
if ||errk || < eps then break
9: end
1:
2:
3:
4:
5:
6:
Wanneer we een verschuiving σ introduceren, convergeert de fout helaas niet sneller. De optimale verschuiving is wanneer σ = 0.
5.5
De methode van Jacobi
Zoals eerder gevonden in vergelijking 1 kon het eigenwaardeprobleem T r = r herschreven worden tot:
(I − dCD−1 )r = (1 − d)ex,
(2)
een eigenwaardeprobleem. Dit is een lineaire vergelijking die met iteratieve methodes kan worden opgelost.
Een iteratieve methode om Ar = b op te lossen is de methode van Jacobi. Er zal worden geconcludeerd
dat deze methode voor het GeneRank probleem hetzelfde systeem oplevert als de zojuist behandelde
machtsmethode.
Schrijf A = L + I + U , met de matrix L de elementen van A onder diagonaal, matrix I de diagonaalelementen van A, en de matrix U de elementen van A boven de diagonaal. Dan geldt dat: (L + I + U )r = b.
⇒ r = b − (L + U )r
De methode van Jacobi werkt met vaste punt iteratie, zoals beschreven in [18]. Het uitvoeren hiervan op
het systeem geeft:
ri+1 = (b − (L + U )ri )
Verder toepassen van de methode van Jacobi op vergelijking 2 geeft: L + U = −dCD−1 en b = (1 − d)ex.
Hieruit volgt::
ri+1 = (1 − d)ex + dCD−1 ri ,
overeenkomstig met de machtsmethode.
Nu kunnen we concluderen dat alle resultaten die worden verkregen door toepassing van de machtsmethode gelijk zijn aan de resultaten die worden verkregen door toepassing van de methode van Jacobi.
29
6
Het gebruik van symmetrie
De symmetrie van de connectiematrix C is een voordelige eigenschap op verschillende gebieden. In het
komende hoofdstuk zullen we zien dat wegens symmetrie de geconjugeerde gradi¨enten methode kan
worden toegepast, zodat relatief erg snel een rangschikkingsvector te vinden is. Er zal eerst worden
gekeken naar hoe de symmetrie van de connectiematrix C zorgt voor een mogelijke omschrijving van A =
(I−dCD−1 ) tot een symmetrisch systeem, vervolgens wordt de CG methode uitgelegd, en als laatste wordt
er gekeken naar de PCG methode toegepast op het GeneRank probleem, ofwel met preconditionering
toegepast op het GeneRank probleem. Een ander voordelig gevolg van de symmetrie van C is dat er
genoeg informatie aanwezig is om gemakkelijk een begrenzing voor de eigenwaardes van de matrix T van
het systeem te vinden. Er zal een bewering worden gedaan over deze eigenwaardes, en die zal worden
bewezen.
6.1
Aanpassing systeem tot symmetrisch computatieprobleem
Terwijl de connectiematrix C symmetrisch is, beschrijven de meeste bekende algoritmes de berekening
van het GeneRank probleem met methodes voor niet-noodzakelijk symmetrische matrices. Er wordt vaak
een niet-symmetrisch lineair systeem opgelost, of een niet-symmetrisch eigenwaardeprobleem. Met andere
woorden, geen van de beschreven algoritmes gebruikt de structuur van C effici¨ent. Nu is de vraag: kunnen
we het GeneRank probleem transformeren tot een symmetrisch matrix computatie probleem, zodat meer
effici¨ente numerieke algoritmes kunnen worden toegepast?
In dit hoofdstuk, wordt laten zien dat het niet-symmetrische lineaire systeem A = (I−dCD−1 ) kan worden
herschreven tot symmetrisch positief definiet (SPD) lineair systeem, en wordt het systeem opgelost,
gebruik makend van een geconjugeerde gradi¨enten methode met preconditionering (PCG) toegepast op
het GeneRank algoritme. Dit algoritme is gebaseerd op wat geschreven is in [17].
Het geldt dat T r = r om kan worden geschreven, zoals in vergelijking 1, ofwel tot:
(I − dCD−1 )r = (1 − d)ex
Het symmetrisch maken van A = (I − dCD−1 ) kan als volgt:
(D − dC)D−1 r = (1 − d)ex.
Dan geldt, met b
r = D−1 r, dat:
(D − dC)b
r = (1 − d)ex.
b = (D − dC) is symmetrisch. Nu kunnen methodes voor symmetrische matrices worden
De matrix A
toegepast op het GeneRank probleem, om de vector b
r = D−1 r, en vervolgens de rangschikkingsvector r,
te berekenen.
6.2
De CG methode
6.2.1
Voorbereidende kennis
De geconjugeerde gradi¨entenmethode (CG) is een algoritme voor de numerieke oplossing van een bepaald
systeem van lineaire vergelijkingen met een symmetrisch en positiefdefiniete matrix, voorgesteld door
Hestenes en Stiefel in 1952.
30
b = (D − dC) een positief definiete matrix is.
Om de CG methode toe te mogen passen, is het nodig dat A
b = (D − dC) is een symmetrische positief definiete (SPD) matrix.
Stelling: A
b symmetrisch is. Te bewijzen: A
b is
Bewijs: In de vorige paragraaf is door omschrijvingen verkregen dat A
b is positief definiet als rT Ar
b > 0 ∀ vectoren 0 6= r ∈ Rn , en voor 0 ≤ d < 1.
positief definiet. A
Wegens de stelling van Gershgorin geldt dat de diagonaalelementen van een matrix br = b, met b
Het doel is om de CG methode toe te passen op Ab
r = D−1 r. De CG methode minimaliseert
daarmee de benadering: ||b
rk − b
r||Ab.
Een groot voordeel van CG is de bescheiden geheugenruimte die het inneemt, namelijk slechts 4 vectoren
met lengte n, en daarmee de convergentie-eigenschappen.
b als uT Av
b = 0. Ofwel, als ze zijn orthoTwee vectoren u,v 6= 0 zijn geconjugeerd met betrekking tot A
gonaal ten overstaande van dit inproduct. Geconjugeerd zijn is een symmetrische relatie.
De conjugatievoorwaarde is een orthonormaal-type voorwaarde en draagt overeenkomsten met GramSchmidt orthonormalisatie. Door de conjugatie en slimme constructie van getallen en vectoren in het
algoritme, is er een kleine geheugenruimte nodig om het algoritme uit te voeren, en convergeert het
algoritme snel. Hier meer over in de komende paragrafen.
6.2.2
Beschrijving CG methode
Hier volgt een bondige beschrijving van de CG methode uit [18], om een beeld te krijgen van hoe de
methode in elkaar steekt.
De volgende vergelijking lossen we op:
rk+1 = rk + αk pk ,
met αk zodanig dat ||b
rk − b
r||Ab minimaal.
Het residu is:
zk+1 = zk − αk Apk .
De zoekdirectie p nemen we in eerste instantie gelijk aan het residu z, ofwel: p0 = z0 . Om te zorgen dat
op´e´envolgende zoekvectoren, pk en pk+1 onderling geconjugeerd zijn, stellen we pi+1 gelijk aan:
pk+1 = rk+1 + βk pk ,
met βk zodanig dat pk Apk+1 = 0.
6.2.3
Implementatie
De invoer voor het algoritme: beginschattingsvector r0 , waarbij we zvva kunnen aannemen dat r0 = 0.
Beginnend met een beginschattingsvector, mogen we de eerste basisvector willekeurig kiezen, en kiezen
we: p0 = b − Ar0 . De andere vectoren in de basis zullen geconjugeerd zijn met de gradi¨ent: Ar0 − b.
Het algoritme is beschreven in algoritme 5.
31
Algorithm 5 Geconjugeerde gradi¨enten methode
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
procedure
= 1e − 12
z0 = b − Ar0
p0 = z0
k=0
loop:
zT z
αk = pTkApk
k
k
rk+1 = rk + αk pk
zk+1 = zk − αk Apk
if ||zk+1 || < then break
zT
z
k+1
βk+1 = k+1
zT
k zk
12:
pk+1 = zk+1 + βk pk
13:
k =k+1
14: end
11:
6.3
De PCG methode toegepast op het GeneRank probleem
6.3.1
Beschrijving PCG toegepast op het GeneRank probleem
Gegeven een lineair systeem Ar = b, betekent preconditioneren het vervangen van het originele systeem
door het systeem: M −1 Ar = M −1 b, met M de preconditionering. Een precondtionering is in dit geval
voordelig, aangezien de matrix M −1 A dichter bij de identiteitsmatrix I ligt. Dit impliceert dat de methode
sneller zal convergeren naar de gewenste rangschikkingsvector r.
b met de
Voor een SPD lineair systeem, is een goede preconditionering, zie [17], een benadering van A
volgende eigenschappen:
• M is symmetrische en positief definiet;
b ≈ I;
• M −1 A
• operaties met M −1 hebben niet veel geheugenruimte nodig.
We kunnen PCG niet direct toepassen op het systeem M −1 Ar = M −1 b, aangezien M −1 A in het algemeen
niet-symmetrisch is. Om de PCG methode te kunnen toepassen, preconditioneren we het systeem op de
volgende manier:
1
b − 12 )(M 12 r) = M − 21 b.
(M − 2 AM
b en M − 21 AM
b − 21 delen dezelfde eigenwaardes. Nu, is het voorstel om de digonaalmatrix
Merk op: M −1 A
als preconditionering te gebruiken om het lineaire systeem (D − dC)b
r = (1 − d)ex, op te kunnen lossen,
ofwel:
M = D = diag(d1 , d2 , . . . , dn ).
Hieruit volgt dat het lineaire systeem met preconditionering wordt:
1
1
1
1
r) = (1 − d)(D− 2 ex).
(I − dD− 2 CD− 2 )(D 2 e
6.3.2
Implementatie
Invoer voor het algoritme is: tolerantiegetal , dempingsfactor d, deg = [deg1 , deg2 , . . . , degn ]T , en k = 0.
ex
ex
Neem verder als invoer de vectoren: r0 = ||ex||
, z0 = ||ex||
, p0 = ex./deg, en y0 = p0 . Het algoritme
1
1
is beschreven in algoritme 6.
32
Algorithm 6 PCG methode toegepast op het GeneRank probleem
1:
2:
3:
4:
5:
6:
procedure
= 1e − 12
k=0
loop:
k =k+1
qk = deg. ∗ pk − d ∗ (C ∗ pk )
yT
z
11:
k−1
αk = k−1
pT qk
rk = rk−1 + αk ∗ pk
zk = zk−1 − αk ∗ qk
if zk+1 < then break
yk = zk ./deg
12:
βk+1 =
7:
8:
9:
10:
13:
14:
15:
16:
yT
k zk
yT
k−1 rk−1
pk+1 = yk + βk+1 ∗ pk
end
r∗ = deg. ∗ rk
rCG = r∗ /norm(r∗ , 1)
Als uitvoer van het algoritme krijgen we de rangschikkingsvector r.
6.4
Begrenzing eigenwaardes
Stelling: Alle eigenwaarden van de GeneRank matrix zijn re¨eel. Bovendien geldt λmax = 1 en −d ≤ λ ≤ d
voor alle andere eigenwaarden, met d de dempingsfactor.
Bewijs: Immers, er geldt dat:
Tr = r
−1
T
dCD + (1 − d)exe r = r.
Er geldt dat ex =
ex
||ex||1 ,
en CD−1 een stochastische matrix. Dan eT CD−1 = eT en eT ex = 1. Nu volgt:
T
de + (1 − d)eT r = eT r
⇔ eT r = eT r = 1
We weten nu wat de linker eigenvector in dit geval is. Een linker eigenvector is een rij-vector u met
de eigenschap: uA = λu. Zo is de rechter eigenvector een kolomvector v, die genoteerd staat aan de
rechterkant van de matrix A, ofwel: Av = λv.
In dit geval is de linkereigenvector bij eigenwaarde λ = 1: e.
Nu, vul in: xj , de j-de eigenvector, j > 1, ⊥ e. Dan:
dCD−1 xj + 0 = λj xj
.
Dus voor (λj , xj ), met j = 2, . . . , n geldt:
dCD−1 xj = λj xj .
33
Nu gebruiken we de stelling van Gershgorin, die het gewenste resultaat zal geven. De stelling van Gershb een symmetrisch n × n-matrix, met elementen aij , en laat
gorin P
heeft de volgende aannames: Laat A
Ri = j6=i |aij | voor i ∈ {1, 2, . . . , n} de som zijn van de absolute waarde van de niet-diagonaal elementen in de i-de rij. Laat D(aii , Ri ) de gesloten schijf zijn met centrum aii en straal Ri . Voldoenend aan de
aannames, geeft de stelling: Iedere eigenwaarde van A ligt in tenminste ´e´en van de Gershgorin schijven
D(aii , Ri ).
Nu kunnnen we met behulp van de stelling van Gershgorin concluderen dat met straal 1 en centrum 0,
dat λ2 = 1 ∗ d en λn = −1 ∗ d met λ1 > λ2 ≥ . . . ≥ λn , en −d ≤ λ ≤ d, ∀λ 6= λmax = λ1 . 34
7
Numerieke experimenten
In de vorige hoofdstukken is er geschreven over het GeneRank probleem en de theorie van verschillende
methodes die kunnen worden toegepast om de rangschikkingsvector van genen te vinden. De methodes
die achtereenvolgens zijn behandeld: de stelselmethode, de (verschoven) machtsmethode, de methode
van Jacobi, en de geconjugeerde gradi¨enten methode (met preconditionering) (PCG). In het komende
hoofdstuk zullen de numerieke experimenten die uitgevoerd zijn op de data worden getoond. Er wordt
gefocust op wat er gebeurt wanneer de machtsmethode en de PCG methode worden toegepast. Er zal
informatie worden gegeven over zowel de uitvoertijd van het algoritme, als het aantal iteraties dat nodig
is om tot de rangschikkingsvector te komen.
7.1
Dataset met recente data
Er wordt gebruikt gemaakt van drie matrices uit [6], ofwel de meest recente data van de GO website in
april 2005.
1. de matrix behorend bij de expressievector ex met zowel positieve als negatieve expressie getallen;
2. de matrix behorend bij de expressievector ex met strikt negatieve expressie getallen;
3. de matrix behorend bij de expressievector ex met strikt positieve expressie getallen.
De hoeveelheid genen en verbindingen die de drie matrices bevatten, is te vinden in tabel 1.
Naam
Matrix 1
Matrix 2
Matrix 3
Aantal genen
Aantal verbindingen tussen genen
4047
1625
2392
169 798
33 626
128 468
Tabel 1: Informatie over de drie matrices met data
7.2
Convergentie in de machtsmethode
7.2.1
Het residu
De convergentie factor is de ratio: c = | λλ12 | = |λ2 | = d, met d de scalar die ervoor zorgt dat de transitiematrix T = d ∗ T1 + (1 − d) ∗ T2 een convex geheel is, aangezien we te maken hebben met een Markov
matrix, met grootste eigenwaarde gelijk aan 1. De asymptotische convergentiefactor is dus gelijk aan de
grootte van de op ´e´en na grootste eigenwaarde van T , ofwel λ2 = d. De machtsmethode convergeert:
rk → r. Convergentie ratio: |λ2 | < 1, aangezien λ1 = 1 de grootste eigenwaarde is, die uniek is. Er is
|λ2 |
langzame convergentie als |λ2 | ≈ 1. De machtsmethode zal snel convergeren als |λ1|
= |λ2 | klein is, en
langzamer hoe dichter |λ2 | bij 1 in de buurt komt. In de methodes is aangenomen dat d = 0.5, zoals in
artikel [6] is voorgesteld.
7.2.2
Implementatie van het residu
Het residu van de machtsmethode plotten in matlab vindt plaats door eerst een residu-vector aan te
maken, en vervolgens iedere iteratie het residu in de residu-vector op te slaan. De residu-vector heeft
35
ex
.
daarmee lengte gelijk aan het aantal iteraties. We gebruiken de expliciete machtsmethode, met r0 = ||ex||
1
Het algoritme, dat slechts de residu berekening als toevoeging heeft op algoritme 2, ziet er uit als in
algoritme 7.
Algorithm 7 Machtsmethode
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
procedure
resvec0 = ones(N, 1)
eps = 1e − 12
loop:
rk+1 = T ∗ rk
normr = norm(rk+1 , 1)
resveck = norm(rk+1 − rk )
rk+1 = rk+1 /normr
if ||resveck || < eps then break
end
7.2.3
Convergentiefactor machtsmethode
Bekijken we de machtsmethode, dan hangt convergentiesnelheid, en daarmee de convergentiefactor, af
van het residu:
k !
λ2 k
k
k
||T r − r || = O = O |λ2 |
λ1
Bij de machtsmethode heeft de volgende vector cmacht de convergentiefactor berekend:
cmachti =
36
resveci
resveci−1
De resulterende vector, berekend met matrix 1, cmacht =
0.00784931856616076
0.126148343420699
0.249107106630327
0.365224879931140
0.394825912045956
0.437841452728645
0.452171381228828
0.473379069080085
0.479556623359947
0.488102925456741
0.490804094380037
0.494171183387977
0.495487260828022
0.496913902494186
0.497610870410932
0.498275752257954
0.498662880494112
0.499000333131726
0.499221684857274
0.499404530954619
0.499533797480900
0.499637635676406
0.499714478702855
0.499775458935961
0.499821883383959
0.499858588491153
0.499887061416148
0.499909720663980
0.499927309401110
0.499941113832902
De factor is hier ongeveer: 0.5 = d. De tijd die erover wordt gedaan om de loop uit te voeren is ongeveer
0.95 seconde.
7.3
Convergentie in de verschoven machtsmethode
7.3.1
Het residu
Het residu wordt voor de verschoven machtsmethode op analoge wijze verwerkt als bij de machtsmethode
zonder verschuiving.
7.3.2
Implementatie van het residu
Het plotten van het residu van de verschoven machtsmethode vindt plaats door een residu vector aan
te maken, en in iedere iteratie het residu in de residu vector op te slaan. Het gebruikte algoritme, met
ex
, is weergegeven in algoritme 8. Het algoritme is een uitbreiding op algoritme 4.
r0 = ||ex||
1
37
Algorithm 8 Verschoven machtsmethode
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
procedure
T = d ∗ (C ∗ (inv(D))) + (1 − d) ∗ (norm ex ∗ ones(1, n))
eps = 1e − 12
loop
rk+1 = T ∗ pk − s ∗ pk
normr = norm(rk+1 , 1)
resveck = norm(T ∗ pk − pk )
lambda = normr + s
pk+1 = rk+1 /normr
if ||resveck || < eps then break
end
7.3.3
Resultaten
De resultaten zijn op basis van de methode, gebruik maken van matrix 1. In figuur 7, 8 en 9, zijn de
residuen weergegeven voor de verschoven machtsmethode met verschuiving σ = 0.2 σ = 0.1, en σ = 0
respectievelijk.
Figuur 7: Convergentiecurve verschoven machtsmethode met σ = 0.2
38
Figuur 8: Convergentiecurve verschoven machtsmethode met σ = 0.1
Figuur 9: Convergentiecurve verschoven machtsmethode met σ = 0
De resultaten weergegeven in een tabel zijn te vinden in tabel2 met verschuiving σ = 0.2, σ = 0.1,
σ = 0.01 en σ = 0.
39
Matrix dimensie
σ = 0.2
σ = 0.1
σ = 0.01
σ=0
= 10−2
= 10−4
= 10−6
= 10−8
= 10−10
= 10−12
0.0497 sec.
3 it.
0.0307 sec.
3 it.
0.0443 sec.
3 it.
0.0335 sec.
3 it.
0.2381 sec.
19 it.
0.0919 sec.
8 it.
0.0691 sec.
6 it.
0.0604 sec.
6 it.
0.6210 sec.
53 it.
0.2886 sec.
19 it.
0.2184 sec.
13 it.
0.1538 sec.
12 it.
1.0960 sec.
88 it.
0.4049 sec.
31 it.
0.2282 sec.
20 it.
0.2332 sec.
19 it.
1.3816 sec.
122 it.
0.5460 sec.
42 it.
0.3667 sec.
27 it.
0.3298 sec.
26 it.
1.6571 sec.
157 it.
0.6316 sec.
53 it.
0.4166 sec.
34 it.
0.3510 sec.
32 it.
Tabel 2: Optimale verschuiving
Overeenkomstig met dat wat we hebben gezien bij de beschrijving van de vershoven machtsmethode,
laten de figuren 7, 8, 9 en tabel 2 zien dat de optimale verschuiving σ = 0 is.
7.4
Convergentie in de methode van Jacobi
Voor de methode van Jacobi geldt overeenkomstig met de machtsmethode:
ri+1 = (1 − d)ex + dCD−1 ri ,
De volgende twee oplossingen kunnen van elkaar af worden getrokken:
ri+1 = dCD−1 ri + (1 − d)ex en
r = dCD−1 r + (1 − d)ex, zodat overblijft:
r − ri+1 = dCD−1 (r − ri )
waarbij dCD−1 de versterkingsmatrix is.
Zo vinden we:
||r − ri+1 || ≈ λmax (dCD−1 )(r − r0 ) = (dCD−1 )(r − r0 ).
met λmax de spectraalradius, het supremum over de absolute waardes van de elementen in het spectrum,
ofwel ρ(T ) = λmax = maxi (|lambdai |).
De methode van Jacobi heeft dezelfde convergentiefactor als de machtsmethode, namelijk: d.
7.5
Convergentie in de PCG methode toegepast op het GeneRank probleem
7.5.1
Het residu
De convergentie ratio van CG wordt verwacht af te nemen als d afneemt. We bekijken de convergentie van
de PCG methode toegepast op het GeneRank algoritme. Stel ven A ≡ D − dC, en b
r de exacte oplossing
van het systeem. Zie in [18] dat geldt :
||b
r−r
(k)
||A ≤ 2
p
κ2 (D − dC) − 1
p
κ2 (D − dC) + 1
!k
||b
r − r(0) ||A ,
met r(k) is geproduceerd door de PCG methode,
κ2 (D − dC) =
λmax (D − dC)
,
λmin (D − dC)
is het 2-norm conditie getal van A, en ||v||2A ≡ vT (D − dC)v.
40
De convergentie is te berekenen, aangezien alle informatie bekend is, namelijk:
κ2 (D − dC) =
zodat:
q
||b
r − r(k) ||A ≤ q
3
2
−1
3
2
+1
3
,
2
≈
0.7
≈ 0.25.
0.2
Dit vertelt dat het aantal CG iteraties nodig om de fout door een vaste factor kleiner dan 1 te reduceren,
proportioneel is met de wortel van het twee-conditie getal κ2 (D − dC). Hier is het twee-conditie getal de
ratio van de grootste eigenwaarde door de kleinste eigenwaarde van (D − dC). De convergentiefactor is
niet alleen door de factor van de grootste en de kleinste eigenwaarde verkregen, maar ook door de gehele
verdeling van eigenwaardes. Het algoritme gaat als in algoritme 9.
7.5.2
Implementatie van het residu
Het residu van de PCG methode wordt in iedere iteratie in de residu-vector opgeslagen. Het algoritme is
ex
ex
, z0 = ||ex||
,
beschreven in algoritme 9, met dempingsfactor d, deg = [deg1 , deg2 , . . . , degn ]T , r0 = ||ex||
1
1
p0 = ex./deg, en y0 = p0 .
Algorithm 9 PCG methode toegepast op het GeneRank probleem
procedure
k=1
eps = 1e − 12
4: loop:
5:
qk = deg. ∗ pk − d ∗ (C ∗ pk )
1:
2:
3:
yT
z
10:
11:
12:
k−1
αk = k−1
pT qk
rk = rk−1 + αk ∗ pk
zk = zk−1 − αk ∗ qk
if zk+1 < then break
normz = norm(z, 1)
resveck+1 = [resveck ; normz]
yk = zk ./deg
13:
βk+1 =
6:
7:
8:
9:
yT
k zk
yT
k−1 rk−1
pk+1 = yk + βk+1 ∗ pk
k =k+1
16: end
14:
15:
7.5.3
Convergentiefactor PCG methode
Bij de PCG methode heeft de volgende vector cpcg de convergentiefactor berekend:
cpcgi =
resveci
resveci−1
41
De resulterende vector, berekend met matrix 1 cpcg =
0.434234349250101
0.260640047854833
0.308150479923512
0.217011918052124
0.250249959369377
0.241909388916971
0.242381422469359
0.209685140833943
0.221689207074938
0.230440411529726
0.198277841625850
0.237966672216885
0.255432972172009
0.231662348035950
0.206759197121159
0.224333871104718
0.222044757461802
0.217144955678176
0.219335415876562
0.223585027269756
De convergentiefactor is ongeveer: 0.25. De tijd die erover wordt gedaan om de loop uit te voeren is
ongeveer 0.06 seconde.
7.6
Vergelijken resultaten machts- en PCG methode
7.6.1
Variatie dempingsfactor d
In artikel [6], pag. 11, is berekend dat een goede waarde voor d zou zijn: d = 0.5. Andere waardes voor
d zijn daarmee niet uitgesloten. Om een beeld te geven van wat het residu is voor verschillende waardes
van d, is de afname van de grootte van het residu weergegeven in de tabellen 3, 4 en 5, matrix 1, 2 en 3
respectievelijk, met alle residuwaardes ∗10−8 . Daar waar een 0 staat is het residu kleiner dan 1 ∗ 10−8 .
d
Iteratie
Machtsmethode
PCGmethode
1
5
10
20
30
40
50
0.85
0.75
0.65
0.55
0.5
1334384.16
1177397.78
1020411.41
863425.04
784931.86
50502.19
27009.83
13206.34
5728.34
3556.85
15730.23
4499.44
1075.67
202.38
78.03
2912.98
238.33
13.62
0.48
0.07
571.21
13.37
0.18
0
0
112.41
0.75
0
0
0
22.13
0.04
0
0
0
0.85
0.75
0.65
0.55
0.5
86793836.89
72921072.32
60314394.13
48808043.73
43423434.93
9162313.58
2920430.21
1008137.32
338522.20
189402.53
231844.14
26240.76
3279.09
377.92
118.96
136.93
1.87
0.03
0
0
0.03
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Tabel 3: Matrix 1, afname van de norm van het residu bij verschillende parameters d
42
d
Iteratie
Machtsmethode
PCGmethode
1
5
10
20
30
40
50
0.85
0.75
0.65
0.55
0.5
1794123.59
1583050.23
1371976.86
1160903.50
1055366.82
66402.27
35513.59
17364.21
7531.85
4676.69
20843.98
5962.17
1425.36
268.18
103.39
3819.98
312.54
17.86
0.63
0.09
746.22
17.46
0.24
0
0
146.71
0.98
0
0
0
28.87
0.06
0
0
0
0.85
0.75
0.65
0.55
0.5
88559731.51
73222250.47
59701302.97
47692212.91
42177527.24
7790332.31
2425437.29
824303.58
273410.27
152146.18
246233.26
27258.54
3365.89
385.34
121.04
177.07
2.43
0.04
0
0
0.11
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Tabel 4: Matrix 2, afname van de norm van het residu bij verschillende parameters d
d
Iteratie
Machtsmethode
PCGmethode
1
5
10
20
30
40
50
0.85
0.75
0.65
0.55
0.5
1822015.15
1607660.43
1393305.71
1178950.98
1071773.62
92078.47
49245.87
24078.55
10444.24
6485.05
36359.80
10400.28
2486.37
467.80
180.36
7029.41
575.13
32.87
1.16
0.17
1382.28
32.35
0.44
0
0
272.10
1.82
0.01
0
0
53.57
0.10
0
0
0
0.85
0.75
0.65
0.55
0.5
87626402.85
73556847.46
60792443.04
49159632.62
43721570.51
8829500.45
2780371.64
954294.23
319246.06
178345.22
204184.96
23801.66
3035.29
355.12
112.47
87.02
1.20
0.02
0
0
0.02
0
0
0
0
0
0
0
0
0
0
0
0
0
0
Tabel 5: Matrix 3, afname van de norm van het residu bij verschillende parameters d
In de tabellen 3, 4 en 5 is af te lezen dat het residu bij meer iteraties afneemt. De machtsmethode
begint met een groter residu dan de PCG methode, en de PCG methode convergeert sneller dan de
machtsmethode, zodat de PCG methode na een kleinere hoeveelheid iteraties afgerond is. Ook is duidelijk
te zien dat hoe verder de parameter d verwijderd van 1, hoe sneller het residu convergeert. Het valt op
dat de drie matrices na genoeg hetzelfde gedrag vertonen.
Nu enkele grafische weergaves voor d = 0.75, anders dan de dempingsfactor die we door het rapport heen
aanhouden, d = 0.5, om te laten zien wat een ander dempingsgetal voor resultaat geeft. De machtsmethode
genereert een stationaire rangschikkingsvector na 73 iteraties. De tijd dat het algoritme wordt uitgevoerd
neemt ongeveer 2.2 seconde in. Dit is te zien in figuur 10.
43
Figuur 10: Convergentiecurve machtsmethode met d = 0.75
De PCG methode genereert voor d = 0.75 een stationaire rangschikkingsvector na 31 iteraties zoals te
zien in figuur 11. De tijd dat het algoritme wordt uitgevoerd neemt ongeveer 0.09 seconde in.
Figuur 11: Convergentiecurve PCG methode met d = 0.75
44
En voor d = 0.5 genereert de machtsmethode een stationaire rangschikkingsvector na 31 iteraties. De tijd
dat het algoritme wordt uitgevoerd neemt ongeveer 0.94 seconde in. Een weergave is in 12.
Figuur 12: Convergentiecurve machtsmethode met d = 0.5
De PCG methode genereert voor d = 0.5 een stationaire rangschikkingsvector na 21 iteraties. De tijd dat
het algoritme wordt uitgevoerd neemt ongeveer 0.06 seconde in. Een weergave hiervan is 13.
45
Figuur 13: Convergentiecurve PCG methode met d = 0.5
7.6.2
Tijdsduur en aantal iteraties
Aangezien de rangschikkingsvector r een stationaire verdeling van een Markov keten is, gebruiken we
de 1-norm van het residu als stopcriterium in dit rapport. Het stopcriterium is dan ||Ark − λk rk || =
||Ark − rk || < , met > 0, en λk = 1 de eigenwaarde bij de k-de iteratie. Deze stopcriteria voor de
machtsmethode en de PCG methode komen overeen. Dit is als volgt in te zien.
Voor de machtsmethode geldt dat:
||T rk − rk || = ||(dCD−1 − (1 − d)exeT )rk − rk || = ||(dCD−1 − I)rk − (1 − d)ex||,
terwijl voor de PCG methode geldt dat:
||T rk − rk || = ||b − Ark || = ||(1 − d)ex − (I − dCD−1 )rk ||,
met b = (1 − d)ex en A = (I − dCD−1 ).
In tabel 6 en 7 zijn de tijdsduur (in secondes) en het aantal iteraties (it.) weergegeven, voor verschillende
. De tolerantie van het residu, , wordt gevarieerd tussen 10−2 en 10−12 , en d = 0.5.
Matrix dimensie
Matrix 1
Matrix 2
Matrix 3
= 10−2
= 10−4
= 10−6
= 10−8
= 10−10
= 10−12
0.0444 sec.
2 it.
0.0212 sec.
3 it.
0.0362 sec.
3 it.
0.1384 sec.
5 it.
0.0402 sec.
6 it.
0.0957 sec.
6 it.
0.3219 sec.
11 it.
0.0687 sec.
12 it.
0.1372 sec.
12 it.
0.5274 sec.
18 it.
0.0977 sec.
18 it.
0.2514 sec.
19 it.
0.7421 sec.
24 it.
0.1354 sec.
25 it.
0.3141 sec.
26 it.
0.9175 sec.
31 it.
0.1766 sec.
31 it.
0.3776 sec.
32 it.
Tabel 6: Machtsmethode, tijdsduur en aantal iteraties
46
Matrix dimensie
Matrix 1
Matrix 2
Matrix 3
= 10−2
= 10−4
= 10−6
= 10−8
= 10−10
= 10−12
0.0141 sec.
5 it.
0.0031 sec.
5 it.
0.0038 sec.
5 it.
0.0226 sec.
9 it.
0.0037 sec.
8 it.
0.0045 sec.
8 it.
0.0358 sec.
12 it.
0.0047 sec.
12 it.
0.0061 sec.
12 it.
0.0376 sec.
15 it.
0.0051 sec.
15 it.
0.0072 sec.
15 it.
0.0491 sec.
18 it.
0.0058 sec.
18 it.
0.0081 sec.
18 it.
0.0676 sec.
21 it.
0.0067 sec.
21 it.
0.0084 sec.
20 it.
Tabel 7: PCG methode, tijdsduur en aantal iteraties
Uit de tabellen 6 en 7 kan het volgende worden geconcludeerd:
• het aantal iteraties dat nodig is, zowel voor de machtsmethode als voor de PCG methode neemt
toe, naarmate kleiner, zoals verwacht;
• de spreiding van het aantal iteraties is bij de machtsmethode groter dan bij de PCG methode;
• wat tijdsduur betreft, is het duidelijk dat deze toeneemt bij het verkleinen van . Het programma
heeft meer tijd nodig te convergeren bij kleine . Ook neemt de tijdsduur toe bij het vergroten van
de matrixdimensie n;
• de tijdsduur en hebben een duidelijk verband; meer iteraties impliceert een langere tijdsduur;
• ook hier valt op dat de drie matrices na genoeg hetzelfde gedrag vertonen.
47
8
Conclusie
Nu alle resultaten zijn beschreven, zetten we de belangrijke zaken die kunnen worden geconcludeerd,
op een rij. Onderzocht in dit rapport is de verbindingen tussen genen, en methodes om het systeem
T r = r op te lossen, met T de GeneRank transitiematrix, om een rangschikkingsvector r te verkrijgen.
De onderzoeksvragen waarmee dit rapport opende, waren de volgende:
1. Wat houdt de verbinding tussen genen in?
2. Welke effici¨ente en snelle methodes zijn er om de rangschikkingsvector voor het GeneRank probleem
te vinden?
3. Zijn er voordelen voor het vinden van de rangschikkingsvector door de symmetrie van de connectiematrix?
We zullen de onderzoeksvragen achtereenvolgens behandelen.
Het GeneRank probleem werkt met twee soorten data van verbindingen tussen genen, de experimentele
data en de biologisch, ofwel verwachte, data. De experimentele data wordt verkregen door het stimuleren
van genen om prote¨ınes te maken, en zo te kunnen meten hoe actief ze zijn, en daarmee wat hun mate
van belang is. Deze data geeft al een rangschikking, maar is door de onnauwkeurigheid van experimenten
zeer gevoelig. De biologische, ofwel verwachte, data in dit rapport bestond uit de drie componenten:
moleculaire functie, biologische proces en cellulaire component.
We hebben ´e´en directe methode bekekeken: de stelselmethode. Verder bekeken we iteratieve methodes
toegepast op het GeneRank probleem: de (verschoven) machtsmethode die een eigenwaardeprobleem
oploste, en de methode van Jacobi die een lineaire stelsel van vergelijkingen oploste. Alles methode
hadden dezelfde rangschikkingsvector r als uitvoer. Het vergelijken van de methodes was gebaseerd op
convergentiesnelheid. We hebben gezien dat voor het GeneRank probleem geldt dat de methode van Jacobi
gelijk is aan de machtsmethode. We hebben geconcludeerd dat de optimale verschuiving bij toepassing
van de machtsmethode, verschuiving σ = 0 is. De machtsmethode convergeert sneller dan de directe
methode.
Wegens de symmetrie van de connectiematrix C konden we gemakkelijk bewijzen dat alle eigenwaarden
van de GeneRank transitiematrix T re¨eel zijn, dat de grootste eigenwaarde gelijk aan 1 is, en dat de overige
eigenwaardes tussen - de dempingsfactor en + de dempingsfactor liggen, ofwel tussen −d en +d. Door
de symmetrische connectiematrix was de niet-symmetrische matrix A = (I − dCD−1 ) om te schrijven
b Deze omschrijving maakte dat we de geconjugeerde gradi¨entenmethode
tot een symmetrisch matrix A.
met preconditionering toegepast op het GeneRank probleem konden toepassen, een methode die erg snel
convergeert, en ons een goede oplossing heeft gegeven. De PCG methode toegepast op het GeneRank
probleem convergeert sneller dan de machts- en directe methode. We concluderen dat de methode die het
snelst een goed resultaat geeft de PCG methode toegepast op het GeneRank probleem is.
Als laatste, concluderen we dat de drie matrices die zijn gebruikt voor de numerieke experimenten,
oplossingen geven die heel dicht bij elkaar liggen, met de reden dat de eigenwaardes λ1 = 1, λ2 = d en
λn = −d voor de matrices hetzelfde zijn. Omdat de eigenwaardes voor de drie matrices gelijk verdeeld
liggen, hebben de methodes oplossingen geven die dicht bij elkaar liggen.
48
9
Discussie
Het rapport over het onderzoek naar het verbeteren van het GeneRank algoritme ligt achter ons. Het
rapport dat nu bijna ten einde is, ging over de verbinding tussen genen, over verschillende algoritmes om
de rangschikking te verkrijgen, en over de symmetrie van de connectiematrix van het GeneRank probleem,
die voordelen met zich meebrengt.
We zagen dat de GeneRank methode de verbinding tussen genen in het algoritme betrouwbaar probeert
te laten zijn, door op twee soorten data te steunen: de experimentele data van mate van belang van genen,
en de verwachte mate van belang van genen, gekeken naar de biologische verbinding tussen genen. We
bekeken verschillende methodes. De directe methode nam erg veel geheugenruimte in, het algoritme deed
er daardoor lang over om tot een oplossing te komen. Verder bekeken we de (verschoven) machtsmethode,
de methode van Jacobi, en de PCG methode toegepast op het GeneRank probleem.
We zagen dat voor het GeneRank probleem geldt dat de machtsmethode gelijk is aan de methode van
Jacobi, en we zagen dat de optimale verschuiving voor de machtsmethode om het proces zo snel mogelijk te laten convergeren, geen verschuiving was. De machtsmethode nam een kortere tijd in om tot de
oplossing te komen dan de directe methode, zoals we verwachtten. We konden door de symmetrische connectiematrix C een het systeem zo omschrijven dat we de PCG methode konden toepassen, een methode
voor symmetrische matrices, die snel een goede rangschikkingsvector berekende. We verwachtten dat de
methode die gebruik maakt van een symmetrische matrix sneller zou convergeren, dan de machtsmethode,
de snelst convergerende methode in dit rapport. Deze verwachting klopte. We zagen dat de PCG methode
dezelfde rangschikkingsvector als de andere methodes gaf, en veel sneller dan de andere methodes een
oplossing vond, overeenkomstig met wat we verwachtten.
De resultaten die we kregen door numerieke experimenten te doen, brachten ons tot de conclusie dat de
drie matrices die waren gebruikt als data, ons na genoeg dezelfde oplossing gaven door de overeenkomstige
ligging van de eigenwaardes van de matrices.
Vervolgonderzoek kan gaan over het vinden van een ander testprobleem van bijvoorbeeld de GO website,
ofwel het omschrijven van een recente GO matrix naar een symmetrische matrix met 0’en en 1’en. Er
is een begin meegemaakt met het onderzoek hiervan, zoals te lezen in de appendix. We verwachten
gelijke resultaten te krijgen na toepassing van de methodes met andere GO matrices, en kunnen dit
bevestigen wanneer de nieuwe testproblemen zijn gereaslieerd. Verder vervolgonderzoek kan gaan over de
dempingsfactor d. In dit rapport namen we naar aanleiding van artikel [6] aan dat d = 0.5. Dit onderzoek
kan nog verder worden ondersteund en uitgebreid.
49
Referenties
[1] http://www.oeiikgroei.nl/baby/medisch/dna-genen-en-chromosomen/
[2] www.bioontology.org
[3] Adi L. Tarca, Roberto Romero en Sorin Draghici. Analysis of microarray experiments of gene
expression profiling. American Journal of Obstetrics and Gynecology, 195(2): 373-388. 2008
[4] HC King, AA Sinha Gene expression profile analysis by DNA microarrays: promise and pitfalls.
JAMA: the journal of the American Medical Association, 286(18):2280-8. 2001.
[5] http://en.wikipedia.org/wiki/Gene expression
[6] Julie L. Morrison, Rainer Breitling, Desmond J. Higham en David R. Gilbert. GeneRank:
Using search engine technology for the analysis of microarray experiments. BMC Bioinformatics,
6:233. Glasgow, 2005.
[7] Gilbert Strang. Linear Algebra and Its Applications (2e editie). Pagina 288. Academic Press. New
York, 1980.
[8] http://www.geneontology.org/GO.format.annotation.shtml.
[9] http://www.geneontology.org/GO.format.obo-1 2.shtml
[10] http://www.geneontology.org/
[11] GO.annotation.conventions.shtml
[12] http://bioportal.bioontology.org/
[13] http://www.intermine.org/
[14] http://ratmine.mcw.edu
[15] http://www.geneontology.org/ontology/gene ontology.obo
[16] M. Ashburner, CA Ball, JA Blake et al. The Gene Ontology Consortium. Gene ontology: tool
for the unification of biology. Nature Genetics, 25(1):25-9. California, 2000.
[17] Gang Wu, Wei Xu, Ying Zhang en Yimin Wei. A preconditioned conjugate gradient algorithm
for GeneRank with application to microarray data mining. Data Min Knowl Disc, 26:27-56. 2011.
[18] Yousef Saad. Iterative Methods for Sparse Linear Systems (1e editie). SIAM. California, 2003.
[19] Schlomo Sternberg. Dynamical Systems. Dover Publications. Harvard, 2010.
[20] BaoJun Yue, Heng Liang, Fengshan Bai. Understanding the GeneRank Model. IEEE 1st Int
Conf. Bioinform. Biomed. Eng., 6-8. Beijing, 2007
[21] http://homes.esat.kuleuven.be/ bioiuser/gpp/index.php
[22] Davod Khojasteh Salkuyeh, Vahid Edalatpour, en Davod hezari. A new preconditioner
for the GeneRank problem. arXiv, 1403.3925 [math.NA]. 2104.
50
10
Appendix
10.1
Methodes
10.1.1
Stelselmethode
clear all; close all;
load(’GO_matrix.mat’)
format long
C=w_All; ex=expr_data; d=0.5;
% De experimentele data, de beginvector, expressie verandering gen i-->
ex = abs(ex);
norm_ex = ex/norm(ex,1);
c = sparse(C);
% Som van de rijen van matrix C-->
deg = sum(C);
% In het pagerank geval dangling node, bij genen komt dit niet voor-->
ind = find(deg == 0);
% Deze actie verandert niets aan de matrix, zie vorige punt-->
deg(ind) = 1;
% D invers, som van de rijen van W op de diagonaal-->
D1 = sparse(diag(1./deg));
% (I-dW’Dinv)-->
A = eye(size(c)) - d * (c’*D1);
% (1-d)ex-->
b = (1-d) * norm_ex;
% Oplossen van r-->
r = A\b;
% Rangschikkingsvector r-->
xStelsel = r/norm(r,1);
10.1.2
Machtsmethode
clear all; close all;
load(’GO_matrix.mat’)
format long
C=w_All; ex=expr_data; d=0.5; eps=1e-12; N=100; n=length(C);
ex = abs(ex);
% Beginschattingsvector-->
r0 = ex/norm(ex,1);
norm_ex = ex/norm(ex,1);
lambda = zeros(N,1);
% Construeren matrix D-->
deg = sum(C); D = diag(deg);
51
% Transitiematrix-->
T = d*(C*(inv(D))) + (1-d)*(norm_ex*ones(1,n));
resvec = ones(N,1);
for i = 2:N
r = T*r0;
%
r = d *(C*(D\r0)) + (1-d)*norm_ex;
% In dit geval is het 1 als ex vanaf het begin niet-negatief-->
normr = norm(r,1);
% Fout in de norm van het residu-->
resvec(i) = norm(T*r0-r0);
lambda(i) = normr;
% Eigenvector maken met normr=lambda=1, dus r0=r-->
r0 = r/normr;
% Tolerantie van het residu-->
if resvec(i)<eps
break
end
end
% Afdrukken van het residu-->
semilogy(resvec(1:i),’*’);
title(’Machtsmethode voor GeneRank’);
xlabel(’Aantral matrix-vector producten’)
ylabel(’1-norm van het residu’)
10.1.3
Verschoven machtsmethode
clear all; close all;
load(’GO_matrix.mat’)
format long
shift = input(’shift = ’);
C=w_All; ex=expr_data; d=0.5; eps=1e-12; N=100; n=length(C);
r0 = ex/norm(ex,1);
ex = abs(ex);
norm_ex = ex/norm(ex,1);
lambda = zeros(N,1);
% Construeren matrix D-->
deg = sum(C); D = diag(deg);
T = d*(C*(inv(D))) + (1-d)*(norm_ex*ones(1,n));
resvec = ones(N,1); %error vector
for i = 1:N
r = T*r0 - shift*r0;
normr = norm(r,1);
52
resvec(i) = norm(T*r0-r0);
lambda(i) = normr + shift;
r0 = r/normr;
if resvec(i)<eps
break
end
end
semilogy(resvec(1:i), ’*’)
title(’Machtsmethode voor GeneRank’);
xlabel(’Aantral matrix-vector producten’)
ylabel(’1-norm van het residu’)
10.1.4
PCG methode toegepast op het GeneRank probleem
De berekeningen voordat het itereren begint:
clear all; close all;
load(’GO_matrix.mat’)
format long
C=w_All; ex=expr_data; d=0.5; eps =1e-12; N=100;
ex=abs(ex);
norm_ex = ex/norm(ex,1);
deg=sum(C); ind = find(deg==0); deg(ind)=1;
n=size(C,1);
% D=sparse(diag(deg))-->
deg=deg’;
Beginschattingsvector-->
r=norm_ex;
% 0’de residu-->
z=norm_ex;
% Norm van het 1-residu-->
normz=norm(z,1);
resvec=[1];
rho=1;
Het itereerproces en afdrukken residu:
for i=2:N
% Constructie van p-->
y=z./deg;
rho1=rho;
rho=z’*y;
% Iteratie 1-->
if(i==2)
p=y;
53
% Iteratie 2,3,..-->
else
beta=rho/rho1;
p=y+beta*p;
end
% (D-dC)p-->
q=deg.*p-d*(C*p);
pq= p’*q;
alpha = rho/pq;
% Rangschikkingsvector berekenen-->
r=r+alpha*p;
% Residu berekenen-->
z=z-alpha*q;
normz = norm(z,1);
% normr wordt in iedere iteratie berekend, en aan de vector toegevoegd-->
resvec = [resvec;normz];
% Tolerantie van het residu-->
if(normz<=eps)
break;
end
end
% Berekening uiteindelijke rangschikkingsvector r-->
r=deg.*r;
rCG = r/norm(r,1);
% Afdrukken van het residu-->
semilogy(resvec, ’*’)
title(’PCG voor GeneRank’)
xlabel(’Aantral matrix-vector producten’)
ylabel(’1-norm van het residu’)
10.2
Gen Ontologie (GO) databank
In matlab zit een “Bioinformatics Toolbox”, die voorzien is van verschillende commando’s om data van
buitenaf, van een website, in te lezen en daar mee te werken.
Om te kijken hoe het werkt om data te downloaden en dan te verwerken in matlab, wordt een random
GO matrix in matlab geconstrueerd. Om een matrix te maken kunnen we het commando “getmatrix()”
gebruiken. Nodig om dit commando te gebruiken is een “GeneontObj”, een object, zoals gecre¨eerd door
de “geneont” constructor functie. De output is een connectiematrix, waarin 0 geen relatie representeert,
1 een “is een”-relatie representeert, en 2 een “deel van”-relatie representeert tussen genen.
De “is een”-relatie en “deel van”-relatie zijn transitieve relaties; als A is een deel van B, en B is een deel
van C, dan is A een deel van C.
54
De meest recente versie van de GO databank, is gemakkelijk in matlab te laden met het commando
geneont. De matrix wordt gecre¨eerd met het commando getmatrix.
De meest recente versie van gene ontology.obo kan van de GO website worden ingelezen, vanaf de pagina
[15]. Een willekeurig netwerk van de GO site downloaden, daar wordt zo meteen op ingegaan. Het geldt
zelfs dat ieder willekeurig obo-file door het commando ’geneont’ kan worden verwerkt, en wat door het
commando ’getmatrix’ tot een matrix wordt gevormd.
De recente matrix gedownload van de GO website bevat de getallen 0, 1 en 2. De connectiematrix C die
gewenst is, bevat slechts de getallen 0 en 1, en is symmetrisch. De twee¨en in de matrix geven een “deel
van”-relatie aan. Omdat dit een duidelijke relatie tussen twee genen is, worden de twee¨en in de matrix
getransformeerd naar enen. Om de matrix klaar voor gebruik te maken, moet worden gezorgd dat de
matrix symmetrisch is.
Nu is het de bedoeling dat er niet slechts de meest recente versie van de GO matrix kan worden gebruikt,
maar dat een netwerk van genen kan worden gedownload van de GO website, om die vervolgens als data
in matlab te kunnen gebruiken. De matrices die je van de GO website kunt downloaden zijn zoals de
meest recente versie niet meteen bruikbare matrices, maar is info waarvan de gewenste connectiematrix
kan worden gemaakt. Dit proces gaat op de volgende manier:
Om te beginnen wordt het gedownloade bestand geunzipt, en wordt er gelezen wat voor info er in het
bestand staat:
Nu zijn de gegevens binnen, en wordt de matrix geconstrueerd. De volgende twee vectoren van gegevens
zijn hiervoor nodig: GOid en DB Object Symbol.
Nu moet het commando “geneont” worden gebruikt om een file van de GO site in te lezen, zodat er staat:
“geneont(’File’, FileValue)”, die een gen object cre¨eert, een string, die de bestandsnaam specificeert van
een Open Biomedische Ontology (OBO)-formaat bestand, die op het MATLAB zoekpad is. Een string
wordt gebruikt om output te labelen en om de namen van m-files in andere functies te gebruiken.
Om relaties, zoals voorouders enzo, te visualiseren, wordt de “biograph” functie gebruikt en de “getmatrix”, allen onderdeel van de “Bioinformatics Toolbox”. De “getmatrix ” methode geeft een vierkante matrix terug van relaties van het gegeven “Gen Ontology” object. Deze graaf wordt soms een ‘ge¨ınduceerde’
graaf genoemd.
10.3
10.3.1
Constructie matrix met meest recente data
Inlezen
GO = geneont(’LIVE’,true);
[MATRIX, ID, REL] = getmatrix(GO);
10.3.2
“Deel van”-relatie ⇒ “is een”-relatie
i=1;
j=1;
for i = 1:length(MATRIX)
for j = 1:length(MATRIX)
if MATRIX(i,j) == 2
MATRIX(i,j) = 1;
55
end
j=j+1;
end
i = i+1;
j = 1;
end
10.3.3
Symmetrisch maken van de matrix
i=1;
j=1;
for i = 1:length(MATRIX)
for j = 1:length(MATRIX)
if MATRIX(i,j) == 1
MATRIX(j,i) = 1;
end
j=j+1;
end
i = i+1;
j = 1;
end
10.4
Inlezen willekeurige datasets van de GO website
\subsubsection{Unzippen}
\begin{verbatim}
gunzip(’gene_association3.gz’);
SGDGenes = goannotread(’gene_association3’);
10.4.1
Constructie matrix uit gedownloade informatie
Twee bruikbare vectoren:
genen = [SGDGenes.GOid];
mensgenen = {SGDGenes.DB_Object_Symbol};
Voorouders en weergave in netwerk met het commando biograph:
ancestors = getancestors(GO,5840)
riboanc = GO(ancestors)
cm = getmatrix(riboanc);
BG = biograph(cm, get(riboanc.Terms.’name’))
view(BG)
56

Innovatie Matrix

Brochure Matrix F2

¡ ¡ 5 22 ве ¤ ¡ 2 3 5 9 вз ¦ 120 213 ве ий

MA T RIX F 12 - Frama Nederland

Lees verder - ME/cvs Vereniging

Profeten van lang geleden.mus

MA T RIX F 32 - Frama Nederland

ENDURANCE CUP - MExT Racing Team

Nu ga ik mijn eigen gang

Proefstamboom van een nakomeling uit Synouk Gipsy v.h.