Problem B
Labb F3: Evolutionära träd och I/O

Det här är en laboration som dels testar din förmåga att arbeta funktionellt genom att överföra en abstrakt beskriven algoritm till ett funktionellt program, och dels låter dig skriva ett fullständigt Haskell-program med inläsning och utskrifter – operationer som har sido-effekter. Vi fortsätter arbeta med problem från molekylärbiologin och det är meningen att det du skriver ska bygga på laboration F2 – du kan antagligen med fördel använda din lösning från F2 utan ändringar, och skriva din lösning på denna labb i en ny modul F3 som importerar F2.

1 Bakgrund

Det finns många anledningar till att vara intresserad av hur arter och/eller gener har uppkommit, bland annat är det för många viktigt att helt enkelt förstå hur olika arter har uppstått. Inom medicin kan kunskap om geners utveckling ge kunskap om hur de fungerar och vilken funktion de har. En grundläggande fråga är då hur man rekonstruerar det evolutionära träd, en fylogeni, som gav upphov till sekvenserna? Figur 1 ger ett exempel på en fylogeni. Indata är alltså en mängd sekvenser, DNA eller protein, och utdata är ett träd där alla inre hörn har grad 3.

Figure 1: Exempelfylogeni för en hypotetisk gen funnen i fem arter. Lägg märke till att det här trädet inte är avsett att påstå något om var evolutionen började, dvs vilken punkt i trädet som är äldst. Man säger att trädet är orotat. (Med vår övriga kunskap om de olika arterna kan vi dock vara ganska säkra på att en eventuell rot skulle ligga på kanten mellan insekterna och fisken i den här specifika fylogenin).

2 Algoritmen Neighbor Joining

Den vanligaste och mest kända algoritmen för att återskapa träd givet avståndsdata är Neighbor Joining (NJ). Det är en iterativ algoritm som successivt väljer ut par av löv och slår ihop dem till ett nytt hörn. Vilket par man väljer är avgörande för att resultatet ska bli bra, och NJ garanterar faktiskt inte att det är det bästa trädet (i meningen: passar bäst med avstånden) som returneras: NJ är en girig heuristik.

Låt $F_1$ vara den mängd hörn som ska vara löv i det träd $T$ vi bygger. Låt $D_1$ vara avståndsmatrisen över $F_1$. I varje iteration $i$ av NJ kommer vi att välja ut två element $a$, $b$ i $F_ i$ och kombinera ihop dessa till ett träd. Detta skapar en ny mängd $F_{i+1}$ och avståndsmatris $D_{i+1}$.

Urvalsfunktionen $S$

Vi väljer de två element $a,b \in F_ i$ som minimerar följande urvalsfunktionen $S$:

Funktionen ser vid en första anblick en smula underlig ut, men man kan göra en tolkning av den. Den andra termen mäter hur långt bort övriga hörn ligger från $x$ och $y$, och den första termen mäter hur nära $x$ och $y$ ligger varandra, skalat med en faktor för att göra de två termerna jämförbara. Det $S$ kan sägas välja ut är alltså de två hörn som ligger längst ifrån de andra.

Pseudokod

Vi använder oss av en förenklad version av NJ. Den som tittar på andra beskrivningar av NJ kommer att finna att steg 4 är lite hårigare än vad som ges här. Låt $F_1$ och $D_1$ vara indata.

1. $i \leftarrow 1$

2. Så länge $|F_ i| > 3$:

   1. Hitta det par $a, b \in F_ i$ som minimerar $S_ i(a, b)$

   2. Skapa ett nytt träd $T_ i$ där träden $a$ och $b$ är barn till en nyskapad nod.

   3. $F_{i+1} \leftarrow F_ i \cup \{ T_ i \} \setminus \{ a, b \} $ – Lägg till det nya trädet till $F$ och ta bort de gamla

   4. Skapa $D_{i+1}$ från $D_ i$ enligt

      \begin{align*} D_{i+1}(x, y) & = D_ i(x, y) & & \text {för $x, y \in F_{i+1} \setminus \{ T_ i\} $} \\ D_{i+1}(x, T_ i) = D_{i+1}(T_ i, x) & = \frac{D_ i(x, a) + D_ i(x, b)}{2} & & \text {för $x \in F_{i+1} \setminus \{ T_ i\} $} \end{align*}

   5. $i \leftarrow i+1$

3. Skapa ett nytt träd $T$ där de tre kvarvarande träden i $F_ i$ är barn till en nyskapad nod

4. Returnera $T$

Exempelkörning

Antag att vi har de fem sekvenserna $a, b, \ldots , e$ som beskrivs i exempel-data 1 nedan. Använder vi vår kod från F2 för att beräkna avståndsmatrisen $D_1$ för dessa får vi:

I algoritmens första iteration har vi “träden” $F_1 = \{ a,b,c,d,e\} $, alla bestående av en enda nod. Urvalsfunktionen $S_1$ för första iterationen får följande värden:

Den minimeras alltså av paret $(a,b)$, så vi skapar ett nytt träd $T_1$ som består av en ny nod med $a$ och $b$ som barn. Den nya avståndsmatrisen $D_2$ över de aktiva träden $F_2 = \{ T_1, c,d,e\} $ och nya urvalsfunktionen $S_2$ kommer se ut som följer:

Urvalsfunktionen $S_2$ minimeras av endera paret $(c,d)$ eller paret $(T_1, e)$, och vi kan välja vilket som helst av dem. Låt oss säga att vi väljer det senare, dvs $(T_1, e)$. Vi bildar då ett träd $T_2$ som består av en ny nod med $T_1$ och $e$ som barn. Vi har nu bara tre aktiva träd kvar $F_3 = \{ T_2, c, d\} $, och vi går till steg 3 i algoritmen och skapar trädet $T$ bestående av en ny nod med $T_2$, $c$, och $d$ som barn. Figur 2 illustrerar de olika träden som byggs upp under körningen.

Figure 2: Träden som byggs upp när algoritmen kör på första exempel-fallet

Hade vi i iteration 2 istället valt att para ihop $(c,d)$ istället för $(T_1,e)$ hade vi fått resultatet som visas i figur 3. Eftersom vi betraktar träden som orotade så är slutresultatet i själva verket samma träd som $T$ från figur 2 även om det ritats annorlunda pga att vi kopplade ihop noderna i en annan ordning.

Figure 3: Alternativ körning på första exempel-fallet

3 Indata

Ditt program ska läsa en uppsättning DNA-sekvenser på standard input (till end of file). Vi kommer använda ett enkelt format där varje sekvens beskrivs av en rad med namnet på sekvensen följt av själva sekvensen, åtskiljda av mellanslag (varken namnet eller sekvensen kommer att innehålla några mellanslag).

Du kan anta följande om indata-strömmen:

• Alla sekvenser kommer att ha samma längd.

• Den innehåller minst $3$ och högst $32$ sekvenser.

• Varje sekvens är minst $1$ och högst $1000$ tecken lång.

• Namnen består bara av tecknen ’a’-’z’, ’A’-’Z’, och ’0’-’9’, och är mellan 1 och 15 tecken långa.

• Sekvenserna är DNA-sekvenser dvs består bara av tecknen ACGT.

4 Utdata

Ditt program ska skriva ut ett evolutionärt träd på standard output. Ett evolutionärt träd (eller fylogeni), skrivs ofta på ett format som kallas Newick. Det träd som visas i figur 2 kan då se ut så här:

(((a,b),e),c,d)

Så länge trädet saknar rot kan man skriva trädet på flera ekvivalenta sätt. Detta träd kan alltså även skrivas på följande sätt:

((a,b),e,(c,d))
(a,b,((c,d),e))
(b,a,((c,d),e))

och många fler. Ett löv i trädet representeras alltså av en sträng med lövets/sekvensens namn. Ett inre hörn i trädet representeras med hjälp av parenteser runt om två komma-åtskilda delträd. På den översta nivån använder vi parenteser runt tre delträd.

Alla ekvivalenta sätt att formatera trädet kommer att godkännas.

Tips

Börja med I/O-delen: skriva ett program som läser in indatasekvenserna, beräknar deras avståndsmatris ($D_1$) med hjälp av din lösning från F2, och sedan skriver ut avståndsmatrisen. Håll delarna av programmet som är inkapslade i IO-monaden minimala – du ska inte behöva gå in i din kod från F2 och ändra den!

I den här uppgiften är det bra att använda ghc för att kompilera Haskell-koden till ett körbart program (utöver att använda ghci för att testa olika funktioner i ert program). I terminalen i Unix kan man sedan använda omdirigering av standard input för att skicka en fil till programmet. Om ni t.ex. har sparat ned första exempel-fallet till filen sample1.in och kompilerat ert program till en körbar binär Main så ska ni i terminalen kunna skriva

> ./Main < sample1.in

för att köra ert program på det aktuella indatat.

Implementera sedan något sätt att representera träden som används i algoritmen, och utskrift av dessa. Modulerna Data.Either eller Data.Maybe kan vara behändiga här – kolla upp dessa!

Med kringarbetet avklarat kan man gripa sig an själva algoritmen. Du vill nog ha något sätt att representera aktuellt tillstånd i algoritmen (mängden $F_ i$ av träd vi har kvar att knyta ihop och avståndsmatrisen $D_ i$), och en funktion som givet avståndsmatrisen $D_ i$ och två element $x, y \in F_ i$ beräknar urvalsfunktionen $S_ i(x, y)$.

Det är antagligen bra att följa den ovan givna algoritmens struktur ganska nära. Indata kommer vara relativt litet, högst $32$ sekvenser, vilket innebär att din implementation inte behöver vara speciellt effektiv och du kan använda naiva lösningar till de olika delproblemen. Om du vill är det såklart tillåtet att göra en smartare, snabbare implementation – modulerna Data.Set och Data.Map kan vara en bra början till detta.

a AACCGGTT
b AACCGGGG
c AATTTTTT
d AACTTTTT
e AACTTTTG

(((a,b),e),c,d)

C CCCCCCCCCCCC
Go GGGGGGGGGGGG
Python GAGGCACCGGGG
Java ACACAACCCCCC
Prolog TTTCATCTTTTT
Haskell CGCGCGTTTGTT

((Prolog,Haskell),(Go,Python),(C,Java))

1 AAA
2 CCC
3 GGG

(1,2,3)