Paper再现：MD+AI自动编码机探测蛋白变构(四):DIO的生成和聚类

四、 DIO的生成和聚类

这一部分比较简单，就是一个sklearn的聚类分析而已，所以没必要说什么，直接上代码。这里一个主要问题是，sklearn里面的聚类方法，是单核运行的，于是运行的速度很慢，这个蛋白只有500左右的残基，已经需要3个小时。所以，应该找其他方法代替。

4.1导入相关模块

#加载模块
import pandas as pd
import numpy as np
import time

import torch 
import torch.nn as nn
import torch.nn.functional as F
from torch import optim

#将torch的数据类型设置为64位，默认是32位的
torch.set_default_dtype(torch.float32)

device = torch.device('cuda:0' if torch.cuda.is_available() else 'cpu')
device

4.2 导入之前的模型

即：

class Aotuencoder_for_md(nn.Module):
    
    def __init__(self):
        super(Aotuencoder_for_md, self).__init__()
        
        self.encoder_layer1 = nn.Linear(1000, 800)
        self.encoder_layer2 = nn.Linear(800, 500)
        self.encoder_layer3 = nn.Linear(500, 200)
        
        self.decoder_layer1 = nn.Linear(200, 500)
        self.decoder_layer2 = nn.Linear(500, 800)
        self.decoder_layer3 = nn.Linear(800, 1000)
    
    def forward(self, x):
        y = F.relu(self.encoder_layer1(x))
        y = F.relu(self.encoder_layer2(y))
        y = F.relu(self.encoder_layer3(y))
        
        y = F.relu(self.decoder_layer1(y))
        y = F.relu(self.decoder_layer2(y))
        y = F.relu(self.decoder_layer3(y))
        return y

4.3 加载模型并训练

#加载模型 opo1_model，如果是高cuda版本训练出来的模型是不能在低cuda版本加载
opo_model= torch.load('./model.pth')
opo_model = opo_model.to(device)
opo_model

记录函数：

def run_log(t, file='run_lof.txt', mode = 'a+'):
    f = open (file, mode)
    time_ = time.strftime("%Y-%m-%d %H:%M:%S", time.localtime()) 
    f.write(t)
    f.write(time_)
    f.write('\n')

加载数据

#opo1-holo1 DIOs
#加载数据集
#请注意很吃很迟内存和显存！！！建议批量加载，批量运行哈
#保存成csv格式很浪费时间和空间，建议保存成数组哈
holo_data = pd.read_csv('dataset/dataset/holo_csv.csv',index_col=0)
index = holo_data.index
run_log('index:\n', file='opo_holo_index.txt', mode='w')
run_log(' '.join(index.values), file='opo_holo_index.txt', mode='a+')

DIO生成：

def get_DIOs(model, data_path, dios_path):
    data = pd.read_csv(data_path,index_col=0)
    index = data.index.values
    data = data.values
    mask = np.all(np.isnan(data), axis=1) | np.all(data == 0, axis=1)
    data = data[~mask]
    index = index[~mask]
    data = torch.from_numpy(data)
    data = data.to(device)
    dios = data-model(data)
    dios = dios.cpu().detach().numpy()
    np.save(dios_path, dios)
    return dios, index

opo_holo, index = get_DIOs(model=opo_model, data_path='dataset/dataset/holo_csv.csv', 
                      dios_path='dataset/dataset/opo_holo_DIOs.npy')
opo_holo, index

4.4 DIO聚类分析

下面进行聚类，根据文章的说明，分为7类：

from sklearn.cluster import AgglomerativeClustering

def DIOs_cluster(X):
    clustering = AgglomerativeClustering(n_clusters=3, 
                                         affinity='euclidean', #不能在X中包含0向量的时候使用，。。。opo1_holo2中有！！！！
                                         linkage='single', compute_distances=True).fit(X) #单核运行很占时间！！
    # cosine was provided as affinity. Ward can only work with euclidean distances.
    # compute_distances=True 这个参数只有在0.24版本才能用的哈，
    #请升级至最高版本，没有这个参数，是没办法实现可视化的。
    #由于使用该参数，运行时间更长，内存占用非常厉害！！！
    return clustering

输出每一个标签的分类label,运行时间好久。

opo_holo_agg = DIOs_cluster(opo_holo)
opo_holo_lables = opo_holo_agg.labels_
opo_holo_lables

4.5 保存一下每个残基对的分类情况

data = pd.DataFrame(index = index, columns=['label'])
data['label'] = opo_holo_lables
data.to_csv('dataset/dataset/opo_holo_lables.csv')

结合之前的一、二、三部分，MD+AI分析变构网络的任务基本就完成了。剩下的就是可视化，这个太简单了，就不在这里介绍了，所以pass。

另外，由于蛋白是我们的研究对象，是保密的，所以原始数据不能分享，大家可以自己做一个MD然后试一下。其实，这个方法不是很稳定，一方面，MD要做的出差别，然后就是聚类分析的时候要好好的处理，还有就是模型这一块，这类简单的自动编码机其实很不适合这个任务的。总之，多尝试~

总的说，文献复制完毕。