基于python的非加权分组平均法构造简单系统发生树（DNA）

---

能实现什么

1.完成用户自定义输入DNA序列个数及序列中碱基排列
2.根据用户输入的序列构造系统发生树，该树结构存储于列表中
3.使用matplotlib将树结构可视化

一、实现步骤

1.算法思想描述

找出所有序列中距离（这里的距离是由两序列中碱基不同的个数决定的，有多少个不同的碱基距离就为多少）最小的两条合并成为一个新的节点，重新计算其与其他序列的距离，形成新的距离矩阵，重复上述步骤直到所有序列都合并完毕。下图能更加直观的演示树的构造过程

2.代码实现过程

a.首先计算距离矩阵，考虑到碱基存在缺失和增添的情况，所以先对不等长度（只对不等长度的）的序列进行全局序列比对之后再计算其不等的碱基数，计算距离函数单独写在globalAlignment.py文件中，函数代码如下（其中使用到的全局比对算法参考https://blog.csdn.net/mmqwqf/article/details/108922579）

def getDistance(s,t):
    distance = 0  # 两个序列间的距离
    if len(s) == len(t):              #如果长度相等，不进行全局比对
        for i in range(len(s)):
            if s[i] != t[i]:
                distance += 1
    else:
        createList()
        gap = -5
        a = len(s)                      #获取s的长度
        b = len(t)                      #获取t的长度
        mark,direction = createMark(s,t,gap)
        remount(mark,direction,a,b,s,t) #调用回溯函数 进行全局比对
        c = a if a > b else b          #判断有多少种比对结果得到最终比对序列的长度
        for i in range(0,c):
            if finalOrder1[i] != finalOrder2[i]:
                distance += 1
    return distance

b.根据距离矩阵选取出最小值，将两点合并，将这两点及其需要在图上画出的距离存储到List树中（树结构存储在列表中），根据上述步骤依次递归完成树的构建，并输出List。
c.根据列表树结构使用matplotlib完成树的绘制。从根节点开始，在坐标图上将点画出然后连接起来，并将距离使用注释标上，然后递归画字树，完成整棵树的绘制。
b，c过程代码如下：

import globalAlignment as ga          #导入全局序列比对文件
import matplotlib.pyplot as plt
import numpy as np
import copy

fig = plt.figure()
ax = fig.add_subplot(111)
ax.set_title('Phylogenetic tree')        #给图添加标识
plt.ylabel('Distance')                   #给y轴添加标识
ax.set_xticks([])                        #不显示x轴上的坐标，让图更好看

def drawArrow(auxiliaryList, treeMark, n, m, upperX, upperY):
    if isinstance(treeMark[0], float):
        x = (m + n) / 2                 #将最上面的点显示在画布的中央位置让图更好看
        if upperX != 0:                 #根节点是没有父节点的，以0来表示没有父节点，不需要向上连接
            ax.plot((x, upperX), (treeMark[0], upperY), color='k')
        else:
            ax.set_xlim(0, n)            #规定画布坐标范围
            ax.set_ylim(0, treeMark[0] + 1)
        ax.scatter(x, treeMark[0], c='k',marker='o')
        ax.text(x + 0.05, treeMark[0]+0.05 , treeMark[0])                       #标识出距离
        ax.scatter((x + n) / 2, treeMark[0], c='k',marker='o',alpha= 0 )
        ax.scatter((x + m) / 2, treeMark[0], c='k',marker='o',alpha= 0)
        ax.plot(((x + n) / 2, x), (treeMark[0], treeMark[0]), color='k')
        ax.plot(((x + m) / 2, x), (treeMark[0], treeMark[0]), color='k')
        drawArrow(auxiliaryList[2], treeMark[2], n, x, (x + n) / 2, treeMark[0])     #递归画右子树
        drawArrow(auxiliaryList[1], treeMark[1], x, m, (x + m) / 2, treeMark[0])     #递归画左子树
    else:                                                   #到达叶子节点，将叶子节点画出
        x = (m + n) / 2
        ax.plot((x, upperX), (0.05, upperY), color='k')
        ax.scatter(x, 0.02, c='k')
        ax.text(x+0.05, 0.05, treeMark[0])


#创建全局变量存储最终结果
def createList():
    global treeFinalList
    global treeFinalMark
    treeFinalList = []
    treeFinalMark = []

#返回矩阵中最小值及其位置的函数
def getMatrixMin(matrix):
    m,n = np.shape(matrix)
    matrixMin = matrix[0][0]
    y = 0                    #记录的是在字典中对应的值而不是矩阵下标
    x = 1                    #第一行表示的是字典中的第二个
    for i in range(m):
        for j in range(i,n):
            if matrix[j][i] < matrixMin :
                matrixMin = matrix[j][i]
                x = j+1
                y = i
                             #x表示的是最小值在矩阵中的行+1，y表示的是最小值在矩阵中的列
    return x,y,matrixMin

#构造下一个距离矩阵  Nextdistancematrix
def createNdm(dic,odm,auxiliaryList,treeMark):
    if len(odm) == 0:
        print('The phylogenetic tree is stored as follows:')
        print(auxiliaryList)
        print(treeMark)
        return
    x,y,matrixMin = getMatrixMin(odm)   #得到距离最小的两个序列在字典中的值
    #dic1 = dic      注意这种属于浅拷贝 改变dic1的值也会同时改变dic的值          #先复制一个字典，再对其进行后续操作
    dic1 = copy.deepcopy(dic)                                #需要使用深拷贝
    dic1.pop(list(dic.keys())[list(dic.values()).index(x)])     #后续字典中需要删除距离最小的两个序列
    dic1.pop(list(dic.keys())[list(dic.values()).index(y)])
    dic1len = len(dic1)
    ndm = np.zeros((dic1len,dic1len))                          #因为字典中还需要加一个键值对，所以创建矩阵是长度不需要减一
    index = 0
    for s in dic1.keys():
        if x < dic1[s]:                                           #如果删去的值是排在这次要求的前面那么就是自己的序号减一作为行
            d1 = odm[dic1[s]-1][x]
        else :
            d1 = odm[x-1][dic1[s]]
        if y < dic1[s]:                                           #如果删去的值是排在这次要求的前面那么就是自己的序号减一作为行
            d2 = odm[dic1[s]-1][y]
        else :
            d2 = odm[y-1][dic1[s]]
        ndm[dic1len-1][index] = (d1+d2)/2
        dic1[s] = index                                           #对字典中的值进行更新，方便对矩阵赋值，之所以在这里才赋值是因为上面的步骤还需要用到之前的值
        index  += 1
    for i in range(dic1len):
        for j in range(i+1,dic1len):
            ndm[j-1][i] = odm[dic[list(dic1.keys())[list(dic1.values()).index(j)]]-1][dic[list(dic1.keys())[list(dic1.values()).index(i)]]]
    newSquence = list(dic.keys())[list(dic.values()).index(y)] + list(dic.keys())[list(dic.values()).index(x)]                    #获得新节点的名字
    dic1[newSquence] = index
    list1 = [newSquence,auxiliaryList[y], auxiliaryList[x]]
    list2 = [matrixMin/2,treeMark[y], treeMark[x]]
    #不能这么删，因为删掉一个之后下标就变了
    #del auxiliaryList[y]
    #del auxiliaryList[x]
    #要先删掉后面那个，才不会影响前面的删除
    del auxiliaryList[x]
    del auxiliaryList[y]
    del treeMark[x]
    del treeMark[y]
    auxiliaryList.append(list1)
    treeMark.append(list2)
    createNdm(dic1,ndm,auxiliaryList,treeMark)

#计算距离初始矩阵函数
def createOdm(sequenceList):
    l = len(sequenceList)                   #输入序列的个数
    odm = np.zeros((l-1,l-1))               #初始化距离矩阵,距离矩阵只需要l-1的行列
    for i in range(l-1):
        for j in range(i+1,l):              #调用全局序列比对中的获得距离函数，进行序列比对是因为碱基可能会发生缺失或增加
            odm[j-1][i] = ga.getDistance(sequenceList[i],sequenceList[j])
    return odm

#根据距离矩阵和对应关系字典计算下一步的距离矩阵和字典的函数


if __name__ == '__main__':
    createList()
    print('Please enter the sequence number(no less than 2):')
    n = int(input())
    treeMark = []
    auxiliaryList = []                 #辅助列表
    dic = {}                           #用于存储序列名称与矩阵下标的对应关系
    for i in range(n):                 #按照输入的顺序 第一个命名为A，以此类推
        dic[chr(65+i)] = i
        auxiliaryList.append([chr(65+i)])
        treeMark.append([chr(65+i)])
    sequenceList = []                  #用来存储序列的列表
    for i in range(1,n+1):
        print('Please enter sequence '+str(list(dic.keys())[i-1]))
        enter = input()
        sequenceList.append(enter)
    odm = createOdm(sequenceList)      #创建初始距离矩阵
    getMatrixMin(odm)
    createNdm(dic,odm,auxiliaryList,treeMark)
    drawArrow(auxiliaryList[0],treeMark[0],n+1,0,0,0)
    plt.show()

二、实验结果

结果一

结果二与结果三中都有序列不等长的情况，计算距离时都使用了全局序列比对。

结果二

结果三

---

总结

1.实验中遇到的问题及困难

a.因为上课的时候算距离讲的都是序列长度相等的，然后就被如果序列有增添或缺失导致长度不相等的时候怎么算距离困扰，然后就想到了先用序列比对，然后再计算距离。
b.树怎么存？开始想的是使用字典来存储，但是始终不能满足后续的画图要求，然后就改用了列表。叶子节点存的是序列名称，其他节点存储的是距离，方便画图，同时为了更清楚的显示给用户，也用了一个列表存储所有节点的名称并显示出来了。
c.因为对于python还不熟悉，在进行树的构建和画图的时候都遇到了字典赋值和深拷贝的问题。
d.画图的时候也是绝了，刚开始是想像决策树一样画的，结果还是太丑也不好画，然后还是使用散点图将所有节点都画出来，还有转角那里本来也是有画点的，但是用了alpha参数使其为零，该点就变透明了，然后也把x坐标轴上的坐标去掉了，会好看一点。每次画图都会学到一些新东西。

2.还能改进的地方

a.距离有多个最小值怎么办，后面的想法是如果是AB，CD都是最小值的话可以不用管，因为无论先合并哪两个，剩下的两个就会在下一次被合并，但如果出现的是AB与AC都为最小值的话就会是不同的情况，画的图会不一样，暂时想法是有几个这样的就画几个不一样的图，但是具体还没有实现，只是一个想法
b.本来想把画布逆时针旋转90度让图更好看的，但是没找到解决办法。