1.读取文件完成随机概率矩阵构建，并找到已知相互作用数最多的蛋白质P 2.使用重启随机游走算法找到与P相关的top20蛋白质 3.使用五折交叉验证方法评判算法性能 4.使用matplotlib画出roc曲线并计算出auc值
本人参考了下面的链接博文，从随机游走算法开始介绍了重启随机游走 重启随机游走算法
a.逐行读取文本文件，求出概率转移矩阵W，求出其中度最大的位置对应的蛋白质P，根据P的位置初始化r和e，带入公式进行五十次迭代计算出r的最终值，根据r得到与P相关性最大的top20蛋白质并输出。（重启概率自定义） b.取出与P相关的蛋白质位置重新生成一个List，打乱其顺序，将其分成五份，每次取其中的一份做测试集，重新计算r，将r中作为训练集的蛋白质的相关性值取出，对剩下的进行排序，完成TPR与FPR的计算，该过程需重复五次。 c.根据计算出的TPR与FPR画出Roc曲线。 代码如下
import numpy as np import re import copy as cp import numpy as np import random import matplotlib.pyplot as plt from sklearn.metrics import roc_curve, auc a = 0.85 #重启概率 #根据文件得到转移概率矩阵 def createW(f): dic = {} # 用来存储蛋白质名称与其在w中对应位置的关系 interrelation = [] #存储蛋白质之间的关系 一个元组表示一个关系 line = f.readline() i = 0 while line: line = re.split(' |\t', line) if line[0] not in dic.keys(): dic[line[0]] = i i += 1 if line[3] not in dic.keys(): dic[line[3]] = i i += 1 relation = (line[0], line[3]) interrelation.append(relation) line = f.readline() W = np.zeros((len(dic),len(dic))) #概率转移矩阵 for k in interrelation: W[dic[k[0]]][dic[k[1]]] = 1 W[dic[k[1]]][dic[k[0]]] = 1 s = 0 Wt = np.transpose(W) #转置之后更好算 max1 = sum(Wt[0]) max_index = 0 # 存储度最大的节点 for n in range(len(dic)): num = sum(Wt[n]) Wt[n] = Wt[n] / num if max1 <= num: #完成矩阵标准化 max_index = n max1 = num W = np.transpose(Wt) #返回的分别是 概率转移矩阵，蛋白质对应矩阵位置字典，关系元组列表，度最大的蛋白质位置，度的最大值 return W,dic,interrelation,max_index,max1 #得到Top20蛋白质的位置 def getTop_20Index(r): #m表示所要求最大值的列表，n表示要求排名前多少的 r = list(r) top_20Index = [] for i in range(20): a = max(r) # 最大值 b = r.index(a) # 最大值的位置 top_20Index.append(b) #添加最大值位置 r[b] = 0 #将此次的最大值删除 return top_20Index #得到真实标记列表函数 def getTureValueList(W,relatedIndex,max1,max_index,dic,dataSet): surplusNum = max1 - len(dataSet) newCon = 1.0 / surplusNum wP = cp.deepcopy(W) relatedIndex1 = cp.deepcopy(relatedIndex) for j in relatedIndex: if j not in dataSet: wP[max_index][j] = newCon wP[j][max_index] = newCon else: relatedIndex1.remove(j) wP[max_index][j] = 0 wP[j][max_index] = 0 e = np.zeros((len(dic), 1)) e[int(max_index)][0] = 1 r = np.zeros((len(dic), 1)) r[int(max_index)][0] = 1 for i in range(50): # 50次迭代将r算出 if i != 0: r = r1 r1 = a * np.dot(wP, r) + (1 - a) * e r2 = [] # 存储与p没有连接的相关度得分 r2_index = [] for k in range(len(dic)): if k not in relatedIndex1: r2.append(r1[k][0]) r2_index.append(k) tureValue = [] # 存储是否与P有关联 r3 = sorted(r2) for k in range(len(r3)): if r2_index[r2.index(r3[k])] in dataSet: tureValue.append(1) else: tureValue.append(0) return tureValue #得到TPR和FPR的函数 def getTPRFPR(tureValueList): tureValueList1 = list(reversed(tureValueList)) #反转列表 tureNum = sum(tureValueList) faulseNum = len(tureValueList) - tureNum TPR = [] FPR = [] tureSum = 0 faulseSum = 0 for i in tureValueList1 : if i == 1: tureSum += 1 else: faulseSum += 1 TPR.append(tureSum / tureNum) FPR.append(faulseSum/faulseNum) return TPR,FPR #画ROC曲线函数 def drawROC(TPRList,FPRList): plt.figure(figsize=(5, 5)) #循环画出五条Roc曲线 tprtotal = np.zeros((1,len(TPRList[0])))[0] fprtotal = np.zeros((1,len(TPRList[0])))[0] for i in range(5): fpr = FPRList[i] tpr = TPRList[i] tprtotal += np.array(tpr) fprtotal += np.array(fpr) roc_auc = auc(fpr, tpr) plt.plot(fpr, tpr, lw=2,label='ROC curve (area = %0.3f)' % roc_auc) #假正率为横坐标，真正率为纵坐标做曲线 plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--') #画出平均的roc曲线 mean_tpr = list(tprtotal/5) # 求数组的平均值 mean_fpr = list(fprtotal/5) mean_auc = auc(mean_fpr, mean_tpr) plt.plot(mean_fpr, mean_tpr, 'k--', label='Mean ROC (area = {0:.2f})'.format(mean_auc), lw=2) #设置x，y plt.xlim([-0.05, 1.0]) plt.ylim([0.0, 1.05]) plt.xlabel('False Positive Rate') plt.ylabel('True Positive Rate') plt.title('ROC') plt.legend(loc="lower right") plt.show() #参数 概率转移矩阵 度最大的节点位置 最大的度的值 节点与位置对应的字典 def crossValidation(W,max_index,max1,dic): relatedIndex = [] deletNum = int(int(max1)/5) #每个测试集中的节点数 这里应该是54 for i in range(len(dic)): #得到所有与p相关点的下标 存入relatedIndex中 if W[max_index][i] != 0: relatedIndex.append(i) random.shuffle(relatedIndex) #将relatedIndex中的元素顺序打乱 完成数据集的随机划分 #完成数据集的划分 tureValueList = [] TPRList = [] FPRList = [] for i in range(5): #分别计算五个数据集的TPR和FPR dataSet = relatedIndex[i * deletNum : (i + 1) * deletNum] tureValueList.append(getTureValueList(W,relatedIndex,max1,max_index,dic,dataSet)) TPR,FPR = getTPRFPR(tureValueList[-1]) TPRList.append(TPR) FPRList.append(FPR) drawROC(TPRList, FPRList) if __name__ == '__main__': f = open('F:\BINARY_PROTEIN_PROTEIN_INTERACTIONS.txt') W,dic,interrection,max_index,max1 = createW(f) print('The protein with the largest number of interactions is known to be：') print(list(dic.keys())[max_index]) f.close() e = np.zeros((len(dic), 1)) e[int(max_index)][0] = 1 r = np.zeros((len(dic), 1)) r[int(max_index)][0] = 1 for i in range(50): # 50次迭代将r算出 if i != 0: r = r1 r1 = a * np.dot(W, r) + (1 - a) * e top_20Index = getTop_20Index(r) #与p相关性最大蛋白质top20的位置 print('The top 20 proteins associated with it:') print('=========================') for i in range(20): name = list(dic.keys())[top_20Index[i]] print('\t'+str(i+1)+'\t'+name) print('=========================') crossValidation(W, max_index, max1, dic) #进行五折交叉验证top20的蛋白质是不会变的，所以两个实验结果左边都是一样的，但是roc会不一样，因为每一次的数据集划分都是随机的。 结果一 结果二
排名第17的蛋白质是-，我开始以为是读取文件的时候出现了问题，但是文件中确实有这个命名
本次实验主要遇到的问题还是在于对算法的理解，通过b站视频的学习和后面老师的多次指导才弄清白整个过程是如何实现的。实现过程还算顺利，通过博客查找都能找到解决办法。
a.我觉得程序最大的问题在于运行时间比较久，因为在程序中多次使用到了for循环而且矩阵有比较大，但是自己也还没有找到方法优化，后期应该更注重时间性能，对程序进行修改。 b.在划分数据集的时候不能平均分就舍弃了1个蛋白质节点，这个应该是可以再想想解决办法的。
