|
5
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 ##interrogate umelt service at Univ of Utah
|
|
|
4 import sys
|
|
|
5 import urllib
|
|
|
6 import urllib2
|
|
|
7 import xml.etree.ElementTree as ET
|
|
|
8 import numpy as np
|
|
|
9 import scipy as sp
|
|
|
10 from scipy import interpolate
|
|
|
11 from scipy.interpolate import interp1d
|
|
|
12
|
|
|
13
|
|
|
14 url='https://www.dna.utah.edu/db/services/cgi-bin/udesign.cgi'
|
|
|
15 timeout_sec=500 ## default for timeout
|
|
|
16
|
|
|
17 # Query UW melt prediction service for a single sequence, returning default array for helicity, assuming within temperature range of 65-95
|
|
|
18 def getmelt(input_seq):
|
|
|
19 values = {'seq' : input_seq, 'rs':0, 'dmso':0,'cation': 20 ,'mg': 2} # Note the buffer conditions
|
|
|
20 data = urllib.urlencode(values)
|
|
|
21 req = urllib2.Request(url, data)
|
|
|
22 try:
|
|
|
23 response = urllib2.urlopen(req,timeout=timeout_sec)
|
|
|
24 except urllib2.HTTPError, e:
|
|
6
|
25 print >> sys.stderr, 'The server couldn\'t fulfill the request.'
|
|
|
26 print >> sys.stderr, 'Error code: ', e.code
|
|
5
|
27 except urllib2.URLError, e:
|
|
6
|
28 print >> sys.stderr, 'We failed to reach a server.'
|
|
|
29 print >> sys.stderr, 'Reason: ', e.reason
|
|
5
|
30 else:
|
|
|
31 melt_data = response.read()
|
|
|
32 tree = ET.fromstring(melt_data)
|
|
|
33 helicity = [amp.find('helicity').text.split() for amp in tree.findall('amplicon')]
|
|
|
34 hels = np.array(helicity[0], dtype=np.float32).transpose()
|
|
|
35 return hels
|
|
|
36
|
|
|
37 # helicity[0] used because the default retreived data is a list of 3 lists (for WT, mu and hets). The 3 lists are identical for our data (no IUPAC) so only [0] is used.
|
|
|
38
|
|
|
39 def getTm(hel_array):
|
|
|
40 temps=np.arange(65,100.5,0.5) # Temperature range of 65-100.5. Step of 0.5 NEEDS TO BE CHECKED..REcent change?
|
|
|
41 tck = interpolate.splrep(temps,hel_array,s=0)
|
|
|
42 xnew =np.arange(65,95.5,0.05)
|
|
|
43 yder = interpolate.splev(xnew,tck,der=1) # der=1, first derivative
|
|
|
44 return xnew[yder.argmin()] # Returns the x value corresponding to the minimum (peak) y value -> Tm
|
