comparison Process/Clustering.py @ 2:46fb79167a61 tip

Main Code

+'''
+Created on 10/11/2014
+
+@organization: Lancaster University & University of Leeds
+@version: 1.0
+Created on 11/12/2014
+
+@author: Victor Padilla
+@contact: [e-mail]
+
+Class related to Phylogenetic trees
+
+'''
+import copy
+import numpy
+from Bio import Phylo
+from cStringIO import StringIO
+
+class cluster:
+    '''
+    cluster object for Clustering class
+    '''
+    pass
+    
+class Clustering:
+    '''
+    This class allows to transform a triangular matrix
+    in a philogenetic tree cutting it at some point 
+    
+    
+    usage:
+    ###  The omr name array
+    species = [ "omr1", "omr2", "omr3", "omr4", "omr5","omr6" ]
+    matr = [ [ 0., 6., 8., 1., 2. ,6. ],
+             [ 0., 0., 8., 6., 6. ,4. ],
+             [ 0., 0., 0., 8., 8. ,8. ],
+             [ 0., 0., 0., 0., 2. ,6. ],
+             [ 0., 0., 0., 0., 0. ,6. ],
+             [ 0., 0., 0., 0., 0. ,0. ], ]
+    
+    clustering=Clustering()
+    ### Transforms the triangular matrix in a complete matrix
+    matr= clustering.getCompleteMatrix(matr)
+    
+    ### establishes the clusters
+    clu = clustering.make_clusters(species)  
+    
+    ### Regroup the cluster 
+    tree = clustering.regroup(clu, matr)
+    
+    ### get the string in newick format
+    strTree=clustering.getStringTree(tree,tree.height,"")
+    
+    print(strTree)
+    clustering.showTree(strTree)
+    
+    ### Takes the main tree with 5 branches
+    maintree=clustering.getBetterTree(tree,5)
+    strTree=clustering.getStringTree(maintree,maintree.height,"")
+    print(strTree)
+    clustering.showTree(strTree)
+        
+    '''
+    
+    def make_clusters(self,species):
+        '''
+        Organises the cluster based on the species array
+        '''
+        clusters = {}
+        Id = 1
+        for s in species:
+            c = cluster()
+            c.id = Id
+            c.data = s
+            c.size = 1
+            c.height = 0
+            clusters[c.id] = c
+            Id = Id + 1
+        return clusters
+    
+    def __find_min(self,clu, d):
+        '''
+        finding the minimum distance
+        private function
+        '''
+        mini = None
+        i_mini = 0
+        j_mini = 0
+        for i in clu:
+            for j in clu:
+                if j>i:
+                    tmp = d[j -1 ][i -1 ]
+                    if not mini:
+                        mini = tmp
+                    if tmp <= mini:
+                        i_mini = i
+                        j_mini = j
+                        mini = tmp
+        return (i_mini, j_mini, mini)
+    
+    def regroup(self,clusters, dist):
+        '''
+        organizing the cluster
+        '''
+        i, j, dij = self.__find_min(clusters, dist)
+        ci = clusters[i]
+        cj = clusters[j]
+        # create new cluster
+        k = cluster()
+        k.id = max(clusters) + 1
+        k.data = (ci, cj)
+        k.size = ci.size + cj.size
+        k.height = dij / 2.
+        # remove clusters
+        del clusters[i]
+        del clusters[j]
+        # compute new distance values and insert them
+        dist.append([])
+        for l in range(0, k.id -1):
+            dist[k.id-1].append(0)
+        for l in clusters:
+            dil = dist[max(i, l) -1][min(i, l) -1]
+            djl = dist[max(j, l) -1][min(j, l) -1]
+            dkl = (dil * ci.size + djl * cj.size) / float (ci.size + cj.size)
+            dist[k.id -1][l-1] = dkl
+        # insert the new cluster
+        clusters[k.id] = k
+    
+        if len(clusters) == 1:
+            # we're through !
+            return clusters.values()[0]
+        else:
+            return self.regroup(clusters, dist)
+    
+    def __pprint(self,tree, length):
+        '''
+        print the tree in newick format 
+        (A:0.1,B:0.2,(C:0.3,D:0.4):0.5);       distances and leaf names
+        '''
+        if tree.size > 1:
+            # it's an internal node
+            print "(",
+            self.__pprint(tree.data[0], tree.height)
+            print ",",
+            self.__pprint(tree.data[1], tree.height)
+            print ("):%2.2f" % (length - tree.height)),
+        else :
+            # it's a leaf
+            print ("%s:%2.2f" % (tree.data, length)),
+    
+    def getStringTree(self,tree, length,strOut):
+        '''
+        returns the string of the tree in newick format
+        (A:0.1,B:0.2,(C:0.3,D:0.4):0.5);       distances and leaf names
+        '''
+        if tree.size > 1:
+            # it's an internal node
+            strOut+="("
+            strOut=self.getStringTree(tree.data[0], tree.height,strOut)
+            strOut+=","
+            strOut=self.getStringTree(tree.data[1], tree.height,strOut)
+            strOut+="):"+str(length- tree.height)
+        else :
+            # it's a leaf
+            strOut+=str(tree.data)+":"+str(length)
+        return strOut
+    
+    
+    def createClusters(self,tree,length):
+        '''
+        separates the tree at the point determined by length
+        '''
+        clusters=[]
+        self.__separateClusters(clusters,tree,length)
+        self.clusters=self.__filterClusters(clusters,length)
+        return clusters
+    
+    def getLeafs(self,tree):
+        '''
+        returns the number of leafs from a tree
+        '''
+        leafs=[]
+        self.getTreeLeafs(leafs,tree,tree.height)
+        return leafs
+    
+    def __separateClusters(self,clusters,tree,length):
+        '''
+        divides the main tree in multiples trees based on the length
+        '''
+        if tree.size > 1:
+            if tree.height>=length:
+                clusters.append(tree.data[0])
+                clusters.append(tree.data[1])
+            else:
+                pass
+            self.__separateClusters(clusters,tree.data[0], length) 
+            self.__separateClusters(clusters,tree.data[1], length)
+        
+    def __filterClusters(self,clusters,length): 
+        '''
+        removes trees higher than length
+        ''' 
+        clustersOut=[]  
+        for c in clusters:
+            if c.height<length:
+                clustersOut.append(c)
+        return clustersOut
+    
+    def getTreeLeafs(self,leafs,tree,length):
+        '''
+        returns the number of leafs from a tree. Recursive function
+        '''
+        if tree.size > 1:
+            # it's an internal node
+            self.getTreeLeafs(leafs,tree.data[0], tree.height)
+            self.getTreeLeafs(leafs,tree.data[1], tree.height)
+        else :
+            # it's a leaf
+            leafs.append(tree.data)
+            
+    def getBetterTree(self,tree,numberLeafs):
+        '''
+        takes the smallest tree that has <= numberLeafs
+        '''
+        while True:
+            leafs=self.getLeafs(tree)
+            if len(leafs)<=numberLeafs:
+                return tree
+            clusters=self.createClusters(tree,tree.height)
+            tree=self.getMainTree(clusters)
+    
+    
+    def getAverageDistance(self,matr):
+        '''
+        returns the average distance from the matrix
+        '''
+        sumat=0
+        lengthMatrix=len(matr[0])
+        length=lengthMatrix*lengthMatrix-lengthMatrix
+        for i in range(lengthMatrix):
+            for j in range(lengthMatrix):
+                    sumat=sumat+matr[i][j]
+        return (sumat/length)/2
+    
+    def getMaximumDistance(self,matr):
+        '''
+        returns the maximum distance from the matrix
+        '''
+        maxValue=0
+        lengthMatrix=len(matr[0])
+        for i in range(lengthMatrix):
+            for j in range(lengthMatrix):
+                    if matr[i][j]>maxValue:
+                        maxValue=matr[i][j]
+        return maxValue/2
+    
+    def getMainTree(self,clusters):
+        '''
+        Takes the tree with the higher number of leafs
+        '''
+        maxim=0
+        indexmax=0
+        for i in range(len(clusters)):
+            leafs=self.getLeafs(clusters[i])
+            if len(leafs)>maxim:
+                maxim=len(leafs)
+                indexmax=i
+                
+        return clusters[indexmax]
+            
+    def getCompleteMatrix(self,matr):
+        '''
+        return the complete matrix from a 
+        triangular matrix
+        '''
+        newMatrix=copy.copy(matr)
+        if isinstance(matr,numpy.ndarray):
+            newMatrix=newMatrix.tolist()
+        lengthMatrix=len(matr[0])
+        for i in range(lengthMatrix):
+            for j in range(i,lengthMatrix): 
+                newMatrix[j][i]=matr[i][j]
+        return newMatrix
+    
+    def showTree(self,strTree): 
+        '''
+        Shows a graphical representation from a newick tree format
+        '''  
+        handle = StringIO(strTree)
+        tree = Phylo.read(handle, 'newick')
+        tree.ladderize()   
+    #     Phylo.draw(tree) 
+        Phylo.draw_ascii(tree) 
+         
+