#! /usr/bin/env python

# rnaSecStrNussJacAndVienna.py
# P.Clote
# This code computes the Nussinov Jacobson as well as Vienna
# optimal secondary structure of an
# RNA sequence together with the list of base pairing positions.

# Functions defined in this program:
#  computeViennaSecStr(rna)
#  viennaPartitionFun(rna)
#  evalViennaEnergy(rna,secStr)
#  computeViennaSecStrAndMFE(rna)

#WARNING: I added -d2 to Vienna
import sys,os,tempfile,math

from misc import DEBUG, GUbond, AUbond, CGbond, INF
from misc import THRESHOLD, error, guPair, cgPair, auPair, basePair
from misc import watsonCrick, a, list2string, basePairList
from misc import RNAeval,RNAfold,RNAsubopt,NORMALIZE
from misc import RgasConstant, TEMP, _RT_

def runRNAsubopt(rna,energy=5,printFlag=0):
  #printFlag is boolean flag whether output of RNAsubopt should be printed
  tmpfileName = tempfile.mktemp()
  tmpfile     = open(tmpfileName,"w")
  tmpfile.write(rna)
  tmpfile.close()
  cmd  = "%s -e %f < %s" % (RNAsubopt,energy,tmpfileName)
  file = os.popen(cmd)
  L    = []  #L is list of suboptimal secStr
  line = file.readline()
  if printFlag: print line[:-1]
  rna0 = line.split()[0]
  assert( rna == rna0 )
  line = file.readline()
  while line:
    if printFlag: print line[:-1]
    secStr = line.split()[0]
    L.append(secStr) 
    line = file.readline()
  return L

def computeViennaSecStr(seq):
  tmpfileName = tempfile.mktemp()
  tmpfile     = open(tmpfileName,"w")
  tmpfile.write(seq)
  tmpfile.close()
  cmd  = "%s  < %s" % (RNAfold,tmpfileName)
  file = os.popen(cmd)
  seq0 = file.readline().strip()
  assert( seq == seq0 )
  secStr = file.readline()
  file.close()
  positionOfBlank = secStr.find(" ")
  secStr = secStr[:positionOfBlank]
  return secStr

def viennaPartitionFun(rna):
  tmpfileName = tempfile.mktemp()
  tmpfile     = open(tmpfileName,"w")
  tmpfile.write(rna)
  tmpfile.close()
  cmd  = "%s  -p0 < %s" % (RNAfold,tmpfileName)
  file = os.popen(cmd)
  file.readline().strip()  #first line is the RNA which we already have
  secStr = file.readline()
  positionOfBlank = secStr.find(" ")
  secStr = secStr[:positionOfBlank]
  freeEnergy = float(file.readline().split()[-2])
  Z          = math.exp(-freeEnergy/_RT_)
  mfeFreqInEns = float(file.readline().split()[-1])
  file.close()
  return (secStr,freeEnergy,Z,mfeFreqInEns)

def evalViennaEnergy(rna,secStr,T=37):
  #compute the minimum free energy using RNAeval
#  tmpFileName = tempfile.mktemp()
  tmpFileName = 'tmp'
  tmpFile     = open(tmpFileName,"w")
  tmpFile.write("%s\n" % rna)
  tmpFile.write("%s\n" % secStr)
  tmpFile.close()
  cmd = "%s -T %d  < %s" % (RNAeval,T,tmpFileName)
  tmpOutFile  = os.popen(cmd)
  line        = tmpOutFile.readline()  
                #discard first line, which is RNA sequence
  line        = tmpOutFile.readline()
  words       = line.split()
  energy      = words[-1]
  indexLeftParen = energy.find("(")
  energy   = energy[indexLeftParen+1:]
  indexRightParen = energy.find(")")
  if indexRightParen >= 0:
    energy   = energy[:indexRightParen]
  energy   = float(energy)
  #print energy
  if NORMALIZE: energy = energy/len(rna)
  tmpOutFile.close()
  return energy

def computeViennaSecStrAndMFE(rna):
  #compute the minimum free energy using RNAfold
  tmpFileName = tempfile.mktemp()
  tmpFile     = open(tmpFileName,"w")
  tmpFile.write("%s\n" % rna)
  tmpFile.close()
  cmd = "%s   < %s" % (RNAfold,tmpFileName)
  tmpOutFile  = os.popen(cmd)
  line        = tmpOutFile.readline()  
                #discard first line, which is RNA sequence
  line        = tmpOutFile.readline()
  words       = line.split()
  energy      = words[-1]
  secStr      = words[0]
  indexLeftParen = energy.find("(")
  energy   = energy[indexLeftParen+1:]
  indexRightParen = energy.find(")")
  if indexRightParen >= 0:
    energy   = energy[:indexRightParen]
  energy   = float(energy)
  if NORMALIZE: energy = energy/len(rna)
  print energy
  tmpOutFile.close()
  return secStr,energy