jenn_scripts.tar

scripts/0000700370132236761610000000000013562321510012662 5ustar  jdawson3babisgrpscripts/pHtempmelt_pKaconst.py0000777370132236761610000000637113033041775017252 0ustar  jdawson3babisgrp#!/usr/bin/env python
#Bootstrap fit for two-state temperature melt data

import os
import sys
import numpy as np
import matplotlib.pyplot as plt
import scipy

from scipy.optimize import curve_fit
from global_fitting import bootstrap

Tref = 293.15 #reference temperature in Kelvin

try:
    Listfile = sys.argv[1]       #Input file name 
    residue = sys.argv[2]        #Specify residue name
    pHval =sys.argv[3]		 #Specify pH
except:
    print "Usage: %s 'Input filename list Residue name'" %sys.argv[0]
    sys.exit()

Directory = ""
OutputFile = "test.txt"
#OutputFile = Directory + residue + pHval + "fit.txt"
#OutputFig = Directory + residue + pHval + "fit.png"
f1 = open(OutputFile,"w")

#loading data
data = np.loadtxt(Listfile,delimiter="\t")
temp = data[0]
w1H = data[2]
w15N = data[1]
w = w1H+0.154*w15N
pKa = 5.09
pKa_err = 0.08
pH = float(pHval) - 0.0028*(temp - 298.15) #temperature dependence value from BioBuffer info, referenced to 25C

iter = 10 #number of fits
iter_boot = 10 #number of bootstrap rounds

y = np.zeros((iter,len(w)))
para_mean = np.zeros((iter,6))
pKai = pKa + pKa_err*np.random.randn(iter)
print pKai
for i in range(0,iter):
    pKa_set = pKai[i]
    def threestate(x,w1,w2,w3,G,S,C):  #Define function: three state pHtitration, then temperature melt 
	K1 = 10**-(pH-pKa_set)
    	DGcalc = G - S*(x-Tref) + C*(x-Tref-x*np.log(x/Tref))
    	K2 = np.exp(-DGcalc/(1.987*10**-3*x))
    	return (w1+w2*K1+w3*K1*K2)/(1+K1+K1*K2)

    #Bootstrap
    p = [w[0],w[0],w[len(w)-1],2.,0.,0.]  #initial guess of fitting parameters
    global_fit = False
    lmfit = False
    para = bootstrap(p,w,temp,threestate,global_fit,lmfit,iterations=iter_boot)
    para_mean[i,:] = np.mean(para,axis=0)
    print "Blah"
    print para_mean[i,:]
    y[i,:] = threestate(temp,*para_mean) 

median_y = np.median(y,axis=0)  
median_para = np.median(para_mean,axis=0)
std_w1 = np.median(abs(para[:,0]-median_para[0]))/0.67449 #median absolute deviation converted to stdev
std_w2 = np.median(abs(para[:,1]-median_para[1]))/0.67449
std_w3 = np.median(abs(para[:,2]-median_para[2]))/0.67449
std_DG0 = np.median(abs(para[:,3]-median_para[3]))/0.67449
std_DS0 = np.median(abs(para[:,4]-median_para[4]))/0.67449
std_DCp = np.median(abs(para[:,5]-median_para[5]))/0.67449

print
print median_para,std_w1,std_w2,std_DG0,std_DS0
print median_DCp, std_DCp

resid = median_y -w
chisq = np.dot(resid,resid)
print chisq

#f1.write("%s %s Two state thermal denaturation fit:\n" %(residue,pHval))
#f1.write("w1 = %8.6f error = %8.6f ppm\nw2 = %8.6f error = %8.6f ppm\n" %(median_w[0],std_w1,median_w[1],std_w2,))
#f1.write("DG0 = %8.6f error = %8.6f kcal/mol\nDS0 = %8.6f error = %8.6f kcal/mol*K\n" %(median_para[0],std_DG0,median_para[1],std_DS0,))
#f1.write("DCp = %8.6f error = %8.6f kcal/mol*K\n" %(median_para[2],std_DCp,))
#f1.write("chi sq = %8.6f\n\n" %(chisq))
#f1.write("\n\n#Data and fit\nTemp (deg C)\tw (ppm)\tw fit (ppm)\n")
#for i in range(0,len(w)):
#    f1.write("%5.2f\t%8.6f\t%8.6f\n" %(temp[i]-273.15,w[i],median_y[i]))

#Plotting data
p1, = plt.plot(temp,w,"ko")
p2, = plt.plot(temp,median_y,"g")
plt.xlabel('Temperature (K)')
plt.ylabel('w weighted (ppm)')
plt.title(residue+" "+pHval+' melt')
plt.legend([p1, p2], ["Data", "Two state fit"],loc='best')
#plt.savefig(OutputFig)
plt.show()


f1.close()
scripts/global_fitting.py0000777370132236761610000001773113031032671016250 0ustar  jdawson3babisgrpimport random

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import leastsq, curve_fit
from lmfit import minimize


""" For debugging """
#np.random.seed(1987)

def model(x,params):
    return None

def gaussian(x,mu,sigma):
    a = 1./(sigma*np.sqrt(2*np.pi))
    return a*np.exp(-np.square(x-mu)/(2.*sigma**2.))

def residual(params,y,x,function,global_fit=False,sigma=None,lmfit=False):
    
    """ Function for calculating residuals for leastsq fit.
        
        Args:

            params (list): your list of parameters in same order as arguments for your function/fitting model e.g. [mu,sigma]
            y: your y values. This can be a list of y datasets if global_fit is True
            x: your x values. This can be a list of x datasets if global_fit is True
            function: your fitting model. Called like so, function(x,*params)
            global_fit: If set to True you can supply the function with a list of x's,y's and functions for global fitting.

        Returns:

            y - model(x,*params)

        Example:

            for a global fit:
                x = [x1,x2,x3]
                y = [y1,y2,y3]
                f = [f1,f2,f3]
                global_fit=True
                pO = [list of starting params] # same order as for *params in fitting function
                result, pcov, infodict, errmsg, success = leastsq(residual,p0,args=(y,x,f,global_fit),full_output=True)


    """

    if sigma is None:

        sigma = [1. for _ in range(len(x))]

    else:

        sigma = sigma


    if global_fit:

        errs = []

        for y_i, x_i, f_i, s_i in zip(y,x,function,sigma):
            if lmfit:
                err = (y_i - f_i(x_i, params)) / s_i
		errs.extend(err)

            else:
                err = (y_i - f_i(x_i, *params)) / s_i
		errs.extend(err)

        errs = np.array(errs).flatten()
        return errs
# This commented code does not work with pandas data Series ...
#        err = [y_i - f_i(x, *params) for y_i, x_i, f_i in zip(y,x,function)]
#        err = np.array(err).flatten()
#        return err

    else:
        if lmfit:
            model = function(x,params)
        else:
            model = function(x,*params)
        return (y-model)/sigma

#def resample(x,y,yerr=None,max_replacement=100):
def resample(x,y,yerr=None,max_replacement=100):
    """ Resample dataset with replacement using numpy.random.choice for bootstrap fit.
        
        Args:
            x: set of x values
            y: set of y values
            yerr: set of y errors
            max_replacement: maximum percentage of data points that can be replaced. Important for small datasets.

        Returns:
            if yerrs:
                resample_x, resample_y
            else:
                resample_x, resample_y, resample_yerr
                
        Example:
            
    """ 

    if len(x) == len(y):

        length = len(x)
        #inds = sorted([random.choice(range(length)) for _ in range(length) ])
        #print inds
        #resample_x = x[inds]
        #resample_y = y[inds]

        # Code to use numpy.random.choice
        data_indices = np.arange(length)
        #choice = np.random.choice(data_indices,int(length*max_replacement/100.))
        #replace = np.random.choice(data_indices,int(length*max_replacement/100.))
        choice = np.random.choice(data_indices,int(length*max_replacement/100.))
        replace = np.random.choice(data_indices,int(length*max_replacement/100.))
        data_indices[choice] = data_indices[replace]
        inds = np.sort(data_indices)

        #inds = np.sort(np.random.choice(data_indices,length))
        print inds
        resample_x = x[inds]
        resample_y = y[inds]
        if yerr is None:
            return resample_x, resample_y

        elif len(yerr) == len(y):
            resample_yerr = yerr[inds]
            return resample_x, resample_y, resample_yerr

        else:
            raise TypeError("yerr should be a list of length y containing the errors of y.")
    else:
        print("len(x) != len(y)!")
        raise TypeError("len(x) should be the same as len(y).")


def bootstrap(params, y, x, function, global_fit=False, lmfit=False, yerr=None, iterations=100, montecarloX=False, mc_kwargs={"std_x":0.1,"iterations":200}, **kwargs):
    if montecarloX:
        from errors import monte_carlox
    else:
        pass
    results = []
    
    for _ in range(iterations):

        if global_fit:
            if yerr is None:    
                data = np.array([resample(x_i,y_i,**kwargs) for x_i, y_i in zip(x,y)])
                x_rs = data[:,0]
                y_rs = data[:,1]
                print x_rs, y_rs
            else:
                data = np.array([resample(x_i,y_i,yerr_i,**kwargs) for x_i, y_i, yerr_i in zip(x,y,yerr)])
                x_rs = data[:,0]
                y_rs = data[:,1]
                yerr_rs = data[:,2]
                print x_rs, y_rs, yerr_rs
             
        else:
            if yerr is None:
                x_rs, y_rs = resample(x,y,**kwargs)
                print x_rs,y_rs
                yerr_rs = None
            else:
                x_rs, y_rs, yerr_rs = resample(x,y,yerr,**kwargs)
                print x_rs,y_rs,yerr_rs
             
#        print data
        #if yerrs is None:
        #    yerrs = [1. for i in y_rs]
        #else
        if montecarloX:
            """ Do montecarlo simulation of x for every y bootstrap """
            std_x = mc_kwargs["std_x"]
            mc_iterations = mc_kwargs["iterations"]
            result = monte_carlox(function, x_rs, params, y_rs, std_x, yerr_rs, global_fit, lmfit, iterations=mc_iterations)
            results.append(result)
        else:
            if lmfit:     
                result = minimize(residual,params,args=(y_rs, x_rs, function, global_fit, yerr_rs, lmfit))
            else:
                result, pcov, infodict, errmsg, success = leastsq(residual, params, args=(y_rs, x_rs, function, global_fit, yerr_rs),full_output=True)
            print result
            results.append(result)
#    print results
    if lmfit:
        return np.array(results).ravel()
    else:
        return np.array(results)
if __name__ == "__main__":

    mu = 10.
    sigma = 1.

    d1 = np.random.normal(mu, sigma, 100)
    d2 = np.random.normal(mu, sigma, 100)
    d3 = np.random.normal(mu, sigma, 100)
    d4 = np.random.normal(mu, sigma, 100)
    d5 = np.random.normal(mu, sigma, 100)
    d6 = np.random.normal(mu, sigma, 100)
    data = np.array([d1,d2,d3,d4,d5,d6])

    hist, bin_edges = np.histogram(d1,20,density=True)
    p0 = [mu,sigma]
    x = bin_edges[:-1] 
    y = hist
    result, pcov, infodict, errmsg, success = leastsq(residual,p0,args=(y,x,gaussian),full_output=True)
    print result
    bs = bootstrap(p0, y, x, gaussian, global_fit=False,iterations=1000)
    av_mu = np.mean(bs[:,0])
    std_mu = np.std(bs[:,0])
    print av_mu,std_mu
    

    """ Global fit """
#    bins = np.min(data),np.max(data)
#    print bins
#    data = np.array([np.histogram(d,20,density=True) for d in data])
#    x = data[:,1]
#    x = [i[:-1] for i in x]
#    y = data[:,0]
#    funcs = [gaussian for _ in x]
#    global_fit=True
#
#    result, pcov, infodict, errmsg, success = leastsq(residual,p0,args=(y,x,funcs,global_fit),full_output=True)
#    print result
#    bs = bootstrap(p0, y, x, funcs, global_fit=True,iterations=100)
#    av_mu = np.mean(bs[:,0])
#    std_mu = np.std(bs[:,0])
#    av_sigma = np.mean(bs[:,1])
#    std_sigma = np.std(bs[:,1])
#    print result
#    print av_mu,std_mu,av_sigma,std_sigma
#
#    c = ["r","g","b","y","orange","pink"]
#    #plt.clf()
#    for i,j,k in zip(x,y,c):
#        plt.bar(i,j,width=0.75,color=k,alpha=0.5)#plt.plot(bins,gaussian(bins,mu,sigma))
#    sim_x = np.linspace(min(x[0]),max(x[0]),100)
#    plt.plot(sim_x,gaussian(sim_x,av_mu,av_sigma),label=r"$ \mu = %.3f \pm %.3f, \sigma = %.3f \pm %.3f$" % (av_mu,std_mu,av_sigma,std_sigma) )
#    #plt.plot(x,y)
#    plt.legend(frameon = False)
#    plt.title("Test global bootstrap fit for Gaussian data.")
#    #print result
#    plt.show()
scripts/pHtitr_2state_refine_indivnucl.py0000777370132236761610000000736213044765034021434 0ustar  jdawson3babisgrp#!/usr/bin/env python
#Bootstrap fit for two-state temperature melt data

import os
import sys
import numpy as np
import matplotlib.pyplot as plt

from global_fitting import bootstrap

try:
    Listfile = sys.argv[1]       #Input file name
    residue = sys.argv[2]        #Specify residue name
    tempval = sys.argv[3]          #Specify pH value
except:
    print "Usage: %s 'Input filename list Residue name'" %sys.argv[0]
    sys.exit()
Directory = residue + "/pHtitr_bytemp/"
#OutputFile = "test.txt"
OutputFile = Directory + residue + "_"+ tempval + "C_refinefit.txt"
OutputFig = Directory + residue + "_"+ tempval + "C_refinefit.png"
f1 = open(OutputFile,"w")

#constants
iter = 1000 #Number of iterations
Tref = 293.15 #reference temperature in Kelvin
pKa = 4.6
pKa_err = 0.3

#loading data
temp=float(tempval)
data_pH = np.loadtxt(Listfile,delimiter="\t")
pH = data_pH[0]-0.0028*(temp-25.) #adjusting for temperature dependence of buffer
w1H = data_pH[2]
w15N = data_pH[1]
#w = w1H+0.154*w15N
xdata = np.concatenate((pH+20.,pH),axis=0)
ydata = np.concatenate((w1H,w15N),axis=0)

def twostate(x,w1h,w1n,w2h,w2n,pKa1):  #Define function: two state temperature melt 
    wcalc = np.zeros(len(x))
    for i in range(0,len(x)):
	if x[i]>20:  #pH titration
    	    K1 = 10**-(x[i]-20.-pKa1)
	    wcalc[i] = (w1h+w2h*K1)/(1+K1)
        else:
    	    K1 = 10**-(x[i]-pKa1)
            wcalc[i] = (w1n+w2n*K1)/(1+K1)
    return wcalc

def twostate_weighted(x,w1,w2,pKa1):
    K1=10**-(x-pKa1)
    return (w1+w2*K1)/(1+K1)

#Bootstrap: 2state model
#w1H,w15N fit
p = [w1H[len(w1H)-1],w1H[0],w15N[len(w15N)-1],w15N[0],5.]  #initial guess of fitting parameters
global_fit = False
lmfit = False
bs = bootstrap(p,ydata,xdata,twostate,global_fit,lmfit,iterations=iter)
para = np.median(bs,axis=0)
y = np.zeros((iter,len(xdata)))
for i in range(0,iter):
    pi = [bs[i,0],bs[i,1],bs[i,2],bs[i,3],bs[i,4]]
    y[i,:] = twostate(xdata,*pi)
median_y = np.median(y,axis=0)
median_w1H = median_y[0:len(pH)]
median_w15N = median_y[len(pH):len(ydata)]
resid = median_y -ydata
#resid = median_y -w
chisq = np.dot(resid,resid)
pKai = pKa + pKa_err*np.random.randn(iter)
K2 = 10**(bs[:,4]-pKai)-1
E2 = np.array([0.])
for i in range(0,iter): 
    if K2[i] > 0: 
    	energy = -1.987*10**-3*293.15*np.log(K2[i]) 
    	E2 = np.vstack((E2,energy))
E2 = np.delete(E2,0,0)
median_energy = np.median(E2)

#Median of fit parameters: 2state
std_w1h = np.median(abs(bs[:,0]-para[0]))/0.67449
std_w1n = np.median(abs(bs[:,1]-para[1]))/0.67449
std_w2h = np.median(abs(bs[:,2]-para[2]))/0.67449
std_w2n = np.median(abs(bs[:,3]-para[3]))/0.67449
std_pKa = np.median(abs(bs[:,4]-para[4]))/0.67449
std_energy = np.median(abs(E2-median_energy))/0.67449
#print "w"
#print para[0], std_w1h, para[1], std_w2h
#print para[2], std_w1n, para[3], std_w2n
#print "pKa"
#print para[4], std_pKa
#print "chisq"
#print chisq
f1.write("Two state pH titration refined fit:\n")
f1.write("pKa = %8.6f error = %8.6f energy = %8.6f error = %8.6f kcal/mol\n" %(para[4],std_pKa,median_energy,std_energy))
f1.write("w1h = %8.6f error = %8.6f ppm   w1n = %8.6f error %8.6f ppm\n" %(para[0],std_w1h,para[1],std_w1n))
f1.write("w2h = %8.6f error = %8.6f ppm   w2n = %8.6f error %8.6f ppm\n" %(para[2],std_w2h,para[3],std_w2n))
f1.write("chi sq = %8.6f\n\n" %(chisq))

#Plotting data
plt.figure(1)
plt.subplot(211)
plt.plot(pH,w1H,"ko")
plt.plot(pH,median_w1H,"k-")
plt.xlabel('pH')
plt.ylabel('w1H (ppm)')
plt.gca().invert_yaxis()
#plt.gca().invert_xaxis()
plt.title(residue+" "+tempval+'C pH titration')
plt.subplot(212)
p1,=plt.plot(pH,w15N,"ko")
p2,=plt.plot(pH,median_w15N,"k-")
plt.xlabel('pH')
plt.ylabel('w15N (ppm)')
plt.gca().invert_yaxis()
#plt.gca().invert_xaxis()
plt.legend([p1, p2], ["Data", "Two state fit"],loc='best')
plt.savefig(OutputFig)
#plt.show()

f1.close()
scripts/tempmelt_bootstrap.py0000777370132236761610000001030513040030103017162 0ustar  jdawson3babisgrp#!/usr/bin/env python
#Bootstrap fit for two-state temperature melt data

import os
import sys
import numpy as np
import matplotlib.pyplot as plt

from global_fitting import bootstrap

iter = 1000 #Number of iterations
Tref = 293.15 #reference temperature in Kelvin
w1_init = 28.6
w2_init = 25.5
DG_init = 0.
DS_init = 0.
DCp_init = 0.

try:
    Listfile = sys.argv[1]       #Input file name
    residue = sys.argv[2]        #Specify residue name
    pHval = sys.argv[3]          #Specify pH value
except:
    print "Usage: %s 'Input filename list Residue name'" %sys.argv[0]
    sys.exit()


def twostate(x,w1,w2,G,S,C):  #Define function: two state temperature melt 
    #DGcalc = G - S*(x-Tref) + C*(x-Tref-x*np.log(x/Tref))
    DGcalc = G - S*(x-Tref) - C*(x-Tref)**2/(2*Tref)
    A = np.exp(-DGcalc/(1.987*10**-3*x))
    return (w1+w2*A)/(1+A)

Directory = residue + "/temp_2state/"
OutputFile = Directory + residue + "_pH"+ pHval + "fit.txt"
OutputFig = Directory + residue + "_pH"+ pHval + "fit.png"
f1 = open(OutputFile,"w")

#loading data
data = np.loadtxt(Listfile,delimiter="\t")
temp = data[0]
w1H = data[2]
w15N = data[1]
w = w1H+0.154*w15N

#Bootstrap: 2state model
p = [w[0],22.6,DG_init,DS_init,DCp_init]  #initial guess of fitting parameters
global_fit = False
lmfit = False
bs = bootstrap(p,w,temp,twostate,global_fit,lmfit,iterations=iter)
#bs = bootstrap(p,y,x,Exp,global_fit,lmfit,yerr=y_err,iterations=iter)


#Median of fit parameters: 2state
#med_w1 = np.median(bs[:,0]) #using median and median absolute deviation because they are less sensitive to outliers (here, bad fits)
#std_w1 = np.median(abs(bs[:,0]-med_w1))/0.67449 #median absolute deviation converted to stdev
#med_w2 = np.median(bs[:,1])
#std_w2 = np.median(abs(bs[:,1]-med_w2))/0.67449
#med_DG0 = np.median(bs[:,2])  
#std_DG0 = np.median(abs(bs[:,2]-med_DG0))/0.67449 
#med_DS0 = np.median(bs[:,3])  
#std_DS0 = np.median(abs(bs[:,3]-med_DS0))/0.67449
#med_DCp = np.median(bs[:,4])  
#std_DCp = np.median(abs(bs[:,4]-med_DCp))/0.67449
#print
#print med_w1, std_w1
#print med_w2, std_w2
#print med_DG0, std_DG0
#print med_DS0, std_DS0
#print med_DCp, std_DCp
#y = np.zeros((len(bs),len(w)))
#for i in range (0,len(bs)):
#   p_i = [bs[i,0],bs[i,1],bs[i,2],bs[i,3], bs[i,4]]
#   y[i,:] = twostate(temp,*p_i)
#median_y = np.median(y,axis=0)
#resid = median_y - w
#chisq = np.dot(resid,resid)


#Filtering endpts: Keeping w1 and w2 within reality
#para = np.array([0.,0.,0.,0.,0.])
#for i in range(0,iter):
#    if (w1_init-4.) <= bs[i,0] <= (w1_init+4.):
#	if (w2_init-4.) <= bs[i,1] <= (w2_init+4.):
#	    a = np.array(bs[i,:])
#	    para = np.vstack([para,a])
#para = np.delete(para,0,0)
para = bs

#Median of fit parameters: 2state
med_w1 = np.median(para[:,0]) #using median and median absolute deviation because they are less sensitive to outliers (here, bad fits)
std_w1 = np.median(abs(para[:,0]-med_w1))/0.67449 #median absolute deviation converted to stdev
med_w2 = np.median(para[:,1])
std_w2 = np.median(abs(para[:,1]-med_w2))/0.67449
med_DG0 = np.median(para[:,2])  
std_DG0 = np.median(abs(para[:,2]-med_DG0))/0.67449 
med_DS0 = np.median(para[:,3])  
std_DS0 = np.median(abs(para[:,3]-med_DS0))/0.67449
med_DCp = np.median(para[:,4])  
std_DCp = np.median(abs(para[:,4]-med_DCp))/0.67449
print
print med_w1, std_w1
print med_w2, std_w2
print med_DG0, std_DG0
print med_DS0, std_DS0
print med_DCp, std_DCp
y = np.zeros((len(para),len(w)))
for i in range (0,len(para)):
   p_i = [bs[i,0],bs[i,1],bs[i,2],bs[i,3],bs[i,4]]
   y[i,:] = twostate(temp,*p_i)
median_y = np.median(y,axis=0)
resid = median_y - w
chisq = np.dot(resid,resid)

f1.write("Two state pH titration fit:\n")
f1.write("w1 = %8.6f error = %8.6f ppm\nw2 = %8.6f error = %8.6f ppm\n" %(med_w1,std_w1,med_w2,std_w2,))
f1.write("DG0 = %8.6f error = %8.6f kcal/mol\nDS0 = %8.6f error = %8.6f kcal/mol*K\n" %(med_DG0,std_DG0,med_DS0,std_DS0,))
f1.write("DCp = %8.6f error = %8.6f kcal/mol*K\n" %(med_DCp,std_DCp,))
f1.write("chi sq = %8.6f\n\n" %(chisq))

#Plotting data
p1, = plt.plot(temp,w,"ko")
p2, = plt.plot(temp,median_y,"g")
plt.xlabel('Temperature (K)')
plt.ylabel('w weighted (ppm)')
plt.title(residue+" "+pHval+' melt')
plt.legend([p1, p2], ["Data", "Two state fit"],loc='best')
plt.savefig(OutputFig)
#plt.show()

f1.close()
scripts/Data_Colorgrd_noiseadj_ThreeColor.py0000777370132236761610000000603713145667174022011 0ustar  jdawson3babisgrp#!/usr/bin/env python
#
#Data_Colorgrd.py: 1. Makes string (for use in pymol selection) specifying the analyzed residues in input file
#              2. Outputs command setting color of each analyzed residue in a color gradient based on data values
# 
#JED June 13,2011

import os
import sys

outputFile = "./AnalyzedResidues_pHdepE_test.pml"

#Color settings=[R,G,B].  Gradient from color1 to color2
color1 = [0.02 , 0.50 , 0.72]   #blue
color2 =  [1.00 , 1.00 , 1.00]   #white
color3 = [1.00 , 0.00 , 0.00]   #red
try:
    DataFile = sys.argv[1]       #Specify input data file
except:
    print "Usage: %s 'Input file'" %sys.argv[0]
    sys.exit()

outputFile = open(outputFile, "w")
data = []
Residue = []
select = "select AnalyzedResidues, resi "
if (os.path.isfile(DataFile)):
    f = open(DataFile, "r")
    for line in (f.readlines()):
	lineStriped = line.strip()
	if (lineStriped[0] != '#' and lineStriped[0].isdigit() != 0):
            line = lineStriped.split()
            name = int(line[0])
            select = select + str(name) + "+"
            data = data[:] + [float(line[1])]
            Residue = Residue[:] + [name]      
    f.close()
select = select.rstrip("+")
outputFile.write("%s\n" %(select))
outputFile.write("select N, name n\n")
outputFile.write("select AnalyzedResidues, AnalyzedResidues&N\n")
outputFile.write("show spheres, AnalyzedResidues\n")
#MaxData = max(data)
MaxData=4.0
#MinData = min(data)
MinData=-4.0

print data[0]
for x in range(len(data)):
    if data[x]<= 0.:
    	Diff = -1*MinData
    	if data[x]<MinData:
            scale=0.00
        else:
            scale = (data[x]-MinData)/Diff
        c1 = color3[0] + (color2[0]-color3[0])*scale
        c2 = color3[1] + (color2[1]-color3[1])*scale
        c3 = color3[2] + (color2[2]-color3[2])*scale
    else:
    	Diff = MaxData
    	if data[x]>MaxData:
            scale=0.00
        else:
            scale = (MaxData-data[x])/Diff
        c1 = color1[0] + (color2[0]-color1[0])*scale
        c2 = color1[1] + (color2[1]-color1[1])*scale
        c3 = color1[2] + (color2[2]-color1[2])*scale
    #outputFile.write("select Res%s, resi %s\n" %(Residue[x],Residue[x]))
    #outputFile.write("select Res%s=Res%s&N\n" %(Residue[x],Residue[x]))
    outputFile.write("set_color Res%s= [%1.3f , %1.3f , %1.3f]\n" %(Residue[x],c1,c2,c3))
    outputFile.write("color Res%s, resi %s\n" %(Residue[x],Residue[x]))

outputFile.close()

#[R,G,B] values for different colors in pymol (From pymol manual)
#optimized rgb values for cmyk output:
#dblue= [0.05 , 0.19 , 0.57]
#blue=  [0.02 , 0.50 , 0.72]
#mblue= [0.5  , 0.7  , 0.9 ]
#lblue= [0.86 , 1.00 , 1.00]
 
#green= [0.00 , 0.53 , 0.22]
#lgreen=[0.50 , 0.78 , 0.50]
#yellow=[0.95 , 0.78 , 0.00]
#orange=[1.00 , 0.40 , 0.0 ]
 
# these are trivial
#red=   [1.00 , 0.00 , 0.00]
#mred=  [1.00 , 0.40 , 0.40]
#lred=  [1.00 , 0.80 , 0.80]
#vlred= [1.00 , 0.90 , 0.90]
#white= [1.00 , 1.00 , 1.00]
#vlgray=[0.95 , 0.95 , 0.95]
#lgray= [0.90 , 0.90 , 0.90]
#gray=  [0.70 , 0.70 , 0.70]
#dgray= [0.50 , 0.50 , 0.50]
#dgray=[0.30 , 0.30 , 0.30]
#black= [0.00 , 0.00 , 0.00]
##
scripts/Data_Colorgrd_noiseadj_TwoColor.py0000777370132236761610000000531713045232145021474 0ustar  jdawson3babisgrp#!/usr/bin/env python
#
#Data_Colorgrd.py: 1. Makes string (for use in pymol selection) specifying the analyzed residues in input file
#              2. Outputs command setting color of each analyzed residue in a color gradient based on data values
# 
#JED June 13,2011

import os
import sys

outputFile = "./AnalyzedResidues_triple.pml"

#Color settings=[R,G,B].  Gradient from color1 to color2
color1 = [1.00 , 1.00 , 1.00]   #white
color2 =  [0.02 , 0.50 , 0.72]   #blue

try:
    DataFile = sys.argv[1]       #Specify input data file
except:
    print "Usage: %s 'Input file'" %sys.argv[0]
    sys.exit()

outputFile = open(outputFile, "w")
data = []
Residue = []
select = "select AnalyzedResidues, resi "
if (os.path.isfile(DataFile)):
    f = open(DataFile, "r")
    for line in (f.readlines()):
	lineStriped = line.strip()
	if (lineStriped[0] != '#' and lineStriped[0].isdigit() != 0):
            line = lineStriped.split()
            name = int(line[0])
            select = select + str(name) + "+"
            data = data[:] + [float(line[1])]
            Residue = Residue[:] + [name]      
    f.close()
select = select.rstrip("+")
outputFile.write("%s\n" %(select))
outputFile.write("select N, name n\n")
outputFile.write("select AnalyzedResidues, AnalyzedResidues&N\n")
outputFile.write("show spheres, AnalyzedResidues\n")
#MaxData = max(data)
MaxData=0.03
#MinData = min(data)
MinData=0.01
Diff = MaxData-MinData
print data[0]
for x in range(len(data)):
    if data[x]>MaxData:
        scale=1.00
    elif data[x]<MinData:
        scale=0.00
    else:
        scale = (data[x]-MinData)/Diff
    c1 = color1[0] + (color2[0]-color1[0])*scale
    c2 = color1[1] + (color2[1]-color1[1])*scale
    c3 = color1[2] + (color2[2]-color1[2])*scale
    #outputFile.write("select Res%s, resi %s\n" %(Residue[x],Residue[x]))
    #outputFile.write("select Res%s=Res%s&N\n" %(Residue[x],Residue[x]))
    outputFile.write("set_color Res%s= [%1.3f , %1.3f , %1.3f]\n" %(Residue[x],c1,c2,c3))
    outputFile.write("color Res%s, resi %s\n" %(Residue[x],Residue[x]))

outputFile.close()

#[R,G,B] values for different colors in pymol (From pymol manual)
#optimized rgb values for cmyk output:
#dblue= [0.05 , 0.19 , 0.57]
#blue=  [0.02 , 0.50 , 0.72]
#mblue= [0.5  , 0.7  , 0.9 ]
#lblue= [0.86 , 1.00 , 1.00]
 
#green= [0.00 , 0.53 , 0.22]
#lgreen=[0.50 , 0.78 , 0.50]
#yellow=[0.95 , 0.78 , 0.00]
#orange=[1.00 , 0.40 , 0.0 ]
 
# these are trivial
#red=   [1.00 , 0.00 , 0.00]
#mred=  [1.00 , 0.40 , 0.40]
#lred=  [1.00 , 0.80 , 0.80]
#vlred= [1.00 , 0.90 , 0.90]
#white= [1.00 , 1.00 , 1.00]
#vlgray=[0.95 , 0.95 , 0.95]
#lgray= [0.90 , 0.90 , 0.90]
#gray=  [0.70 , 0.70 , 0.70]
#dgray= [0.50 , 0.50 , 0.50]
#dgray=[0.30 , 0.30 , 0.30]
#black= [0.00 , 0.00 , 0.00]
##