Mercurial > repos > john-mccallum > pcr_markers
comparison umelt_service.py @ 5:b321e0517be3 draft
Uploaded
| author | ben-warren |
|---|---|
| date | Thu, 22 May 2014 20:30:19 -0400 |
| parents | |
| children | f201e8c6e004 |
comparison
equal
deleted
inserted
replaced
| 4:be070a68521e | 5:b321e0517be3 |
|---|---|
| 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: | |
| 25 print 'The server couldn\'t fulfill the request.' | |
| 26 print 'Error code: ', e.code | |
| 27 except urllib2.URLError, e: | |
| 28 print 'We failed to reach a server.' | |
| 29 print 'Reason: ', e.reason | |
| 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 |