Additive Phylogeny (no project to choose as WPB project)

Trimming the Tree

When we look at Figure 2 we can recognize that next we need to come to the tree D_trimmed. This Tree is different to D_bald, because we reduce the matrix and remove the changed columns. After we again construct Dbald and remove again column / row until we have a 2*2 matrix. The function to remove columns / rows can be found in the Fsharp.Stats library and is shown in the next Codeblock

//Hier wird die veränderte Reihe entfernt, da diese dem SpecIndex entspricht und die Matrix damit verringert
//Bsp. delete_species (dBaldMatrix (myMatrixAdditive.DistanceM) 1) führt zu 
//       0 21 22
//       21 0  13
//       22 13 0
  
let delete_species (distMat:matrix) (specIndex: int)  =
    distMat
    |> Matrix.removeColAt specIndex
    |> Matrix.removeRowAt specIndex 

Finding the Attachmentpoint of a leaf

Because we create the tree by Adding leaves backwards, we need a further function, we can prepare. This function tries to identify the point of attachment of leaf y. To find the attachment point in the trimmed tree, you should "consider that Tree D^bald, is the same as Tree (D), except that Limblength(j) = 0" (Compeau and Pevnzer, Bioinformatics Algorithms: An Active Learning Approach). because we know from the Limblength theorom, that is integrated in dBald that leaves i and k muist be so that limblength(j)=0. This leads to the following assumption: The attachmentpoint of a leaf j has to be located at distance_i,j^bald on a path connecting node i with node k in the trimmed tree.
When the attachmentpoint is at an existing node, we connect j to this node. otherwise if it occurs at an edge, we insert a new internal node at distance and connect j to it. This explanation can now be translated in coding language.

//Hier wird der Attachmentpoint einer bestimmten leaf distinctLeaf bestimmt. 
//Dafür wird eine dBaldMatrix bzw. trimmed matrix, sowie die leaf von Interesse als Input benötigt sowie die Anzahl der SpecIndices bestimmt

let attachmentPoint (dBald: matrix) (distinctLeaf: int) =
    // entspricht der Anzahl an Specindices z.B. in einer 29 * 29 Matrix = 28
    let NumberOfSpecies = dBald.NumRows-1

      
    //in dieser verschachtelten For loop werden alle möglichen Attachmentpoints für eine bestimmte leaf bestimmt und in einem Array gespeichert 
    //z.B. für 29*29 matrix 0,26 = 8142 bis 26,27 = 868
     
    [|for i = 0 to NumberOfSpecies-1 do 
        for j = i+1 to NumberOfSpecies-1 do
            if dBald.[i,j] = (dBald.[i,distinctLeaf] + dBald.[distinctLeaf,j]) then 
                   (i,j),dBald.[i,(distinctLeaf)]
    |]
     
    // Anschließend wird das Array in kleinere Arrays groupiert, wobei die Gruppierung nach gleicher Menge (i) erfolgt --> alle mit 0 in eine Gruppe usw.
    |> Array.groupBy (fun ((i,j),d) -> i)
    // Es wird das maximale Element herasugefiltert und dabei das erste Element desa Tupels genutzt, in diesem Fall z.B. 26
    |> Array.maxBy fst
    // nun wird auf alle Elemente dieser Gruppe zugegriffen und das erste Element behalten z.B hier (26,27), 868 --> 
    // diese sagt nun das auf der Strecke zwischen Knoten 26 und Knoten 27 auf der Distanz 868 ein weiterer Knoten liegen muss
    |> snd
    |> Array.head
    
// Test of Attachmentpoint

attachmentPoint (dBaldMatrix myMatrixAdditive.DistanceM 3) 3 // (1,2 -->6)

((attachmentPoint ( dBaldMatrix twentyNineMatrix.DistanceM 28)) 28) // (26,27--->868)

Creating a algorithm for distance-based phylogeny construction

After we have designed all essential functions for the additive phylogeny mechanism we can create an Algorithm for the finding of a simple tree fitting the n*n distance matrix D. For the creation of the Algorithm you can use either Figure 2 or the Pseudocode in Fig 5.

But before we start to design the mechanism, we first create a type, where we can store the pathes and have a list which describes the created tree. This type is called PhylTreeConnection and contains the members SourceToTargetIndex and Distance. SourcetotargetIndex describes, which nodes are connected and distance how the distance is between the nodes. For example Fig.5. which shows a list how it should look like in the end

// Hier wird ein Typ deklariert, welche die Darstellung des phylogenetischen Baum am Ende ermöglicht 
type PhylTreeConnection = { 
    SourceToTargetIndex: (int*int)
    Distance: float
    } with 
        static member Create sourceToTargetIndex distance = { SourceToTargetIndex= sourceToTargetIndex; Distance=distance} 

Now we are able to create the phylogenetic tree with help of combinating the functions. The steps of the additive phylogeny Mechanism are described in the following steps:

• addPhyl is a function, which gets only the distancematrix as input and you proof first if the matrix is additive, because the function works only for additive matrices

• Now we define some variables of the algorithm, that you need during the process:

  • distmat is a mutable variable, that you need only as storage for the actual version of distanceMatrix,because you have to reduce this matrix until its a 2*2 Matrix
  • n is again a mutable variable, that has to be always the last species (Row) of the Matrix and needs to be updated during algorithm
  • lengthOfMatrix is a nonmutable variable, that stores the size of the input matrix
  • resultList is a mutable variable, that is in the beginning empty and is needed as List to store the pathes at the end
  • graph and attachmentsStore are two dictionaries, that stores all limbLength and distancematrices until we reach a length of 3*3 (graph) as well as the computed attachmentpoints for each matrixtype --> both habe the same length
  • internalNodes is a mutable variable and describes the index of inserted nodes. It starts with the highest point and ends with the smallest

• After we have defined all variables we can start with the algorithm, that in general consist of an if else clause. The if Part first proofs if we habe only a 2*2 matrix.

  • If the size of distancematrix is 2*2, then the maximal Specindex n is 1 and we have directly the simple tree consisting of an single edge connecting both species and having size equal to the distance shown in the distancematrix

• if we have no 2*2 Matrix then we have a quite more complex algorithm and we need the else part

  • first we need to trimm the Matrix until its a 2*2 Matrix but store the updated matrix in each step. For each matrix (n), we also need to store the attachmentpoint as well. Here we use the defined dictionaries attachmentpoint and graph.

  • These dictionaries are filled in a simple while loop that has as termination condition size of matrix 2*2 (n < 2). The procedure of filling the dictionary is done in different steps

    • First we fill for every round during loop graph by computing the corresponding limblength and the dBald matrix
    • Then we can fill another dictionary that contains the attacghmentspoints as well as source and target
    • Now we trimm the distancematrix in the same while loop and update n. The trimminmg is done be remove every time the last row / column in the matrix and n has always to be the specindex of the last row / column

  • After we have succesful reduced the the distancematrix iterative and stored the corresponding limblength and attachmentspoints in graph / attachmentsstore, we can start now with the Traceback, where the initial graph consisting of node 0 and 1 is updated by insering all corresponding nodes from attachmentsstore together with limblength stored in graph

    • For the Traceback we first update n, because we need n equal to 2 and we need to insert the single path of the 2*2 matrix being the starting point in resultlist, representing the tree. For this source is 0 , target is 1 and the distance is the distance represented in distancematrix.[0,1]. Because of symmetry source and target can also be changed and have the same distance. In this case target is 0 and source 1.

    • Then we have to choose our first dBaldMatrix and the corresponding limblength and store this in a tuple (limblength,dMatix). For this we use n and assign the entry in graph to n
    • In another tuple we store the first attachmentpoint as well as the corresponding src and target stored in attachmentsStore
    • Next, we need to find the path in the actual tree that has to be modified. This path that has to be modified is assigned to edgeToModify. Attention: This Step is actually errorprone. Becuase we search after the first step in resultList for a path 1,2. But there exists no such Path. For this reason we need to find pathes were 1 is the src and we need to find a path were 2 is the target. This path are in this case (1,55) and (55,2). Then we must look at the distances from this pathes 2291 and 4032. Because the attachmentpoint between 1 and 2 is at 2678 we know that the attachmentpoint now has to be inserted on the path between 55,2 at distance size 387. This step is not working actually, becuase it actually searchs only after pathes being completely actual
    • if we have found the path to modify in the resultList, we also need to assign the inverse path, because both pathes are part in the description of the simple tree. For that we only need to change the value of src in the value of target and target in src. The Distance is because of symmetry the same

    • Next we need to define, which pathes have to be added to the path. In general this are every time 2*3 different edges called A as well as their inverse path

      • edgeAdditionA and edgeAdditionAinverse describe the Path from the src to internalNode and the length between them e.g. (0,55 with Distance 745)
      • edgeAdditionB and edgeAdditionBinverse describe the path from the internal Node to the target and the Distance which is equal to the edgetomodify Distance - the computed Distance from attachmentpoint
      • edgeAdditionC and edgeAdditionCinverse describe the Step where we add the next species (n) back to the tree (e.g.55 -2). This Path has the length limblength
    • After Definition we update internalNodes and the resultList. While internalnodes has to be reduced by one, we have to remove the modified edges from resultList and add the computed edges, also n has to be increased by one
    • To improve clarity we also have to sort the list by the Srcindex and then we can return the unique tree fitting the additive distancematrix

let addPhyl (distMat:matrix) =
    //Hier wird getestet ob die in die Funktion gegebene distanzmatrix die Anforderungen einer additiven Matrix erfüllt, 
    //da der Algorithmus nur für eine additive Matrix funktionsfähig ist. Dafür wird die Funktion testingAdditivity aufgerufen.
    
    testingAdditivity distMat

    // Hier werden verschiedene Variablen (zum Teil veränderbar) deklariert. Diese werden im Laufe des Algorithmus benötigt und müssen zum Teil veränderbar sein

    // veränderbare Variable, welche die aktuelle Version der Matrix beinhaltet, da die Matrix verändert werden muss (Stichwort Dbald)
    let mutable dist_Mat = distMat

    // veränderliche Varialble, die immer auf die letzte Reihe / Spalte der aktuellen Matrix verweist  
    let mutable n = dist_Mat.NumRows-1 

    // Unveränderliche Variable, welche auf die Länge der InputMatrix verweist
    let lengthOfMatrix = distMat.NumRows-1

    // Hier wird eine leere Liste vom Typ PhylTreeConnection deklariert in der die Pfade am Ende gelistet sind und der Tree beschrieben wird 
    let mutable resultList: PhylTreeConnection list = []

    //Hier werden twei mutable dictionaries deklariert. graph ist dafür notwendig, alle Matrixen zu speichern und damit rückwärts auflösen zu können  
    //und enthält den Wert von n als Schlüssel ebensp wie der attachmentsstore der src,target und den attachmentPoint speichert. 
    
    let graph = new Dictionary<int,(float*matrix) >()
    let attachmentsStore = new Dictionary<int,((int*int)*float)>()

    // internalNodes ist die variable, welche den Index des einzufügenden Knoten deklariert
    let mutable internalNodes = 2*dist_Mat.NumRows-3
    
    //  Hier wird geprüft, ob eine 2*2 Matrix vorliegt, oder eine größere Matrix.
    //     Wenn eine 2*2 Matrix vorliegt, dann muss der Algorithmus nicht weiter durchlaufen werden 
    //     und es kann direkt der phylogenetische Baum, bestehend aus einer einfachen Verbindung der beiden Spezies erstellt werden.
    //     Die Länge der Verbindung entspricht der Distanz zwischen den beiden Spezies. Folglich wird  der Baum in der Liste 
    //     dargestellt durch den Eintrag von Source to Target (0,1) bzw. (1,0) und der Distanz zwischen beiden Spezies und die Liste wird zurückgegeben   
    
    if n = 1 then
        let distance = dist_Mat.[0,1]
        resultList <- {SourceToTargetIndex = (0,1); Distance = distance}::resultList
        resultList <- {SourceToTargetIndex = (1,0); Distance = distance}::resultList
        resultList

    // Diese Verzweigung wird ausgeführt falls keine 2*2 Matrix vorliegt
    
    else

     
        // Hier werden die in Zeile 24 und 25 deklarieren Dictionaries graph und attachmentsStore solange gefüllt bis eine Matrixgröße 2+2 erreicht ist.
        // Das Dictionary graph beinhaltet dabei als Schlüssel n und ordnet diesem die jeweilige limbLength sowie die Distanzmatrix zu aus welcher die aktuelle Pfad Länge entnommen werden kann
        // Im Dictionary attachmentsStore werden die berechneten attachmentpoints von zwei Spezies gespeichert sowie src und target
    
        while n >= 2 do
            let limblength = LimbLength (dist_Mat) n
            dist_Mat <- dBaldMatrix dist_Mat n 
            graph.Add (n,(limblength,dist_Mat))

            //    printfn "n: %i" n
            //     printfn "limblength:  %f %A" limblength dist_Mat
            //     Beispielhafte Ausgabe wenn n = 3 erreicht ist, zur Ansicht was in dictionary graph gespeichert wird
             
            //  n: 3
            //  limblength: 409.000000 
            //         0        1        2        3
                                                
            // 0 ->    0.000 3036.000 4777.000 1132.000
            // 1 -> 3036.000    0.000 6323.000 2678.000
            // 2 -> 4777.000 6323.000    0.000 3645.000
            // 3 -> 1132.000 2678.000 3645.000    0.000
             

            let attachment = attachmentPoint dist_Mat n
            let x = snd attachment 
            let (src,target) = fst attachment
            
            attachmentsStore.Add(n,((src,target),x))
            //  printfn " n: %i x: %f src: %i target:%i" n x src target 
            // Beispielshafte Ausgabe wenn n = 3 erreicht ist in der 29*29 Matrix aus Rosalind
            // n: 3 x: 2678.000000 src: 1 target:2
            
                
            // Hier wird die Matrix getrimmt und n geupdatet sodass n wieder dem Index der letzten Reihe / Spalte entspricht

            dist_Mat <- delete_species dist_Mat n
            n <- n-1
        
        
        // Bis hier wurde die distanzmatrix iterativ reduziert und die jeweiligen Attachmentpoints und limblength in graph und attachmentsStore gespeichert
        // Nun muss im TraceBack-Schritt in jeder Iteration ein initial graph bnestehen aus node 0 und 1 mit allen fehleden nodes vervollständigt werden. Diese werden 
        // an den punkten aus attachmentstore angefügt mit der limblength, die in graph gespeichert ist. 
        // Dafür wird zunächst der aktuelle Pfad der in der Matrix als Schnittpunkt 0,1 dargestellt ist in resultList eingefügt 
        

        resultList <- {SourceToTargetIndex = (0,1); Distance = distMat.[0,1]}::resultList
        resultList <- {SourceToTargetIndex = (1,0); Distance = distMat.[1,0]}::resultList

        // Hier wird n wieder angepasst sodass n = 2 ist und damit dem kleinsten Schlüssel in den Dictionaries entspricht
        n <- n+1

         
            // Im folgenden Schritt wird der traceback durchgeführt und nach und nach die bestimmten attachmentpoints eingefügt. 
            // Dafür wird eine zweite While loop erstellt die solange wie n nicht dem maximalen Specindex entspricht den Traceback durchführt.
        
        let traceback =
            while n < lengthOfMatrix do (

            
            //    Hier wird die aktuelle resultList für jeden Durchgang geprintet:  printfn "resultList: %A" resultList, Bsp. zu Beginn und nach dem ersten Durchlauf  
            //     resultList zu Beginn: [{ SourceToTargetIndex = (1, 0) Distance = 3036.0 }; { SourceToTargetIndex = (0, 1) Distance = 3036.0 }]
            //     resultList nach Durchgang 1:[{ SourceToTargetIndex = (0, 55) Distance = 745.0 }; { SourceToTargetIndex = (55, 0) Distance = 745.0 }; 
            //                                  { SourceToTargetIndex = (55, 1) Distance = 2291.0 };
            //                                  { SourceToTargetIndex = (1, 55) Distance = 2291.0 }; { SourceToTargetIndex = (55, 2) Distance = 4032.0 }; 
            //                                  { SourceToTargetIndex = (2, 55) Distance = 4032.0 }]
              


                // Take the currentGraphElement
                let (limblength,dMatrix) = graph.[n]
                // Take the currentAttachmentPoint
                let ((src,trg),distance) = attachmentsStore.[n]

        
                //  Hier kann Beispielshaft durch den Print Befehl überprüft werden; ob der  Attachmentpoint richtig übertragen wurde - Ausgabe für twentyNineMatrix

                //     printfn "index: %i, src: %i, trg: %i, distance: %f"n src trg distance
                //     index: 2, src: 0, trg: 1, distance: 745.000000
                //     index: 3, src: 1, trg: 2, distance: 2678.000000

                
                //  Ab hier ist der Code noch nicht vollständig funcktionsfähig, da aktuell nur nach edgestomodify gesucht wird die genau dem src, target entsprechen, dies ist allerdings nicht ausreichend, da vorher bereits Pfade verändert wurden und im Pfad zwischen 1 und 2 z.B. schon eine internalNode eingefügt wurde
                //     Problem: Es soll nodeindex 3 auf dem path zwischen nodes 1 und 2 angefügt werden, aber dieser path ist nicht verfügbar weil im vorigsen schritt schon 
                //     ein artificail node (55) zwischen 1 und 2 eingefügt wurde. Deswegen müsste man resultlist jetzt untersuchen nach dem Path 1 -> 2 und die edge identifizieren an der 
                //     der zusätzliche node (54) angefügt werden soll (aufgrund der distanzen 1-55: 2291 und 55-2: 4032 müsste node 54 demnach zwischen 55 und 2 in der Distanz 387 eingesetzt werden

                let edgeToModify = 
                    resultList |> List.find (fun element -> 
                        element.SourceToTargetIndex = (src,trg) ||
                        element.SourceToTargetIndex = (trg,src) 
                        )

                // Hier soll der inverse Path von edge to Modify (Bsp. von 0,1 => 1,0) bestimmt und zugeordnet werden, da dieser auch aus der Liste entnommen werden muss
                let edgeToModifyInnverse = {SourceToTargetIndex = (snd (edgeToModify.SourceToTargetIndex),fst (edgeToModify.SourceToTargetIndex)); Distance = edgeToModify.Distance}
               
                //  Hier wird das im vorherigen Codeabschnitt bestimmte edgeToModify und der inverse Pfad bestimmt
                //     printfn " edgetomodify: %A , edgetomodify inverse: %A"     edgeToModify edgeToModifyInnverse 
                //     edgetomodify: { SourceToTargetIndex = (1, 0) Distance = 3036.0 } , edgetomodify inverse: { SourceToTargetIndex = (0, 1) Distance = 3036.0 }
                

                // Hier werden die Pfade bestimmt, die in resultList hinzugefügt werden sollen und an welcher Stelle ein Knoten eingefügt werden soll
                //     edgeAdditionA: Beschreibt den Pfad zwischen der Quelle und dem einzufügnden internal Node (z.B. 0 -55 ==> Distance 745 oder 55-54 ==>  distance 387 )
                //     edgeAdditionB: Beschreibt den Pfad zwischen InternalNodes und dem target also von z.B. 55-1 ==> 2291 oder 54 -2 ==> 1354. 
                //     Dafür muss die edgeToModiy Distanz von den einzelnen Distanzen ( Distanzen der beteiligten Teilpfade / Attachmentpoindistanz) abgezogen werden.
                //     edgeAdditionC: Hier wird die aktuelle Spezies wieder zum Baum hinzugefügt im Abstand der limblength zum aktuellen InternalNode ( z.B. im ersten Schritt wird Spezies 2 wieder eingefüht)
                
                
                let edgeAdditionA = {SourceToTargetIndex = (src,internalNodes); Distance = distance}
                let edgeAdditionA_inverse = {SourceToTargetIndex = (internalNodes,src); Distance = distance}
                let edgeAdditionB = {SourceToTargetIndex = (internalNodes,trg); Distance = edgeToModify.Distance - distance}
                let edgeAdditionB_inverse = {SourceToTargetIndex = (trg, internalNodes); Distance = edgeToModify.Distance - distance}
                let edgeAdditionC = {SourceToTargetIndex = (internalNodes,n); Distance = limblength}
                let edgeAdditionC_inverse = {SourceToTargetIndex = (n,internalNodes); Distance = limblength}

                // Beispielshafte Ausgabe nach dem ersten Durchlauf von edgeAddition A , edgeAddition B und edgeAddition C 
                //     Edgeaddition A: { SourceToTargetIndex = (0, 55) Distance = 745.0 } 
                //     edgeAdditionB: { SourceToTargetIndex = (55, 1) Distance = 2291.0 } 
                //     edgeAdditionC: { SourceToTargetIndex = (55, 2) Distance = 4032.0 } 
                                                                                     
                // Aktualisieren von internal Nodes 

                internalNodes <- internalNodes-1

                //  Hier wird der phylogenetische Baum  geupdated. 
                //      Dafür werden die zu verändernden Pfade zunächst aus der Liste entferent und anschließend werden die vorher bestimmten Pfade der Liste hinzugefügt.
                //      Zuletzt wird der übersichtlichkeithalber die Liste nach den Src sotiert und n geupdatet. Nach Beendigung der While Loop wird zuletzt die fertige resultList, 
                //      welche den Baum beschreibt ausgegeben 
              

                let newResultList = 
                    resultList
                    |> List.where  (fun x -> x <> edgeToModify && x <> edgeToModifyInnverse)
                    |> List.append [edgeAdditionA;edgeAdditionA_inverse;edgeAdditionB;edgeAdditionB_inverse;edgeAdditionC;edgeAdditionC_inverse]
                    |>List.sortBy (fun x -> fst (x.SourceToTargetIndex))
                resultList <- newResultList
                n <- n+1
                )

        resultList
      
        
// Beispielshafte Ausgabe zur Darstellung des Problems
// printfn "resultList: %A" resultList
// printfn "index: %i, src: %i, trg: %i, distance: %f"n src trg distance

// resultList: [{ SourceToTargetIndex = (1, 0)
//    Distance = 3036.0 }; { SourceToTargetIndex = (0, 1)
//                           Distance = 3036.0 }]
// index: 2, src: 0, trg: 1, distance: 745.000000
// resultList: [{ SourceToTargetIndex = (0, 55)
//    Distance = 745.0 }; { SourceToTargetIndex = (55, 0)
//                          Distance = 745.0 }; { SourceToTargetIndex = (55, 1)
//                                                Distance = 2291.0 };
//  { SourceToTargetIndex = (1, 55)
//    Distance = 2291.0 }; { SourceToTargetIndex = (55, 2)
//                           Distance = 4032.0 }; { SourceToTargetIndex = (2, 55)
//                                                  Distance = 4032.0 }]
// index: 3, src: 1, trg: 2, distance: 2678.000000



// test off addPhyl

addPhyl twentyNineMatrix.DistanceM

No value returned by any evaluator