atom_bond_ring_and_pharmacophore_search.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 分子基础操作与药效团查找"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1.原子操作"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "在rdkit中，分子中的每一个原子都是对象，可以通过原子对象的属性和函数来获取各种信息。\n",
    "- 对原子进行遍历：m.GetAtoms()\n",
    "- 获取原子索引：GetIdx()\n",
    "- 获取原子序号：GetAtomicNum()\n",
    "- 获取原子符号：GetSymbol()\n",
    "- 获取原子连接数（受H是否隐藏影响）：GetDegree()\n",
    "- 获取原子总连接数（与H是否隐藏无关）：GetTotalDegree()\n",
    "- 获取原子形式电荷：GetFormalCharge()\n",
    "- 获取原子杂化方式：GetHybridization()\n",
    "- 获取原子显式化合价：GetExplicitValence()\n",
    "- 获取原子隐式化合价：GetImplicitValence()\n",
    "- 获取原子总的化合价：GetTotalValence()\n",
    "- ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "id\tnum\tsymbol\tdegree\tcharge\thybrid\n",
      "0\t6\tC\t2\t0\tSP3\n",
      "1\t8\tO\t2\t0\tSP3\n",
      "2\t6\tC\t2\t0\tSP3\n"
     ]
    }
   ],
   "source": [
    "from rdkit import Chem\n",
    "m = Chem.MolFromSmiles('C1OC1')\n",
    "print('\\t'.join(['id', 'num', 'symbol', 'degree', 'charge', 'hybrid']))\n",
    "for atom in m.GetAtoms():\n",
    "    print(atom.GetIdx(), end='\\t')\n",
    "    print(atom.GetAtomicNum(), end='\\t')\n",
    "    print(atom.GetSymbol(), end='\\t')\n",
    "    print(atom.GetDegree(), end='\\t')\n",
    "    print(atom.GetFormalCharge(), end='\\t')\n",
    "    print(atom.GetHybridization())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 也可以通过索引获取原子：GetAtomWithIdx()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "C\n"
     ]
    }
   ],
   "source": [
    "print(m.GetAtomWithIdx(0).GetSymbol())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 获取相连的原子：GetNeighbors()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[6, 6]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "atom = m.GetAtomWithIdx(1)\n",
    "[x.GetAtomicNum() for x in atom.GetNeighbors()]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2.键操作\n",
    "同样，每一个键也都是对象，可以通过属性和函数来获取键的信息。\n",
    "- 对键进行遍历：m.GetBonds()\n",
    "- 获取键的索引：GetIdx()\n",
    "- 获取键的类型：GetGetBondType()\n",
    "- 以数字形式显示键的类型：GetBondTypeAsDouble()\n",
    "- 是否为芳香键：GetIsAromatic()\n",
    "- 是否为共轭键：GetIsConjugated()\n",
    "- 是否在环中：IsInRing()\n",
    "- 是否在n元环中：IsInRingSize(n)\n",
    "- 获取起始原子：GetBeginAtom()\n",
    "- 获取末尾原子：GetEndAtom()\n",
    "- ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "id\ttype\tdouble\taromic\tconjug\tring\tbegin\tend\n",
      "0\tDOUBLE\t2.0\tFalse\tTrue\tFalse\t0\t1\n",
      "1\tSINGLE\t1.0\tFalse\tTrue\tFalse\t1\t2\n",
      "2\tDOUBLE\t2.0\tFalse\tTrue\tFalse\t2\t3\n"
     ]
    }
   ],
   "source": [
    "m = Chem.MolFromSmiles('C=C-C=C')\n",
    "print('\\t'.join(['id', 'type', 'double', 'aromic', 'conjug', 'ring', 'begin', 'end']))\n",
    "for bond in m.GetBonds():\n",
    "    print(bond.GetIdx(), end='\\t')\n",
    "    print(bond.GetBondType(), end='\\t')\n",
    "    print(bond.GetBondTypeAsDouble(), end='\\t')\n",
    "    print(bond.GetIsAromatic(), end='\\t')\n",
    "    print(bond.GetIsConjugated(), end='\\t')\n",
    "    print(bond.IsInRing(), end='\\t')\n",
    "    print(bond.GetBeginAtomIdx(), end='\\t')\n",
    "    print(bond.GetEndAtomIdx())  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 也可以通过索引获取键：GetBondWithIdx()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DOUBLE\n"
     ]
    }
   ],
   "source": [
    "print(m.GetBondWithIdx(0).GetBondType())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.环操作\n",
    "- 查看原子的环信息"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n",
      "True\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcIAAACWCAIAAADCEh9HAAAABmJLR0QA/wD/AP+gvaeTAAATk0lEQVR4nO3da0xT5x8H8KedICLDC3i/Tu5YdIIioGJQyLwUXknfOEu2JeBeSGsQwelG2TQBIbEHM03JsoTimxUT/7EIGhAviApemBa54xXvIiiCCq79v3jwDB0qFtqnl+8nvpDfSvPNlnx3zunp7wgMBgMBAABjCVkHAACwbqhRAIAhQY0CAAwJahQAYEhQowAAQzKCdQAAe9TW1tbd3c06BfRxdXUdM2aM0b+OGgUwqyNHjqjV6t7e3kOHDrHOAn2SkpKys7ON/nXUKICZNDY2JiUlFRYWEkLc3NxYx4F/HTt2rL29fd++fSNHjjTi1wW4/R7A1Do6OtLT03///ffe3l5+6ODgMGnSJIapgLp79y6twc7OThcXFyPeAUejACak1+sPHDiQnJz86NEjOhEKhXq9nhDi4+Oj0+mYpgNCCHFxcenq6hrKO+CTegBTOXnyZGBgYFxcHN+hy5cv/+uvv9imgmGHGgUYfnfu3JFKpREREVeuXKGT6dOn5+XlnThxwsfHh202GHY4qQcYTt3d3bt3787MzHz16hWdODs7Jycnp6amOjk5sc0GJoIaBRgeBoPh4MGDW7ZsuX37Np0IBIJ169ZlZWXNmjWLbTYwKdQowDC4ePGiXC6vqKjgJ0FBQRzHLVmyhGEqMA9cGwUYknv37iUkJCxevJjv0ClTpqhUqqqqKnSoncDRKICRenp69u/f/8svvzx//pxOHB0dN27c+Ntvv7m6urLNBuaEGgUwhlarlcvl169f5ydisZjjuDlz5jBMBUzgpB7g89TV1a1atSomJobvUF9f3+LiYq1Wiw61T6hRgMF6+vSpTCYLCAg4duwYnYwfP16pVOp0ulWrVrHNBgzhpB7g0968efPnn39u3779yZMndDJixIjvv/9+165d7u7ubLMBc6hRgE84fvy4XC6vqanhJytXrlQqlSKRiGEqsBw4qQf4oKamJolEEhkZyXeol5eXRqMpLS1FhwIPR6MAA3jx4kV2dnZGRsbr16/pZPTo0Vu2bNm2bZtxKynBhqFGAd5BV9tt3br14cOHdCIUCtevX5+VlYX1oDAg1CjAvyorK2UyWWVlJT8JDg7mOC4kJIRhKrBwuDYKQAghra2tUqk0NDSU71C62u78+fPoUPg4HI2Cvevu7t67d+/OnTtfvHhBJ87Ozps2bdqxY4dxj5QAe4MaBbum1Wo3bdp069YtfiIWi/fu3Tt79mx2ocDK4KQe7NTly5eXLVsWExPDd2hgYODp06e1Wi06FD4LahTszv379xMSEoKDg8+cOUMn7u7uSqWyqqpq2bJlbLOBNcJJPdiR3t7effv29V9t5+Dg8OOPP2K1HQwFahTshVar3bx5c0tLCz+JjIzMycnx8/NjmApsAE7qwfbV19evWbMmJiaG71AfH58jR46UlJSgQ2HoUKNgy9rb2+lqu+LiYjoZN26cUqmsqalZs2YN22xgM3BSD7aJrrbbsWPH48eP6YSuttu5c+eECRPYZgMbgxoFG1RWViaXy3U6HT9ZsWLFnj175s2bxzAV2Cqc1INNaW5ulkgkK1eu5DvU09NTo9EcP34cHQomgqNRsBFdXV1ZWVn/XW2Xmprq5OTENhvYNtQoWD2DwZCfn5+SkvLgwQM6EQgE33777e7duydPnsw2G9gD1ChYt6qqKrlcfu7cOX6yaNEijuNCQ0MZpgK7gmujYK3u3r0rlUpDQkL4Dp02bVpeXl5lZSU6FMwJR6NgfV6+fJmTk9N/td2oUaMSExO3b9/+5Zdfss0Gdgg1ClZGq9UmJibevHmTn4jF4pycnK+++opdKLBrOKkHq1FdXR0eHh4TE8N36IIFC06dOqXVatGhwBBqFKxAW1ubTCZbtGhReXk5nbi5uSmVygsXLoSHh7PNBoCTerBodLVdWlras2fP6ISutvv111/HjBnDNhsAhRoFy1VaWiqTyWpra/lJZGQkx3H+/v4MUwG8Byf1YIkaGhrEYnFUVBTfod7e3oWFhSUlJehQsDQ4GgXL0t7enpmZuWfPnp6eHjoZO3Zsamrq5s2bHR0d2WYDGBBqFCyFXq8/cOBAcnLyo0eP6EQoFK5fvz47O3vixIlsswF8BGoULMKJEyfkcvnVq1f5SURExJ49e+bPn88wFcBg4NooMHbnzh2pVLpixQq+Q2fMmJGXl1dWVoYOBauAo1Fghq62y8zMfPXqFZ04OzsnJydjtR1YF9QoMGAwGA4ePLhly5bbt2/TiUAgWLduXXZ29syZM9lmA/hcqFEwtwsXLsjl8rNnz/KThQsXchwXFhbGMBWA0XBtFMzn3r17CQkJISEhfIdOnTpVpVJVVlaiQ8F64WgUzKGnp2f//v0///xzZ2cnnTg6Om7cuHHnzp1YbQfWDjUKJqfVamUy2Y0bN/iJWCzmOG7OnDkMUwEMF5zUgwnV1dWtWrUqJiaG71A/P7+jR49qtVp0KNgM1CiYxNOnT2UyWUBAwLFjx+hk/PjxSqVSp9N98803bLMBDC+c1MMwo6vtFApFR0cHnTg4OHz33Xe7du1yd3dnmw3AFIw9Gi0uJgLB+3847p3XcNwAr2luHvgN33vxkiVGBgOmSktLFyxYIJfL+Q6NjIy8fPmySqVCh4KtMrZGV68mBgMpKur7UakkBgORyd55jUxGDAYSH9/3Y1ERMRiIp+f7b9XcTAQCIpeTpiZiMPT9EYkG6GWwYI2NjdHR0VFRUdeuXaMTLy8vjUZTUlIiEonYZgMwqaGd1Ht59f3F2/uDr+G3Q65ePcA/bW7uexOD4Z25SkX8/YlcTgh5v53BwnR0dGRkZCiVytevX9OJi4tLUlLStm3bRo4cyTYbgBmwvjZKO5Q/qu1PJiMaDZHLibf3wBUMrNHVdlu3bn348CGd0NV2WVlZkyZNYpsNwGyYflKfkEAIIWFhH2zJHTsIIWTnTvNFgkE7depUYGBgXFwc36HLly+/dOmSWq1Gh4JdYVejzc0kN5eQt105IFqvZ8+S4mIzpYJBaG1tlUqlERERV65coZPp06fn5eWdOHHi66+/ZpsNwPzY1eiRI31/4S+wDoh+QoUDUsvQ3d2tUCi8vb3z8/MNBgMhxNnZOS0trbGxUSqVCgQC1gEBGGB3bVSjIYSQsLABPrvvj35CdfYsaW7+xCvBlOhqu+Tk5Fu3btEJXW2XlZU1a9YsttkA2GJ3NEp3/HzyVhj+HoCmJtPmgQ+7dOnSsmXLJBIJ36FBQUGnT5/WaDToUABGNTr4a538Kf///meiLPAR9+/fT0hICA4OrqiooJMpU6aoVKqqqqqlS5eyzQZgIVjf8DT4Z47X1JCWFpKfb8o0plLg6nrt+XPWKT7blStXSkpKurq66I8jR47cvHnzTz/9hNV2AP0xqtHGRmN+q6WFpKcPdxRzKAgJKTh/nnWKz+bk5MQ/JUksFiuVSg8PD7aRACzQMNXomjXD8z7/xX+s1O+ZEwAAlmOYarSo6IO30HNc33c6h87Dg6SlDc9bmVesq6u/FW6H639SX1hYWFJSgpN6gP9ifW30s3h4EIWCdQhjxBISyzqDce7fv69QKP744w+9Xv/69euMjIy8vDyFQvHDDz988cUXrNMBWASLX9vML9bDI89Y4D+XX/J2dSH97H7x4sVnzpxhmw3AQjCq0Y9shALLExQUVF5e3v8u0UuXLoWHh/e/kxTAbjGqUf5u0Nrawf4KdlYyJRAIYmNja2tr09LSRo0aRQgxGAwFBQX+/v4KheLly5esAwIww6hG+c/fa2o+8Ur+y0uDv8MUTMbZ2VmhUDQ2Nm7YsIF+g767uzs9Pd3b21utVhveWxoLYB/YXRul1zo/eRsTf4fp2rWmzQODNn36dLVa3X+fU2tra1xcXGhoaGVlJdtsAObHrkYlkr6/fOjpTBQ96//kBhMwO7pdNC8vj98uWllZGRYWJpVK+Q2kAPaAXY3yjwbhN+YNiO4k5TsXLIlQKJRKpfX19SkpKfR5IXq9Pj8/39PTU6FQ8M8UAbBtTG94UioJebsxb0D8BhM8jsmCjR07NiMjQ6fTxcb23R374sWL9PT0gICAgoICttkAzIBpjcpkJCzsY8vt6bbmAZ/UBBaGfw7o3Llz6aSpqUkikURFRdV88oNEAGvG+vZ7un5tzZoBrpByHDl7liiVeJ6dFYmMjKyuru7/VHr65PqEhIQnT56wzQZgIkOrUf5upI9sbOLvDP3QISe9S8bL650XLFlC5HKiVOJ03uo4ODjEx8c3NDQkJibSL4y+efMmNzfXx8eH47h//vmHdUCA4WYwTlGRgZD3/yiV77xGqRzgNU1Ng3rDsDAjg4Elqa2t/ebdnSx+fn5Hjx5lnYulq1ev0n8VIpGIdRYwGAyG0aNH0/8inZ2dxr2DwIBbpsHEtFqtXC6/fv06PxGLxRzHzZkzh2EqVnQ63bx58wghIpFIp9OxjgPExcWFrjHr7Ox0cXEx4h1YXxsFOxAdHV1XV6dUKvkNe4WFhX5+fjKZrLOzk202gKFDjYI5ODo6ymSy+vr6+Ph4oVBICOnp6cnJyfH19c3NzdXr9awDAhgPNQrmM3XqVJVKRb/sRCf37t2ja/fO4ukGYLVQo2BuCxcuPHPmjEajmTlzJp1cvHhx6dKlEonk9u3bbLMBGAE1Cgz0X7vn5ORE3q7d8/PzUygU/HP0AKwCahSYGT16NL92j076r91jmw1g8FCjwNiMGTPUanVZWRm9DYgQcufOnbi4uIiIiCtXrrDNBjAYqFGwCBEREdXV1Xl5eRMnTqSTkydPBgYGSqXSR48esc0G8HGoUbAUdO1eQ0NDSkqKo6Mjebt2z8fHJzMzs6enh3VAgIGhRsGy8Gv31r593kFHR0dqampAQMCRj6+mBWAENQqWyNvbu7CwsKSkxP/tM7gaGxvFYnFUVNS1a9fYZgN4D2oULFdkZOTff/+tVCrHjBlDJ3Ttnkwme/bsGdtsADzUKFg0BwcHmUzW0tLCr93r7e3Nycnx8PDA2j2wEKhRsAJubm4cx124cCE8PJxO2tra5HL5okWLTp8+zTYbAGoUrMaCBQtOnTp1+PDh2bNn00l1dfXy5cujo6Nv3LjBNBrYNdQoWJno6Oja2tqMjIz+a/fmzp2bmpqKtXvABGoUrM+oUaNSUlLq6uo2bNggEAgIIS9fvszMzPTz81Or1dhEDmaGGgVrNW3aNLVaXVlZGRoaSid3796Ni4tbvHjxuXPn2GYDu4IaBeu2aNGiioqKvLy8yZMn08mFCxeWLFkilUofPHjANhvYCdQoWD2BQCCVSpubm/uv3cvPz/f09MTaPTAD1CjYCLp2T6fTxcbG0klXV1d6erpIJCooKGCbDWwbahRsiqenp0ajOX78eEBAAJ20tLRIJJKVK1fyTzYGGF6oUbBBK1asuHz5skqlmjBhAp2UlZUtWLBAKpU+fvyYbTawPahRsE0jRoyIj49vaGhITEwcMWIEwdo9MBnUKNiycePGcRxXU1OzevVqOmlvb09NTZ03b15RURHbbGAzUKNg+3x8fIqKig4fPuzh4UEnDQ0Na9eujYqKqqurY5sNbABqFOxFdHR0XV2dUql0dXWlk9LS0vnz52PtHgwRahTsCNbugSmgRsHuuLu7cxxXVVW1dOlSOqFr94KDg8vLy9lmA2uEGgU7FRgYWF5efvjw4VmzZtHJ5cuXw8PDo6Ojb968yTQaWBnUKNg1fu2ei4sLnfBr9168eME2G1gL1CjYO2dn5/fW7nV3d2dmZvr6+mLtHgwGahSAEEKmT5+uVqvPnTsXEhJCJ3TtXkhIyPnz59lmAwuHGgX41+LFi+navUmTJtFJVVUV1u7Bx6FGAd4hFAr5tXsjR44kb79FStfuvX79mnVAsDioUYABuLi4YO0eDBJqFOCDvLy8NBpNaWmpSCSik+bmZolEEhkZqdPp2GYDy4EaBfiElStXVldXq1Qqd3d3Ojl+/HhgYGBCQsKTJ0/YZgNLgBoF+LT/rt178+ZNbm6uj48Px3Fv3rxhHRBYQo0CDNb48eM5jtPpdKtWraKTp0+fyuXygICAo0ePss0GDKFGAT6Pr69vcXHx4cOH58yZQyf19fWrV6+Ojo6+fv0622zABGoUwBj/XbtXWFjo5+cnk8meP3/ONhuYGWoUwEiOjo4ymay+vj4+Pl4oFBJCenp6cnJyfH19c3Nz9Xo964BgJqhRgCGZMmWKSqWiX3aik/v37yckJAQHB1dUVLDNBuaBGgUYBkFBQeXl5RqNZubMmXRy6dKlZcuWSSSSW7dusc0GpoYaBRgeAoEgNja2rq4uLS3NycmJEGIwGAoKCvz9/RUKxcuXL1kHBFNBjQIMJ2dnZ4VC0dTU1H/tXnp6ure3N9bu2SrUKMDwo2v3ysrK5s+fTyetra1xcXERERH19fVss8GwE+B/jwCmo9frDxw4kJyc/OjRIzoRCoX0Q3yRSIQv5lsCFxeXrq4uQkhnZyf/EITPghoFMLmOjo709PTff/+9t7eXHzo4OPBbTYGhu3fv0hpEjQJYusbGxqSkpMLCQkKIm5tbW1sb60TQRyQSBQcH79u3j26Y/VyoUQCzOnLkiFqt7u3tPXToEOss0CcpKSk7O9voX0eNAjDQ1tbW3d3NOgX0cXV1HTNmjNG/jhoFABgS3PAEADAkqFEAgCFBjQIADAlqFABgSFCjAABD8n/mYj3DXceLkQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<rdkit.Chem.rdchem.Mol at 0x4f7d8a0>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "m = Chem.MolFromSmiles('OC1C2C1CC2')\n",
    "print(m.GetAtomWithIdx(0).IsInRing())\n",
    "print(m.GetAtomWithIdx(2).IsInRing())\n",
    "print(m.GetAtomWithIdx(2).IsInRingSize(3))\n",
    "m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "Chem.Draw.MolToFile(m, 'data/atom_bond_ring_1.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 但是要注意，IsInRingSize()函数只能判断是否在最小的环中"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True\n",
      "False\n"
     ]
    }
   ],
   "source": [
    "print(m.GetAtomWithIdx(2).IsInRingSize(3))\n",
    "print(m.GetAtomWithIdx(2).IsInRingSize(5))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 可以查看所有最小环(smallest set of smallest rings, SSSR）的信息：GetSymmSSSR()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "[1, 2, 3]\n",
      "[4, 5, 2, 3]\n"
     ]
    }
   ],
   "source": [
    "ssr = Chem.GetSymmSSSR(m)\n",
    "print(len(ssr))\n",
    "print(list(ssr[0]))\n",
    "print(list(ssr[1]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- 直接获取环的信息：GetRingInfo()\n",
    "- 查看一共有几个环：NumRings()\n",
    "- 查看原子在几个环中：NumAtomRings()\n",
    "- 查看id为n的原子是否在n1元环中.IsAtomInRingOfSize(n, n1)\n",
    "- 查看id为n的键是否在n1元环中.IsAtomInRingOfSize(n , n1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2\n",
      "0\n",
      "True\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "ri = m.GetRingInfo()\n",
    "print(ri.NumRings())\n",
    "print(ri.NumAtomRings(6))\n",
    "print(ri.IsAtomInRingOfSize(1,3))\n",
    "print(ri.IsBondInRingOfSize(1,3))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4.手动实现氧族药效团查找\n",
    "这里只做一个简单的演示，通过原子属性查找目标药效团，更复杂的操作类似。\n",
    "- 假设要查找的氧族氢供体标准：氧或硫原子，不带电荷，含有1个氢。\n",
    "- 氧族氢受体标准（部分）：不带氢原子，化合价为2，且不与氮原子相连。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{1: [8, 'HAcceptor'], 3: [8, 'HAcceptor'], 4: [8, 'HDonor']}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pharmacophore(m):\n",
    "    atomPharma = {}\n",
    "    # 定义氧族原子\n",
    "    chalcogen = [8, 16]\n",
    "    mol = Chem.MolFromSmiles(m)\n",
    "    mol = Chem.AddHs(mol)\n",
    "    # 开始查找\n",
    "    for atom in mol.GetAtoms():\n",
    "        p = []\n",
    "        if atom.GetAtomicNum() == 1 or atom.GetAtomicNum() not in chalcogen:\n",
    "            continue\n",
    "        else:\n",
    "            # 氢供体\n",
    "            if atom.GetFormalCharge() == 0:\n",
    "                nbrs = [x for x in atom.GetNeighbors()]\n",
    "                HDflag = False\n",
    "                for nbr in nbrs:\n",
    "                    if nbr.GetAtomicNum() == 1:\n",
    "                        HDflag = True\n",
    "                if HDflag:\n",
    "                    p.append('HDonor')\n",
    "                    \n",
    "            # 氢受体\n",
    "            if atom.GetTotalValence() == 2:\n",
    "                nbrs = [x for x in atom.GetNeighbors()]\n",
    "                HAflag_1 = True\n",
    "                HAflag_2 = True\n",
    "                if len(nbrs) == 1:\n",
    "                    nbr = nbrs[0]\n",
    "                    if nbr.GetAtomicNum() == 7:\n",
    "                        HAflag_1 = False\n",
    "                else:\n",
    "                    for nbr in nbrs:\n",
    "                        if nbr.GetAtomicNum() == 1:\n",
    "                            HAflag_2 = False\n",
    "                if HAflag_1 and HAflag_2:\n",
    "                    p.append('HAcceptor')\n",
    "        atomPharma[atom.GetIdx()] = [atom.GetAtomicNum(), ' '.join(p)]\n",
    "    return atomPharma\n",
    "\n",
    "m = 'COC(=O)O'\n",
    "res = pharmacophore(m)\n",
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "更简单的操作，直接写成SMARTS查找就可以了\n",
    "- 氢供体：[O,S;H1;+0]\n",
    "- 氢受体（部分）：[O;H0;v2;!$(O=N-*)] "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{4: [8, 'HDonor'], 1: [8, 'HAcceptor'], 3: [8, 'HAcceptor']}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def pharmacophore_smarts(m):\n",
    "    # 定义smarts\n",
    "    HDonorSmarts = '[O,S;H1;+0]'\n",
    "    HAcceptorSmarts = '[O;H0;v2;!$(O=N-*)]'\n",
    "    HDonor = Chem.MolFromSmarts(HDonorSmarts)\n",
    "    HAcceptor = Chem.MolFromSmarts(HAcceptorSmarts)\n",
    "    atomPharma = {}\n",
    "    mol = Chem.MolFromSmiles(m)\n",
    "    # 氢供体\n",
    "    HDonors = mol.GetSubstructMatches(HDonor)\n",
    "    for i in HDonors:\n",
    "        atom = mol.GetAtomWithIdx(i[0])\n",
    "        atomPharma[atom.GetIdx()] = [atom.GetAtomicNum(), 'HDonor']\n",
    "    # 氢受体\n",
    "    HAcceptors = mol.GetSubstructMatches(HAcceptor)\n",
    "    for i in HAcceptors:\n",
    "        atom = mol.GetAtomWithIdx(i[0])\n",
    "        atom_prop = atomPharma.get(atom.GetIdx(), [])\n",
    "        if atom_prop:\n",
    "            atom_prop[1] = atom_prop[1] + ' HAcceptor'\n",
    "            atomPharma.update({atom.GetIdx():atom_prop})\n",
    "        else:\n",
    "            atomPharma[atom.GetIdx()] = [atom.GetAtomicNum(), 'HAcceptor']\n",
    "    return atomPharma\n",
    "\n",
    "res = pharmacophore_smarts(m)\n",
    "res"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}