others/Attribute Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ongoing-liabilities",
   "metadata": {},
   "outputs": [],
   "source": [
    "import tensorflow as tf\n",
    "\n",
    "from tensorflow.keras import Model\n",
    "from tensorflow.keras.layers import Dense, Input, concatenate, BatchNormalization\n",
    "from tensorflow.keras.losses import MeanSquaredError,MSE\n",
    "from tensorflow.keras.metrics import MeanAbsoluteError\n",
    "from tensorflow.keras.optimizers import Adam,RMSprop\n",
    "from tensorflow.keras.regularizers import l2\n",
    "from tensorflow.keras.models import load_model\n",
    "\n",
    "from spektral.data import MixedLoader,Dataset,DisjointLoader,Graph,BatchLoader\n",
    "from spektral.datasets.mnist import MNIST\n",
    "from spektral.layers import GCNConv,GlobalSumPool,ECCConv,CrystalConv,GlobalMaxPool\n",
    "from spektral.layers import GlobalAvgPool,GlobalAttnSumPool,GlobalAttentionPool\n",
    "from spektral.layers.ops import sp_matrix_to_sp_tensor\n",
    "from spektral.data import Graph\n",
    "from spektral.data import Dataset\n",
    "\n",
    "import rdkit.Chem as Chem\n",
    "from rdkit.Chem import AllChem\n",
    "import csv\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "violent-kuwait",
   "metadata": {},
   "outputs": [],
   "source": [
    "Index, X_smiles, M_adducts, CCS = [],[],[],[]\n",
    "f = csv.reader(open('data/Attribute importance data/data.csv','r', encoding='gbk',errors='ignore'))\n",
    "\n",
    "for i in f:\n",
    "    Index.append(i[0])\n",
    "    X_smiles.append(i[1])\n",
    "    M_adducts.append(i[2])\n",
    "    CCS.append(float(i[3]))\n",
    "\n",
    "GCN_smiles, GCN_adducts, GCN_Index, GCN_CCS = [],[],[],[]\n",
    "for i in range(len(X_smiles)):\n",
    "    try:\n",
    "        GCN_smiles.append(X_smiles[i])\n",
    "        GCN_adducts.append(M_adducts[i])\n",
    "        GCN_Index.append(Index[i])\n",
    "        GCN_CCS.append(CCS[i])\n",
    "    except:\n",
    "        ;\n",
    "\n",
    "smiles, ccs, adduct, Coordinate = [], [], [], []\n",
    "for i in range(len(GCN_Index)):\n",
    "    MOL = Chem.MolFromSmiles(GCN_smiles[i])\n",
    "    atoms = [atom.GetSymbol() for atom in MOL.GetAtoms()]\n",
    "    one_atom = []\n",
    "    f = csv.reader(open('data/Attribute importance data/Coordinate data/'+str(GCN_Index[i])+'.csv','r'))\n",
    "    files = [i for i in f]\n",
    "    for j in range(len(atoms)):\n",
    "        one_atom.append([float(iii) for iii in files[j+1][2:]])\n",
    "    Coordinate.append(one_atom)\n",
    "    smiles.append([GCN_smiles[i]])\n",
    "    ccs.append([GCN_CCS[i]])\n",
    "    adduct.append(GCN_adducts[i])\n",
    "###########################################################################################\n",
    "Max_Coor =  15.615155868453662 \n",
    "Min_Coor = -15.475082312818216\n",
    "for i in range(len(Coordinate)):\n",
    "    Coordinate[i] = (np.array((Coordinate[i])) - Min_Coor) / (Max_Coor - Min_Coor)\n",
    "\n",
    "Atom_radius = {'N' :71, 'Se':116, 'F':64, 'Co':111, 'O':63,'As':121,\n",
    "               'Br':114,'Cl':99,  'S':103,'C' :75, 'P':111, 'I':133,'H':32}\n",
    "Atom_radius_list = [Atom_radius[i] for i in Atom_radius]\n",
    "Max_radius, Min_radius = np.max(Atom_radius_list), np.min(Atom_radius_list)\n",
    "for i in Atom_radius:\n",
    "    Atom_radius[i] = (Atom_radius[i] - Min_radius) / (Max_radius-Min_radius)\n",
    "    \n",
    "Atom_mass = {'N':14.00674,'Se':78.96,'F':18.9984032,'Co':58.933195,'As':74.92160,\n",
    "             'O':15.9994,'Br':79.904,'Cl':35.453,'S':32.065,'C':12.0107,\n",
    "             'P':30.973762,'I':126.90447,'H':1.00794}\n",
    "Atom_mass_list = [Atom_mass[i] for i in Atom_mass]\n",
    "\n",
    "Max_mass, Min_mass = np.max(Atom_mass_list), np.min(Atom_mass_list)\n",
    "for i in Atom_mass:\n",
    "    Atom_mass[i] = (Atom_mass[i] - Min_mass) / (Max_mass-Min_mass)\n",
    "\n",
    "All_Atoms = ['As', 'Br', 'C', 'Cl', 'F', 'I', 'N', 'O', 'P', 'S', 'Se']\n",
    "###########################################################################################\n",
    "\n",
    "def convertToGraph(smi_lst):\n",
    "    adj,adj_norm, features, edge_features = [], [], [], []\n",
    "    maxNumAtoms = 50 \n",
    "    NodeNumFeatures, EdgeNumFeatures, INDEX = 0, 4, -1\n",
    "    for smi in smi_lst:\n",
    "        INDEX += 1\n",
    "        iMol = Chem.MolFromSmiles(smi[0]) \n",
    "        maxNumAtoms = iMol.GetNumAtoms() \n",
    "        iAdjTmp = Chem.rdmolops.GetAdjacencyMatrix(iMol)\n",
    "        \n",
    "        one_edge_features = edge_feature(iMol)\n",
    "        edge_features.append(one_edge_features)\n",
    "        \n",
    "        iFeature = np.zeros((maxNumAtoms, NodeNumFeatures))\n",
    "        iFeatureTmp = []\n",
    "        for atom in iMol.GetAtoms():\n",
    "            iFeatureTmp.append(atom_feature(atom,INDEX))\n",
    "        features.append(np.array(iFeatureTmp))\n",
    "        adj.append(iAdjTmp)\n",
    "            \n",
    "    features = np.asarray(features)\n",
    "    edge_features = np.asarray(edge_features)\n",
    "    return adj, features, edge_features\n",
    "\n",
    "def atom_feature(atom,INDEX):\n",
    "    return np.array(\n",
    "        one_of_k_encoding_unk(atom.GetSymbol() ,All_Atoms) +\n",
    "        one_of_k_encoding_unk(atom.GetDegree(), [0, 1, 2, 3, 4]) +\n",
    "        [Atom_radius[atom.GetSymbol()],Atom_mass[atom.GetSymbol()]] +\n",
    "        one_of_k_encoding_unk(atom.IsInRing(), [0, 1]) +\n",
    "        list(Coordinate[INDEX][atom.GetIdx()])\n",
    "    )\n",
    "\n",
    "def one_of_k_encoding_unk(x, allowable_set):\n",
    "    if x not in allowable_set:\n",
    "        x = allowable_set[-1]\n",
    "    return list(map(lambda s: x == s, allowable_set))\n",
    "\n",
    "def edge_feature(iMol):\n",
    "    iAdjTmp = Chem.rdmolops.GetAdjacencyMatrix(iMol)\n",
    "    Edge_feature = []\n",
    "    count = 0\n",
    "    for bond in iMol.GetBonds():\n",
    "        count += 1\n",
    "        bond_feature = np.array(\n",
    "            one_of_k_encoding_unk(bond.GetBondTypeAsDouble(),[1,1.5,2,3])\n",
    "        )\n",
    "        Edge_feature.append(bond_feature)\n",
    "        Edge_feature.append(bond_feature)\n",
    "    Edge_feature = np.array(Edge_feature)\n",
    "    Edge_feature = Edge_feature.astype(np.float)\n",
    "    return Edge_feature\n",
    "\n",
    "adj, features, edge_features = convertToGraph(smiles)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "similar-progressive",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MyDataset(n_graphs=10)\n",
      "['[M+H]+', '[M+Na]+', '[M-H]-']\n"
     ]
    }
   ],
   "source": [
    "class MyDataset(Dataset):\n",
    "    def __init__(self, features, adj, edge_features, ccs, **kwargs):\n",
    "        self.features = features\n",
    "        self.adj = adj\n",
    "        self.edge_features = edge_features\n",
    "        self.ccs = ccs\n",
    "        super().__init__(**kwargs)\n",
    "        \n",
    "    def read(self):\n",
    "        return [Graph(x = self.features[i], \n",
    "                      a = self.adj[i], \n",
    "                      e = self.edge_features[i],\n",
    "                      y = float(self.ccs[i][0])) for i in range(len(self.adj))]\n",
    "    \n",
    "DataSet = MyDataset(features, adj, edge_features, ccs)\n",
    "print(DataSet)\n",
    "adduct_SET = ['[M+H]+', '[M+Na]+', '[M-H]-']\n",
    "adduct_SET.sort()\n",
    "print(adduct_SET)\n",
    "\n",
    "dataset_te = DataSet\n",
    "adduct_te  = adduct"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "dangerous-premises",
   "metadata": {},
   "outputs": [],
   "source": [
    "import spektral\n",
    "import umap\n",
    "ECC_model = load_model('model/model.h5',\n",
    "                       custom_objects = {\"ECCConv\": spektral.layers.ECCConv,\n",
    "                                         \"GlobalSumPool\": spektral.layers.GlobalSumPool})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "sorted-province",
   "metadata": {},
   "outputs": [],
   "source": [
    "# 2 : Coordinates\n",
    "# 3 : Elemental symbols\n",
    "# 4 : Degree\n",
    "# 5 : Atomic radius\n",
    "# 6 : Atomic volume\n",
    "# 7 : Is on the ring\n",
    "\n",
    "def convertToGraph_2(smi_lst, SWITCH):\n",
    "    adj,adj_norm, features, edge_features = [], [], [], []\n",
    "    NodeNumFeatures, EdgeNumFeatures, INDEX = 0, 4, -1\n",
    "    for smi in smi_lst:\n",
    "        INDEX += 1\n",
    "        iMol = Chem.MolFromSmiles(smi[0]) # Convert Smiles strings to mol objects    \n",
    "        maxNumAtoms = iMol.GetNumAtoms()\n",
    "        iAdjTmp = Chem.rdmolops.GetAdjacencyMatrix(iMol) # Obtain the adjacency matrix of mol\n",
    "        # Characterization of structural chemical bonds\n",
    "        one_edge_features = edge_feature(iMol)\n",
    "        edge_features.append(one_edge_features)\n",
    "        # Constructing node feature data\n",
    "        iFeature = np.zeros((maxNumAtoms, NodeNumFeatures))\n",
    "        iFeatureTmp = []\n",
    "        for atom in iMol.GetAtoms():\n",
    "            if SWITCH ==   1:\n",
    "                iFeatureTmp.append(atom_feature_2(atom,INDEX))\n",
    "            elif SWITCH == 2:\n",
    "                iFeatureTmp.append(atom_feature_3(atom,INDEX))\n",
    "            elif SWITCH == 3:\n",
    "                iFeatureTmp.append(atom_feature_4(atom,INDEX))\n",
    "            elif SWITCH == 4:\n",
    "                iFeatureTmp.append(atom_feature_5(atom,INDEX))\n",
    "            elif SWITCH == 5:\n",
    "                iFeatureTmp.append(atom_feature_6(atom,INDEX))\n",
    "            elif SWITCH == 6:\n",
    "                iFeatureTmp.append(atom_feature_7(atom,INDEX))\n",
    "        features.append(np.array(iFeatureTmp))\n",
    "        adj.append(iAdjTmp)\n",
    "            \n",
    "    features = np.asarray(features)\n",
    "    edge_features = np.asarray(edge_features)\n",
    "    return adj, features, edge_features\n",
    "\n",
    "'''Coordinates'''\n",
    "def atom_feature_2(atom,INDEX):\n",
    "    return np.array(\n",
    "        one_of_k_encoding_unk(atom.GetSymbol() ,All_Atoms) +\n",
    "        one_of_k_encoding_unk(atom.GetDegree(), [0, 1, 2, 3, 4]) +\n",
    "        [Atom_radius[atom.GetSymbol()],Atom_mass[atom.GetSymbol()]] +\n",
    "        one_of_k_encoding_unk(atom.IsInRing(), [0, 1]) +\n",
    "        [0,0,0]\n",
    "    )\n",
    "\n",
    "'''Elemental symbols'''\n",
    "def atom_feature_3(atom,INDEX):\n",
    "    return np.array(\n",
    "        [0,0,0,0,0,0,0,0,0,0,0] + \n",
    "        one_of_k_encoding_unk(atom.GetDegree(), [0, 1, 2, 3, 4]) +\n",
    "        [Atom_radius[atom.GetSymbol()],Atom_mass[atom.GetSymbol()]] +\n",
    "        one_of_k_encoding_unk(atom.IsInRing(), [0, 1]) +\n",
    "        list(Coordinate[INDEX][atom.GetIdx()])\n",
    "    )\n",
    "'''Degree'''\n",
    "def atom_feature_4(atom,INDEX):\n",
    "    return np.array(\n",
    "        one_of_k_encoding_unk(atom.GetSymbol() ,All_Atoms) +\n",
    "        [0,0,0,0,0] + \n",
    "        [Atom_radius[atom.GetSymbol()],Atom_mass[atom.GetSymbol()]] +\n",
    "        one_of_k_encoding_unk(atom.IsInRing(), [0, 1]) +\n",
    "        list(Coordinate[INDEX][atom.GetIdx()])\n",
    "    )\n",
    "'''Atomic radius'''\n",
    "def atom_feature_5(atom,INDEX):\n",
    "    return np.array(\n",
    "        one_of_k_encoding_unk(atom.GetSymbol() ,All_Atoms) +\n",
    "        one_of_k_encoding_unk(atom.GetDegree(), [0, 1, 2, 3, 4]) +\n",
    "        [0,Atom_mass[atom.GetSymbol()]] +\n",
    "        one_of_k_encoding_unk(atom.IsInRing(), [0, 1]) +\n",
    "        list(Coordinate[INDEX][atom.GetIdx()])\n",
    "    )\n",
    "'''Atomic volume'''\n",
    "def atom_feature_6(atom,INDEX):\n",
    "    return np.array(\n",
    "        one_of_k_encoding_unk(atom.GetSymbol() ,All_Atoms) +\n",
    "        one_of_k_encoding_unk(atom.GetDegree(), [0, 1, 2, 3, 4]) +\n",
    "        [Atom_radius[atom.GetSymbol()],0] +\n",
    "        one_of_k_encoding_unk(atom.IsInRing(), [0, 1]) +\n",
    "        list(Coordinate[INDEX][atom.GetIdx()])\n",
    "    )\n",
    "'''Is on the ring'''\n",
    "def atom_feature_7(atom,INDEX):\n",
    "    return np.array(\n",
    "        one_of_k_encoding_unk(atom.GetSymbol() ,All_Atoms) +\n",
    "        one_of_k_encoding_unk(atom.GetDegree(), [0, 1, 2, 3, 4]) +\n",
    "        [Atom_radius[atom.GetSymbol()],Atom_mass[atom.GetSymbol()]] +\n",
    "        list([0,0]) +\n",
    "        list(Coordinate[INDEX][atom.GetIdx()])\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "compound-humidity",
   "metadata": {},
   "outputs": [],
   "source": [
    "adj_2, features_2, edge_features_2 = convertToGraph_2(smiles, 1)\n",
    "adj_3, features_3, edge_features_3 = convertToGraph_2(smiles, 2)\n",
    "adj_4, features_4, edge_features_4 = convertToGraph_2(smiles, 3)\n",
    "adj_5, features_5, edge_features_5 = convertToGraph_2(smiles, 4)\n",
    "adj_6, features_6, edge_features_6 = convertToGraph_2(smiles, 5)\n",
    "adj_7, features_7, edge_features_7 = convertToGraph_2(smiles, 6)\n",
    "\n",
    "DataSet_2 = MyDataset(features_2, adj_2, edge_features_2, ccs)\n",
    "DataSet_3 = MyDataset(features_3, adj_3, edge_features_3, ccs)\n",
    "DataSet_4 = MyDataset(features_4, adj_4, edge_features_4, ccs)\n",
    "DataSet_5 = MyDataset(features_5, adj_5, edge_features_5, ccs)\n",
    "DataSet_6 = MyDataset(features_6, adj_6, edge_features_6, ccs)\n",
    "DataSet_7 = MyDataset(features_7, adj_7, edge_features_7, ccs)\n",
    "\n",
    "dataset_te_2 = DataSet_2\n",
    "dataset_te_3 = DataSet_3\n",
    "dataset_te_4 = DataSet_4\n",
    "dataset_te_5 = DataSet_5\n",
    "dataset_te_6 = DataSet_6\n",
    "dataset_te_7 = DataSet_7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "soviet-atlanta",
   "metadata": {},
   "outputs": [],
   "source": [
    "def Fun(dataset_te,X):\n",
    "    loader_te = BatchLoader(dataset_te,batch_size=1,epochs=1,shuffle=False); loader_te_data = (); ltd_index = 0;\n",
    "    for i in loader_te.load():\n",
    "        adduct_one_hot = [one_of_k_encoding_unk(adduct_te[ltd_index+ltd_index_i],adduct_SET) for ltd_index_i in range(len(i[1]))]\n",
    "        adduct_one_hot = np.array(adduct_one_hot)\n",
    "        one_sample = ((adduct_one_hot,i[0][0],i[0][1],i[0][2]),i[1])\n",
    "        loader_te_data += (one_sample,)\n",
    "        ltd_index += len(i[1])\n",
    "    loader_te_data = (i for i in loader_te_data)\n",
    "    for batch in loader_te_data:\n",
    "        inputs, target = batch\n",
    "        predictions = ECC_model(inputs, training=False)\n",
    "        predictions = np.array(predictions[0])\n",
    "        X.append(predictions[0])\n",
    "    return np.array(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "taken-channels",
   "metadata": {},
   "outputs": [],
   "source": [
    "Target = np.array([i[0] for i in ccs])\n",
    "A = Fun(dataset_te,[])\n",
    "B = Fun(dataset_te_2,[])\n",
    "C = Fun(dataset_te_3,[])\n",
    "D = Fun(dataset_te_4,[])\n",
    "E = Fun(dataset_te_5,[])\n",
    "F = Fun(dataset_te_6,[])\n",
    "G = Fun(dataset_te_7,[])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "turned-bottom",
   "metadata": {},
   "outputs": [],
   "source": [
    "A2 = abs(Target-A)/Target*100.\n",
    "B2 = abs(Target-B)/Target*100.\n",
    "C2 = abs(Target-C)/Target*100.\n",
    "D2 = abs(Target-D)/Target*100.\n",
    "E2 = abs(Target-E)/Target*100.\n",
    "F2 = abs(Target-F)/Target*100.\n",
    "G2 = abs(Target-G)/Target*100."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "loaded-copying",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.9966394827944038 0.4945458276011555 \n",
      "\n",
      "[0.2286734655417619, 0.2120133110483818, 0.17139201078671926, 0.08007338174632164, 0.011328149511868692, 0.29651968136494666]\n",
      "[7.629376953609119, 7.073533719090296, 5.718259676657221, 2.6715387019022674, 0.3779481917944938, 9.892973012577597]\n"
     ]
    }
   ],
   "source": [
    "from sklearn.metrics import r2_score\n",
    "R2_Score = r2_score(A,Target)\n",
    "print(R2_Score,np.median(A2),'\\n')\n",
    "\n",
    "ALL_ARE = [np.mean(B2)-np.mean(A2),\n",
    "           np.mean(C2)-np.mean(A2),\n",
    "           np.mean(D2)-np.mean(A2),\n",
    "           np.mean(E2)-np.mean(A2),\n",
    "           np.mean(F2)-np.mean(A2),\n",
    "           np.mean(G2)-np.mean(A2)]\n",
    "\n",
    "ALL_ARE_P = [i/np.sum(ALL_ARE) for i in ALL_ARE]\n",
    "print(ALL_ARE_P)\n",
    "print(ALL_ARE)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "ECC",
   "language": "python",
   "name": "ecc"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}