example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# RNA-Puzzles Assessment example"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> This notebook provides a tutorial on RNA 3D structure comparison in RNA-Puzzles"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-01T15:11:49.527225Z",
     "iopub.status.busy": "2024-02-01T15:11:49.526475Z",
     "iopub.status.idle": "2024-02-01T15:11:49.722136Z",
     "shell.execute_reply": "2024-02-01T15:11:49.721120Z",
     "shell.execute_reply.started": "2024-02-01T15:11:49.527157Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "import sys,os\n",
    "# import warnings\n",
    "# warnings.filterwarnings(\"ignore\")\n",
    "\n",
    "import RNA_normalizer\n",
    "\n",
    "from operator import attrgetter\n",
    "\n",
    "RESIDUES_LIST = \"data/residues.list\"\n",
    "ATOMS_LIST = \"data/atoms.list\"\n",
    "\n",
    "def CleanFormat(f):\n",
    "\tos.system( \"mac2unix -q %s\" %f )\n",
    "\tos.system( \"dos2unix -q %s\" %f )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> CleanFormat is a function used to format different platform formats to unix. Users need to install dos2unix"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Normalize structures\n",
    "> Different prediction methods may give different PDB formats. The standard PDB format considered is the \n",
    "[1992](https://www.rcsb.org/pdb/file_formats/pdb/pdbguide2.2/PDB_format_1992.pdf) format. To normalize the structure format, we need to:   \n",
    "* Correct the residue names and the atom names: The name mapping dictionaries are: defined as \"RESIDUES_LIST\" and \"ATOMS_LIST\" in the data folder. \n",
    "* select the part of structure we are interested. The part of structure is assigned in \"index_file\". The selection grammar is: [Chain id]:[Beg of Residue number]:[length of fragment]\n",
    "    + [Chain id] is single letter character\n",
    "    + [Beg of Residue number] is the number of the beginning residue in the PDB file. This number is exact the \"res. seq no.\" field of the PDB file without any renumbering. \n",
    "    + [length of fragment] is the length of the fragment. numbering of the residues in the fragment is not considered.\n",
    "    + for example, A:1:31,A:33:29 includes two fragments of 31 residue and 29 residues. The first fragment starts from residue 1 and until residue 31. leaving out residue 32, the second fragment starts from residue 33 and ends till 29 residues after. The number of these 29 residues are not considered except the first one. \n",
    "* extract the interesting part. To do this, we need to set the output file for \"extract_file\" parameter. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-01T15:11:51.085220Z",
     "iopub.status.busy": "2024-02-01T15:11:51.084338Z",
     "iopub.status.idle": "2024-02-01T15:11:51.096119Z",
     "shell.execute_reply": "2024-02-01T15:11:51.094569Z",
     "shell.execute_reply.started": "2024-02-01T15:11:51.085161Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def normalize_structure(struct, out_file = None, index_file=None, extract_file = None):\n",
    "\tpdb_normalizer = RNA_normalizer.PDBNormalizer( RESIDUES_LIST, ATOMS_LIST )\n",
    "\tok = pdb_normalizer.parse( struct, out_file )\n",
    "\tif not ok:\n",
    "\t\tsys.stderr.write(\"ERROR: structure not normalized!\\n\")\n",
    "\telse:\n",
    "\t\tsys.stderr.write(\"INFO: Normalization succeded!\\n\")\n",
    "\tif not extract_file is None:\n",
    "\t\tcoords=open(index_file).read()\n",
    "\t\tRNA_normalizer.extract_PDB(SOLUTION_NORMAL,coords, extract_file)\n",
    "\t\tsys.stderr.write(\"INFO:\tstructure extracted\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculate RMSD and P-value\n",
    "> This function calculate the all atom RMSD and the P-value of prediction. \n",
    "* Before calculation, we need to first make sure that the two structure have the same sequence. \n",
    "If we need to use the whole structure instead of the selected part assigned by index_file, we can give an empty index_file.\n",
    "* The RMSD calculated is all atom RMSD.\n",
    "* the P value set is according to ***Hajdin et al., RNA (7) 16, 2010***. By default, it is set to \"-\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-01T15:11:53.677579Z",
     "iopub.status.busy": "2024-02-01T15:11:53.676922Z",
     "iopub.status.idle": "2024-02-01T15:11:53.692028Z",
     "shell.execute_reply": "2024-02-01T15:11:53.690482Z",
     "shell.execute_reply.started": "2024-02-01T15:11:53.677521Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "# PVALUE set according to Hajdin et al., RNA (7) 16, 2010, either \"+\" or \"-\"\n",
    "def calc_RMSD(native_file, native_index, prediction_file, prediction_index, PVALUE = \"-\"):\n",
    "\tres_struct = RNA_normalizer.PDBStruct()\n",
    "\tres_struct.load( native_file, native_index )\n",
    "\tres_raw_seq = res_struct.raw_sequence()\n",
    "\t\n",
    "\tsol_struct = RNA_normalizer.PDBStruct()\n",
    "\tsol_struct.load( prediction_file, prediction_index )\n",
    "\tsol_raw_seq = sol_struct.raw_sequence()\n",
    "\t\n",
    "\tif( sol_raw_seq != res_raw_seq ):\n",
    "\t\tsys.stderr.write(\"ERROR Result sequence != Solution sequence!\\n\")\n",
    "\t\tsys.stderr.write(\"DATA Solution sequence --> '%s'\\n\" %sol_raw_seq )\n",
    "\t\tsys.stderr.write(\"DATA Result sequence   --> '%s'\\n\" %res_raw_seq )\n",
    "\t\treturn(-1)\n",
    "\t# computes the RMSD\n",
    "\tcomparer = RNA_normalizer.PDBComparer()\n",
    "\trmsd = comparer.rmsd( sol_struct, res_struct )\n",
    "\tsys.stderr.write(\"INFO Partial RMSD --> %f\\n\" %rmsd )\n",
    "\tpvalue = comparer.pvalue( rmsd, len(sol_raw_seq), PVALUE )\n",
    "\tsys.stderr.write(\"INFO Partial P-Value --> %e\\n\" %pvalue )\n",
    "\treturn(rmsd, pvalue)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interaction Network Fidelity\n",
    "> Interaction Network Fidelity is defined by ***RNA. 2009 Oct; 15(10): 1875–1885.***.\n",
    "* It checks the percentage of interactions that are found in the predicted structure. \n",
    "* The interaction types are assigned according to [MC-annotate](http://www-lbit.iro.umontreal.ca/mcannotate-simple/). \n",
    "* The interaction types can be sorted to: All interactions (INF_ALL), Watson-Crick interactin (INF_WC), non-Watson-Crick interaction \n",
    "   (INF_NWC) and Stacking (INF_STACK). \n",
    "* Deformation Index is defined as RMSD/INF_ALL (***RNA. 2009 Oct; 15(10): 1875–1885.***). \n",
    "* This function returns a tuple of the metrices [RMSD, Deformation Index, INF_ALL, INF_WC, INF_NWC,INF_STACK]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-01T15:11:56.465366Z",
     "iopub.status.busy": "2024-02-01T15:11:56.464348Z",
     "iopub.status.idle": "2024-02-01T15:11:56.480362Z",
     "shell.execute_reply": "2024-02-01T15:11:56.478768Z",
     "shell.execute_reply.started": "2024-02-01T15:11:56.465297Z"
    },
    "tags": []
   },
   "outputs": [],
   "source": [
    "def InteractionNetworkFidelity(native_file, native_index, prediction_file, prediction_index):\n",
    "\tres_struct = RNA_normalizer.PDBStruct()\n",
    "\tres_struct.load( native_file, native_index )\n",
    "\tres_raw_seq = res_struct.raw_sequence()\n",
    "\t\n",
    "\tsol_struct = RNA_normalizer.PDBStruct()\n",
    "\tsol_struct.load( prediction_file, prediction_index )\n",
    "\tsol_raw_seq = sol_struct.raw_sequence()\n",
    "\t\n",
    "\tif( sol_raw_seq != res_raw_seq ):\n",
    "\t\tsys.stderr.write(\"ERROR Result sequence != Solution sequence!\\n\")\n",
    "\t\tsys.stderr.write(\"DATA Solution sequence --> '%s'\\n\" %sol_raw_seq )\n",
    "\t\tsys.stderr.write(\"DATA Result sequence   --> '%s'\\n\" %res_raw_seq )\n",
    "\t\treturn(-1)\n",
    "\t# computes the RMSD\n",
    "\tcomparer = RNA_normalizer.PDBComparer()\n",
    "\trmsd = comparer.rmsd( sol_struct, res_struct )\n",
    "\tINF_ALL = comparer.INF( sol_struct, res_struct, type=\"ALL\" )\n",
    "\tDI_ALL = rmsd / INF_ALL\n",
    "\tINF_WC = comparer.INF( sol_struct, res_struct, type=\"PAIR_2D\" )\n",
    "\tINF_NWC = comparer.INF( sol_struct, res_struct, type=\"PAIR_3D\" )\n",
    "\tINF_STACK = comparer.INF( sol_struct, res_struct, type=\"STACK\" )\n",
    "\treturn (rmsd,DI_ALL, INF_ALL, INF_WC, INF_NWC,INF_STACK)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Local Distance Difference Test (lDDT)\n",
    "> The Local Distance Difference Test (lDDT) is defined by ***Bioinformatics. 2013 Nov 1;29(21):2722-8.***.\n",
    "* It checks the percentage of interactions that are found in the predicted structure. \n",
    "* [OpenStructure](https://openstructure.org/download) for lDDT calculations.\n",
    "* This function returns the lDDT score"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-01T15:11:58.056499Z",
     "iopub.status.busy": "2024-02-01T15:11:58.055784Z",
     "iopub.status.idle": "2024-02-01T15:11:58.083725Z",
     "shell.execute_reply": "2024-02-01T15:11:58.082651Z",
     "shell.execute_reply.started": "2024-02-01T15:11:58.056437Z"
    }
   },
   "outputs": [],
   "source": [
    "import os\n",
    "#import extract\n",
    "import sys\n",
    "sys.path.append('./RNA_normalizer')\n",
    "sys.path.append('./lddt/bin')\n",
    "import extract\n",
    "\n",
    "from structure_processed import preprocess\n",
    "from metrics_lddt import LDDT\n",
    "\n",
    "def extract_and_preprocess(parent_dir, pdb_list):\n",
    "    # extract and process PDB\n",
    "    for pdb_name in pdb_list:\n",
    "        coords = \"\"\n",
    "        index_file_path = os.path.join(parent_dir, f\"{pdb_name.split('.pdb')[0]}.index\")\n",
    "        with open(index_file_path, \"r\") as index_file:\n",
    "            for row in index_file:\n",
    "                row = row.strip()\n",
    "                if not row.startswith(\"#\") and row:\n",
    "                    coords += row + \",\"  \n",
    "        coords = coords.rstrip(\",\") \n",
    "        \n",
    "        pdb_path = os.path.join(parent_dir, pdb_name)\n",
    "        new_pdb_path = os.path.join(parent_dir, f\"{pdb_name.split('.pdb')[0]}_extract.pdb\")\n",
    "        extract.extract_PDB(pdb_path, coords, new_pdb_path)\n",
    "        preprocess(new_pdb_path)  # run on the new pdb structure\n",
    "\n",
    "def calculate_lddt(reference, model,chain_mapping,force=False):\n",
    "    # calculate the lDDT score\n",
    "    lddt = LDDT()\n",
    "    lddt_output=lddt.calculate(reference, model,chain_mapping,force=force)\n",
    "    return lddt_output\n",
    "    \n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# call the functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-01T15:12:01.943022Z",
     "iopub.status.busy": "2024-02-01T15:12:01.942291Z",
     "iopub.status.idle": "2024-02-01T15:12:02.325249Z",
     "shell.execute_reply": "2024-02-01T15:12:02.324402Z",
     "shell.execute_reply.started": "2024-02-01T15:12:01.942960Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "INFO: Normalization succeded!\n",
      "WARNING    >> Atom P from residue (' ', 1, ' ') not found in target atom list\n",
      "WARNING    >> Atom OP1 from residue (' ', 1, ' ') not found in target atom list\n",
      "WARNING    >> Atom OP2 from residue (' ', 1, ' ') not found in target atom list\n",
      "INFO Partial RMSD --> 7.640878\n",
      "INFO Partial P-Value --> 4.551914e-15\n",
      "WARNING    >> Atom P from residue (' ', 1, ' ') not found in target atom list\n",
      "WARNING    >> Atom OP1 from residue (' ', 1, ' ') not found in target atom list\n",
      "WARNING    >> Atom OP2 from residue (' ', 1, ' ') not found in target atom list\n"
     ]
    }
   ],
   "source": [
    "# %%capture\n",
    "\n",
    "# Normalize PDB format, correct residue names and atom names. \n",
    "normalize_structure('example/PZ14Bound_solution_4.pdb','example/PZ14Bound_solution_normalized.pdb')\n",
    "\n",
    "# calculate RMSD for RNA structures\n",
    "# require biopython\n",
    "rmsd,pval = calc_RMSD(\"example/PZ14Bound_solution_4.pdb\",\n",
    "          \"example/PZ14Bound_solution_4.index\",\n",
    "          \"example/14_ChenPostExp_2.pdb\",\n",
    "          \"example/14_ChenPostExp_2.index\")\n",
    "\n",
    "# calculate InteractionNetworkFidelity and Deformation Index for RNA structures\n",
    "# need to have MA-annotate in the directory or set in mcannotate.py\n",
    "rmsd2,DI_ALL, INF_ALL, INF_WC, INF_NWC,INF_STACK = InteractionNetworkFidelity(\"example/PZ14Bound_solution_4.pdb\",\n",
    "          \"example/PZ14Bound_solution_4.index\",\n",
    "          \"example/14_ChenPostExp_2.pdb\",\n",
    "          \"example/14_ChenPostExp_2.index\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-01T15:12:06.359044Z",
     "iopub.status.busy": "2024-02-01T15:12:06.358294Z",
     "iopub.status.idle": "2024-02-01T15:12:08.142657Z",
     "shell.execute_reply": "2024-02-01T15:12:08.140625Z",
     "shell.execute_reply.started": "2024-02-01T15:12:06.358981Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Processing /home/bufan/RNA_assessment/example/14_ChenPostExp_2_extract.pdb\n",
      "Did not create new folder.\n",
      "/home/bufan/RNA_assessment/example/processed already exists.\n",
      "Processing /home/bufan/RNA_assessment/example/14_solution_0_extract.pdb\n",
      "Did not create new folder.\n",
      "/home/bufan/RNA_assessment/example/processed already exists.\n",
      "command: /home/bufan/RNA_assessment/lddt/ema --mount /home/bufan/RNA_assessment/example/processed /home/bufan/RNA_assessment/lddt/bin/complex_lddt_no_stereocheck.py /home/bufan/RNA_assessment/example/processed/14_ChenPostExp_2_extract_processed.pdb /home/bufan/RNA_assessment/example/processed/14_solution_0_extract_processed.pdb '{\"A\": \"A\"}'\n"
     ]
    }
   ],
   "source": [
    "parent_dir = '/home/xxx/RNA_assessment/example/'\n",
    "pdb_list = ['14_ChenPostExp_2.pdb', '14_solution_0.pdb']\n",
    "\n",
    "# first step: extract and process\n",
    "extract_and_preprocess(parent_dir, pdb_list)\n",
    "\n",
    "# construct the name of the new file(_processed.pdb) \n",
    "reference = os.path.join(parent_dir, 'processed', '14_solution_0_extract_processed.pdb')\n",
    "model = os.path.join(parent_dir, 'processed', '14_ChenPostExp_2_extract_processed.pdb')\n",
    "lDDT=calculate_lddt(reference, model,chain_mapping='{\"A\": \"A\"}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Show the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "execution": {
     "iopub.execute_input": "2024-02-01T15:12:10.746661Z",
     "iopub.status.busy": "2024-02-01T15:12:10.745937Z",
     "iopub.status.idle": "2024-02-01T15:12:10.762951Z",
     "shell.execute_reply": "2024-02-01T15:12:10.761782Z",
     "shell.execute_reply.started": "2024-02-01T15:12:10.746588Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'RMSD: 7.641; P-value: 4.552e-15; Deformation Index: 10.567; INF_all: 0.723; INF_wc: 0.938; INF_nwc: 0.250; INF_stack: 0.701,lDDT: 0.548'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'RMSD: %.3f; P-value: %.3e; Deformation Index: %.3f; \\\n",
    "INF_all: %.3f; INF_wc: %.3f; INF_nwc: %.3f; INF_stack: %.3f,lDDT: %.3f'%(rmsd,pval,DI_ALL, INF_ALL, INF_WC, INF_NWC,INF_STACK,lDDT)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:UTRLM]",
   "language": "python",
   "name": "conda-env-UTRLM-py"
  },
  "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.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}