第5章AlphaFold与蛋白质结构预测
长期以来，研究人员使用诸如X射线晶体衍射、核磁共振和冷冻电镜等实验方法确定蛋白质的结构，每种方法都需要大量的实验、错误试探和纠错，每一个蛋白质结构的成功解析，背后都是无数艰辛的科学探索与实验支撑。近年来，随着AI技术的崛起，借助AI方法，根据氨基酸残基序列直


图5.1基于残基序列直接预测蛋白质结构

接计算预测蛋白质结构成为计算生物学的热点领域。如图5.1所示，直接根据残基序列推断蛋白质三维结构的方法，相当于用超强的计算过程取代了科学家们繁重的实验过程，是一条极富科学前景的探索之路。
5.1什么是AlphaFold
蛋白质结构预测技术的关键评估（Critical Assessment of protein Structure Prediction，CASP），是基于残基序列直接预测蛋白质三维结构的全球科学盛会，目标是对来自世界各地的科学家团队提交的蛋白质结构预测模型的关键指标做出科学评估。
参加CASP比赛的团队一般由科学家、工程师和人工智能专家等组成，CASP每两年举办一次。从1994年的第1届CASP1到2018年的第13届CASP13，从未间断。在撰写本教材的过程中，第14届CASP14已经于2020年4月拉开序幕。有人将两年一度的CASP誉为结构生物学领域的“奥林匹克”。
Google旗下DeepMind公司的一个人工智能研究团队参加了CASP13，团队名称为A7D，其开发的蛋白质结构预测系统称为AlphaFold（又称A7D系统）。相关研究成果发表在Nature上(Senior，Evans，et al.，2020)，深度学习模型的代码开源在GitHub（https://github.com/deepmind/deepmindresearch/tree/master/alphafold_casp13）。AlphaFold系统的逻辑结构如图5.2所示。






图5.2AlphaFold系统的逻辑结构



图5.2源自Andrew W. Senior等人发表于Nature的论文。图5.2（a）表示AlphaFold的系统逻辑，包括以下几个关键步骤。
（1） 特征提取阶段。使用HHblits和BLAST软件，通过对已有的蛋白质数据库的搜索、构建和计算，完成给定残基序列的多序列比对对齐（Multiple Sequence Alignment，MSA）分析，实现MSA特征提取。
（2） 模型训练与结构预测。定义了一个包含220个残差块的ResNet深度卷积神经网络，基于序列和MSA特征，训练和预测蛋白质残基对之间的距离，以及每个残基的二面角Phi（）和Psi（ψ）。
单个残差块的结构如图5.2（b）所示。图5.2（b）是一个ResNet的残差结构，输入层的通道数为128，残差块包含3个卷积层，第一个卷积层采用了1×1的卷积降维到64通道，第2个卷积层采用3×3的卷积，第3个卷积层采用1×1的卷积上采样到128通道，卷积层采用Batch norm方法和Elu激励函数。
（3） 二次特征提取与整合。整合深度学习得到的二面角、残基对距离，结合原子的范德华半径，构建生成蛋白质结构的特征集。
（4） 生成蛋白质结构。根据模型预测的距离分布、二面角Phi（）和Psi（ψ）的分布以及范德华半径构建蛋白质的候选潜在结构，通过一个深度卷积网络估算候选结构的准确性，最终确定蛋白质结构。
AlphaFold在全球约一百个研究团队提交的数万个模型中脱颖而出，取得综合排名第一的成绩，取得了用计算方法预测蛋白质结构前所未有的进步。图5.3所示为AlphaFold








图5.3AlphaFold的蛋白质结构建模过程解析


对CASP13中标签为T0986s2的蛋白质结构建模过程的解析。
图5.3源自Andrew W. Senior等人发表于Nature的论文。T0986s2的序列长度为L=155，PDB数据库中的蛋白质ID为6N9V。
图5.3中a表示结构预测步骤，依次是特征提取，深度神经网络进行距离和二面角预测，二次特征提取，深度神经网络对模型结构进行评估与筛选。
图5.3中b表示ResNet神经网络根据MSA特征和序列特征预测得到L×L的距离分布图，预测过程以64×64作为残基窗口进行累积预测。
图5.3中c显示的是经过1200步梯度下降迭代，TM分数和均方根误差（r.m.s.d.）的变化图形，TM分数逐步提升，均方根误差逐步下降，并且给出了第0步、第280步、第500步、第600步、第1200步预测的蛋白质二级结构形状。
TM是CASP采用的蛋白质结构评分方法，表示预测的结构总体上与真实结构的匹配程度，取值范围为［0，1］。AlphaFold在蛋白质结构自由建模领域遥遥领先于其他方法。
图5.3中d是对模型的预测结构与真实结构的匹配程度的直观展示，灰色表示预测的结构，彩色表示真实的结构。
图5.3中e表示整个测试集（377个蛋白质样本）的TM平均分数随着迭代次数的变化趋势。
5.2肽键、多肽与肽链
蛋白质几乎参与了所有的生命活动过程，是生命活动的物质基础，是构成生物体最基本的功能物质。蛋白质的主要组成元素占比如表5.1所示。


表5.1蛋白质组成元素占比



元 素 名 称占比元 素 名 称占比

碳50％~55%氮16%
氢6.5%~7.3％硫0%~3％
氧19%~24％其他微量剩余占比



图5.4α氨基酸的基本结构

蛋白质的基本构成单位为氨基酸，常见氨基酸有20种。用酸、碱或蛋白酶可以分别将蛋白质水解，生成游离的氨基酸。组成蛋白质的20种常见氨基酸中除脯氨酸外，均为α氨基酸，其基本结构如图5.4所示，除了侧链R以外，其他部分是固定不变的。
肽是由一个氨基酸的羧基（αCOOH）与另一个氨基酸的氨基(αNH2）脱水缩合而成的化合物，氨基酸间脱水后生成的共价键称为肽键（Piptide bond），其中的氨基酸单位称为氨基酸残基。由两个氨基酸缩合而成的肽称为二肽，少于10个氨基酸残基的肽称为寡肽，多于10个氨基酸残基的肽称为多肽。肽链上的各个侧链由不同氨基酸的Ｒ侧链构成。如图5.5所示，两个氨基酸脱水生成二肽。


图5.5两个氨基酸组成二肽


由多个氨基酸通过肽键相互连接形成的链状结构称为肽链。通常在肽链的一端含有一个游离的α氨基，称为氨基端或N端； 在另一端含有一个游离的α羧基，称为羧基端或C端。如图5.6所示，左端是N端，右端是C端。


图5.6肽链的基本结构


肽链中的氨基酸残基按一定的顺序排列，这种排列顺序称为氨基酸顺序。
氨基酸顺序是以N端残基为起点、以C端残基为终点所形成的排列顺序。如图5.6所示的五肽按照残基顺序可表示为SerValTyrAspGln或者SVTAG。
5.3蛋白质的四级结构


图5.7蛋白质的结构层次

每一种天然蛋白质都有自己特有的空间结构，这种空间结构称为蛋白质的（天然）构象。蛋白质的结构层次分为四级，如图5.7所示。
如图5.8所示为蛋白质四级结构示意图。


图5.8蛋白质四级结构示意图


蛋白质的一级结构是指氨基酸残基的排列顺序，这是蛋白质生物功能的基础。一级结构(残基序列)中含有形成高级结构全部必需的信息，一级结构决定高级结构及其功能。所以，根据一级机构（残基序列）可以推断蛋白质的高级空间结构，根据高级空间结构可以进一步推断蛋白质的功能性质。
蛋白质的二级结构是多肽链中各个肽段借助氢键形成有规则的构象，主要包括α螺旋、β折叠、β转角、无规卷曲等。
超二级结构指蛋白质中相邻的二级结构单位（α螺旋或β折叠或β转角）组合在一起，形成有规则的、在空间上能够辨认的二级结构组合体。超二级结构的基本类型有： αα、βαβ、βββ。
结构域是多肽链在二级结构或超二级结构基础上形成三级结构的局部折叠区。结构域是蛋白质的独立折叠单位，一般由100~200个氨基酸残基构成。
蛋白质的三级结构是多肽链借助各种非共价键的作用力，通过弯曲、折叠，形成具有一定走向的紧密球状构象。球状构象的表面积最低，使蛋白质与周围环境的相互作用力最小。
蛋白质的四级结构是寡聚蛋白中各亚基之间在空间上的相互关系和结合方式。
决定蛋白质二级构象的主要因素是二面角和残基之间的距离（接触）。目前，以AlphaFold为代表的机器学习软件，主要是基于二面角和距离（接触）对蛋白质的空间结构做出预测。
5.4数据集
AlphaFold开源在GitHub上的模型，数据集采用已经处理好的TFRecords格式，模型的计算需求较大，初学者的入门门槛较高，因此，本章案例以作者Eric Alcaide发布的MinFold（https://github.com/EricAlcaide/MiniFold）开源项目为起点进行探索，并做了局部修改和优化。
MiniFold可以看作是AlphaFold的一个简易版本，特别是在距离预测与角度预测方面，与AlphaFold的理念与方法保持一致，采用的网络结构也是源自AlphaFold，在此向Eric Alcaide的杰出工作表示致敬！
数据集下载页面（https://github.com/aqlaboratory/proteinnet）包含CASP7～CASP12的比赛数据集。每个数据集包含Textbased与TFRecords两种格式，本章案例选择的是CASP7中的Textbased数据集。
CASP7数据集的下载文件为casp7.tar.gz压缩包，压缩包大小为3.18GB。下载完成后，解压到chapter5＼dataset目录下，相关文件描述如表5.2所示。


表5.2CASP7蛋白质结构数据集



文件名数 据 规 模大小功能

training_3010333个蛋白质样本675MB包含训练集30%的数据
training_5013024个蛋白质样本898MB包含训练集50%的数据
training_7015207个蛋白质样本1.03GB包含训练集70%的数据
training_9017611个蛋白质样本1.2GB包含训练集90%的数据
training_9517938个蛋白质样本1.23GB包含训练集95%的数据
training_10034557个蛋白质样本2.42GB包含训练集100%的数据
validation224个蛋白质样本15.2MB验证集数据
testing93个蛋白质样本6.5MB测试集数据

以training_30文件为例，文件结构组织如表5.3所示。


表5.3training_30文件的组织结构



字 段 域 名行数功能

［ID］单独占1行蛋白质的ID编号
［PRIMARY］单独占1行氨基酸残基序列（单字母表示）
［EVOLUTIONARY］单独占21行 残基序列的PSSM评分矩阵
［TERTIARY］单独占3行主链上N─Calpha─C原子的x，y，z坐标
［MASK］单独占1行残基掩码（‘+’： 有实验数据，‘-’： 数据缺失）

5.5筛选蛋白质序列
从training_30数据集中筛选出残基序列长度不超过136的蛋白质作为训练集。之所以选择长度为136的残基序列，一方面是考虑计算力的问题，另一方面是因为后面计算二面角时采用了宽度为34的剪辑窗口，136是34的整数倍。
为了组织项目文件结构，按照如下步骤创建子目录。
（1） 在chapter5目录下新建两个子目录： preprocessing和models。preprocessing用于数据预处理和特征提取，models用于存放距离预测模型和角度预测模型。
（2） 在models目录下创建angles和distance两个子目录，分别用于存放角度预测和距离预测的相关模型程序。
在preprocessing目录下创建程序calculate_aa_distance.ipynb。执行程序段P5.1导入库。



P5.1#导入库

001import numpy as np

002import matplotlib.pyplot as plt

003import os




执行程序段P5.2，读取数据集。



P5.2#读入文件

004filename = '../dataset/training_30'

005with open(filename, 'r') as f:

006lines = f.readlines()#读取所有行

007print('文件中共有 {0} 行'.format(len(lines)))




运行结果显示： 文件中共有340989行数据。
由于［EVOLUTIONARY］字段域和［TERTIARY］字段域都包含多行数据，这些数据需要解析为浮点型数据列表，所以用程序段P5.3定义一个行数据解析函数。



P5.3#行数据析取函数

008def lines_split(lines, splice) :

   '''

   功能： 按照分隔符splice，将字符行lines中的数值拆分为列表，数据转为float类型

   参数： 

   lines:蛋白质结构文件中的若干行

   splice:分隔符

   返回值： 

   return： 浮点型数据列表

   '''

009data = ［］

010for line in lines :

011data.append(［float(x) for x in line.split(splice)］)

012return data




执行程序段P5.4，读取数据集中所有蛋白质的结构信息，分类存放。



P5.4#读取数据集中所有的蛋白质结构信息

013names = ［］#蛋白质ID

014seqs = ［］   #氨基酸序列

015coords = ［］  #主链上原子的坐标

016pssms = ［］  #PSSM评分矩阵

017masks = ［］  #掩码序列

018for i in range(len(lines)):

019if lines［i］ == '［ID］＼n' :

020names.append(lines［i+1］) 

021elif lines［i］ == '［PRIMARY］＼n' :

022seqs.append(lines［i+1］)

023elif lines［i］ == '［TERTIARY］＼n' :

024coords.append(lines_split(lines［i+1:i+4］, "＼t"))

025elif lines［i］ == '［EVOLUTIONARY］＼n' :

026pssms.append(lines_split(lines［i+1:i+22］, "＼t"))

027elif lines［i］ == '［MASK］＼n' :

028masks.append(lines［i+1］)

029print('数据集中包含{0}个蛋白质序列'.format(len(seqs)))

030print('第1个蛋白质的序列长度为： {0}'.format(len(seqs［0］)-1))

031print('第1个蛋白质主链原子坐标矩阵的维度为： {0}'.format(np.array(coords［0］).shape))

032print('第1个蛋白质的Mask掩码长度为： {0}'.format(len(masks［0］)-1))




程序运行结果如下。



数据集中包含10333个蛋白质序列

第1个蛋白质的序列长度为： 307

第1个蛋白质主链原子坐标矩阵的维度为： (3, 921)

第1个蛋白质的Mask掩码长度为： 307




执行程序段P5.5，按照残基序列长度统计training_30数据集中的蛋白质分布，如图5.9所示，该图给出的是残基序列长度在6个范围段上的汇总情况。



P5.5#统计蛋白质分布

033under64,under128,under136,under200,under300,above300 = 0,0,0,0,0,0

034for i in range(len(seqs)):

035if (len(seqs［i］) -1<= 64) :

036under64 +=1

037elif (len(seqs［i］) -1<= 128) :

038under128 +=1

039elif (len(seqs［i］) -1<= 136) :

040under136 +=1

041elif (len(seqs［i］) -1<= 200) :

042under200 +=1

043elif (len(seqs［i］) -1<= 300) :

044under300 +=1

045else :

046above300 +=1

047x = ［'<=64','(64, 128］','(128, 136］','(136, 200］','(200, 300］','>300'］

048y = ［under64, under128, under136, under200, under300, above300］

049plt.ylabel('Counts',fontsize=15)

050plt.title('Amino Acid Sequence Length Counts',fontsize=15)

051plt.bar(x, y, alpha=0.7)






图5.9training_30数据集的蛋白质序列分布


执行程序段P5.6，删除连续两个原子坐标为0的蛋白质，因为连续坐标为0，说明相邻原子数据缺失，无法计算角度。



P5.6#删除连续两个原子坐标为0的蛋白质

052flag = 1

053index = ［］

054for i in range(len(seqs)):#遍历所有蛋白质  

055for j in range(len(seqs［i］)):#第j个原子与第j+1个原子坐标均为0

056if ((coords［i］［0］［j］ == 0.0 and coords［i］［0］［j+1］ == 0.0) and 

057(coords［i］［1］［j］ == 0.0 and coords［i］［1］［j+1］ == 0.0) and 

058(coords［i］［2］［j］ == 0.0 and coords［i］［2］［j+1］ == 0.0)):

059flag = 0

060break

061if (flag == 1): 

062index.append(i)

063flag = 1   

064print('筛选前包含的蛋白质数量： {0}'.format(len(names)))

065print('筛选后包含的蛋白质数量： {0}'.format(len(index)))




运行结果如下。



筛选前包含的蛋白质数量： 10333

筛选后包含的蛋白质数量： 6241




执行程序段P5.7，根据残基序列长度和MASK筛选蛋白质。



P5.7#根据氨基酸序列长度筛选序列，长度小于或等于L，并且Mask全部为"+"

066print("原有的蛋白质结构总数:{0} ".format(len(seqs)))

067L = 136 #指定长度

068under = ［］

069for i in range(len(index)):

070k = index［i］

071if (len(seqs［k］-1) <= L) and (masks［k］.find('-') < 0):

072under.append(k)#记录满足条件的蛋白质序号

073print('氨基酸序列长度小于或等于 {0} 且 Mask全部为 ＼"+ ＼" 的蛋白质总数： 

{1}'.format( L, len(under)))




运行结果如下。



原有的蛋白质结构总数:10333

氨基酸序列长度小于或等于 136 且 Mask全部为 "+ " 的蛋白质总数： 2919




如果计算力不够，可以将L设置为更小的值，例如L=68，正好是剪辑窗口宽度34的两倍。
5.6计算残基之间的距离
Ni、Ciα、Ci表示第i个残基的三个原子，Ni-1、Ci-1α、Ci-1表示第i-1个残基的三个原子，Ni+1、Ci+1α、Ci+1表示第i+1个残基的三个原子，这些原子都分布在肽链的主链上。第i个残基与第i-1个残基之间的距离，定义为Ciα与Ci-1α两个原子之间的距离，第i个残基与第i+1个残基之间的距离，定义为Ciα与Ci+1α两个原子之间的距离，如图5.10所示，虚线表示相邻两个残基之间的距离。


图5.10两个残基之间的距离


如图5.10所示，Ni-Ci-1是连接两个氨基酸的肽键，Ci-Ni+1也是肽键，肽键不能自由旋转，所以肽键连接的6个原子(Ci-1α、Ci-1、Oi-1、Hi、Ni、Ciα)处于同一平面内。执行程序段P5.8，完成残基对之间的距离计算。



P5.8#根据氨基酸C-alpha原子的坐标计算氨基酸残基对之间的距离

074dists = ［］

075protein_number = len(under)

076for k in range(protein_number) :

077dist = ［］

078key = under［k］

079coords_length = np.array(coords［key］).shape［1］

080for i in range(coords_length) :   

081if i%3 == 1 :  #每个氨基酸在主链上有N-Calpha-C三个原子，只处理Calpha
#的坐标 

082aad = ［］ #存储第i个AA到所有AA之间的距离

083for j in range( coords_length) :

084if j%3 == 1 :   #只处理Calpha的坐标

#第i个AA的Calpha的坐标 

085coordi = np.array(［coords［key］［0］［i］,coords［key］［1］［i］,coords［key］［2］［i］］)

#第j个AA的Calpha的坐标 

086coordj = np.array(［coords［key］［0］［j］,coords［key］［1］［j］,coords［key］［2］［j］］)

087aa_dist = np.sqrt(np.sum(np.square(coordi - coordj)))#计算距离

088aad.append(aa_dist)  #第i个AA到第j个AA的距离加入列表

089dist.append(aad) #第i个AA到所有AA的距离加入列表

090dists.append(dist)#当前序列的AA距离加入总列表




执行程序段P5.9，对前面计算的数据做一个检测。



P5.9#数据测试与检验

091n = 0

092key = under［n］

093print(key)

094print("蛋白质ID: ", names［key］)

095print("残基序列: ", seqs［key］) 

096print("序列长度", len(seqs［key］)-1)

097print("当前序列第1个残基N原子的坐标(x,y,z)=［{0}, {1}, {2}］".

format(coords［key］［0］［0］, coords［key］［1］［0］, coords［key］［1］［0］))

098print("当前序列第1个残基Calpha原子的坐标(x,y,z)=［{0}, {1}, {2}］".

format(coords［key］［0］［1］, coords［key］［1］［1］, coords［key］［1］［1］))

099print("当前序列第1个残基C原子的坐标(x,y,z)=［{0}, {1}, {2}］".

format(coords［key］［0］［2］, coords［key］［1］［2］, coords［key］［1］［2］))

100print("当前序列第1个残基与第9个残基之间的距离: ", dists［0］［0］［9］)

101print("当前序列生成的距离矩阵的维度： ", np.array(dists［0］).shape)

102print("当前序列的最大距离： ",np.array(dists［0］).max())

103print("当前序列的最小距离： ",np.array(dists［0］).min())

104print("当前序列的平均距离： ",np.array(dists［0］).mean())




运行结果如下。



key= 1

蛋白质ID:  2EUL_d2euld1

残基序列:  MAREVKLTKAGYERLMQQLERERERLQEATKILQELMESSDDY

DDSGLEAAKQEKARIEARIDSLEDILSRAVILEE

序列长度 77

当前序列第1个残基N原子的坐标(x,y,z)=［981.8, 4076.1, 4076.1］

当前序列第1个残基Calpha原子的坐标(x,y,z)=［1093.7, 4174.2, 4174.2］

当前序列第1个残基C原子的坐标(x,y,z)=［1219.7, 4106.1, 4106.1］

当前序列第1个残基与第9个残基之间的距离:  2526.233815386058

当前序列生成的距离矩阵的维度：  (77, 77)

当前序列的最大距离：  5536.502383274119

当前序列的最小距离：  0.0

当前序列的平均距离：  2063.4028892563697




执行程序段P5.10，构建训练集文件distance_aa_under136.txt，存放到dataset/distances/目录下。



P5.10#将处理完成的数据保存到新文件中

105out_path = '../dataset/distances/'

106if not os.path.exists(out_path):

107os.makedirs(out_path)

108filename = out_path + 'distance_aa_under136.txt'

109with open(filename, 'w') as f:

110for k in range(len(under)) :

111key = under［k］

#保存ID






112f.write("＼n［ID］＼n")

113f.write(names［key］)

#保存序列

114f.write("＼n［PRIMARY］＼n")

115f.write(seqs［key］)

#保存PSSM矩阵，PSSM有21行

116f.write("＼n［EVOLUTIONARY］＼n")

117for line in range(21):

118for item in pssms［key］［line］:

119f.write(str(item)+'＼t')

120f.write('＼n')

#保存坐标矩阵，坐标有三行，分别代表x,y,z

121f.write("＼n［TERTIARY］＼n")

122for line in range(3):  

123for item in coords［key］［line］:

124f.write(str(item)+'＼t')

125f.write('＼n')

#保存距离矩阵，距离矩阵的维度与序列长度L有关： (L，L)

126f.write("＼n［DIST］＼n")

127for line in range(len(dists［k］)):

128for item in dists［k］［line］:

129f.write(str(item)+'＼t')

130f.write('＼n')




5.7二面角与拉氏构象图
肽键连接的6个原子处于同一平面，如图5.11所示，Ca、C、N、H、O、Ca六个原子共处同一平面。肽键中的CN键具有部分双键性质，不能自由旋转，结果使肽键处在一个刚性的平面上，此平面被称为肽键平面(酰胺平面)。
两个肽键平面之间的α碳原子，可以作为一个旋转点形成二面角。二面角的变化，决定着多肽主键在三维空间的排列方式，是形成不同蛋白质构象的基础。
组成肽键的6个原子处于同一平面，每个残基的CαN键和CαC键是单键，可以相对自由旋转，围绕CαN键轴旋转产生的角度称为Phi（Φ），围绕CαC键轴旋转产生的角度称为Psi（Ψ）。
根据拉氏构象图（Ramachandran plot）原理，二面角（Φ、Ψ）的变化范围为-180°～180°，但是Φ、Ψ不能任意取值，有部分（Φ、Ψ）构象是不存在的，根据原子的范德华半径，在旋转过程中，原子间会发生碰撞冲突，无法实现（Φ、Ψ）的自由变化，如图5.12所示的拉氏构象图，以Φ为横坐标，以Ψ为纵坐标，显示了（Φ、Ψ）的分布规律。
图5.12显示了100000个（Φ、Ψ）数据点的分布，每个数据点代表出现在单个氨基酸中的Φ和Ψ角度的组合。这些数据点不包括甘氨酸、脯氨酸，因为这两种氨基酸的拉氏构象分布具有自己的特点。



图5.11二面角Phi（Φ）和Psi（Ψ）




图5.12二面角拉氏构象图



图5.12中α螺旋构象的二面角区域标记为α，β折叠构象的二面角区域标记为β。右上象限中的数据簇主要表示残基之间的转角。
5.8计算二面角Phi（Φ）和Psi（Ψ）
在preprocessing目录下创建程序calculate_phi_psi.ipynb。执行程序段P5.11导入库。



P5.11#导入库

131import numpy as np

132import matplotlib.pyplot as plt




执行程序段P5.12，定义函数，解析行数据。



P5.12#定义函数，将行数据解析为float型数据列表

133def parse_line(raw):

134return np.array(［［float(x) for x in line.split("＼t") if x != ""］ for line in raw］)




执行程序段P5.13，读取数据集中的蛋白质ID、残基序列、PSSM评分矩阵和原子坐标。



P5.13#读取数据集

135names = ［］  #ID

136seqs = ［］   #序列

137psis = ［］   #Psi

138phis = ［］   #Phi

139pssms = ［］  #PSSM

140coords = ［］ #坐标






141path = "../dataset/distances/distance_aa_under136.txt"

142with open(path, "r") as f:

143lines = f.read().split('＼n')#读取文件内容，按行拆分

#从文本文件析取蛋白质的ID、序列、PSSM和坐标

144for i,line in enumerate(lines):

145if line == "［ID］":

146names.append(lines［i+1］)

147elif line == "［PRIMARY］":

148seqs.append(lines［i+1］)

149elif line == "［EVOLUTIONARY］":

150pssms.append(parse_line(lines［i+1:i+22］))  #PSSM

151elif line == "［TERTIARY］":

152coords.append(parse_line(lines［i+1:i+4］))  #坐标

153print (len(names))




运行结果显示数据集中有2919个序列长度不超过136的蛋白质。
执行程序段P5.14，定义函数获取单种类型原子的坐标。



P5.14#定义函数获取单种类型原子的坐标

154def separate_coords(full_coords, pos):

'''

full_coords: 主链上所有原子的坐标

pos:  0： n_term原子, 1： calpha原子, 2： cterm原子

return: 单种类型原子的坐标

'''

155res = ［］

156for i in range(len(full_coords［1］)):

157if i%3 == pos:

158res.append(［full_coords［j］［i］ for j in range(3)］)  #坐标(x,y,z)

159return np.array(res)

#根据原子类型分离坐标

160coords_nterm = ［separate_coords(full_coords, 0) for full_coords in coords］ #n_term原子

161coords_calpha = ［separate_coords(full_coords, 1) for full_coords in coords］ #calpha原子

162coords_cterm = ［separate_coords(full_coords, 2) for full_coords in coords］  #cterm原子

#坐标检查

163print("第一个蛋白质第1个残基的三个原子N-Ca-C的坐标： ")

164print(coords_nterm［0］［0］)

165print(coords_calpha［0］［0］)

166print(coords_cterm［0］［0］)




运行结果如下。



第一个蛋白质第1个残基的三个原子N-Ca-C的坐标： 

［  981.8  4076.1 -7423.1］

［ 1093.7  4174.2 -7419.3］

［ 1219.7  4106.1 -7366.7］




执行程序段P5.15，定义函数，用向量法计算二面角。



P5.15#定义函数，用向量法计算二面角

167def get_dihedral(coords1, coords2, coords3, coords4):

"""

coords1, coords2, coords3, coords4： 对应主链上四个原子的坐标

计算phi角时： C-N-Calpha-C

计算psi角时： N-Calpha-C-N

Returns： 返回phi角或psi角

"""

168a1 = coords2 - coords1

169a2 = coords3 - coords2

170a3 = coords4 - coords3

171v1 = np.cross(a1, a2)

172v1 = v1 / (v1 * v1).sum(-1)**0.5

173v2 = np.cross(a2, a3)

174v2 = v2 / (v2 * v2).sum(-1)**0.5

175porm = np.sign((v1 * a3).sum(-1))

176rad = np.arccos((v1*v2).sum(-1) / ((v1**2).sum(-1) * (v2**2).sum(-1))**0.5)

177if not porm == 0:

178rad = rad * porm

179return rad




执行程序段P5.16，完成所有蛋白质序列主链上二面角的计算。



P5.16#完成所有蛋白质序列主链上二面角的计算

180phis, psis = ［］, ［］ 

181ph_angle_dists, ps_angle_dists = ［］, ［］

#遍历每一个蛋白质

182for k in range(len(coords)):

183phi, psi = ［0.0］, ［］   # phi从0开始，psi以0结束

#遍历每一个氨基酸

184for i in range(len(coords_calpha［k］)):

#计算 phi，第1个残基不计算phi

185if i>0:

186phi.append(get_dihedral(coords_cterm［k］［i-1］, coords_nterm［k］［i］, 

coords_calpha［k］［i］, coords_cterm［k］［i］))   

#计算psi，最后1个残基不计算psi

187if i<len(coords_calpha［k］)-1: 

188psi.append(get_dihedral(coords_nterm［k］［i］, coords_calpha［k］［i］, 

coords_cterm［k］［i］, coords_nterm［k］［i+1］))

#最后一个残基的psi为0

189psi.append(0)

#添加到列表

190phis.append(phi)

191psis.append(psi)




执行程序段P5.17，绘制Ramachandran图形，观察100个蛋白质的Phi和Psi分布。



P5.17#绘制Ramachandran图形，观察100个蛋白质的Phi和Psi分布

192n = 100#蛋白质的数量

193test_phi = ［］

194for i in range(n):

195for test in phis［i］:

196test_phi.append(test)

197test_phi = np.array(test_phi)

198test_psi = ［］

199for i in range(n):

200for test in psis［i］:

201test_psi.append(test)

202test_psi = np.array(test_psi)

#按照象限统计Phi和Psi角度分布

203quads = ［0,0,0,0］

204for i in range(len(test_phi)):

205if test_phi［i］ >= 0 and test_psi［i］ >= 0:

206quads［0］ += 1

207elif test_phi［i］ < 0 and test_psi［i］ >= 0:

208quads［1］ += 1

209elif test_phi［i］ < 0 and test_psi［i］ < 0:

210quads［2］ += 1

211else:

212quads［3］ += 1

213print("象限分布: ", quads, " 总数： ", len(test_phi))

#绘制 Ramachandran 分布图

214plt.scatter(test_phi, test_psi, marker=".")

215plt.xlim(-np.pi, np.pi)

216plt.ylim(-np.pi, np.pi)

217plt.xlabel("Phi")

218plt.ylabel("Psi")

219plt.show()




运行结果显示这100个蛋白质共包含3980对Phi和Psi角度，一、二、三、四象限分布数量分别为294、1585、1965、136，对应的拉氏构象如图5.13所示。


图5.13100个蛋白质的拉氏构象分布


执行程序段P5.18，定义函数将向量转换为字符串。



P5.18#定义函数，将向量格式转换为字符串行

220def stringify(vec):

221line = ""

222for v in vec:

223line = line+str(v)+" "

224return line

225print(［stringify(［1,2,3,4,5,6］)］)




执行程序段P5.19，保存文件，包含蛋白质ID、氨基酸序列、评分矩阵以及二面角Phi和Psi。



P5.19#保存文件，包含蛋白质ID、氨基酸序列、评分矩阵以及二面角Phi和Psi

226with open("../dataset/angles/angles_phi_psi_under64.txt", "w") as f:

227for k in range(len(names)):

228f.write("＼n［ID］＼n")

229f.write(names［k］)#ID

230f.write("＼n［PRIMARY］＼n") 

231f.write(seqs［k］)#序列

232f.write("＼n［EVOLUTIONARY］＼n")

233for j in range(len(pssms［k］)): 

234f.write(stringify(pssms［k］［j］)+"＼n")#PSSM

235f.write("［PHI］＼n")

236f.write(stringify(phis［k］))#PHI

237f.write("＼n［PSI］＼n")

238f.write(stringify(psis［k］)+'＼n')#PSI




5.9裁剪残基序列的OneHot矩阵
在preprocessing目录下创建程序prepare_angle_data.ipynb。执行程序段P5.20，将行数据解析为列表。



P5.20#将行数据解析为列表

239import numpy as np

240def parse_lines(raw):

241return np.array(［［float(x) for x in line.split(" ") if x != ""］ for line in raw］)

242def parse_line(line):

243return np.array(［float(x) for x in line.split(" ") if x != ""］)




执行程序段P5.21，读取5.8节完成的Phi和Psi数据集，包含2919个蛋白质的相关数据。



P5.21#打开文件，读取Phi和Psi数据集

244path = "../dataset/angles/angles_phi_psi_under136.txt"

245with open(path, "r") as f:

246lines = f.read().split('＼n')






247names = ［］

248seqs = ［］

249psis = ［］

250phis = ［］

251pssms = ［］

#读取每一个蛋白质的结构信息

252for i,line in enumerate(lines):  

253if line == "［ID］":

254names.append(lines［i+1］)

255elif line == "［PRIMARY］":

256seqs.append(lines［i+1］)

257elif line == "［EVOLUTIONARY］":

258pssms.append(parse_lines(lines［i+1:i+22］))

259elif lines［i］ == "［PHI］":

260phis.append(parse_line(lines［i+1］))

261elif lines［i］ == "［PSI］":

262psis.append(parse_line(lines［i+1］))

263print(len(names))




DNA翻译为蛋白质的编码图谱如图5.14所示，根据这个图谱不难看出，自然界已知氨基酸（不考虑最新发现的）只有20种，将这20种氨基酸按照一定顺序排列，例如HRKDENQSYTCPAVLIGFWM。根据这个字母排列可以定义残基序列的OneHot编码。


图5.14DNA翻译为蛋白质的编码图谱

氨基酸序列通过肽键构成了肽链，肽链折叠形成蛋白质。根据氨基酸序列预测二面角，需要首先完成残基序列的OneHot编码。以数据集中ID为1KX6_1_A的蛋白质序列（长度为29）为例，对应的OneHot编码如图5.15所示。


图5.15残基序列的OneHot编码


执行程序段P5.22，完成残基序列片段的OneHot编码和范德华半径编码。单个残基的OneHot编码为（20，1）的向量，后面加上范德华半径和范德华表面半径，所以维度为（22，1），并对给定的残基序列seq，从pos位置向左向右按照（34，22）的窗口裁剪。



P5.22#残基序列片段的One-Hot编码和范德华半径编码，并裁剪

264def onehotter_aa(seq, pos):

'''

功能： 每次截取长度为34(17×2个残基)的氨基酸序列片段

seq： 氨基酸序列

pos:截取片段的中心位置，截取起点为pos-17，截取终点为pos+17

return： 截取片段的One-Hot编码

'''

265pad = 17

#定义20种氨基酸序列表

266key = "HRKDENQSYTCPAVLIGFWM"

#20种氨基酸的范德华半径

267vdw_radius = {"H": 118, "R": 148, "K": 135, "D": 91, "E": 109, "N": 96, "Q": 114,

"S": 73, "Y": 141, "T": 93, "C": 86, "P": 90, "A": 67, "V": 105,

"L": 124, "I": 124, "G": 48, "F": 135, "W": 163, "M": 124}

268radius_rel = vdw_radius.values()

269basis = min(radius_rel)/max(radius_rel)

#20种氨基酸的范德华表面半径 

270surface = {"H": 151, "R": 196, "K": 167, "D": 106, "E": 138, "N": 113, "Q": 144,

"S": 80, "Y": 187, "T": 102, "C": 104, "P": 105, "A": 67, "V": 117,

"L": 137, "I": 140, "G": 0, "F": 175, "W": 217, "M": 160}






271surface_rel = surface.values()

272surface_basis = min(surface_rel)/max(surface_rel)

#One-Hot 编码

273one_hot = ［］

274for i in range(pos-pad, pos+pad): #遍历截取的片段，片段长度为 pad×2

275vec = ［0 for i in range(22)］#先将One-Hot编码置为0

#将当前氨基酸与20种氨基酸序列表一一比对

276for j in range(len(key)):

277if seq［i］ == key［j］:

278vec［j］ = 1 #此处标定为1，表示氨基酸类型

#在One-Hot编码的末尾添加范德华半径和范德华表面半径，并归一化

279vec［-2］ = vdw_radius［key［j］］/max(radius_rel)-basis

280vec［-1］ = surface［key［j］］/max(surface_rel)-surface_basis

281one_hot.append(vec)#将当前氨基酸的维度为(1,22)的编码向量加入列表one_hot

282return np.array(one_hot)#One_Hot的维度为(34,22)





5.10裁剪评分矩阵和二面角标签
本节对氨基酸序列片段进行裁剪，完成模型训练的数据准备工作，主要包括： 
（1） 氨基酸序列的裁剪矩阵，单个片段维度为（34，22），如图5.16（a）所示。
（2） 评分裁剪矩阵，单个片段裁剪维度为（34，21），如图5.16（b）所示。
（3） 定义二面角标签，单个标签维度为（1，2），如图5.16（c）所示。


图5.16二面角训练集的裁剪方法

执行程序段P5.23，裁剪PSSM评分矩阵。



P5.23#裁剪PSSM评分矩阵

283def pssm_cropper(pssm, pos):

'''

功能： 每次截取长度为34(17×2个残基)的PSSM评分序列片段

pssm： 评分矩阵

pos:截取片段的中心位置，截取起点为pos-17，截取终点为pos+17

return： 截取的片段矩阵，维度为(21,34)

'''

284pssm_out = ［］

285pad = 17

286for i,row in enumerate(pssm):

287pssm_out.append(row［pos-pad:pos+pad］)

288return np.array(pssm_out)

#检查所有的矩阵是否包含相同的蛋白质结构数

289print("Names: ", len(names))

290print("Seqs: ", len(seqs))

291print("PSSMs: ", len(pssms))

292print("Phis: ", len(phis))

293print("Psis: ", len(psis))




运行结果如下。



Names:  2919

Seqs:  2919

PSSMs:  2919

Phis:  2919

Psis:  2919




执行程序段P5.24，裁剪氨基酸序列、评分矩阵和二面角矩阵。



P5.24#裁剪氨基酸序列、评分矩阵和二面角矩阵

294input_aa = ［］

295input_pssm = ［］

296outputs = ［］

297count = 0 #计数需要裁剪多少次

#遍历每一个蛋白质结构

298for i in range(len(seqs)): 

299if len(seqs［i］)>17*2:  #序列长度大于34

300count += len(seqs［i］)-17*2

301for j in range(17,len(seqs［i］)-17):  #设定pos的位置j

#根据pos的位置j向左向右裁剪序列

302input_aa.append(onehotter_aa(seqs［i］, j))

303input_pssm.append(pssm_cropper(pssms［i］, j))

304outputs.append(［phis［i］［j］, psis［i］［j］］)

305print("全部的蛋白序列，共裁剪 {0} 次，这是样本的数量".format(count))

#得到的矩阵长度

306print("标签矩阵outputs的长度: ", len(outputs))

307print("残基片段矩阵Inputs_aa的长度: ", len(input_aa))

308print("评分片段矩阵input_pssm的长度: ", len(input_pssm))




运行结果如下。



全部的蛋白序列，共裁剪 140756 次，这是样本的数量

标签矩阵outputs的长度:  140756

残基片段矩阵Inputs_aa的长度:  140756

评分片段矩阵input_pssm的长度:  140756




执行程序段P5.25，调整特征矩阵的维度。



P5.25#调整特征矩阵的维度

309input_aa = np.array(input_aa).reshape(len(input_aa), 17*2, 22)

310print('氨基酸序列片段矩阵的维度： {0}'.format(input_aa.shape))

311input_pssm = np.array(input_pssm).reshape(len(input_pssm), 17*2, 21)

312print('评分片段矩阵的维度： {0}'.format(input_pssm.shape))




运行结果如下。



氨基酸序列片段矩阵的维度： (140756, 34, 22)

评分片段矩阵的维度： (140756, 34, 21)




执行程序段P5.26，将向量转换为字符串，以空格分隔数据项。



P5.26#将向量转换为字符串，以空格分隔数据项

313def stringify(vec):

314return "".join(str(v)+" " for v in vec)




执行程序段P5.27，保存二面角训练集的标签，即将Phi和Psi保存到文件outputs.txt。



P5.27#保存二面角训练集的标签文件

315out_path = '../dataset/angles/'

316with open(out_path + "outputs.txt", "w") as f:

317for o in outputs:

318f.write(stringify(o)+"＼n")




执行程序段P5.28，将残基片段序列的矩阵保存到文件input_aa.txt。



P5.28#将残基片段序列的矩阵保存到文件input_aa.txt

319with open(out_path + "input_aa.txt", "w") as f:

320for aas in input_aa:

321f.write("＼nNEW＼n")

322for j in range(len(aas)):

323f.write(stringify(aas［j］)+"＼n")




执行程序段P5.29，将PSSM片段序列的矩阵保存到文件input_pssm.txt。



P5.29#将PSSM片段序列的矩阵保存到文件input_pssm.txt

324with open(out_path + "input_pssm.txt", "w") as f:

325for k in range(len(input_pssm)):






326f.write("＼nNEW＼n")

327for j in range(len(input_pssm［k］)):

328f.write(stringify(input_pssm［k］［j］)+"＼n")




5.11定义二面角预测模型
二面角预测模型采用ResNet结构，卷积层采用一维卷积，残差块结构如图5.17所示。


图5.17残差块结构


二面角预测模型结构如图5.18所示。样本的特征矩阵维度为（34，42），标签维度为（1，4）。标签由原来的（Φ，Ψ）调整为（sinΦ，cosΦ，sinΨ，cosΨ），可以借助角度的正弦值和余弦值，将标签值归一化为［-1，1］的值。


图5.18二面角预测模型结构


模型定义存放于models/angles/resnet_1d_angles.py文件，程序源码如程序段P5.30所示，将在二面角模型训练和预测程序models/angles/predicting_angles.ipynb中对其引用。



P5.30#ResNet 1D卷积模型定义

329import keras

330from keras.models import Model

331from keras.regularizers import l2

332from keras.losses import mean_squared_error, mean_absolute_error

333from keras.layers.convolutional import Conv1D

334from keras.layers import Dense, Flatten, Input, BatchNormalization, Activation

335from keras.layers.pooling import AveragePooling1D

336def custom_mse_mae(y_true, y_pred):

""" 自定义损失函数 - MSE + MAE """

337return mean_squared_error(y_true, y_pred)+mean_absolute_error(y_true, y_pred)

#定义1D卷积层






338def resnet_layer(inputs,

num_filters=16,

kernel_size=3,

strides=1,

activation='relu',

batch_normalization=True,

conv_first=False):

"""卷积层有两种设计顺序，即： BN-Relu-Conv 或者 Conv-BN-Relu



#参数

inputs (tensor): 输入层或来自上一层的输入

num_filters (int): 过滤器数量

kernel_size (int): 过滤器尺寸

strides (int): 步长

activation (string): 激励函数名称

batch_normalization (bool): 是否包含BN层

conv_first (bool): conv-bn-activation (True) 或者 bn-activation-conv (False)

#Returns

x (tensor): 输出向量x，作为下一层的输入

"""

#定义1D卷积







339conv = Conv1D(num_filters,

 kernel_size=kernel_size,

 strides=strides,

 padding='same',

 kernel_initializer='he_normal',

 kernel_regularizer=l2(1e-4))

340x = inputs

341if conv_first:

342x = conv(x)

343if batch_normalization:

344x = BatchNormalization()(x)

345if activation is not None:

346x = Activation(activation)(x)

347else:

348if batch_normalization:

349x = BatchNormalization()(x)

350if activation is not None:

351x = Activation(activation)(x)

352x = conv(x)

353return x

#定义resnet_v2的1D卷积网络

354def resnet_v2(input_shape, depth, num_classes=4, conv_first=True):

"""ResNet Version 2 卷积网络的结构： 

瓶颈残差块采用 (1×1)-(3×3)-(1×1)卷积堆叠

每个stage第1个残差块的直连采用 1×1 卷积做维度变换

后续残差块采用恒等变换

在stage的第1层，采用strides=2的步长下采样将特征图尺寸减半，过滤器数量翻倍

同一stage内部，过滤器数量和特征图尺寸不变

#参数

input_shape (tensor): 输入层的维度

depth (int): 核心卷积层的数量

num_classes (int): 类别的数量

#Returns

model (Model): Keras model实例

"""

355if (depth - 2) % 9 != 0:

356raise ValueError('网络深度depth按照公式 depth = 9n+2 计算 例如depth=56或110')

#初始化模型参数

357num_filters_in = 16

358num_res_blocks = int((depth - 2) / 9)

359inputs = Input(shape=input_shape)

#resnet_v2 第一层采用 conv-bn-activation 卷积

360x = resnet_layer(inputs=inputs,

num_filters=num_filters_in,

conv_first=True)

#一个resnet分为多个stage,每个stage中有多个block，每个block中包含多个层

361for stage in range(3):

#残差块block







362for res_block in range(num_res_blocks):   

363activation = 'relu'

364batch_normalization = True

365strides = 1

366if stage == 0:

367num_filters_out = num_filters_in * 4#第1个stage的过滤器数量

368if res_block == 0:#第1个stage的第1层

369activation = None

370batch_normalization = False

371else:

372num_filters_out = num_filters_in * 2

373if res_block == 0:  #非第1个stage的第1层

374strides = 2#下采样

#定义瓶颈残差块

#残差块第1层： 1×1的1D卷积

375y = resnet_layer(inputs=x,

num_filters=num_filters_in,

kernel_size=1,

strides=strides,

activation=activation,

batch_normalization=batch_normalization,

conv_first=conv_first)

#残差块第2层： 3×3的1D卷积

376y = resnet_layer(inputs=y,

num_filters=num_filters_in,

conv_first=conv_first)

#残差块第3层： 1×1的1D卷积

377y = resnet_layer(inputs=y,

num_filters=num_filters_out,

kernel_size=1,

conv_first=conv_first)

378if res_block == 0:

#不同的stage之间需要做维度变换，用1×1卷积完成

379x = resnet_layer(inputs=x,

num_filters=num_filters_out,

kernel_size=1,

strides=strides,

activation=None,

batch_normalization=False)

380x = keras.layers.add(［x, y］) #完成残差块直连

381num_filters_in = num_filters_out #后续每个stage，过滤器数量翻倍

382x = BatchNormalization()(x)

383x = Activation('relu')(x)

384x = AveragePooling1D(pool_size=3)(x)  #1D平均池化

385y = Flatten()(x)

#定义输出层

386outputs = Dense(num_classes,

activation='linear',

kernel_initializer='he_normal')(y)

#实例化模型

387model = Model(inputs=inputs, outputs=outputs)

388return model




5.12二面角模型参数设定与训练
执行程序段P5.31，导入库。



P5.31#导入库

389import numpy as np

390import matplotlib.pyplot as plt

391from keras.callbacks import EarlyStopping

#导入ResNet 1D模型结构

392from resnet_1d_angles import *




执行程序段P5.32，读取Phi和Psi的标签文件，显示标签数据集的维度为(140756，2)。



P5.32#读取Phi和Psi的标签文件

393outputs = np.genfromtxt("../../dataset/angles/outputs.txt")

394outputs［np.isnan(outputs)］ = 0.0

395outputs.shape






图5.19训练集标签扩展为

（4，1）向量

计算角度的sin值和cos值，用Phi和Psi的正余弦作为标签，如图5.19所示，单个残基对应的标签向量由（2，1）扩展为（4，1）。
执行程序段P5.33，整个训练集的标签维度变为(140756，4)。



P5.33#计算角度的sin值和cos值，用Phi和Psi的正余弦作为标签

396out = ［］

397out.append(np.sin(outputs［:,0］))

398out.append(np.cos(outputs［:,0］))

399out.append(np.sin(outputs［:,1］))

400out.append(np.cos(outputs［:,1］))

401out = np.array(out).T

402print(out.shape)




执行程序段P5.34，定义函数，读取残基序列的特征集或评分矩阵的特征集。



P5.34#定义函数，读取残基序列的特征集或评分矩阵的特征集

403def get_ins(path = "../../dataset/angles/input_aa.txt", pssm=None):

404#pssm文件路径

405if pssm: path = "../../dataset/angles/input_pssm.txt"

406with open(path, "r") as f:

407lines = f.read().split('＼n')

408pre_ins = ［］

#遍历每一个样本的特征矩阵

409for i,line in enumerate(lines):






410if line == "NEW":

411prot = ［］

412raw = lines［i+1:i+(17*2+1)］

#对每一个特征行做解析

413for r in raw:

414prot.append(np.array(［float(x) for x in r.split(" ") if x != ""］))

#汇总每一个样本的特征矩阵

415pre_ins.append(np.array(prot))  

416return np.array(pre_ins)




二面角特征矩阵的构建方法如图5.20所示。


图5.20二面角训练集特征矩阵的构建方法


执行程序段P5.35，堆叠残基序列和PSSM评分矩阵。



P5.35#堆叠残基序列和PSSM评分矩阵

417aas = get_ins() #读取残基序列矩阵

418pssms = get_ins(pssm=True)#读取PSSM评分矩阵

#检查矩阵维度

419print('残基序列维度： ', aas.shape, '＼n评分矩阵维度： ',pssms.shape)

#特征矩阵堆叠

420inputs = np.concatenate((aas［:, :, :20］, pssms［:, :, :20］, aas［:, :, 20:］), axis=2) 

421print('堆叠矩阵维度： ', inputs.shape)




运行结果如下。



残基序列维度： (140756, 34, 22)

评分矩阵维度： (140756, 34, 21)

堆叠矩阵维度： (140756, 34, 42)




执行程序段P5.36，模型定义和编译，设定优化算法、学习率等超参数。



P5.36#模型定义和编译，设定超参数

422adam = keras.optimizers.Adam(lr=0.001, beta_1=0.9, beta_2=0.999, 

epsilon=1e-8, decay=0.0, amsgrad=True)#采用AMSGrad优化算法

#创建模型






423model = resnet_v2(input_shape=(17*2,42), depth=29, num_classes=4, conv_first=True)

424model.compile(optimizer=adam, loss=custom_mse_mae, metrics=［'accuracy'］)

425model.summary()




模型结构摘要显示，设定模型深度为29时，模型共有31个卷积层，需要学习训练的参数总量为473108个。
执行程序段P5.37，划分训练集与验证集，由于样本数量为14万条左右，所有设定训练集与验证集的划分比例为8∶2。



P5.37#训练集与验证集划分

426split = int(len(inputs)*0.8)

427x_train, x_val = inputs［:split］, inputs［split:］

428y_train, y_val = out［:split］, out［split:］




执行程序段P5.38，模型训练，设定epoch数量为20，minibatch大小为16，如果连续3次迭代的损失不下降，则提前停止模型训练。



P5.38#模型训练，如果连续3次迭代的损失不下降，则提前停止模型训练

429early_stopping = EarlyStopping(monitor='val_loss', patience=3)

430his = model.fit(x_train, y_train, epochs=20, batch_size=16, verbose=1, shuffle=True, 

validation_data=(x_val, y_val), callbacks=［early_stopping］)




训练集规模为112604个样本，验证集包含28152个样本。训练主机CPU配置为Intel CoreTM i76700 CPU@3.40Hz，内存配置为16GB。模型的训练提前终止于第20代，完成20代训练，用时1h11min31s。读者可根据自己的计算力，修正数据集规模。一般情况下，数据集越大，取得的效果会越好。
执行程序段P5.39，绘制准确率曲线如图5.21所示，绘制损失函数曲线如图5.22所示。



P5.39#绘制准确率和损失函数曲线

431x = range(1, len(his.history［'accuracy'］)+1)

432plt.plot(x, his.history［'accuracy'］)

433plt.plot(x, his.history［'val_accuracy'］)

434plt.title('Model Accuracy')

435plt.ylabel('Accuracy')

436plt.xlabel('Epoch')

437plt.xticks(x)

438plt.legend(［'Train', 'Val'］, loc='upper left')

439plt.show()

440plt.plot(x, his.history［'loss'］)

441plt.plot(x, his.history［'val_loss'］)

442plt.title('Model Loss')

443plt.ylabel('Loss')

444plt.xlabel('Epoch')

445plt.xticks(x)

446plt.legend(［'Train', 'Val'］, loc='lower left')

447plt.show()






图5.21训练集和验证集准确率曲线对比




图5.22训练集和验证集损失函数下降曲线对比



准确率曲线显示模型的准确率起点较高，随着迭代次数增加，呈缓慢上升趋势，过拟合现象不明显。
损失函数曲线显示模型在验证集和训练集上的趋势基本一致，方差较小，模型的泛化能力较好。
5.13二面角模型预测与评价
执行程序段P5.40，在验证集上做预测，并将预测后得到的sin(Φ)、cos(Φ)、sin(Ψ)、cos(Ψ)还原为角度值Φ和Ψ，维度变为(28152，2)。



P5.40#在验证集上做预测

448preds = model.predict(x_val)

#将验证集的预测结果sin值和cos值还原为二面角Phi和Psi

449refactor = ［］

450for pred in preds:

451angles = ［］

452phi_sin, phi_cos, psi_sin, psi_cos = pred［0］, pred［1］, pred［2］, pred［3］

#Phi分布在一、四象限

453if (phi_sin>=0 and phi_cos>=0) or (phi_cos>=0 and phi_sin<=0):

454angles.append(np.arctan(phi_sin/phi_cos))

455elif (phi_cos<=0 and phi_sin>=0) :

456angles.append(np.pi + np.arctan(phi_sin/phi_cos)) # Phi在第二象限

457else: 

458angles.append(-np.pi + np.arctan(phi_sin/phi_cos)) # Phi在第三象限

#Psi分布在一、四象限

459if (psi_sin>=0 and psi_cos>=0) or (psi_cos>=0 and psi_sin<=0):

460angles.append(np.arctan(psi_sin/psi_cos))

461elif (psi_cos<=0 and psi_sin>=0) :

462angles.append(np.pi + np.arctan(psi_sin/psi_cos))  # Psi在第二象限

463else:

464angles.append(-np.pi + np.arctan(psi_sin/psi_cos))  # Psi在第三象限

465refactor.append(angles)

#将预测的角度限定在(-pi, pi)

466refactor = np.array(refactor)

467refactor［refactor>np.pi］ = np.pi

468refactor［refactor<-np.pi］ = -np.pi

469print(refactor.shape)




执行程序段P5.41，绘制真实值与预测值的拉氏分布对比图，如图5.23所示。



P5.41#绘制真实值与预测值的拉氏分布对比图

470plt.scatter(outputs［split:,0］, outputs［split:,1］, marker=".")

471plt.scatter(refactor［:,0］, refactor［:,1］, marker=".")

472plt.legend(［"Truth distribution", "Predictions distribution"］, loc="lower right")

473plt.xlim(-np.pi, np.pi)

474plt.ylim(-np.pi, np.pi)

475plt.xlabel("Phi")

476plt.ylabel("Psi")

477plt.show()






图5.23验证集上的预测值与真实值对比


执行程序段P5.42，判断预测值与真实值的相关性。



P5.42#判断预测值与真实值的相关性

478cos_phi = np.corrcoef(np.cos(refactor［:,0］), np.cos(outputs［split:,0］))

479cos_psi = np.corrcoef(np.cos(refactor［:,1］), np.cos(outputs［split:,1］))

480print("理想的相关系数应该是：  Phi: 0.65 ，Psi: 0.7")

481print("模型的Phi相关系数: ", cos_phi［0,1］)

482print("模型的Psi相关系数: ", cos_psi［0,1］)




运行结果如下。



理想的相关系数应该是：  Phi: 0.65 ，Psi: 0.7

模型的Phi相关系数:  0.45145041419891035

模型的Psi相关系数:  0.5169855298544158




显然在2919个蛋白质上的训练结果还是比较理想的，增大蛋白质数据集规模，可望取得更好的效果。
执行程序段P5.43，保存训练模型。



P5.43#保存训练模型

483model.save("resnet_1d_angles.h5")




根据需要，可以执行程序段P5.44，加载预训练模型。



P5.44#加载预训练模型

484from keras.models import load_model

485model = load_model("resnet_1d_angles.h5", custom_objects=

{'custom_mse_mae': custom_mse_mae})




5.14定义距离预测模型
距离预测模型仍然采用ResNet结构，与二面角预测模型采用1D卷积不同，距离预测模型采用2D卷积，结构如图5.24所示。


图5.24距离预测模型结构


模型编码定义如程序段P5.45所示。



P5.45#定义距离预测模型

486import numpy as np

487import keras

488import keras.backend as K






489from keras.models import Model

490from keras.regularizers import l2

491from keras.activations import softmax

492from keras.layers.convolutional import Conv2D, Conv2DTranspose

493from keras.layers import Input, BatchNormalization, Activation

#softmax回归函数

494def softMaxAxis2(x):

""" 在axis 2维度上计算softmax """

495return softmax(x,axis=2)

#带权重的损失函数定义

496def weighted_categorical_crossentropy(weights):

"""

功能:根据类别的权重计算损失

参数:

weights: 权重系数向量，维度为(C,) ，C是类别的数量

"""

497weights = K.variable(weights)

498def loss(y_true, y_pred):

#预测值归一化

499y_pred /= K.sum(y_pred, axis=-1, keepdims=True)

#数据范围剪辑






500y_pred = K.clip(y_pred, K.epsilon(), 1 - K.epsilon())

#计算损失

501loss = y_true * K.log(y_pred) * weights

502loss = -K.sum(loss, -1)

503return loss

504return loss

#残差块定义

505def resnet_block(inputs,

num_filters=64,

kernel_size=3,

strides=1,

activation='elu',

batch_normalization=True,

conv_first=False):

"""#参数： 

inputs (tensor): 输入层或上一层的输入

num_filters (int): Conv2D过滤器数量

kernel_size (int): Conv2D过滤器尺寸

strides (int): Conv2D步长

activation (string): 激励函数名称

batch_normalization (bool):是否进行BN计算

conv_first (bool): conv-bn-activation (True) 或者 bn-activation-conv (False)

"""

506x = inputs

#下采样

507x = BatchNormalization()(x)

508x = Activation(activation)(x)

509x = Conv2D(num_filters//2, kernel_size=1, strides=1, padding='same',

kernel_initializer='he_normal', kernel_regularizer=l2(1e-4))(x)

#卷积计算

510x = BatchNormalization()(x)

511x = Activation(activation)(x)

512x = Conv2D(num_filters//2, kernel_size=3, strides=strides, padding='same',

kernel_initializer='he_normal', kernel_regularizer=l2(1e-4))(x)

#上采样

513x = BatchNormalization()(x)

514x = Activation(activation)(x)

515x = Conv2DTranspose(num_filters, kernel_size=1, strides=strides, padding='same',

kernel_initializer='he_normal', kernel_regularizer=l2(1e-4))(x)

516return x

#定义ResNet二维卷积层

517def resnet_layer(inputs,

num_filters=32,

kernel_size=3,

strides=1,

activation='relu',

batch_normalization=True,

conv_first=False):







"""

功能： 二维卷积： 按照 BN-Activation-Conv 顺序

Returns： x (tensor):输入下一层的向量

"""

518conv = Conv2D(num_filters,

 kernel_size=kernel_size,

 strides=strides,

 padding='same',

 kernel_initializer='he_normal',

 kernel_regularizer=l2(1e-4))

519x = inputs

520if conv_first:

521x = conv(x)

522if batch_normalization:

523x = BatchNormalization()(x)

524if activation is not None:

525x = Activation(activation)(x)

526else:

527if batch_normalization:

528x = BatchNormalization()(x)

529if activation is not None:

530x = Activation(activation)(x)

531x = conv(x)

532return x

#RestNet残差网络定义

533def resnet_v2(input_shape, depth, num_classes=4, conv_first=True):

""" 

采用ELU作为激励函数的ResNet，网络深度参数 depth 应是4的倍数

"""

534if depth % 4 != 0:

535raise ValueError('depth should be 4n (eg 8 or 16)')

#模型定义

536num_filters_in = 128

537inputs = Input(shape=input_shape)

#定义第一个卷积层，Conv-BN-ReLU

538x = resnet_layer(inputs=inputs,

num_filters=num_filters_in,

conv_first=True)

#定义残差块的卷积步长

539striding = ［1,2,4,8］

#遍历每一个stage

540for stage in range(depth):  

541activation = 'elu'

542batch_normalization = True

#瓶颈残差块

543y = resnet_block(inputs=x,

num_filters=128,

kernel_size=3,

strides=striding［stage%4］)   







544x = keras.layers.add(［x, y］)  #直连求和   

#顶层添加一个线性卷积分类器

545y = Conv2D(num_classes, kernel_size=1, strides=1, padding='same',

kernel_initializer='he_normal', kernel_regularizer=l2(1e-4))(x)

546outputs = Activation(softMaxAxis2)(y)

#初始化模型

547model = Model(inputs=inputs, outputs=outputs)

548return model




5.15构建残基序列3D特征矩阵
前面筛选蛋白质序列时，设定的蛋白质序列长度为136，考虑到距离模型的计算需求较大，因此本节将蛋白质序列的长度范围限定为不超过64。
重新打开preprocessing目录下的calculate_aa_distance.ipynb程序，将程序段P5.7中L的值由136修改为64，同时将程序段P5.10中的文件名称修改为distance_aa_under64.txt，重新运行程序calculate_aa_distance.ipynb，生成新的距离数据集。
在model/distances目录下新建程序文档predicting_distances.ipynb，执行程序段P5.46，导入库。



P5.46#导入库

549import numpy as np

550import matplotlib.pyplot as plt

551from keras.callbacks import EarlyStopping

#导入距离预测模型结构

552from elu_resnet_2d_distances import *




执行程序段P5.47，定义行解析函数，将行数据转换为float型数据向量。



P5.47#行解析函数，将行数据转换为float型数据向量

553def parse_lines(raw):

554return ［［float(x) for x in line.split("＼t") if x != ""］ for line in raw］

555def parse_line(line):

556return ［float(x) for x in line.split("＼t") if x != ""］




执行程序段P5.48，读取文件，解析为行的集合。为了降低计算需求，程序段P5.48指定从distance_aa_under64.txt中读取数据。



P5.48#读取文件，解析为行的集合

557path = "../../dataset/distances/distance_aa_under64.txt"

558with open(path, "r") as f:

559lines = f.read().split('＼n')




执行程序段P5.49，读取所有蛋白质的结构信息，结果显示数据集共包含890个蛋白质结构信息。



P5.49#读取所有蛋白质的结构信息

560names = ［］

561seqs = ［］

562dists = ［］

563pssms = ［］

#遍历数据集文件的每一行

564for i,line in enumerate(lines):

565if line == "［ID］":

566names.append(lines［i+1］)

567elif line == "［PRIMARY］":

568seqs.append(lines［i+1］)

569elif line == "［EVOLUTIONARY］":

570pssms.append(parse_lines(lines［i+1:i+21］))

571elif line == "［DIST］":

572dists.append(parse_lines(lines［i+1:i+len(seqs［-1］)+1］))

573print('蛋白质数量： {0}'.format(len(seqs)))




构建氨基酸序列的3D特征矩阵，即行列均为氨基酸序列，高和宽均为L，不足L的部分用0补齐。深度为氨基酸的OneHot编码向量，如图5.25所示。


图5.25残基序列的3D特征矩阵


执行程序段P5.50，将单个氨基酸序列转换为L×L×N的特征矩阵。



P5.50#将单个氨基酸序列转换为L×L×N的特征矩阵

574def wider(seq, L=64, N=20):

""" 将氨基酸序列转换为OneHot编码，维度为： L×L×N，不是 L×N  

seq： 氨基酸序列

L： 氨基酸序列的最大长度

N： 20种氨基酸序列长度

return： L×L×N特征矩阵

"""







575key = "HRKDENQSYTCPAVLIGFWM"#20种氨基酸

576tensor = ［］

577for i in range(L):

578d2 = ［］

579for j in range(L):

580d1 = ［1 if (j<len(seq) and i<len(seq) and key［x］ == seq［i］ and key［x］ == seq［j］) else 0 

for x in range(N)］  

581d2.append(d1) #d1是维度为(1,20)的One-Hot特征向量

582tensor.append(d2) #d2是维度为(L,20)的特征矩阵

583return np.array(tensor) #tensor是维度为(L,L,20)的特征矩阵




执行程序段P5.51，将所有蛋白质的残基序列编码为L×L×N的特征矩阵。



P5.51#将所有蛋白质的残基序列转换为L×L×N的特征矩阵

584inputs_aa = np.array(［wider(seq) for seq in seqs］)  

585print('所有蛋白质残基序列构成的特征矩阵维度为:{0}'.format(inputs_aa.shape))




结果显示所有蛋白质残基序列构成的特征矩阵维度为(890，64，64，20)。
5.16构建3D评分矩阵
将单个蛋白质的评分矩阵PSSM转换为L×L×N的3D特征矩阵，L表示序列的长度，N表示氨基酸种类，固定为20，如图5.26所示。


图5.263D评分矩阵的定义方法


实践中，在图5.26的基础上，根据两个残基间的相对关系，增加两个平衡值，扩展评分向量，深度方向的维度调整为22。执行程序段P5.52，将单个蛋白质的评分矩阵PSSM转换为L×L×(N+2)的特征矩阵。



P5.52#将单个蛋白质的评分矩阵PSSM转换为L×L×(N+2)的特征矩阵

586def wider_pssm(pssm, seq, L=64, N=20):

""" 将PSSM转换为 L×L×N 矩阵，而不是 L×N """

587key_alpha = "ACDEFGHIKLMNPQRSTVWY"

588tensor = ［］

589for i in range(L):

590d2 = ［］

591for j in range(L):

#残基序列范围内，添加一个维度为(1，20)的评分向量

592if j<len(seq) and i<len(seq): 

593d1 = ［(aa［i］+aa［j］)/2 for aa in pssm］  #取i、j列的平均值

594else:  #范围外向量(1,20)的元素取值均为0

595d1 = ［0 for i in range(N)］

#残基序列范围内，扩展评分向量，取值 pssm［i］*pssm［j］

596if j<len(seq) and i<len(seq):

597d1.append(pssm［key_alpha.index(seq［i］)］［i］ *

598  pssm［key_alpha.index(seq［j］)］［j］)   # i、j列评分乘积

599else:# 残基序列范围外，取值 0

600d1.append(0)

#根据残基序列i、j位置的相对关系扩展评分向量

601if j<len(seq) and i<len(seq):

602d1.append(1-abs(i-j)/L)

603else:

604d1.append(0)

605d2.append(d1)#d1是一个(1,22)的向量

606tensor.append(d2)#d2是一个(L，22)的矩阵

607return np.array(tensor)#tensor是(L, L，22)的矩阵




执行程序段P5.53，将所有蛋白质的评分矩阵转换为L×L×(N+2)的特征矩阵。



P5.53#将所有蛋白质的评分矩阵转换为L×L×(N+2)的特征矩阵

608inputs_pssm = np.array(［wider_pssm(pssms［i］, seqs［i］) for i in range(len(pssms))］)

609print('所有蛋白质构成的评分特征矩阵维度为:{0}'.format(inputs_pssm.shape))




结果显示所有蛋白质构成的评分特征矩阵维度为(890，64，64，22)。
执行程序段P5.54，将残基序列矩阵和评分矩阵堆叠在一起，构建训练集特征矩阵（m，L，L，42）。



P5.54#构建训练集特征矩阵(m，L，L，42)

610inputs = np.concatenate((inputs_aa, inputs_pssm), axis=3)

611print('训练集特征矩阵维度为:{0}'.format(inputs.shape))

#删除不再使用的数据集，释放内存

612del inputs_pssm 

613del inputs_aa




结果显示训练集特征矩阵维度为(890，64，64，42)。
5.17定义距离标签的3D矩阵
距离是表征残基对之间关系的一种度量方法，5.6节计算得到的残基对之间的距离是一个维度为（L，L）的2D矩阵，距离是非负实数。本节根据距离的分段范围，将距离的2D矩阵转换为3D矩阵。
首先对单个蛋白质距离矩阵的维度进行扩展，如果残基序列的实际长度小于L，则空白处填充-1，如图5.27所示。


图5.27残基距离矩阵扩展


执行程序段P5.55，将单个蛋白质的距离矩阵的维度扩展为（L，L）。



P5.55#将单个蛋白质的距离矩阵的维度扩展为(L，L)

614def embedding_matrix(matrix, L=64):

#列方向扩展

615for i in range(len(matrix)):

616if len(matrix［i］)<L:

617matrix［i］.extend(［-1 for i in range(L-len(matrix［i］))］)

#行方向扩展

618while len(matrix)<L:

619matrix.append(［-1 for x in range(L)］)

620return np.array(matrix)




执行程序段P5.56，将所有蛋白质的距离矩阵拓展，构成训练集的标签矩阵（m，L，L）。



P5.56#将所有蛋白质的距离矩阵拓展，构成训练集的标签矩阵(m，L，L)

621dists = np.array(［embedding_matrix(matrix) for matrix in dists］)

622print('训练集的距离标签矩阵维度为:{0}'.format(dists.shape))




运行结果显示训练集的距离标签矩阵维度为(890，64，64)。
以［<-0.5，≤500，≤750，≤1000，≤1400，≤1700，>1700］作为距离分段的标准，将距离划分为7个类别，对如图5.27所示的距离矩阵进行OneHot编码，从而将2D距离矩阵拓展为3D标签矩阵。如图5.28所示，以距离为801.23为例，其对应的OneHot向量为［0，0，0，1，0，0，0］，序列范围以外的距离值-1，均转换为向量［1，0，0，0，0，0，0］。


图5.28构建距离标签的3D矩阵


定义程序段P5.57，将残基间的距离编码为7个类别。



P5.57#将残基间的距离编码为7个类别

623def treshold(matrix, cuts=［-0.5, 500, 750, 1000, 1400, 1700］, L=64):

#将(L，L)特征矩阵转换为(L，L，7)特征矩阵 

624trash = (np.array(matrix)<cuts［0］).astype(np.int)#<-0.5 段0

625first = (np.array(matrix)<=cuts［1］).astype(np.int)-trash#0.5~500 段1

626sec = (np.array(matrix)<=cuts［2］).astype(np.int)-trash-first#500~750 段2

627third = (np.array(matrix)<=cuts［3］).astype(np.int)-trash-first-sec #7500~
#1000 段3

628fourth = (np.array(matrix)<=cuts［4］).astype(np.int)-trash-first-sec-third#1000~1400 段4

629fifth = (np.array(matrix)<=cuts［5］).astype(np.int)-trash-first-sec-third-fourth #1400~1700段5

630sixth = np.array(matrix)>cuts［5］#>1700 段6

631return np.concatenate((trash.reshape(L,L,1),

632first.reshape(L,L,1),

633sec.reshape(L,L,1),

634third.reshape(L,L,1),

635fourth.reshape(L,L,1),

636fifth.reshape(L,L,1),

637sixth.reshape(L,L,1)),axis=2) #矩阵堆叠为(L，L，7)




执行程序段P5.58，将所有蛋白质的距离矩阵编码为（m，L，L，7）的标签矩阵。



P5.58#将所有蛋白质的距离矩阵编码为(m，L，L，7)的标签矩阵

638outputs = np.array(［treshold(d) for d in dists］)

639print('距离标签矩阵维度为:{0}'.format(outputs.shape))

640del dists#释放内存




运行结果显示距离标签矩阵维度为(890，64，64，7)。
5.18距离模型参数设定与训练
执行程序段P5.59，完成模型参数定义与编译。



P5.59#模型参数定义与编译

#设定优化算法

641adam = keras.optimizers.Adam(lr=0.001, beta_1=0.9, beta_2=0.999, epsilon=1e-8, 

decay=0.0, amsgrad=True)

#创建模型

642model = resnet_v2(input_shape=(64,64, 42),  depth=16,  num_classes=7)

#［1e-07, 0.45, 1.65, 1.75, 0.73, 0.77, 0.145］为输出标签的权重参数，用于损失值的计算

643model.compile(optimizer=adam, loss=weighted_categorical_crossentropy(

np.array(［1e-07, 0.45, 1.65, 1.75, 0.73, 0.77, 0.145］)), metrics=［"accuracy"］)

644model.summary()




运行结果显示模型包含16个残差块，50个卷积层，需要学习训练的参数为913927个。
执行程序段P5.60，划分训练集与验证集，蛋白质总数为890个，训练集设定为600个，验证集设定为290个。



P5.60#划分训练集与验证集

645split = 600

646x_train, x_val = inputs［:split］, inputs［split:］

647y_train, y_val = outputs［:split］, outputs［split:］




执行程序段P5.61，设定模型训练参数，开启训练过程，设定epochs为50代，early_stopping参数的作用是，如果连续5次在验证集的损失值不下降，则提前停止训练。



P5.61#设定模型训练参数，开启训练过程

#连续5次验证集损失值不下降，停止训练

648early_stopping = EarlyStopping(monitor='val_loss', patience=5) 

649his = model.fit(x_train, y_train, epochs=50, batch_size=4, verbose=1, shuffle=True, 

validation_data=(x_val, y_val), callbacks=［early_stopping］)




训练主机CPU配置为Intel CoreTM i76700 CPU@3.40Hz，内存配置为16GB，训练过程提前终止于第21个epoch，用时1h21min51s。读者可根据自己的计算能力，修正数据集规模。数据集越大，取得的效果越好。
执行程序段P5.62，保存模型，绘制准确率曲线如图5.29所示，损失函数曲线如图5.30所示。



P5.62#保存模型，绘制准确率曲线和损失函数曲线

650model.save("model_under_64.h5")

651plt.figure(figsize=(8,4))

652x=range(1, len(his.history［'accuracy'］)+1)

653plt.plot(x,his.history［"accuracy"］)

654plt.plot(x,his.history［"val_accuracy"］)

655plt.legend(［"accuracy", "val_accuracy"］, loc="lower right")

656plt.xlabel('Epoch')

657plt.xticks(x)

658plt.show()

659plt.figure(figsize=(8,4))

660plt.plot(x,his.history［"loss"］)

661plt.plot(x,his.history［"val_loss"］)

662plt.legend(［"loss", "val_loss"］, loc="upper right")

663plt.xlabel('Epoch')

664plt.xticks(x)

665plt.show()






图5.29训练集与验证集的准确率曲线对比


如图5.29所示，验证集的准确率波动较大，应该是数据集样本的个体间的差异造成的，随着迭代次数的增加，准确率曲线在训练集与验证集均有下降趋势，再次证明数据集的样本不够充分。准确率在训练集与验证集的趋势保持一致，证明模型的方差控制在合理范围内。
如图5.30所示，训练集与验证集的损失函数下降趋势保持了高度一致，证明模型的方差小，泛化能力好。


图5.30训练集与验证集的损失函数曲线对比


5.19距离模型预测与评价
执行程序段P5.63，对给定范围的蛋白质做出距离预测。



P5.63#对给定范围的蛋白质做出距离预测

666start, num = 0, 5  #预测范围： start表示起点，num表示数量

667sample_pred = model.predict(［inputs［start : start+num］］)




执行程序段P5.64，根据预测结果与真实标签中的最大值索引重新编码，将预测结果与真实标签中的trash类别的索引值0设为7，从而用1、2、3、4、5、6、7表示7种距离值。



P5.64#预测结果与真实标签的再次编码，用1、2、3、4、5、6、7表示7种距离值

668preds_matrix = np.argmax(sample_pred, axis=3)  

669preds_matrix［preds_matrix==0］ = 7 #将 trash 类别设置为 7，表示无穷远

670outs_matrix = np.argmax(outputs［start : start+num］, axis=3)

671outs_matrix［outs_matrix==0］ = 7 #将标签中的trash类别设为7，表示无穷远




执行程序段P5.65，根据准确率选择最好的5个预测结果。



P5.65#根据准确率选择最好的5个预测结果

672results = ［np.sum(np.equal(pred［:len(seqs［start+j］),:len(seqs［start+j］)］,

outs_matrix［j,:len(seqs［start+j］),:len(seqs［start+j］),］),

axis=(0,1))/len(seqs［start+j］)**2 for j,pred in enumerate(preds_matrix)］

673best_score = max(results)

674print("准确率最高值为: ", best_score)

675sorted_scores = ［acc for acc in sorted(results, key=lambda x: x, reverse=True)］

676print("准确率最好的5个值为: ", sorted_scores［:5］)

677print("准确率最好的蛋白质序号为： ", ［results.index(x) for x in sorted_scores［:5］］)

678best_score_index = results.index(best_score)

679print("准确率最高的蛋白质序号为： ", best_score_index)  




运行结果如下。



准确率最高值为:  0.7552083333333334

准确率最好的5个值为:  ［0.7552083333333334, 0.48404542996214167, 0.4822485207100592, 

0.4359438660027162, 0.4351961950059453］

准确率最好的蛋白质序号为：  ［1, 2, 4, 3, 0］

准确率最高的蛋白质序号为：  1




执行程序段P5.66，分别根据预测值与标签值绘制蛋白质距离图，真实距离分布如图5.31所示，预测距离分布如图5.32所示。



P5.66#绘制蛋白质距离图(真实图和预测图对比)

680best_score_index=3   #显示序号为3的蛋白质

681plt.title('Ground Truth of ' + names［best_score_index］)#真实图  

682norm = plt.Normalize(1, 7)

683plt.imshow(outs_matrix［best_score_index, :len(seqs［start+best_score_index］),

:len(seqs［start+best_score_index］)］, cmap='viridis_r', interpolation='nearest' , norm=norm)

684plt.colorbar()

685plt.show()

686plt.title("Prediction by model of " + names［best_score_index］)#预测图

687plt.imshow(preds_matrix［best_score_index, :len(seqs［start+best_score_index］),

:len(seqs［start+best_score_index］)］, cmap='viridis_r', interpolation='nearest', norm=norm)

688plt.colorbar()

689plt.show()







图5.31蛋白质1GPT_1_A的距离真实图




图5.32蛋白质1GPT_1_A的距离预测图



图5.31和图5.32中亮色区域（黄色区域）表示残基对的距离较近，深色区域（蓝色区域）表示距离较远。虽然预测图形比较模糊，与真实图形有较大误差，但是仍然可以看到轮廓骨架还是具有较高的相似性。
执行程序段P5.67，通过混淆矩阵，进一步观察预测结果与真实值的偏差，得到混淆矩阵如图5.33所示。



P5.67#用混淆矩阵观察预测值与真实值的偏差

690from sklearn.metrics import confusion_matrix

691from sklearn.utils.multiclass import unique_labels

692preds_crop = np.concatenate( ［pred［:len(seqs［start+j］), :len(seqs［start+j］)］.flatten() 

for j,pred in enumerate(preds_matrix)］ )

693outs_crop = np.concatenate( ［outs_matrix［j, :len(seqs［start+j］), :len(seqs［start+j］)］.flatten() 

for j,pred in enumerate(preds_matrix)］ )

694matrix = cm = confusion_matrix(outs_crop, preds_crop)

695classes = ［i+1 for i in range(7)］

696title = "Comfusion matrix"

697cmap = "YlOrRd"

698normalize = True

699if normalize:

700cm = cm.astype('float') / cm.sum(axis=1)［:, np.newaxis］

701print("归一化的混淆矩阵")

702else:

703print('非归一化的混淆矩阵')

704fig, ax = plt.subplots()

705im = ax.imshow(cm, interpolation='nearest', cmap=cmap)

706ax.figure.colorbar(im, ax=ax)

#设置刻度，表示距离类型

707ax.set(xticks=np.arange(cm.shape［1］),

708yticks=np.arange(cm.shape［0］),

709xticklabels=classes, yticklabels=classes,

710title=title,

711ylabel='True label',

712xlabel='Predicted label')

713plt.setp(ax.get_xticklabels(), rotation=45, ha="right",

714 rotation_mode="anchor")

#设置显示格式

715fmt = '.2f' if normalize else 'd'

716thresh = cm.max() / 2.

717for i in range(cm.shape［0］):

718for j in range(cm.shape［1］):

719ax.text(j, i, format(cm［i, j］, fmt),

720ha="center", va="center",

721color="white" if cm［i, j］ > thresh else "black")

722fig.tight_layout()

723print("总均方误差: ", np.linalg.norm(outs_crop-preds_crop))

724print("平均均方误差: ", np.linalg.norm(outs_crop-preds_crop)/len(preds_matrix))






图5.33预测值与真实值的归一化混淆矩阵


图5.33的行方向为真实值，列方向为预测值。不难看出，模型对1、2、3、5、6这五种距离的准确率都超过了50%。对距离1、2的预测效果最好。对距离4的预测效果较差。
以准确率最高的距离标签1为例，有10%的情况预测为标签2，2%的情况预测为标签3。程序段P5.67显示模型对五个蛋白质结构预测的总均方误差为76.27，平均均方误差为15.25。
修改第698行程序代码，设置normalize=False，可以得到用数量表示偏差的混淆矩阵。
目前，AlphaFold还没有做到100%开源，按照DeepMind团队给出的解释，其采用的特征工程方法严重依赖Google的内部计算设施以及一些工具方法。在GitHub上，DeepMind团队只共享了AlphaFold在距离和角度预测方面的深度学习模型，没有开源全部细节，特别是特征生成和蛋白质3D结构的生成。
二面角预测模型和距离预测模型是AlphaFold进行蛋白质结构预测中的两个关键部分，本章基于Eric Alcaide在GitHub上开源的MinFold项目做了优化设计，基于CASP7的蛋白质文本数据集完成了二面角和距离的建模过程与训练过程，取得了较好的训练效果。
小结
本章以蛋白质结构预测的基本理论与方法为基础，梳理了AlphaFold预测蛋白质结构的算法逻辑，基于CASP7数据集构建了面向距离预测与面向角度预测的残基序列的特征矩阵，自定义实现了ResNet_V2一维卷积和二维卷积网络模型，在此基础上实践了蛋白质结构预测中关键的两个环节，即二面角模型和残基距离模型的构建、训练、预测与评估。
习题
一、 判断题
1. X射线晶体衍射、核磁共振和冷冻电镜等实验方法可以确定蛋白质的结构。
2. AlphaFold（又称A7D系统）是一套蛋白质结构预测系统。
3. 蛋白质的基本构成单位为氨基酸，常见氨基酸有20种。
4. 组成蛋白质的20种常见氨基酸中除脯氨酸外，均为α氨基酸。
5. 肽链的一端含有一个游离的α氨基，称为氨基端或N端； 在另一端含有一个游离的α羧基，称为羧基端或C端。
6. 氨基酸间脱水后生成的共价键称为肽键，其中的氨基酸单位称为氨基酸残基。
7. 每一种天然蛋白质都有自己特有的空间结构，这种空间结构称为蛋白质的（天然）构象。
8. 蛋白质的一级结构指蛋白质多肽链中氨基酸残基的排列顺序。
9. 一级结构(残基序列)中含有形成高级结构全部必需的信息，一级结构决定高级结构及其功能。
10. 蛋白质的二级结构是多肽链中各个肽段借助氢键形成有规则的构象。主要包括α螺旋、β折叠、β转角等。
11. 蛋白质的三级结构是多肽链借助各种非共价键的作用力，通过弯曲、折叠，形成具有一定走向的紧密球状构象。
12. 蛋白质的四级结构是寡聚蛋白中各亚基之间在空间上的相互关系和结合方式。
13. 第i个残基与第i+1个残基之间的距离，定义为Ciα与Ci+1α两个原子之间的距离。
14. 围绕Cα-N键轴旋转产生的角度称为Phi（Φ），围绕Cα-C键轴旋转产生的角度称为Psi（Ψ）。
15. 拉氏构象图表明二面角（Φ、Ψ）的变化范围为-180°～180°，但是Φ、Ψ不能任意取值。
16. 因为氨基酸有20种，所以单个残基的OneHot编码的维度为（20，1）。
17. 二面角预测模型采用了ResNet一维卷积结构，距离预测模型采用了ResNet二维卷积结构。
二、 编程题
根据本章案例采用的ResNet_V2网络结构，自定义ResNet_V2模型，基于CIFAR10数据集实现图像分类识别。
程序的主体框架已经在习题5文件夹中的程序文档cifar10_begin.ipynb中给出，请根据程序中的逻辑提示与运行结果提示，补充完成整个程序，并进行测试。
为了增强对数据集结构的理解，本编程项目采用CIFAR10作者发布的原始数据集，不采用Keras等框架提供的在线版本。
CIFAR10数据集和CIFAR100数据集是由Alex Krizhevsky、Vinod Nair和Geoffrey Hinton从8000万个微型图像中收集标记的图像子集。数据集下载地址： 
https://www.cs.toronto.edu/~kriz/cifar10python.tar.gz
CIFAR10数据集包含10个类别的60000个32×32彩色图像，每个类别包含6000个图像。50000个图像用于训练模型，10000个图像用于模型测试。
数据集分为5个训练批次和1个测试批次，每个批次包括10000个图像。测试批次包含每个类别的1000个随机选择的图像。训练批次按随机顺序包含其余图像，某些训练批次并不保证10类图像一样多，但是5个批次加起来，每个类别都是5000个图像。
图5.34显示了数据集包含的10个类别，以及每个类别中随机选出的10个图像。


图5.3410种类别及随机选取的10个图像


数据集解压目录如图5.35所示，包含文件data_batch_1、data_batch_2、data_batch_3、data_batch_4、data_batch_5和test_batch。每个文件都是由cPickle生成的Python“pickled”对象。


图5.35CIFAR10数据集目录


可以采用自定义函数unpickle(file)读取数据集文件。



def unpickle(file):  #读取数据集文件

import pickle

with open(file, 'rb') as fo:

dict = pickle.load(fo, encoding='bytes')

return dict





函数unpickle(file)将读取的数据返回为一个字典对象。字典包含data和labels两个元素。
data： 是10000×3072的numpy数组矩阵，数据类型为uint8。矩阵的每一行都存储一个32×32彩色图像。前1024个像素值为红色通道值，中间1024个像素值为绿色通道值，最后1024个像素值为蓝色通道值。
labels： 是10000个数字构成的列表，数字范围为0~9。索引i处的数字表示矩阵数据中第i个图像的标签。
数据集包含的文件batchs.meta返回的也是Python字典对象。
label_names： 是由10个元素组成的列表，为数字标签提供有意义的名称。例如，label_names［0］=="airplane"，label_names［1］=="automobile"等。