基于长短期记忆神经网络的蛋白质二级结构预测

问题描述

首先我们应该对对蛋白质二级结构有一定的了解，知道蛋白质二级结构预测指的是什么。蛋白质二级结构预测其实就是在蛋白质一级结构已知的情况下去预测该位置对应的二级结构。

数据来源

业界常用的数据集包括CB513数据集、CULLPDB数据集、CASP10数据集。

数据的处理

其实我个人觉得这个预测最为关键的一步就是数据的处理，数据处理的方法很大部分决定了预测的准确性。下面以CB513数据集为例：
首先到上面链接下载下载CB513数据集：
然后打开readme.txt文件：

可用看到该数据集目前是以numpy格式存储，可 以先将其重塑为（514 个蛋白×700 个氨基酸×57 特征）。也就是说它的数据结构是有 N 个矩阵，每个矩阵有 700 行，每行有 57 列，57 个列分别是：[0,22]列：氨基酸残基，顺序 为’A’，‘C’，‘E’，‘D’，‘G’，‘F’，‘I’，‘H’，‘K’ ，‘M’，‘L’，‘N’，‘Q’，‘P’，‘S’，‘R’，‘T’，‘W’，‘V’， ‘Y’，‘X’，‘NoSeq’ 、[22,31]列：二级结构标签，其序列为’L’，‘B’，‘E’，‘G’，‘I’，‘H’，‘S’，‘T’， ‘NoSeq’ 、[31,33]：N-和 C-终端、[33,35]：相对和绝对溶剂可及性，仅用于训练（绝对可 达性定为 15；相对可及性通过蛋白质中最大的可达性值归一化，阈值为 0.15；原始溶剂可 及性由 DSSP 计算）、[35,57]：22 位 PSSM 谱编码。[22,31]就是我们数据集的标签，他对 应的是’L’，‘B’，‘E’，‘G’，‘I’，‘H’，‘S’，‘T’，‘NoSeq’ 九种二级结构所对应的独热编码。[0,22] 就 是 21 种 氨 基 酸 和 一 种 未 确 定 氨 基 酸 对 应 的 正 交 编 码 ， 比 如 [1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]表示 A（丙氨酸）。[35,57]22 位 PSSM 谱编码， 也是作为特征。也就是说特征包括 22 位正交编码和 22 位 PSSM 谱编码，标签是 9 位的二 级结构独热编码。可以将数据集划分80%用来训练，20%用来预测。
对于 CB513 数据集，由于前人已经将数据集处理成 N 蛋白乘上 k 特征的矩阵，所以我 只需要将数据 reshape 成（矩阵的长度链的长度特征属性的个数）的矩阵。在预测时我采用一种蛋白质二级结构预测常用的滑动窗口技术（滑动窗口的大小需要经过反复调整以达到最佳预测效果为止）。例如：

// 读取数据
input_data = np.load('D:/cb513+profile_split1.npy.gz','r')
//将数据reshape成所需要的格式
input_data = np.reshape(input_data, (514, 700, 57))
//[0-22]位是特征，也就是输入
x = input_data[:,:,0:22]
//[22-31]位是标签
y = input_data[:,:,22:31]
dataX = np.reshape(x, (len(x), 700, 22))

长短期记忆神经网络

本文采用的是长短期记忆神经网络来进行预测，构造该网络的方法有很多，我才此用的是keras来构造神经网络。

// 存放输入的数组
X = []
i = 0
j = 0
k = 0
//滑动窗口大小
window = 17
while i<((700-window+1)*514):
	if((j%(700-window+1)==0) and (j!=0)):
		j=0
		k+=1
	x = dataX[k,j:j+window,:]
	X.append(x)
	i+=1
	j+=1
X = np.reshape(X, (len(X), window, 22))
# print(len(X))
i = 0
j = 0
k = 0
//标签存放数组
Y = []
while i<((700-window+1)*514):
	if((j%(700-window+1)==0) and (j!=0)):
		j=0
		k+=1
	Y.append(y[k,j+8,:])
	i+=1
	j+=1
y = np.reshape(Y, (len(Y), 9))
print(len(y))
# # # 数据的划分
Xtrain, Xtest, ytrain, ytest = train_test_split(X, y, test_size=0.2)
模型的构建
model = Sequential()
model.add(Bidirectional(LSTM(32, dropout=0.5, recurrent_dropout=0.2,return_sequences = True),input_shape=(X.shape[1], X.shape[2])))
model.add(LSTM(32,dropout=0.5))
model.add(Dense(32))
model.add(Dropout(0.3))
model.add(Dense(16))
model.add(Dropout(0.3))
model.add(Dense(y.shape[1],use_bias=True))
model.add(Activation("sigmoid"))
model.compile(loss="binary_crossentropy", optimizer="rmsprop",metrics=["accuracy"])
网络训练
model.fit(Xtrain, ytrain, batch_size=32, epochs=80,validation_data=(Xtest, ytest))
# 模型验证
score, acc = model.evaluate(Xtest, ytest, batch_size=32)
print("\nTest score: %.3f, accuracy: %.3f" % (score, acc))

训练结果

可以看到预测的准确率非常的高，这是很不合理的，因为目前世界上最高的准确率也不过是80%出头。通过查看数据集的数据本身，可以发现数据集中有很多是属于 Noseq 的结构，而 NoSeq 是指 CB513 数据集把每一个蛋白质链都补成了长度为 700 个氨基酸，当一个链的氨基酸个数不足 700 时，数据集的创始者通过补 NoSeq 结构填补，也就是说 NoSeq 结构的预测准确率对整体准确率有影响，但是实际却是没有意义的，所以必须把 NoSeq 结构从 CB513 数据集中去掉，总体的准确率不应该包括 NoSeq 结构的准确率。 不包括 Noseq 结构准确率大概为 76.7%。

// 调整去掉Noseq结构，重新计算准确率![在这里插入图片描述](https://img-blog.csdn.net/2018101717233993?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM0MjE4MjU1/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70)
acc1 = 0
acc2 = 0
acc3 = 0
acc4 = 0
sum1 = 0
sum2 = 0
sum3 = 0
sum4 = 0 

for i in range(len(Xtest)):
    # idx = np.random.randint(len(X))
    xtest = Xtest[i,:,:].reshape(1,17,22)
    # print(xtest)
    ylabel = ytest[i,:]
    # print(ylabel)
    ypre = model.predict(xtest)[0,:]
    max = 0
    i = -1
    j = 0
    for y4 in ypre:
        if(y4>max):
            max = y4
            i+=1
            j = i
        else:
            i+=1
    y1 = [0,0,0,0,0,0,0,0,0]
    y1[j] = 1
    if((y1==[0,1,0,0,0,0,0,0,0]) or (y1 == [0,0,1,0,0,0,0,0,0])):
        y1 = [1,0,0]
    else:
        if((y1==[0,0,0,1,0,0,0,0,0]) or (y1==[0,0,0,0,0,1,0,0,0])):
            y1 = [0,1,0]
        else:
            if(y1==[0,0,0,0,0,0,0,0,1]):
                y1 = [0,0,0]
            else:
                y1 = [0,0,1]
    if(y1!=[0,0,0]):
        print("预测：")
        print(y1)
    y2 = []
    for y3 in ylabel:
    	y2.append(int(y3))
    if((y2==[0,1,0,0,0,0,0,0,0]) or (y2 == [0,0,1,0,0,0,0,0,0])):
        y2 = [1,0,0]
        sum2+=1
        # QE
    else:
        if((y2==[0,0,0,1,0,0,0,0,0]) or (y2==[0,0,0,0,0,1,0,0,0])):
            y2 = [0,1,0]
            sum3+=1
            # QH
        else:
            if(y2==[0,0,0,0,0,0,0,0,1]):
                y2 = [0,0,0]
                sum1+=1
                # NoSeq
            else:
                y2 = [0,0,1]
                sum4+=1
                # QC 
    if(y2!=[0,0,0]):
        print("标签：")
        print(y2)
    if((y1==y2==[1,0,0])):
        acc2+=1
    if((y1==y2==[0,1,0])):
        acc3+=1
    if((y1==y2==[0,0,1])):
        acc4+=1


print((acc2+acc3+acc4)/(len(Xtest)-sum1))
print(acc2/sum2)
print(acc3/sum3)
print(acc4/sum4)

预测结果：