methods/rna-state-inf.ipynb

{"cells":[{"cell_type":"markdown","metadata":{"id":"TDPiVD2lK5sU"},"source":["# Standard packages"]},{"cell_type":"code","execution_count":1,"metadata":{"executionInfo":{"elapsed":5046,"status":"ok","timestamp":1645064033916,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"RDB7DJGMoEYK"},"outputs":[],"source":["import os\n","import pandas as pd\n","import time\n","\n","import numpy as np\n","import tensorflow as tf\n","\n","# disable warning TF logs\n","tf.get_logger().setLevel('ERROR')"]},{"cell_type":"markdown","metadata":{"id":"I6bOkFD6K-8K"},"source":["# Specific installs\n","\n","rna-state-inf is a pre-trained TensorFlow model."]},{"cell_type":"code","execution_count":12,"metadata":{"executionInfo":{"elapsed":308,"status":"ok","timestamp":1645064584015,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"xwSM-JQ5yKF4"},"outputs":[],"source":["method_name  = \"rna-state-inf\"\n","\n","model_fname  = 'rna-state-inf_rnnepoch50.h5'\n","gtfold_fname = 'gtfold-3.0_x86-64_ubuntu'\n","rs_fname     = 'RNAstructureLinuxTextInterfaces64bit.tgz'\n","\n","url1 = f\"https://raw.githubusercontent.com/sinc-lab/lncRNA-folding/main/methods/{model_fname}\"\n","url2 = f\"https://master.dl.sourceforge.net/project/gtfold/{gtfold_fname}.tar.gz?viasf=1\"\n","url3 = f\"http://rna.urmc.rochester.edu/Releases/current/{rs_fname}\""]},{"cell_type":"code","execution_count":13,"metadata":{"executionInfo":{"elapsed":344,"status":"ok","timestamp":1645064585437,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"LQoHDUbwqxoa"},"outputs":[],"source":["# download trained model used for estimating SHAPE data. It was trained following the instructions in this repository:\n","# https://github.com/dwillmott/rna-state-inf\n","!wget -q {url1}"]},{"cell_type":"code","execution_count":14,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":1243,"status":"ok","timestamp":1645064587004,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"1ieT2dHxrY9p","outputId":"4443bb43-14fc-4949-ec13-2b41a7a76c03"},"outputs":[{"name":"stdout","output_type":"stream","text":["GTfold: A Scalable Multicore Code for RNA Secondary Structure Prediction\n","(c) 2007-2011  D.A. Bader, C.E. Heitsch, S.C. Harvey\n","Georgia Institute of Technology\n"]}],"source":["# download method to predict secondary structure using SHAPE data\n","!wget -q {url2} -O {gtfold_fname}.tar.gz\n","!tar xfz {gtfold_fname}.tar.gz\n","\n","gfold_data = os.path.join(gtfold_fname, 'data')\n","gfold_bin  = os.path.join(gtfold_fname, 'gtmfe')\n","\n","!{gfold_bin} | head -3"]},{"cell_type":"code","execution_count":15,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":3990,"status":"ok","timestamp":1645064590992,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"h8lcipJt1Q0D","outputId":"d2cf7b84-e13e-4eb9-ff35-8970db916347"},"outputs":[{"name":"stdout","output_type":"stream","text":["ct2dot: Version 6.4 (December 8, 2021).\n","Copyright Mathews Lab, University of Rochester.\n","Invalid structure number given.\n"]}],"source":["# download utility for converting ct files to dot files\n","!wget -q {url3}\n","!tar xfz {rs_fname}\n","\n","source_path = \"/content/RNAstructure/\"\n","data_table  = os.path.join(source_path, 'data_tables')\n","ct2dot_bin  = os.path.join(source_path, 'exe', 'ct2dot')\n","\n","!{ct2dot_bin} --version"]},{"cell_type":"markdown","metadata":{"id":"bCdjGI0grQpm"},"source":["# S. cerevisiae (sce) 18 long non-coding RNA dataset\n","Data source: https://genie.weizmann.ac.il/pubs/PARS10/data/sce_genes_folded.tab.gz"]},{"cell_type":"code","execution_count":16,"metadata":{"executionInfo":{"elapsed":519,"status":"ok","timestamp":1645064591507,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"Fmvg_3l4_TBH"},"outputs":[],"source":["gh_path = \"https://raw.githubusercontent.com/sinc-lab/lncRNA-folding/master/data/\"\n","sce = pd.read_csv(gh_path + \"sce_genes_folded.tab\", delimiter='\\t', \n","                  header=None, index_col=0, \n","                  names=(\"Gene ID\", \"sequence\", \"PARS-assisted folding\"))"]},{"cell_type":"code","execution_count":17,"metadata":{"executionInfo":{"elapsed":6,"status":"ok","timestamp":1645064591508,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"eeup_dp2AHrN"},"outputs":[],"source":["# Sequences to process\n","sce18 = [\"snR81\", \"snR34\", \"snR43\", \"snR44\",  \"snR31\",  \"snR10\",\n","         \"snR63\", \"snR11\", \"snR82\", \"snR17b\", \"snR17a\", \"snR37\",\n","         \"SCR1\",  \"SRG1\",  \"snR19\", \"snR30\",  \"LSR1\",   \"TLC1\"]"]},{"cell_type":"markdown","metadata":{"id":"Vz9-OyOlshXJ"},"source":["# Compute structures"]},{"cell_type":"code","execution_count":18,"metadata":{"executionInfo":{"elapsed":6,"status":"ok","timestamp":1645064591509,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"5R8U_xEcodwa"},"outputs":[],"source":["#@title Utils\n","sequencedict = {'A' : 0, 'C' : 1, 'G' : 2, 'U' : 3}\n","for letter in 'BDEFHIJKLMNOPQRSTVWXYZ':\n","    sequencedict.update({letter : 4})\n","\n","def process_sequence(fin):\n","    with open(fin) as f:\n","        for a in f:\n","            _, seq, _ = a.strip().split('\\t')\n","\n","    seq = [sequencedict[s] for s in list(seq)]\n","    seq = tf.one_hot(seq, 5)\n","    seq = tf.convert_to_tensor([seq])\n","\n","    return seq\n","\n","def prob2shape(lstm_probs):\n","    \"\"\"This method was taken from the following link: \n","    https://github.com/dwillmott/rna-state-inf/blob/master/method.py\n","    \"\"\"\n","    a = 0.214\n","    b = 0.6624\n","    s = 0.3603\n","\n","    outstr = \"\"\n","    #for each position, write a SHAPE value depending on the probability\n","    for i, lstm_prob in enumerate(lstm_probs):\n","        if lstm_prob >= 0.5:\n","            shapevalue = ((a-s)/0.5)*(lstm_prob-1)+a\n","        else:\n","            shapevalue = ((s-b)/0.5)*lstm_prob+b\n","\n","        outstr += '%d %.3f ' % (i+1, shapevalue)\n","        outstr += '\\n'\n","\n","    return outstr\n","\n","def make_prediction(seq_fin, model_fin, outfname):\n","    seq   = process_sequence(seq_fin)\n","    model = tf.keras.models.load_model(model_fin)\n","\n","    pred  = model.predict(seq)\n","    #np.savetxt(sys.stdout.buffer, pred[0,:], delimiter='\\t')\n","    strout = prob2shape(pred[0,:,1])\n","\n","    fout = open(outfname, 'w')\n","    fout.write(strout)\n","    fout.close()\n","\n","def fas2tsv(fin):\n","  seq = \"\"\n","  with open(fin) as f:\n","    for a in f:\n","      b = a.strip()\n","      if b[0] == \">\":\n","        desc = b[1:]\n","      else:\n","        seq += b\n","  \n","  # write files\n","  foutname = fin + '.tsv'\n","  out = open(foutname, 'w')\n","  out.write('{}\\t{}\\t{}\\n'.format(desc, seq, desc))\n","  out.close()\n","\n","  foutseq = fin + '.seq'\n","  out = open(foutseq, 'w')\n","  out.write('{}\\n'.format(seq))\n","  out.close()\n","\n","  return foutname, foutseq"]},{"cell_type":"code","execution_count":19,"metadata":{"executionInfo":{"elapsed":7,"status":"ok","timestamp":1645064591510,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"6WbzBYt7uvJO"},"outputs":[],"source":["def run_folding(fasta_name):\n","  out_file_name = fasta_name + '.dot'\n","  tmp_file1 = out_file_name + '.tmp1'\n","  tmp_file2 = out_file_name + '.tmp2'\n","\n","  tsv_file, seq_file = fas2tsv(fasta_name)\n","\n","  # estimate SHAPE-like data\n","  make_prediction(tsv_file, model_fname, tmp_file1)\n","  \n","  # predict structure from estimated SHAPE data\n","  os.system(f\"export GTFOLDDATADIR={gfold_data}; {gfold_bin} {seq_file} --useSHAPE {tmp_file1} --output {out_file_name} > /dev/null 2>&1\")\n","  \n","  # convert ct to dot-bracket format\n","  os.system(f\"export DATAPATH={data_table}; {ct2dot_bin} {out_file_name + '.ct'} -1 {tmp_file2} >/dev/null 2>&1\")\n","\n","  # fix description line\n","  os.system(f\"head -1 {fasta_name}   > {out_file_name} # original description line\")\n","  os.system(f\"tail -n+2 {tmp_file2} >> {out_file_name} # sequence + predicted structure\")\n","\n","  return out_file_name"]},{"cell_type":"code","execution_count":20,"metadata":{"colab":{"base_uri":"https://localhost:8080/"},"executionInfo":{"elapsed":64738,"status":"ok","timestamp":1645064657670,"user":{"displayName":"Leandro Bugnon","photoUrl":"https://lh3.googleusercontent.com/a-/AOh14GgzVvcUV49jcdycZ0hiYEpk15dSatLCOzmAxByhTR8=s64","userId":"04594583387142575984"},"user_tz":180},"id":"EvIlcy1Vsd8t","outputId":"08ebc212-d520-4a00-e30a-717fc39186e3"},"outputs":[{"name":"stdout","output_type":"stream","text":["   \t lnc \t len \t time\n","1/18\tsnR81 \t 201\t 3.5 s\n","2/18\tsnR34 \t 203\t 3.8 s\n","3/18\tsnR43 \t 209\t 3.0 s\n","4/18\tsnR44 \t 211\t 3.3 s\n","5/18\tsnR31 \t 225\t 3.0 s\n","6/18\tsnR10 \t 245\t 3.4 s\n","7/18\tsnR63 \t 255\t 3.1 s\n","8/18\tsnR11 \t 258\t 3.0 s\n","9/18\tsnR82 \t 268\t 3.5 s\n","10/18\tsnR17b \t 332\t 3.1 s\n","11/18\tsnR17a \t 333\t 3.4 s\n","12/18\tsnR37 \t 386\t 3.2 s\n","13/18\tSCR1 \t 522\t 3.3 s\n","14/18\tSRG1 \t 551\t 3.8 s\n","15/18\tsnR19 \t 568\t 3.4 s\n","16/18\tsnR30 \t 606\t 3.8 s\n","17/18\tLSR1 \t 1175\t 5.1 s\n","18/18\tTLC1 \t 1301\t 5.9 s\n"]}],"source":["out_fasta_name = method_name + \"_yeast18\"\n","if os.path.exists(out_fasta_name + \".fasta\"): os.remove(out_fasta_name + \".fasta\")\n","\n","print(\"   \\t lnc \\t len \\t time\")\n","for i, lnc in enumerate(sce18): \n","\n","  start_time = time.time()\n","  seq = sce.loc[lnc][\"sequence\"]\n","  print(f\"{i+1}/{len(sce18)}\\t{lnc} \\t {len(seq)}\", end='\\t')\n","\n","  # Write a one-sequence fasta\n","  with open(\"tmp.fasta\", \"w\") as ofile: \n","    ofile.write(f\">{lnc}\\n{seq}\\n\")\n","  \n","  dot_file_name = run_folding(\"tmp.fasta\")\n","\n","  # Concatenate outputs\n","  os.system(\"cat \" + dot_file_name + \" >> \" + out_fasta_name + \".fasta\") \n","\n","  print(f\"{time.time() - start_time: .1f} s\")"]}],"metadata":{"colab":{"collapsed_sections":[],"name":"rna-state-inf.ipynb","provenance":[{"file_id":"1Pri-xC9r064Y0XfYRGYGkQfAXalNXH6B","timestamp":1643658668660},{"file_id":"1dhERI_Orc5Y4tAv_GaeCNvZRyR99RjUQ","timestamp":1643394834915},{"file_id":"1AsNf4NrlSwDMqmjDiSbx9YrTNeKJYYNo","timestamp":1643288018739},{"file_id":"1GwdNuzwwDx7oru5KT-sB7K8AGtgEO94P","timestamp":1639659140840},{"file_id":"1nLXu1g0Ws9Dhm4-p88DpJDHwsGfnZLWZ","timestamp":1639618378257}]},"kernelspec":{"display_name":"Python 3","name":"python3"},"language_info":{"name":"python"}},"nbformat":4,"nbformat_minor":0}