文章目录

算法描述
设计思想
源代码及运行成果
- 1.运行成果
- 2.源代码
遇到的问题及总结
- 问题
- 总结

算法描述

分子进化与系统发生：计算分子进化——利用算法在分子水平上构建物种的进化树。这里说的分子水平是指DNA、RNA、以及蛋白质序列。

1. 非加权分组平均法：
对于某几条不同的序列，找出距离最小的一对序列（两条序列之间的距离为两条序列相同位置上不同的碱基总数），将其合并聚集，形成一个分支，聚集后将两条序列看成一个整体，再分别计算与其他序列之间的距离，直到所有序列归为一类，系统发生树构建完成，树枝的长短直接反应了它们与共同祖先的距离。

2. 邻接法：

![在这里插入图片描述](https://img-blog.csdnimg.cn/20201027021904765.png#pic_center)

设计思想

- 非加权分组平均法：

读取输入的序列并将进行命名
构建距离矩阵（使用pandas库的dataframe进行存储）
使用循环逐步更新距离矩阵，直到系统发生树构建完成（即距离矩阵行只剩一行），更新距离矩阵分为以下四种情况：

	A	B	D
B	①
C		③
D
E	④		②

①最小距离在AB处：列A与列B进行相加计算得聚类AB后的距离，赋值给列A，删除行B及列B得到最新距离矩阵如下

	A B	C	D
C
D
E

②最小距离在DE处：行C与行D进行相加计算得聚类CD后的距离，赋值给行E，删除行D及列D得到最新距离矩阵如下

	A	B	C
B
C
DE

③最小距离在BC处：列B与列C进行相加计算得聚类BC后的距离，赋值给列B，删除行C；行B与行C进行相加计算得聚类BC后的距离，赋值给行B，删除行C，得到最新距离矩阵如下

	A	BC	D
BC
D
E

④最小距离在AE处：A列B、C、D行分别与E行B、C、D列相加后减半，删除行E，删除最后一列D，得到最新距离矩阵如下

	AE	B	C
B
C
D

画出系统发生树，根据各分支的关系及聚类之间最小距离设置节点坐标，并使用matplotlib的plot函数从根节点往分支方向依次画出系统发生树，各节点的纵坐标为距离值

- 邻接法：

读取输入的序列并将进行命名
构建距离矩阵（使用pandas库的dataframe进行存储）
计算初始初始净分歧度和最小速率校正距离
使用循环逐步更新距离矩阵，情况与非加权分组平均法类似，但计算方法略有不同，还需计算净分歧度r及最小速率校正距离

源代码及运行成果

1.运行成果

非加权分组平均法

运行截图如下：

邻接法

2.源代码

非加权分组平均法 代码如下：


```python
#!/usr/bin/env python
# -*- coding:utf-8 -*-
from pandas import DataFrame
from pylab import *
from matplotlib import pyplot as plt
mpl.rcParams['font.sans-serif'] = ['SimHei']#解决中文乱码#计算两序列之间的距离
def CalDistance(str1,str2):dis = 0for i in range(len(str1)):if str1[i] != str2[i]:dis += 1return dis#获取dataframe中最小值的行名列名
def getMinName(df):#最小值所在列minColNameumns = df.min().idxmin()#最小值Min = df.min().min()#最小值所在行minIndexName = ''.join(df[(df[minColNameumns] == Min)].index.tolist()[0])#转换为list再转换为stringreturn minIndexName, minColNameumns,Min#输入序列命名
LetterIndex = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
#输入序列
Sequences = []
print("请输入序列（每条序列以换行符相隔）：")
for line in iter(input, ''):Sequences.append(line.replace(',',''))
length = len(Sequences)
#创建列名序列
columns = list(LetterIndex)[:length-1]
#创建行名序列
index = list(LetterIndex)[1:length]
#创建距离矩阵
distance = DataFrame(np.zeros((length-1,length-1)),index = index,columns = columns)
#存储系统发生树的节点坐标
PointList = {}
#存储父节点与子节点的关系
LeafDict = {}
#画树所需节点集合
xPointList = []
yPointList = []
MinDistance = {}
order = []print("输入序列如下   ")
for i in range(length):print("序列",LetterIndex[i],":",Sequences[i])
#计算距离矩阵
for i in range(length - 1):for j in range(length - 1):if columns[i] != index[j] :# 序列名字转换并计算各序列的距离distance.loc[index[j],columns[i]] = CalDistance(Sequences[LetterIndex.index(columns[i])],Sequences[LetterIndex.index(index[j])])else:distance.loc[index[j], columns[i]] = np.NaN
#构建距离矩阵
while len(distance) >= 1:#判断是否构建为最终的距离矩阵len(distance) >= 1print("============================")print(distance)#获取表格中最小值的行名、列名minIndexName,minColName,Min  = getMinName(distance)#连接行名及列名newSeName = minColName + minIndexName#存储系统发生树中父子节点关系LeafDict[newSeName] = [minColName,minIndexName]#存储最小距离MinDistance[newSeName] = Min/2#存储父节点构建顺序order.append(newSeName)#更新序列合并后的距离矩阵if len(distance) > 1:if minIndexName in columns and minColName not in index :distance.drop(minIndexName,axis=0,inplace=True)# 删除多余行distance[minColName] = (distance[minColName] + distance[minIndexName])/2# 求得两序列连接后的距离矩阵distance.rename(columns={minColName: newSeName}, inplace=True)# 改列名distance.drop(minIndexName, axis=1, inplace=True)# 删除多余列columns = list(distance)# 更新列名elif minColName in index and minIndexName not in columns:distance.drop(minColName, axis=1, inplace=True)  # 删除多余列distance.loc[minIndexName] = (distance.loc[minIndexName] + distance.loc[minColName]) / 2 # 算新距离矩阵distance.rename(index={minIndexName: newSeName}, inplace=True)  # 改行名distance.drop(minColName, axis=0, inplace=True)  # 删除多余行index = distance._stat_axis.values.tolist() # 更新行名elif minIndexName in columns and minColName in index:# 若行名列名都出现在行索引及列索引集合中distance[minColName] = (distance[minColName] + distance[minIndexName]) / 2distance.rename(columns={minColName: newSeName}, inplace=True)  # 改列名distance.drop(minIndexName, axis=1, inplace=True)  # 删除minIndexName列distance.loc[minColName] = (distance.loc[minIndexName] + distance.loc[minColName]) / 2  # 算新距离矩阵distance.rename(index={minColName: newSeName}, inplace=True)  # 改行名distance.drop(minIndexName, axis=0, inplace=True)  # 删除minColName行distance.at[newSeName,newSeName] = np.NaN# 更新行名及列名columns = list(distance)index = distance._stat_axis.values.tolist()else:#获取行名列表newindex = distance._stat_axis.values.tolist()newindex.remove(minIndexName)#计算新距离矩阵for id in newindex:distance.at[id,minColName] = (distance.at[id,minColName] + distance.at[minIndexName,id])/2distance.rename(columns={minColName: newSeName}, inplace=True)  # 改列名distance.drop(minIndexName, axis=0, inplace=True)  # 删除minIndexName行distance.drop(list(distance)[-1], axis=1, inplace=True)  # 删除最后一列#更新行名及列名columns = list(distance)index = distance._stat_axis.values.tolist()else:break# 设置树的节点坐标
for i in range(len(order)-1,-1,-1):root = order[i]if i == len(order) - 1:interval = 0.4maxDistance = MinDistance[root]PointList[root] = (interval,0)interval /= 2leaf1 = LeafDict[root][0]leaf2 = LeafDict[root][1]if leaf1 not in order:PointList[leaf1] = (PointList[root][0] + interval, maxDistance)else:PointList[leaf1] = (PointList[root][0] + interval,maxDistance - MinDistance[leaf1])if leaf2 not in order:PointList[leaf2] = (PointList[root][0] - interval, maxDistance)else:PointList[leaf2] = (PointList[root][0] - interval, maxDistance - MinDistance[leaf2])# 为画树构建节点列表xPointList.append([PointList[leaf1][0], PointList[root][0]])xPointList.append([PointList[leaf2][0], PointList[root][0]])yPointList.append([PointList[leaf1][1], PointList[root][1]])yPointList.append([PointList[leaf2][1], PointList[root][1]])# 画树
fig=plt.figure()
ax=fig.add_subplot(111)
for i in range(len(xPointList)):plt.plot(xPointList[i],yPointList[i],color = 'black',alpha = 0.5)plt.scatter(xPointList[i],yPointList[i],color = 'black',alpha=0.8)
# 标明注释
for k,v in PointList.items():if k in list(LetterIndex)[:length]:strName = k +":" + Sequences[LetterIndex.index(k)]plt.annotate(s = strName,xy=v,xytext=(v[0]-10,v[1]+10) ,rotation = 45,textcoords = 'offset points',alpha = 0.5)
plt.ylabel("距离")
ax = plt.gca()
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['bottom'].set_visible(False)
ax.spines['left'].set_visible(False)
plt.show()

邻接法 代码如下：

#!/usr/bin/env python
# -*- coding:utf-8 -*-
from pandas import DataFrame
from pylab import *
from matplotlib import pyplot as plt
mpl.rcParams['font.sans-serif'] = ['SimHei']#解决中文乱码#计算两序列之间的距离
def CalDistance(str1,str2):dis = 0for i in range(len(str1)):if str1[i] != str2[i]:dis += 1return dis#获取dataframe中最小值的行名列名
def getMinName(df):#最小值所在列minColNameumns = df.min().idxmin()#最小值Min = df.min().min()#最小值所在行minIndexName = ''.join(df[(df[minColNameumns] == Min)].index.tolist()[0])#转换为list再转换为stringreturn minIndexName, minColNameumns#计算r
def createR(columns,index,distance):r = {}for s in columns+index:if s == columns[0]:#第一个序列求和r[s] = distance[s].sum()elif s == index[-1]:#最后一个序列求和r[s] = distance.loc[s].sum()else:#中间序列求和r[s] = distance[s].sum() + distance.at[s,columns[0]]return r
#计算M
def createM(columns,index,M,distance,r,length):for col in columns:for ind in index:if col != ind:M.at[ind,col] = distance.at[ind,col] - (r[col] + r[ind])/(length-2)#输入序列命名
LetterIndex = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ'
#输入序列
Sequences = []
# Sequences = ['tagg','tacg','aagc','agcc']print("请输入序列（每条序列以换行符相隔）：")
for line in iter(input, ''):Sequences.append(line.replace(',',''))
length = len(Sequences)
#创建列名序列
columns = list(LetterIndex)[:length-1]
#创建行名序列
index = list(LetterIndex)[1:length]
#创建距离矩阵
distance = DataFrame(np.zeros((length-1,length-1)),index = index,columns = columns)
#存储父节点与子节点的关系
LeafDict = {}
order = []print("输入序列如下   ")
for i in range(length):print("序列",LetterIndex[i],":",Sequences[i])#计算距离矩阵
for i in range(length - 1):for j in range(length - 1):if columns[i] != index[j] :# 序列名字转换并计算各序列的距离distance.loc[index[j],columns[i]] = CalDistance(Sequences[LetterIndex.index(columns[i])],Sequences[LetterIndex.index(index[j])])else:distance.loc[index[j], columns[i]] = np.NaN#构建距离矩阵
while len(distance) >= 1:#判断是否构建为最终的距离矩阵len(distance) >= 1print("============================")print("distance:")print(distance)# 计算r净化距离r = createR(columns, index, distance)# 计算最小速率校正矩阵M = distance.copy(deep=True)createM(columns, index, M, distance, r, length)# print("M:" )# print(M)#获取表格中最小值的行名、列名minIndexName,minColName  = getMinName(M)#连接行名及列名newSeName = minColName + minIndexNamenum = distance.at[minIndexName,minColName]LeafDict[newSeName] = [minColName,minIndexName]order.append(newSeName)#更新序列合并后的距离矩阵if len(distance) > 1:if minIndexName in columns and minColName not in index :distance.drop(minIndexName,axis=0,inplace=True)# 删除多余行distance[minColName] = (distance[minColName] + distance[minIndexName] - num)/2# 求得两序列连接后的距离矩阵distance.rename(columns={minColName: newSeName}, inplace=True)# 改列名distance.drop(minIndexName, axis=1, inplace=True)# 删除多余列columns = list(distance)# 更新列名index = distance._stat_axis.values.tolist()  # 更新行名elif minColName in index and minIndexName not in columns:distance.drop(minColName, axis=1, inplace=True)  # 删除多余列distance.loc[minIndexName] = (distance.loc[minIndexName] + distance.loc[minColName] - num) / 2 # 算新距离矩阵distance.rename(index={minIndexName: newSeName}, inplace=True)  # 改行名distance.drop(minColName, axis=0, inplace=True)  # 删除多余行columns = list(distance)  # 更新列名index = distance._stat_axis.values.tolist() # 更新行名elif minIndexName in columns and minColName in index:# 若行名列名都出现在行索引及列索引集合中distance[minColName] = (distance[minColName] + distance[minIndexName] -num) / 2distance.rename(columns={minColName: newSeName}, inplace=True)  # 改列名distance.drop(minIndexName, axis=1, inplace=True)  # 删除minIndexName列distance.loc[minColName] = (distance.loc[minIndexName] + distance.loc[minColName] -num) / 2  # 算新距离矩阵distance.rename(index={minColName: newSeName}, inplace=True)  # 改行名distance.drop(minIndexName, axis=0, inplace=True)  # 删除minColName行distance.at[newSeName,newSeName] = np.NaN# 更新行名及列名columns = list(distance)index = distance._stat_axis.values.tolist()else:#获取行名列表newindex = distance._stat_axis.values.tolist()newindex.remove(minIndexName)#计算新距离矩阵for id in newindex:distance.at[id,minColName] = (distance.at[id,minColName] + distance.at[minIndexName,id] - num)/2distance.rename(columns={minColName: newSeName}, inplace=True)  # 改列名distance.drop(minIndexName, axis=0, inplace=True)  # 删除minIndexName行distance.drop(list(distance)[-1], axis=1, inplace=True)  # 删除最后一列#更新行名及列名columns = list(distance)index = distance._stat_axis.values.tolist()else:break
print("============================")
print("聚类如下:")
print(LeafDict)

遇到的问题及总结

问题

输入序列的读取与命名 解决方法-使用含有26个字母的字符串进行依次转换命名（待完善）
思维误区1：更新距离矩阵考虑情况过于简单，只考虑到情况①②，导致调试报错 解决方法-逐步调试完善情况
思维误区2：画树时从叶子节点往根节点方向画导致叶子节点无法改变，造成画树时树枝交叉缠绕 解决方法-另外创建数据结构存储画树所需节点信息，从根节点往叶子方向画（待完善）

总结

系统发生树的距离矩阵我使用了pandas库的dataframe进行存储，相对于numpy的二维数组，dataframe表格在对行列进行操作时更加方便，不需要使用循环而可以直接相加减，实现距离矩阵的更新时简化了代码，但对于是否提升时间复杂度还有待学习，需要了解pandas的dataframe的内部实现，暂时没有查到资料（了解到dataframe比较适用于操作表格）。

最初实现算法时由于头脑简单忽略了很多矩阵更新的情况导致后来代码大改写，严重影响到整个作业的进度和心态。画树的过程中也遇到了很大的问题，因为最初头脑简单考虑不周导致只能画出某个特定情况下的树，由于对matplotlib的不熟悉可能也没有选择最佳的画图函数，第一次画树失败后卡了很久甚至一度想放弃画图，最后牺牲了空间开销，转变思维方向画出了树的结构，但树的形状属实还有待美化。

在实现邻接法构建树的过程中由于时间原因没有完成树的构建，只是输出了距离矩阵和分类。

归根到底还是对算法理解不到位以及对python的pandas和画图工具库不熟悉，导致实现非加权分组平均法时耗费了过多时间。总的来说虽然使用非加权分组平均法构建出了系统发生树，但代码中有很多不成熟的处理方式，且邻接法画树的步骤还有待完成，还十分需要改进和完善。

2.系统发生树构建（非加权分组平均法、邻接法python实现）相关推荐

进化树构建之邻接法(Neighbor-Joining)的介绍
进化树构建进化树构建的问题是推断可能产生给定基因序列数据的进化树的拓扑结构和分支长度.推断树中叶节点的数量应等于给定数据中基因序列的数量. Neighbor-Joining Algorithm Ne ...
数学建模——层次分析法Python代码
数学建模--层次分析法Python代码 import numpy as np class AHP: """ 相关信息的传入和准备 """ d ...
Python训练营2021:构建8个真实世界的Python项目
时长:19h 27m |视频:. MP4,1280×720 30 fps |音频:AAC,44.1 kHz,2ch |大小:9.54 GB 语言:英语+中英文字幕机译从Python Web开发的初 ...
keras构建前馈神经网络（feedforward neural network）进行分类模型构建基于早停法（Early stopping)
keras构建前馈神经网络(feedforward neural network)进行分类模型构建基于早停法(Early stopping) 当我们训练深度学习神经网络的时候通常希望能获得最好的泛化性 ...
灰色关联与TOPSIS法 —— python
目录 1.简介 2.算法详解 2.1 指标正向化及标准化 2.2 找到最大最小参考向量 2.3 计算与参考向量的相关系数 2.4 求评分 3.实例分析 3.1 读取数据 3.2 数据标准化 3.3 得 ...
非线性规划的拉格朗日乘子法python编程python包编程
非线性规划的拉格朗日乘子法&python编程&python包编程一.拉格朗日乘子法 1.1 拉格朗日乘子法定义 1.2 KKT条件定义 1.3 拉格朗日乘子法手工推导例题二.Pyt ...
matlab套利,期现套利-现货组合构建(1)-市值权重法
本帖最后由 faruto 于 2011-12-27 23:58 编辑期现套利-现货组合构建(1)-市值权重法最近抽空做的一点东西,常见的期现套利现货构建的方法之1--市值权重法. 一直想把跟踪指数 ...
大津法Python实现
大津法Python实现 1.简介在计算机视觉和图像处理中,大津法被用于自动获取图像的阈值,或者将灰度图像转换为二值化图像.该算法假设图像包含两个类别的像素(前景像素和背景像素),然后它计算一个最优的 ...
二分插入排序法-Python版
简介: 传统的插入法思想为: 1.从头开始,构建有序数列: 2.再将之后需要排序的数据从头至尾(或从尾至头)进行比较,插入到其相应的位置. 而二分插入排序法第一步与传统插入法相同,第二步而是采用二分查 ...

2.系统发生树构建（非加权分组平均法、邻接法python实现）