diff --git a/bin/hdxworkbench.py b/bin/hdxworkbench.py new file mode 100755 index 0000000..a55f407 --- /dev/null +++ b/bin/hdxworkbench.py @@ -0,0 +1,196 @@ +import argparse + +from scipy.stats import sem + +import analysis +import hdx_models +import input_data +import plots +import sampling +import system_setup + +inseq = "" +benchmark = False +# Arguments: HDXWorkbench file, output_dir, run_type, inseq +parser = argparse.ArgumentParser(description='Analysis of HDXWorkbench DHDX') +parser.add_argument('-w', help='HDXWorkbench CSV file', required=True) +parser.add_argument('--benchmark', help='Run a short benchmark instead of full sampling run', + action="store_true") +parser.add_argument('--inseq', metavar='s', type=str, + help='Input sequence string', required=False) +parser.add_argument('--nsteps', help='Equilibrium steps. 5000 to 10000. ' + 'Default: 1000', + default=1000, + type=int, + required=False) +parser.add_argument('--init', help='How to initialize - either "random" or "enumerate". ' + 'Enumerate is slower but sampling will converge faster. ' + 'Default: enumerate', + default="enumerate", + required=False) +parser.add_argument('--num_exp_bins', help='Number of log(kex) values for sampling. 20 is generally sufficient. ' + 'Default: 20', + default=20, + type=int, + required=False) +parser.add_argument('--mol_name', help='Molecule name', type=str, + required=True) +parser.add_argument('-o', '--outputdir', help='Output directory.', required=True) +parser.add_argument('--annealing_steps', help='Steps per temperature in annealing - 100-200 sufficient. ' + 'Default: 20', + default=20, + type=int, + required=False) +parser.add_argument('--sigma0', help='Estimate for experimental error in %%D Units. ' + 'Default: 5', + default=5, + type=float, + required=False) +parser.add_argument('--saturation', help='Deuterium saturation in experiment. ' + 'Default: 1.0', + default=1.0, + type=float, + required=False) +parser.add_argument('--offset', help='Offset between fragment start/end values and FASTA sequence. ' + 'Default: 0', + default=0, + type=float, + required=False) + +args = parser.parse_args() + +if args.inseq is not None: + inseq = args.inseq + +#################### +### What kind of run do you want to do? + +if args.benchmark: + run_type = "benchmark" # Use this to run a quick simulation to see how long a full sampling run will tak +else: + run_type = "sampling" # Use this to do the full sampling/analysis + +########################################## +### File/Directory Setup +### +### output directory for this simulation. +outputdir = args.outputdir +### +### HDX data file location +workbench_file = args.w + +f = open(workbench_file, "r") +for line in f.readlines(): + if line.split(",")[0] == "Offset": + offset = int(line.split(",")[1].strip()) + if line.split(",")[0] == "Experiment Protein Sequence": + inseq = line.split(",")[1] + if line.split(",")[0] == "Deuterium solution concentration": + saturation = float(line.split(",")[1].strip()) + if line.split(",")[0] == "Experiment name": + name = line.split(",")[1].strip().replace(" ", "_") + +if inseq == "": + raise Exception( + "HDX Workbench file does not contain FASTA sequence. Please manually add the sequence to the command line using the flag -s or --inseq") + +########################################### +### Simulation Parameters + + +# User-controlled variables +if args.benchmark: + benchmark = True + +num_exp_bins = args.num_exp_bins # Number of log(kex) values for sampling. 20 is generally sufficient. +init = args.init # How to initialize - either "random" or "enumerate". Enumerate is slower but sampling will converge faster +nsteps = args.nsteps # equilibrium steps. 5000 to 10000 +num_best_models = 1000000 # Number of best models to consider for analysis +outputdir = args.outputdir + +# Non user controlled vbl - for now. +offset = args.offset +annealing_steps = args.annealing_steps # steps per temperature in annealing - 100-200 sufficient +sigma0 = args.sigma0 # Estimate for experimental error in %D Units +saturation = args.saturation # Deuterium saturation in experiment +percentD = True # Is the data in percent D (True) or Deuterium units? - Always percentD for Workbench. +############################### +### System Setup: +############################### + + +# Initialize model (name, FASTA sequence, offset) +model = system_setup.HDXModel(args.mol_name, + inseq, + offset=offset) + +# Add data to model (model, filename) +input_data.HDXWorkbench(model, workbench_file) + +# Initialize a sampling model for each state (Multiexponential in this case) +for state in model.states: + hdxm = hdx_models.MultiExponentialModel(model=model, + state=state, + sigma=sigma0, + init=init) +############################### +### Sampling: +### + +# If benchmark is set to true, run a short simulation to estimate runtime +if run_type == "benchmark": + sampling.benchmark(model, sample_sigma=True) + exit() + +# Simulated Annealing macro runs high temperature dynamics and relaxes +# to low temperature, followed by an equilibration run of "nsteps" + +if run_type == "sampling": + sampling.simulated_annealing(model, sigma=sigma0, equil_steps=nsteps, sample_sigma=True, + annealing_steps=annealing_steps, + outfile_prefix="Macro", outdir=outputdir) + +bsm, scores = analysis.get_best_scoring_models(model.states[0].modelfile, model.states[0].scorefile, + num_best_models=num_best_models, + prefix=model.states[0].state_name, write_file=False) + +num_best_models_list = [100, 100] +minstderr2 = 1.0E+34 +for i in range(100, len(scores)): + subset = scores[0:i] + stderr = sem(subset) + if stderr < minstderr2: + minstderr2 = stderr + num_best_models_list[0] = i +print('No. model: {} stderr: {}'.format(num_best_models_list[0], minstderr2)) + +bsm, scores = analysis.get_best_scoring_models(model.states[1].modelfile, model.states[1].scorefile, + num_best_models=num_best_models, + prefix=model.states[1].state_name, write_file=False) + +minstderr2 = 1.0E+34 +for i in range(100, len(scores)): + subset = scores[0:i] + stderr = sem(subset) + if stderr < minstderr2: + minstderr2 = stderr + num_best_models_list[1] = i +print('No. model: {} stderr: {}'.format(num_best_models_list[1], minstderr2)) + +############################### +### Analysis: + +plots.plot_2state_fragment_avg_model_fits_num_model_per_state(model.states[0], model.states[1], + sig=5.0, + num_best_models=num_best_models_list, + write_file=False, + outdir=outputdir, + show_plot=False) + +plots.plot_apo_lig_dhdx(model, show_plot=False, save_plot=True, + outfile="dhdx.png", + outdir=outputdir, + noclobber=False) + +plots.plot_fragment_chi_values(model.states[0], sig="model", outdir=outputdir, show_plot=False) +plots.plot_fragment_chi_values(model.states[1], sig="model", outdir=outputdir, show_plot=False) diff --git a/pyext/src/analysis.py b/pyext/src/analysis.py index e19a0ab..5cb74e9 100644 --- a/pyext/src/analysis.py +++ b/pyext/src/analysis.py @@ -2,15 +2,21 @@ Analysis functions for HDX simulations """ from __future__ import print_function -import numpy + import os.path +from matplotlib import * +from pylab import * + +import hxio + class Convergence(object): ''' Given two ParseOutputFile classes, allow testing of various convergence and clustering metrics ''' + def __init__(self, parse_output1, parse_output2): self.sample1 = parse_output1 self.sample2 = parse_output2 @@ -45,9 +51,7 @@ def parse_header(self): # #-symbol meas datasets elif line[0:2] == "# ": - self.datafiles.append( - (line[2:].split("|")[0].strip(), - float(line[2:].split("|")[1].strip()))) + self.datafiles.append((line[2:].split("|")[0].strip(), float(line[2:].split("|")[1].strip()))) # @-symbol means sectors elif line[0:2] == "@ ": @@ -81,11 +85,11 @@ def parse_header(self): def get_datasets(self): if len(self.datasets) == 0: self.generate_datasets() - return datasets + return self.datasets def generate_datasets(self): for f in self.datafiles: - datasets.append(hxio.import_json(f)) + self.datasets.append(hxio.import_json(f)) def get_best_scoring_models(self, N, sigmas=False, return_pf=False): ''' Get the N best scoring models from the output file @@ -104,8 +108,7 @@ def get_best_scoring_models(self, N, sigmas=False, return_pf=False): score = float(line.split("|")[1].strip()) # if the score is better than the last best score - if (len(best_scoring_models) < N - or score < best_scoring_models[-1][0]): + if len(best_scoring_models) < N or score < best_scoring_models[-1][0]: if len(best_scoring_models) >= N: del best_scoring_models[-1] model_string = line[1:].split("|")[0].strip() @@ -113,11 +116,9 @@ def get_best_scoring_models(self, N, sigmas=False, return_pf=False): for m in model_string.split(" "): model_list.append(int(m)) if return_pf: - model_list = self.models_to_protection_factors( - model_list) + model_list = self.models_to_protection_factors(model_list) best_scoring_models.append((score, model_list)) - best_scoring_models = sorted( - best_scoring_models, key=lambda x: x[0]) + best_scoring_models = sorted(best_scoring_models, key=lambda x: x[0]) self.best_scoring_models = best_scoring_models return best_scoring_models @@ -134,7 +135,7 @@ def models_to_protection_factors(self, models): for res in range(len(m)): if res + 1 in self.observed_residues: # print(res+1, m[res-1], self.pf_grids[res+1]) - pf_model.append(float(self.pf_grids[res+1][m[res-1]-1])) + pf_model.append(float(self.pf_grids[res + 1][m[res - 1] - 1])) else: pf_model.append(numpy.nan) @@ -148,6 +149,7 @@ class OutputAnalysis(object): Class that analyzes an output file and produces standard graphs ''' + def __init__(self, output): self.output = output @@ -155,10 +157,8 @@ def parse_output_files(self): pass -def get_best_scoring_models(modelfile, scorefile, num_best_models=100, - prefix=None, write_file=True): - # takes a model and score file and writes a new model file with the - # best X scoring models. +def get_best_scoring_models(modelfile, scorefile, num_best_models=100, prefix=None, write_file=True): + # takes a model and score file and writes a new model file with the best X scoring models. # This new file can then be imported into an HDXModel class for analysis i = 0 if prefix is None: @@ -187,9 +187,10 @@ def get_best_scoring_models(modelfile, scorefile, num_best_models=100, for i in range(num_best_models): top_score_indices.append(scores.index(min(scores))) top_scores.append(min(scores)) - scores[scores.index(min(scores))] = max(scores)+1 + # print(scores, min(scores), scores.index(min(scores)), top_score_indices) + scores[scores.index(min(scores))] = max(scores) + 1 if write_file: - _ = open(outfile, "w") + output_file = open(outfile, "w") return top_models, top_scores else: for i in top_score_indices: @@ -207,7 +208,7 @@ def sector_sort(sectors): if s.start_res > last_first_res: last_first_res = s.start_res - for n in range(last_first_res+1): + for n in range(last_first_res + 1): for s in sectors: if s.start_res == n: sorted_sectors.append(s) @@ -225,8 +226,7 @@ def array_frequency(a): return numpy.vstack((unique, count)).T -def get_residue_rate_probabilities(modelfile, scorefile, sectors, seq, grid, - num_models=5, +def get_residue_rate_probabilities(modelfile, scorefile, sectors, seq, grid, num_models=5, outfile="rate_probabilities.dat", offset=0): # Given a set of models (from a best_models.dat file) # returns the probability of observing each rate @@ -238,8 +238,7 @@ def get_residue_rate_probabilities(modelfile, scorefile, sectors, seq, grid, if hasattr(grid, '__iter__'): grid = len(grid) - best_models, best_scores = get_best_scoring_models( - modelfile, scorefile, num_models, write_file=False) + best_models, best_scores = get_best_scoring_models(modelfile, scorefile, num_models, write_file=False) best_models = numpy.array(best_models) @@ -248,35 +247,34 @@ def get_residue_rate_probabilities(modelfile, scorefile, sectors, seq, grid, # Loop over all sectors and add up instances of each rate bin for s in sorted_sectors: # get all instances of each rate - freq = array_frequency(best_models[:, s.start_res:s.end_res+1]) + freq = array_frequency(best_models[:, s.start_res:s.end_res + 1]) model = numpy.zeros(grid) for i in freq: - model[i[0]] = 1.0*i[1]/s.num_amides/num_models + model[i[0]] = 1.0 * i[1] / s.num_amides / num_models - for n in range(s.start_res, s.end_res+1): - if seq[n+offset] == "P": + for n in range(s.start_res, s.end_res + 1): + if seq[n + offset] == "P": of.write(n, "P", numpy.zeros(len(grid), numpy.int)) else: - of.write(str(n+1)+", "+seq[n+offset]+", " + - str([m for m in model])+"\n") + of.write(str(n + 1) + ", " + seq[n + offset] + ", " + str([m for m in model]) + "\n") def get_convergence(state, num_points=None): """ Takes all models in the exp_model of the given state and - divides them into two halves. The average and SD for each residue - is computed for each ensemble and compared via a Two-sample - Kolmogorov-Smirnov Test""" + divides them into two halves. The average and SD for each residue is computed + for each ensemble and compared via a Two-sample Kolmogorov-Smirnov Test""" these_states = model.states for state in these_states: es = state.exp_model - # sectors = state.sectors - if num_points is None or num_points > len(es.exp_models)/2: - num_points = len(es.exp_models)/2 - exp_model1 = es.exp_models[ - len(es.exp_models)/2-num_points:len(es.exp_models)/2] - exp_model2 = es.exp_models[len(es.exp_models)-num_points:-1] + sectors = state.sectors + if num_points is None or num_points > len(es.exp_models) / 2: + num_points = len(es.exp_models) / 2 + exp_model1 = es.exp_models[len(es.exp_models) / 2 - num_points:len(es.exp_models) / 2] + exp_model2 = es.exp_models[len(es.exp_models) - num_points:-1] zscore = calculate_zscore(exp_model1, exp_model2, state, state) - print(zscore) + # print(state.state_name) + # print zscore + # time.sleep(1) def get_average_sector_values(exp_models, state): @@ -291,48 +289,53 @@ def get_cdf(exp_models): equivalent to the empirical density function for each residue. """ exp_model_edf = numpy.empty((len(exp_models), len(exp_models[0]))) A = numpy.array(exp_models) - # y = numpy.linspace(1./len(exp_models), 1, len(exp_models)) + y = numpy.linspace(1. / len(exp_models), 1, len(exp_models)) print(len(exp_models[0])) for i in range(len(exp_models[0])): - counts, edges = numpy.histogram( - A[:, i], len(A), range=(-6, 0), density=False) + counts, edges = numpy.histogram(A[:, i], len(A), range=(-6, 0), density=False) # print i,A[:,i],counts,numpy.cumsum(counts*1.0/len(A)) - exp_model_edf[:, i] = numpy.cumsum(counts*1.0/len(A)) + exp_model_edf[:, i] = numpy.cumsum(counts * 1.0 / len(A)) return exp_model_edf def get_chisq(exp_models1, exp_models2, nbins): - """Takes two lists of exp_models and returns the chi2 value along - the second axis""" + """ Takes two lists of exp_models and returns the chi2 value along the second axis """ A = numpy.array(exp_models1) B = numpy.array(exp_models2) + # y=numpy.linspace(1./len(exp_models1),1,len(exp_models1)) print(len(exp_models1[0])) for i in range(269, len(exp_models1[0])): meanA = numpy.mean(A[:, i]) - ssd = numpy.std(A[:, i])**2 + numpy.std(B[:, i])**2 + ssd = numpy.std(A[:, i]) ** 2 + numpy.std(B[:, i]) ** 2 sstdev = numpy.sqrt(ssd / 5000) meanB = numpy.mean(B[:, i]) t = 1.96 ci = t * sstdev dm = meanA - meanB - print(i, dm, ci, dm/ci) - # return exp_model_edf + print(i, dm, ci, dm / ci) + # fig=plt.figure() + # ax1 = fig.add_subplot(111) + # ax1.plot(range(nbins), countsA) + # ax1.plot(range(nbins), countsB) + # plt.show() + + # data = [countsA, countsB] + # print(i, chi2_contingency(data)) + return exp_model_edf def calculate_ks_statistic(edf1, edf2): """ Takes a two edfs and returns a vector of the Kolmogorov-Smirnov statistic for each residue""" maxdiff = numpy.zeros(len(edf1[0])) - threshold = 1.98*numpy.sqrt( - 1.0*(len(edf1)+len(edf2))/(1.0*len(edf1)*len(edf2))) + threshold = 1.98 * numpy.sqrt(1.0 * (len(edf1) + len(edf2)) / (1.0 * len(edf1) * len(edf2))) if len(edf1[0]) != len(edf2[0]): - print( - "Different Number of Residues for EDFs in KS calculation: Exiting") + print("Different Number of Residues for EDFs in KS calculation: Exiting") exit() for r in range(len(edf1[0])): maxdiff[r] = 0 for m in range(len(edf1[:, 0])): - diff = abs(edf1[m, r]-edf2[m, r]) + diff = abs(edf1[m, r] - edf2[m, r]) if diff > maxdiff[r]: maxdiff[r] = diff return maxdiff, threshold @@ -341,18 +344,15 @@ def calculate_ks_statistic(edf1, edf2): def get_sector_averaged_models(exp_models, state): sector_avg_models = [] for n in range(len(exp_models)): - sector_avg_models.append( - state.exp_model.get_sector_averaged_protection_values( - exp_models[n], state.sectors)) + sector_avg_models.append(state.exp_model.get_sector_averaged_protection_values(exp_models[n], state.sectors)) return sector_avg_models def calculate_zscore(exp_models1, exp_models2, state1, state2): avg1, sd1 = get_average_sector_values(exp_models1, state1) avg2, sd2 = get_average_sector_values(exp_models2, state2) - zscore = numpy.subtract(avg1, avg2) / \ - numpy.sqrt(numpy.add(numpy.add(numpy.square(sd1), numpy.square(sd2)), - 0.00001)) + zscore = numpy.subtract(avg1, avg2) / numpy.sqrt( + numpy.add(numpy.add(numpy.square(sd1), numpy.square(sd2)), 0.00001)) return zscore diff --git a/pyext/src/hdx_models.py b/pyext/src/hdx_models.py index 29a3777..5b52a18 100644 --- a/pyext/src/hdx_models.py +++ b/pyext/src/hdx_models.py @@ -2,209 +2,208 @@ Classes to initialize various models and calculate likelihoods """ from __future__ import print_function +#import IMP +#import IMP.hdx +import hdx_models +import analysis +import input_data import sampling import numpy import scipy import scipy.special +from numpy import linalg import sys from copy import deepcopy +import os +import os.path +import time from random import randint - class MultiExponentialModel(object): - """ A multi-exponential model for exchange that approximates each - independent sector as a sum of exponentials, corresponding to the - number of exchanging amides in that sector. - The exchange rate constant, kex, for each residue is sampled along - a grid of log(kex) values. - - The model can be initialized to a certain bin in the exponential grid, - or to random values. - The Bayesian likelihood of the model is calculated against a set - of hdx.system_setup.Fragment objects which contain experimental data. + """ A multi-exponential model for exchange that approximates each independent sector + as a sum of exponentials, corresponding to the number of exchanging amides in that sector. + The exchange rate constant, kex, for each residue is sampled along a grid of log(kex) values. + + The model can be initialized to a certain bin in the exponential grid, or to random values. + The Bayesian likelihood of the model is calculated against a set of hdx.system_setup.Fragment objects + which contain experimental data. """ - def __init__(self, model, state, frags=None, sigma=10.0, num_exp_bins=10, - init="enumerate", error_model="gaussian", - output_directory=None, noclobber=True, first_tp=10, - last_tp=3600, marginalize=False): + def __init__(self, model, state, frags=None, sigma=10.0, num_exp_bins=10, init="enumerate", + error_model="gaussian", output_directory=None, noclobber=True, first_tp=10, last_tp=3600, marginalize=False): if frags is None: - self.frags = state.frags + self.frags=state.frags else: - self.frags = frags - # print "using frags" - tpoints = [] - times = [] + self.frags=frags + #print "using frags" + tpoints=[] + times=[] # get first and last time points for f in state.frags: tpoints += [tp for tp in f.timepoints] - times += [t.time for t in tpoints] + times+=[t.time for t in tpoints] first_tp = numpy.min(times) last_tp = numpy.max(times) self.first_tp = first_tp self.last_tp = last_tp - self.state = state - self.state.exp_model = self - self.exp_grid = self.calc_exp_grid(num_exp_bins, first_tp, last_tp) - self.exp_seq = [0]*len(model.seq) - self.name = state.state_name + self.state=state + self.state.exp_model=self + self.exp_grid=self.calc_exp_grid(num_exp_bins, first_tp, last_tp) + self.exp_seq=[0]*len(model.seq) + self.name=state.state_name state.add_exp_model(self) - self.offset = model.offset - self.sigma = sigma - self.oldsigma = -1 - self.scores = [] - self.seq = model.seq - self.flag_outliers = True - self.output_directory = output_directory - self.noclobber = noclobber - self.marginalize = marginalize + self.offset=model.offset + self.sigma=sigma + self.oldsigma=-1 + self.scores=[] + self.seq=model.seq + self.flag_outliers=True + self.output_directory=output_directory + self.noclobber=noclobber + self.marginalize=marginalize # These grids contain the values at which the integral of sigma is # computed - # stdev between 0.01 and 0.5 - self.sigma0_grid = numpy.linspace(0.1, 10, 100) - # stdev between 0.01 and 0.5 - self.omega_grid = numpy.linspace(1, 100, 100) - # fmod between 0.00 and 1.00 - self.fmod_grid = numpy.linspace(0, 100, 101) - # fmod between 0.00 and 1.00 - self.fexp_grid = self.fmod_grid - - # Initialize each log(k) value to one of the grid points. - # If there is no coverage, or there is a proline, make it zero - # so we do not use this position in calculations. Otherwise, use - # a random integer to refer to the exp_grid value to use for log(k). - # Or use a user defined value or just the first value. - - if init == "enumerate": + self.sigma0_grid=numpy.linspace(0.1,10,100) #stdev between 0.01 and 0.5 + self.omega_grid=numpy.linspace(1,100,100) #stdev between 0.01 and 0.5 + self.fmod_grid=numpy.linspace(0,100,101) #fmod between 0.00 and 1.00 + self.fexp_grid=self.fmod_grid #fmod between 0.00 and 1.00 + + #Initialize each log(k) value to one of the grid points. + #If there is no coverage, or there is a proline, make it zero + #so we do not use this position in calculations. Otherwise, use + #a random integer to refer to the exp_grid value to use for log(k). + #Or use a user defined value or just the first value. + + if init=="enumerate": self.guess_init_exp_sequence(5) else: - for n in range(0, model.num_res-1): - if (model.seq[n] == 'P' or model.seq[n] == 'p' - or state.coverage[n] == 0): - self.exp_seq[n] = int(0) - elif init == "random": - self.exp_seq[n] = numpy.random.randint( - len(self.exp_grid)-1)+1 - elif (isinstance(init, int) and init > 0 - and init < len(self.exp_grid)-1): - self.exp_seq[n] = init - elif isinstance(init, int) and len(init) == len(seq): - self.exp_seq = init + for n in range(0,model.num_res-1): + if model.seq[n]=='P' or model.seq[n]=='p' or state.coverage[n]==0: + self.exp_seq[n]=int(0) + elif init=="random": + self.exp_seq[n]=numpy.random.randint(len(self.exp_grid)-1)+1 + elif isinstance(init,int) and init > 0 and init < len(self.exp_grid)-1: + self.exp_seq[n]=init + elif isinstance(init,int) and len(init)==len(seq): + self.exp_seq=init else: - self.exp_seq[n] = 1 + self.exp_seq[n]=1 - state.exp_model = self - self.exp_models = [] + state.exp_model=self + self.exp_models=[] if marginalize: self.calculate_marginal_probabilities() - self.score = self.calculate_bayesian_score( - self.frags, sigma, error_model=error_model) + self.score=self.calculate_bayesian_score(self.frags, sigma, error_model=error_model) - def set_sigma(self, sigma): - self.sigma = sigma + def set_sigma(self,sigma): + self.sigma=sigma - def calculate_bayesian_score(self, frags, sig=None, - error_model="gaussian", report_stats=False): - joint_likelihood = 0 + def calculate_bayesian_score(self, frags, sig=None, error_model="gaussian", report_stats = False): + total_score=0 + joint_likelihood=0 if sig is None: sig = self.sigma # Get the sigma for this timepoint - j = self.get_nearest_value(sig, self.sigma0_grid) - # print("SIG:: CBS:", sig, j, self.sigma0_grid[j]) + j=self.get_nearest_value(sig, self.sigma0_grid) + #print("SIG:: CBS:", sig, j, self.sigma0_grid[j]) for f in frags: frag_likelihood = 0 frag_replicates = 0 + tp_likelihood=0 # Crash if this is not included due to the use outside the loop for tp in f.timepoints: tp.clear_model_deuteration_values() - tp_avg, tp_sd = tp.get_avg_sd() - tp_likelihood = 0 + tp_avg, tp_sd=tp.get_avg_sd() + tp_likelihood=0 - model_deut = 0 + model_deut=0 - sigma = tp.sigma + sigma=tp.sigma # Calculate model deuteration values - for r in range(f.start_res+self.offset+2, - f.end_res+self.offset+1): + #for r in range(f.start_res-self.offset,f.end_res-self.offset-1): + for r in range(f.start_res+self.offset+2,f.end_res+self.offset+1): if self.seq[r] != 'P' and self.seq[r] != 'p': - kr = 10**self.exp_grid[self.exp_seq[r]-1] - model_deut += (1-numpy.exp(-kr*tp.time)) + kr=10**self.exp_grid[self.exp_seq[r]-1] + model_deut+=(1-numpy.exp(-kr*tp.time)) - model_deut = 100.*model_deut/f.num_observable_amides + model_deut=100.*model_deut/f.num_observable_amides tp.add_model_deuteration_value(model_deut) - # fm=self.get_nearest_value(model_deut, self.fmod_grid) + #fm=self.get_nearest_value(model_deut, self.fmod_grid) # Calculate (or lookup) replicate likelihood for n in range(len(tp.replicates)): frag_replicates += 1 if self.marginalize: - fe = tp.replicates[n].get_nearest_gridpoint( - self.fexp_grid) + #print("Marginalized step") + + #delfmod=tp.replicates[n].deut-model_deut + #sig=100.*tp.get_sigma()/f.num_observable_amides + + #fe=self.get_nearest_value(tp.replicates[n].deut, self.fexp_grid) + fe = tp.replicates[n].get_nearest_gridpoint(self.fexp_grid) - replicate_score = self.prob_grid[j, fm, fe] + replicate_score=self.prob_grid[j,fm,fe] else: - if error_model == "gaussian": - replicate_score = self.gaussian_model( - tp.replicates[n].deut, model_deut, sigma) - elif error_model == "truncated_gaussian": - replicate_score = self.gaussian_model( - tp.replicates[n].deut, model_deut, sigma) - replicate_score += self.truncated_gaussian_factor( - tp.replicates[n].deut, -10, 120, sigma) - elif error_model == "lognormal": - replicate_score = self.lognormal_model( - tp.replicates[n].deut, model_deut, sigma) + if error_model=="gaussian": + replicate_score = self.gaussian_model(tp.replicates[n].deut, model_deut, sigma) + elif error_model=="truncated_gaussian": + replicate_score = self.gaussian_model(tp.replicates[n].deut, model_deut, sigma) + replicate_score += self.truncated_gaussian_factor(tp.replicates[n].deut,-10,120,sigma) + elif error_model=="lognormal": + replicate_score = self.lognormal_model(tp.replicates[n].deut, model_deut, sigma) else: print("Unrecognized Error Model:", error_model) sys.exit(1) # Hack to deal with log-based overflows from likelihood=0 if replicate_score != 0: - tp_likelihood = \ - tp_likelihood + -1*numpy.log(replicate_score) + tp_likelihood = tp_likelihood + -1*numpy.log(replicate_score) #numpy.log(replicate_score) #+ -1.0*numpy.log(self.get_sigma_score(self.sigma,tp.sigma)) else: tp_likelihood = tp_likelihood + 10000 - frag_likelihood = frag_likelihood+tp_likelihood + + + + #print(f.seq, tp.time, model_deut, [r.deut for r in tp.replicates], sigma, replicate_score, tp_likelihood) + + #avg_diff=(model_deut-numpy.average([r.deut for r in tp.replicates])) + frag_likelihood=frag_likelihood+tp_likelihood f.likelihood = frag_likelihood joint_likelihood = joint_likelihood + frag_likelihood - # print(f.seq, f.start_res, frag_likelihood, joint_likelihood) + #print(f.seq, f.start_res, frag_likelihood, joint_likelihood) return joint_likelihood def guess_init_exp_sequence(self, size_enum): - ''' Guess initial state by finding most likely exp values for - each fragment + ''' Guess initial state by finding most likely exp values for each fragment Returns a sequence of exp values ''' frags = deepcopy(self.frags) - exp_seq = [0]*len(self.exp_seq) - exp_grid = self.calc_exp_grid(size_enum, self.first_tp, self.last_tp) + exp_seq=[0]*len(self.exp_seq) + exp_grid=self.calc_exp_grid(size_enum, self.first_tp, self.last_tp) print("\n##### Initial Guess: Enumerate #####") print("# Enumerating all fragments with gridsize of", size_enum) print("# State:", self.state.state_name, "\n") print("> Fragment start_res end_res exp_grid") for f in frags: - f.exp_vals = sampling.enumerate_fragment( - f, exp_grid, 2.0, num_models=1) + f.exp_vals = sampling.enumerate_fragment(f, exp_grid, 2.0, num_models = 1) print(f.seq, f.exp_vals) - print("Calculating overlaps and casting to gridsize", - len(self.exp_grid)) + print("Calculating overlaps and casting to gridsize", len(self.exp_grid)) sectors = self.state.get_sectors(frags) # sort sectors by size sectors.sort(key=lambda x: x.num_amides, reverse=True) @@ -212,14 +211,11 @@ def guess_init_exp_sequence(self, size_enum): for s in sectors: sector_exp_vals = numpy.zeros(size_enum) sector_frags = deepcopy(s.fragments) - # if there is only one fragment, just place the residues. - # We are done. - if (len(sector_frags) == 1 - and sum(sector_frags[0].exp_vals) == s.num_amides): + # if there is only one fragment, just place the residues. We are done. + if len(sector_frags) == 1 and sum(sector_frags[0].exp_vals) == s.num_amides: method = "Single fragment" sector_exp_vals = sector_frags[0].exp_vals - sector_frags[0].exp_vals = \ - deepcopy(sector_frags[0].exp_vals) - sector_exp_vals + sector_frags[0].exp_vals = deepcopy(sector_frags[0].exp_vals) - sector_exp_vals else: # Populate sector with overlapping exp values # Is total overlap greater or less than s.num_amides? @@ -227,54 +223,53 @@ def guess_init_exp_sequence(self, size_enum): # Overlap is the minimum value for each field # How to deal with the rest? High to low or low to high? # Three possibilities: - # * Equal - assign overlap to each amide. - # Location doesn't matter - # * sum(overlap) > len(seq) - take (random?) subset - # of overlap and assign to sequence - # * sum(overlap) < len(seq) - Assign overlap to sequence. - # Find closest non-overlapping values - - all_exp_vals = [] + # * Equal - assign overlap to each amide. Location doesn't matter + # * sum(overlap) > len(seq) - take (random?) subset of overlap and assign to sequence + # * sum(overlap) < len(seq) - Assign overlap to sequence. Find closest non-overlapping values + + all_exp_vals = [] #numpy.zeros((len(frags),len(self.exp_grid))) for f in sector_frags: all_exp_vals.append(f.exp_vals) overlap = numpy.amin(numpy.array(all_exp_vals), axis=0) - # ------------------ + #------------------ # For equal number of overlap values and amides: - # ------------------ - if numpy.sum(overlap) == s.num_amides: + #------------------ + if numpy.sum(overlap)==s.num_amides: method = "Equal Overlap" sector_exp_vals = overlap + #print("Overlap = num_amides", s.seq, overlap, s.start_res) for f in sector_frags: f.exp_vals = deepcopy(f.exp_vals) - overlap - # print(exp_seq[s.start_res:s.end_res+1]) - # ------------------ + #print(exp_seq[s.start_res:s.end_res+1]) + #------------------ # For more overlap values than amides: # - randomly pick amides to assign - # ------------------ + #------------------ elif numpy.sum(overlap) > s.num_amides: method = "Overlap > num_amides" overlap_sub = numpy.zeros(len(overlap)).astype(int) - i = 1 + i=1 while i <= s.num_amides: - picked_val = randint(0, len(overlap)-1) + picked_val = randint(0,len(overlap)-1) + #print(i, picked_val, overlap, overlap[picked_val], overlap_sub) if overlap[picked_val] != 0: overlap_sub[picked_val] += 1 overlap[picked_val] -= 1 - i += 1 + i+=1 + #print("Overlap > num_amides", s.seq, overlap, overlap_sub, s.start_res) sector_exp_vals = overlap_sub for f in sector_frags: f.exp_vals = deepcopy(f.exp_vals) - overlap_sub - # ------------------ + #------------------ # For fewer overlap values than amides: # - assign overlap - # - start from smallest fragment and assign values until - # you hit n_remaining - # ------------------ + # - start from smallest fragment and assign values until you hit n_remaining + #------------------ else: method = "Overlap < num_amides" overlap_sub = numpy.zeros(len(overlap)).astype(int) @@ -282,222 +277,228 @@ def guess_init_exp_sequence(self, size_enum): for f in sector_frags: f.exp_vals = deepcopy(f.exp_vals) - overlap - # sort frags from smallest to largest + #sort frags from smallest to largest sector_frags.sort(key=lambda x: x.num_observable_amides) - i = 1 - # randomly pick value from smallest fragment + #print("***Overlap < num_amides", s.seq, s.num_amides - sum(overlap), overlap, sector_frags[0].exp_vals, s.start_res) + i=1 + #randomly pick value from smallest fragment while i <= n_remaining: - picked_val = randint(0, len(overlap)-1) + picked_val = randint(0,len(overlap)-1) if sector_frags[0].exp_vals[picked_val] != 0: overlap_sub[picked_val] += 1 - sector_frags[0].exp_vals[picked_val] = \ - deepcopy( - sector_frags[0].exp_vals[picked_val])-1 + sector_frags[0].exp_vals[picked_val] = deepcopy(sector_frags[0].exp_vals[picked_val])-1 + #print(i, n_remaining, picked_val, overlap_sub, sector_frags[0].exp_vals) i = i+1 sector_exp_vals = overlap + overlap_sub - # For each other frag, delete exp_vals closest - # to overlap_sub + # For each other frag, delete exp_vals closest to overlap_sub for f in sector_frags[1:-1]: - nonzero_elements = f.exp_vals.nonzero()[0] + nonzero_elements=f.exp_vals.nonzero()[0] for i in range(len(exp_grid)): if overlap_sub[i] != 0: - # if we are at i=0, delete first nonzero - if i == 0: + if i==0: # if we are at i=0, delete first nonzero first_nze = nonzero_elements[0] - f.exp_vals[first_nze] = \ - deepcopy(f.exp_vals[first_nze]) - 1 + f.exp_vals[first_nze] = deepcopy(f.exp_vals[first_nze]) - 1 - # if we are at i=N, delete last nonzero - elif i == len(exp_grid)-1: + elif i==len(exp_grid)-1: # if we are at i=N, delete last nonzero last_nze = nonzero_elements[-1] - f.exp_vals[last_nze] = \ - deepcopy(f.exp_vals[last_nze]) - 1 - - # Otherwise a bit more complicated - else: - # decide whether to go + or - - direction = int(2*(randint(0, 1)) - 1) - # print(i, direction, f.exp_vals) - # try 1st direction - if f.exp_vals[i+direction] != 0: - f.exp_vals[i+direction] = \ - deepcopy( - f.exp_vals[i+direction]) - 1 - # if not, try second direction - elif f.exp_vals[i-direction] != 0: - f.exp_vals[i-direction] = \ - deepcopy( - f.exp_vals[i-direction]) - 1 - # cop out for this...just use the - # first one. - else: - f.exp_vals[nonzero_elements[0]] = \ - deepcopy( - f.exp_vals[ - nonzero_elements[0]]) - 1 + f.exp_vals[last_nze] = deepcopy(f.exp_vals[last_nze]) - 1 + + else: #Otherwise a bit more complicated + direction = int(2*(randint(0,1)) - 1) #decide whether to go + or - + #print(i, direction, f.exp_vals) + if f.exp_vals[i+direction] != 0: # try 1st direction + f.exp_vals[i+direction] = deepcopy(f.exp_vals[i+direction]) - 1 + elif f.exp_vals[i-direction] != 0: #if not, try second direction + f.exp_vals[i-direction] = deepcopy(f.exp_vals[i-direction]) - 1 + else: #cop out for this...just use the first one. + f.exp_vals[nonzero_elements[0]] = deepcopy(f.exp_vals[nonzero_elements[0]]) - 1 if numpy.sum(sector_exp_vals) != s.num_amides: - print("ERROR: sector ", s.seq, - "does not have right number of amides,", s.num_amides) - print("Overlap not correctly calculated. Sector exp_vals:", - sector_exp_vals) + print("ERROR: sector ", s.seq, "does not have right number of amides,", s.num_amides) + print("Overlap not correctly calculated. Sector exp_vals:", sector_exp_vals) print("Method:", method) sys.exit(1) + # cast exp_grid to self.exp_grid for each sector - long_sector_exp_vals = numpy.zeros(len(self.exp_grid)).astype(int) + long_sector_exp_vals=numpy.zeros(len(self.exp_grid)).astype(int) # first and last are identical for each. - long_sector_exp_vals[0] = sector_exp_vals[0] - long_sector_exp_vals[-1] = sector_exp_vals[-1] + long_sector_exp_vals[0]=sector_exp_vals[0] + long_sector_exp_vals[-1]=sector_exp_vals[-1] - # Cast middle values + #Cast middle values - # something is wrong with these indicies...need to fix. - # long_sector -1 is getting overwritten - factor = float(len(sector_exp_vals)) \ - / (float(len(long_sector_exp_vals))) + #something is wrong with these indicies...need to fix. long_sector -1 is getting overwritten + factor = float(len(sector_exp_vals))/(float(len(long_sector_exp_vals))) - for i in range(1, len(sector_exp_vals)-1): - for j in range(1, len(long_sector_exp_vals)-1): + for i in range(1,len(sector_exp_vals)-1): + for j in range(1,len(long_sector_exp_vals)-1): if j * factor <= i and (j+1)*factor >= i: - ran_int = randint(0, 1) + ran_int = randint(0,1) long_sector_exp_vals[j+ran_int] += sector_exp_vals[i] break + #print(j*factor, i,long_sector_exp_vals[j], sector_exp_vals[i]) + + #print(s.seq, sector_exp_vals, method) + #print(s.seq, long_sector_exp_vals) if numpy.sum(long_sector_exp_vals) != s.num_amides: - print("ERROR: sector ", s.seq, - "does not have right number of amides,", s.num_amides) - print(sector_exp_vals, "not correctly cast to", - long_sector_exp_vals) + print("ERROR: sector ", s.seq, "does not have right number of amides,", s.num_amides) + print(sector_exp_vals, "not correctly cast to", long_sector_exp_vals) print("Method:", method) sys.exit(1) for n in range(s.start_res, s.end_res+1): - if self.seq[n] == 'P' or self.seq == 'p': - exp_seq[n] = 0 + if self.seq[n]=='P' or self.seq=='p': + exp_seq[n]=0 else: for i in range(len(self.exp_grid)): if long_sector_exp_vals[i] != 0: exp_seq[n] = i+1 - long_sector_exp_vals[i] -= 1 + long_sector_exp_vals[i] = long_sector_exp_vals[i] - 1 break + #print("SSSS", s.seq, sector_exp_vals, exp_seq[s.start_res:s.end_res+1]) + + #for f in frags: + #print(f.seq, f.exp_vals) + + #for n in range(len(self.exp_seq)): + #print(n,self.state.seq[n], exp_seq[n]) + #print(exp_seq) self.exp_seq = exp_seq return exp_seq + + #print("Fragment Amides:",s.num_amides,"Overlap:",overlap, "Diff:" , s.num_amides - numpy.sum(overlap), s.seq) + + ''' + + def calc_overlap_score(self, seq): + coverage = self.state.get_coverage(self.frags) + seclist = self.state.get_sectors().sort + sorted_sec = + + if coverage == 0 or self.state.seq[n] == 'P': + seq[n] = 0 + else + for f in self.frags: + if n > f.start_res + 2 and n < f.end_res: + + ''' + + def get_sigma(self, inval, sig): return inval*sig+0.05 - def gaussian_model(self, exp, model, sig): - # print("hello") - return numpy.exp(-((model-exp)**2)/(2*sig**2)) \ - / (sig*numpy.sqrt(2*numpy.pi)) + def gaussian_model(self,exp,model,sig): + #print("hello") + return numpy.exp(-((model-exp)**2)/(2*sig**2))/(sig*numpy.sqrt(2*numpy.pi)) - def lognormal_model(self, exp, model, sig): - return numpy.exp(-((numpy.log(model)-numpy.log(exp))**2)/(2*sig**2)) \ - / (sig*numpy.sqrt(numpy.pi)*model) + def lognormal_model(self,exp,model,sig): + return numpy.exp(-( (numpy.log(model)-numpy.log(exp) )**2)/(2*sig**2))/(sig*numpy.sqrt(numpy.pi)*model) - def truncated_gaussian_factor(self, exp, a, b, sig): - return 1 / (0.5 * ( - scipy.special.erf((b-exp)/sig * numpy.sqrt(3.1415)) - - scipy.special.erf((a-exp)/sig * numpy.sqrt(3.1415)))) + def truncated_gaussian_factor(self,exp,a,b,sig): + return 1/ ( 0.5 * ( scipy.special.erf( (b-exp)/sig * numpy.sqrt(3.1415) ) - scipy.special.erf( (a-exp)/sig * numpy.sqrt(3.1415) ) ) ) def set_exp_grid(self, grid): - self.exp_grid = grid + self.exp_grid=grid - def calc_exp_grid(self, num_bins, first_time_point=10, - last_time_point=3600): + def calc_exp_grid(self, num_bins, first_time_point=10, last_time_point=3600): if num_bins < 2: - raise Exception("get_exp_bins: num_bins must be greater than 1") - # slow_exp_val should satisfy - # 0.99 = exp(-10**slowest_rate_bin*first_time_point) - # fast_exp_val should satisfy - # 0.01 = exp(-10**fastest_rate_bin*last_time_point) + raise Exception("get_exp_bins: num_bins must be greater than 1" ) + # slow_exp_val should satisfy 0.99 = exp(-10**slowest_rate_bin*first_time_point) + # fast_exp_val should satisfy 0.01 = exp(-10**fastest_rate_bin*last_time_point) - slowest_rate_bin = numpy.log10(-numpy.log(0.01)/first_time_point) - fastest_rate_bin = numpy.log10(-numpy.log(0.99)/last_time_point) + slowest_rate_bin=numpy.log10(-numpy.log(0.01)/first_time_point) + fastest_rate_bin=numpy.log10(-numpy.log(0.99)/last_time_point) + #print("RATE BINS", slowest_rate_bin, fastest_rate_bin, first_time_point, last_time_point) return numpy.linspace(fastest_rate_bin, slowest_rate_bin, num_bins) def get_sector_averaged_protection_values(self, seq, sectors): - sector_averaged_protection_values = [0.]*len(seq) + sector_averaged_protection_values=[0.]*len(seq) for s in sectors: - sum_exp = 0 + sum_exp=0 if s.num_amides != 0: for n in range(s.start_res, s.end_res+1): - sum_exp = sum_exp+self.exp_grid[int(seq[n])-1] - frac_exp = sum_exp/s.num_amides + #print n, seq[n], sum_exp, self.exp_seq[n], self.exp_grid[int(seq[n])-1] + sum_exp=sum_exp+self.exp_grid[int(seq[n])-1] + frac_exp=sum_exp/s.num_amides for n in range(s.start_res, s.end_res+1): - if self.seq[n] != 'p' and self.seq[n] != 'P': - sector_averaged_protection_values[n] = frac_exp + if self.seq[n]!='p' and self.seq[n]!='P': + sector_averaged_protection_values[n]=frac_exp else: - sector_averaged_protection_values[n] = 0 + sector_averaged_protection_values[n]=0 return sector_averaged_protection_values - def get_indiv_sector_averaged_protection_values(self, seq, sectors): - sector_averaged_protection_values = [0.]*len(sectors) + def get_indiv_sector_averaged_protection_values(self,seq, sectors): + sector_averaged_protection_values=[0.]*len(sectors) for snum in range(len(sectors)): s = sectors[snum] - sum_exp = 0 + sum_exp=0 for n in range(s.start_res, s.end_res+1): - sum_exp = sum_exp+self.exp_grid[int(seq[n])-1] - frac_exp = sum_exp/s.num_amides - sector_averaged_protection_values[snum] = frac_exp + #print n, seq[n], sum_exp, self.exp_seq[n], self.exp_grid[int(seq[n])-1] + sum_exp=sum_exp+self.exp_grid[int(seq[n])-1] + frac_exp=sum_exp/s.num_amides + sector_averaged_protection_values[snum]=frac_exp return sector_averaged_protection_values - def get_protection_values(self, seq): - protection_values = [0.]*len(seq) + def get_protection_values(self,seq): + protection_values=[0.]*len(seq) for n in range(len(seq)): - protection_values[n] = self.exp_grid[self.exp_seq[n]-1] + protection_values[n]=self.exp_grid[self.exp_seq[n]-1] return protection_values - def get_protection_bins(self, seq): - protection_values = [0.]*len(seq) + def get_protection_bins(self,seq): + protection_values=[0.]*len(seq) for n in range(len(seq)): - protection_values[n] = self.exp_seq[n] + protection_values[n]=self.exp_seq[n] return protection_values def get_frag_deut_comparison(self, frags, exp_seq="none"): - if exp_seq == "none": - exp_seq = self.exp_seq + if exp_seq=="none": + exp_seq=self.exp_seq for f in frags: for t in f.timepoints: - model_deut = 0 - grid = [] - for r in range(f.start_res+2-self.offset, - f.end_res-self.offset): - kr = 10**self.exp_grid[exp_seq[r]-1] - model_deut = model_deut+(1-numpy.exp(-kr*t.time))*self.sat + model_deut=0 + grid=[] + #print len(grid), f.start_res, f.end_res + for r in range(f.start_res+2-self.offset, f.end_res-self.offset): + kr=10**self.exp_grid[exp_seq[r]-1] + model_deut=model_deut+(1-numpy.exp(-kr*t.time))*self.sat grid.append(exp_seq[r]-1) + #print(f.seq, t.time, grid, model_deut) return model_deut def calc_exp_chiscores_from_file(self, frags, infile, outdir): self.import_model_deuteration_from_file(frags, infile, outdir) for f in frags: - chi = 0 - totd = 0 + chi=0 + totd=0 for t in f.timepoints: t.avg_sd_model() for d in t.deut: - chi = chi+(d-t.model_avg)**2/self.sigma - totd = totd+1 - # print f.seq, chi, t.time, t.model_avg, t.deut - f.chi = chi/totd + chi=chi+(d-t.model_avg)**2/self.sigma + totd=totd+1 + #print f.seq, chi, t.time, t.model_avg, t.deut + f.chi=chi/totd + #print f.seq, chi, [t.model_avg for t in f.timepoints], [t.model_avg for t in f.timepoints] def calc_exp_chiscores_from_gridvals(self, frags, infile, outdir): self.import_model_deuteration_from_gridvals(frags, infile, outdir) for f in frags: - chi = 0 - totd = 0 + chi=0 + totd=0 for t in f.timepoints: t.avg_sd_model() for d in t.deut: - chi = chi+(d-t.model_avg)**2/self.sigma - totd = totd+1 - # print f.seq, chi, t.time, t.model_avg, t.deut - f.chi = chi/totd + chi=chi+(d-t.model_avg)**2/self.sigma + totd=totd+1 + #print f.seq, chi, t.time, t.model_avg, t.deut + f.chi=chi/totd print(f.seq, chi) def get_model_distributions(self, frags): @@ -506,104 +507,112 @@ def get_model_distributions(self, frags): print(t.model) def get_model_deuteration(self, time, frag, exp_seq): - model_deut = 0 + model_deut=0 for r in range(frag.start_res+1, frag.end_res): - # print f.seq, t.time, r, self.seq[r] - if self.seq[r] != 'P' and self.seq[r] != 'p': - kr = 10**self.exp_grid[int(float(exp_seq[r])-1)] - model_deut = model_deut+(1-numpy.exp(-kr*time)) + #print f.seq, t.time, r, self.seq[r] + if self.seq[r] !='P' and self.seq[r] != 'p': + kr=10**self.exp_grid[int(float(exp_seq[r])-1)] + model_deut=model_deut+(1-numpy.exp(-kr*time)) return model_deut/frag.num_observable_amides*100 - def import_model_deuteration_from_file(self, frags, infile, firstline=1, - lastline=-1, append=False): - data = open(infile, "r") - if not append: - self.exp_models = [] + + def import_model_deuteration_from_file(self, frags, infile, firstline=1, lastline=-1, append=False): + data=open(infile,"r") + if append==False: + self.exp_models=[] for line in data: - # print line - exp = map(int, line.split(' ')) + #print line + exp=list(map(int, line.split(' '))) # map is not subscriptable from Python 3 self.exp_models.append(exp) for f in frags: - # print f.seq + #print f.seq for t in f.timepoints: - model_deut = 0 - if not append: + model_deut=0 + if append==False: t.clear_model_deuteration_values() - append = True + append=True for r in range(f.start_res+1, f.end_res): - # print f.seq, t.time, r, self.seq[r] - if self.seq[r] != 'P' and self.seq[r] != 'p': - kr = 10**self.exp_grid[int(float(exp[r])-1)] - model_deut = model_deut+(1-numpy.exp(-kr*t.time)) - model_deut = model_deut/f.num_observable_amides*100 + #print f.seq, t.time, r, self.seq[r] + if self.seq[r] !='P' and self.seq[r] != 'p': + kr=10**self.exp_grid[int(float(exp[r])-1)] + model_deut=model_deut+(1-numpy.exp(-kr*t.time)) + #print f.seq, self.seq[r], r, f.start_res, r-f.start_res-1, (1-numpy.exp(-kr*t.time))*self.sat + model_deut=model_deut/f.num_observable_amides*100 t.add_model_deuteration_value(model_deut) + #print(f.seq, t.time, len(t.models), model_deut, t.get_model_avg(), t.get_model_sd(), [r.deut for r in t.replicates]) + #print f.seq, t.time, t.get_model_avg(), t.get_model_sd(), t.get_avg_sd() def import_model_deuteration_from_gridvals(self, frags, models): for f in frags: - timepoints = f.timepoints + timepoints=f.timepoints for t in timepoints: for m in models: - model_deut = 0 - for r in range(f.start_res+self.offset+2, - f.end_res+self.offset+1): - if self.seq[r] != 'P' and self.seq[r] != 'p': - kr = 10**self.exp_grid[int(float(m[r])-1)] - model_deut += \ - (1-numpy.exp(-kr*t.time))*t.replicates[0].sat - model_deut = model_deut/f.num_observable_amides*100 + model_deut=0 + for r in range(f.start_res+self.offset+2, f.end_res+self.offset+1): + #print r, f.seq, f.start_res+1, f.end_res, self.seq[r] + if self.seq[r] !='P' and self.seq[r] != 'p': + #print(r, self.seq[r], m[r], "|", m) + kr=10**self.exp_grid[int(float(m[r])-1)] + model_deut+=(1-numpy.exp(-kr*t.time))*t.replicates[0].sat + #print(">>", r, m[r], int(float(m[r])), self.exp_grid[int(float(m[r])-1)]) + #print f.seq, self.seq[r], r, f.start_res, r-f.start_res-1, (1-numpy.exp(-kr*t.time))*self.sat + #print(model_deut, t.time, t.sigma, f.seq, f.start_res, f.end_res, m[f.start_res+2+self.offset:f.end_res+1+self.offset], m, self.exp_grid) + model_deut=model_deut/f.num_observable_amides*100 t.add_model_deuteration_value(model_deut) + #t.avg_sd_model() + #print t.time, t.model_avg, t.model_sd, t.deut def import_models_from_gridvals(self, models, append=False): - if append == 'False': - self.exp_models = [] + if append=='False': + self.exp_models=[] for i in models: self.exp_models.append(i) + def import_models_from_file(self, modelfile, append='False'): - data = open(modelfile, "r") - if append == 'False': - self.exp_models = [] + data=open(modelfile,"r") + if append=='False': + self.exp_models=[] for line in data: - self.exp_models.append(map(int, line.split(' '))) + self.exp_models.append(list(map(int, line.split(' ')))) # map is not subscriptable from Python 3 def calc_model_scores(self, frags=None, sig=1.0, error_model="gaussian"): - if self.exp_models == []: + if self.exp_models==[]: print("No models to calculate score") return if frags is None: - frags = self.frags + frags=self.frags - self.model_scores = [] + self.model_scores=[] for m in self.exp_models: - self.exp_seq = deepcopy(m) - score = self.calculate_bayesian_score(frags, sig, error_model) + self.exp_seq=deepcopy(m) + score=self.calculate_bayesian_score(frags,sig,error_model) self.model_scores.append(score) return self.model_scores def import_scores(self, scorefile): try: - data = open(scorefile, "r") + data=open(scorefile,"r") for line in data: self.scores.append(line.split().strip()) - except: # noqa: E722 + except: for i in scorefile: self.scores.append(i) - def get_sigma_score(self, sigma0, sig): - return (2*sigma0 / (numpy.sqrt(numpy.pi) * sig**2)) \ - * numpy.exp(-sigma0**2 / sig**2) + def get_sigma_score(self,sigma0,sig): + return (2*sigma0/ (numpy.sqrt(numpy.pi) * sig**2 ))*numpy.exp(-sigma0**2 /sig**2) def get_model_average(self): - # print(len(self.exp_models)) + #print(len(self.exp_models)) if len(self.exp_models) > 0: - modelavg = numpy.average(numpy.array(self.exp_models), axis=0) + #print(len(self.exp_models), len(self.exp_models[0]), numpy.array(self.exp_models).shape, self.exp_models[0]) + modelavg=numpy.average(numpy.array(self.exp_models), axis=0) return modelavg else: print("No models to get the average of =(") def unimodal_prior(self, omj, sigma0=0.1): - return (2.0 * omj / (numpy.sqrt(numpy.pi) * omj**2)) \ - * numpy.exp(-sigma0**2 / omj**2) + return (2.0 * omj / ( numpy.sqrt(numpy.pi) * omj**2 )) * numpy.exp(-sigma0**2 / omj**2) def get_nearest_value(self, value, array): ''' @@ -611,105 +620,97 @@ def get_nearest_value(self, value, array): ''' return (numpy.abs(array-value)).argmin(), - def calculate_marginal_probabilities(self, - error_model="truncated_gaussian"): - # Calculates marginal probabilities given grids of - # omega, sigma0, fmod, and fexp + def calculate_marginal_probabilities(self, error_model="truncated_gaussian"): + # Calculates marginal probabilities given grids of omega, sigma0, fmod, and fexp - self.prob_grid = numpy.zeros((len(self.sigma0_grid), - len(self.fmod_grid), - len(self.fexp_grid))) + self.prob_grid=numpy.zeros( (len(self.sigma0_grid), len(self.fmod_grid), len(self.fexp_grid)) ) for fm in range(len(self.fmod_grid)): - fmod = self.fmod_grid[fm] + fmod=self.fmod_grid[fm] for f in range(len(self.fexp_grid)): - fexp = self.fexp_grid[f] + fexp=self.fexp_grid[f] for s in range(len(self.sigma0_grid)): - sigma0 = self.sigma0_grid[s] - cumul = 0 - - for j in range(1, len(self.omega_grid)): - # We're going to calculate the likelihoods at - # each point - omj = self.omega_grid[j] - omjm1 = self.omega_grid[j-1] + sigma0=self.sigma0_grid[s] + cumul=0 + + for j in range(1,len(self.omega_grid)): + # We're going to calculate the likelihoods at each point + omj=self.omega_grid[j] + omjm1=self.omega_grid[j-1] priorj = self.unimodal_prior(omj, sigma0) priorjm1 = self.unimodal_prior(omjm1, sigma0) - dom = omj-omjm1 - - if error_model == "truncated_gaussian": - # If we are using the truncated gaussian, - # calculate those pre-factors - # use the truncated gaussian factors - # at -0.1 and 1.2 - factor = self.truncated_gaussian_factor( - fexp, -10, 120, omj) - factorm1 = self.truncated_gaussian_factor( - fexp, -10, 120, omjm1) + dom=omj-omjm1 + + + if error_model=="truncated_gaussian": + # If we are using the truncated gaussian, calculate those pre-factors + # use the truncated gaussian factors at -0.1 and 1.2 + factor=self.truncated_gaussian_factor(fexp,-10,120,omj) + factorm1=self.truncated_gaussian_factor(fexp,-10,120,omjm1) else: - factor = 1.0 - factorm1 = 1.0 + factor=1.0 + factorm1=1.0 - # Calculate the likelihood at each point - pj = self.gaussian_model(fexp, fmod, omj)*factor*priorj - pjm1 = self.gaussian_model( - fexp, fmod, omjm1)*factorm1*priorjm1 - cumul = cumul+(pj+pjm1)/2.0/dom - self.prob_grid[s, fm, f] = -1.0 * numpy.log(cumul) + # Calculate the likelihood at each point + pj=self.gaussian_model(fexp,fmod,omj)*factor*priorj + pjm1=self.gaussian_model(fexp,fmod,omjm1)*factorm1*priorjm1 + + cumul = cumul+(pj+pjm1)/2.0/dom; + #print(pj, pjm1, omj, omjm1, sigma0, "| ", factor, factorm1, "| ", priorj, priorjm1, cumul) + #print(fexp, fmod, fexp-fmod, sigma0, "| ", cumul, -1.0 * numpy.log(cumul)) + #print j,fexp,fmod,omj,(pj+pjm1)/2.0/dom,cumul, cumul2, cumul2/cumul + self.prob_grid[s,fm,f]=-1.0 * numpy.log(cumul) return self.prob_grid - def get_marginalized_sigma_probabilities(self, sigma0, - error_model="gaussian"): + def get_marginalized_sigma_probabilities(self,sigma0,error_model="gaussian"): ''' - Precomputes marginal probabilities for sigma levels using a - cauchy distribution - @param delfmod - difference between forward model and noise model - in delta pct-D units to normalize for fragment length - (num amides) - @param sigma0 - Estimated standard deviation of the instrument - in pct-D units + Precomputes marginal probabilities for sigma levels using a cauchy distribution + @param delfmod - difference between forward model and noise model in delta pct-D + units to normalize for fragment length (num amides) + @param sigma0 - Estimated standard deviation of the instrument in pct-D units + ''' - self.prob_grid = numpy.zeros((len(self.fmod_grid), - len(self.fexp_grid))) + self.prob_grid=numpy.zeros( (len(self.fmod_grid), len(self.fexp_grid)) ) + for fm in range(len(self.fmod_grid)): - fmod = self.fmod_grid[fm] + fmod=self.fmod_grid[fm] # Unimodal Distributions for sigma # for f in range(len(self.fexp_grid)): - fexp = self.fexp_grid[f] - if error_model == "truncated_gaussian": - # If we are using the truncated gaussian, calculate - # those pre-factors + fexp=self.fexp_grid[f] + if error_model=="truncated_gaussian": + # If we are using the truncated gaussian, calculate those pre-factors # use the truncated gaussian factors at -0.1 and 1.2 - factor = self.truncated_gaussian_factor( - fexp, -10, 120, omj) - factorm1 = self.truncated_gaussian_factor( - fexp, -10, 120, omjm1) + factor=self.truncated_gaussian_factor(fexp,-10,120,omj) + factorm1=self.truncated_gaussian_factor(fexp,-10,120,omjm1) + o=1 else: - factor = 1.0 - factorm1 = 1.0 + factor=1.0 + factorm1=1.0 for j in range(len(self.omega_grid)): # We're going to calculate the likelihoods at each point - omj = self.omega_grid[j] - omjm1 = self.omega_grid[j-1] + omj=self.omega_grid[j] + omjm1=self.omega_grid[j-1] priorj = self.unimodal_prior(omj, sigma0) priorjm1 = self.unimodal_prior(omjm1, sigma0) - dom = omj-omjm1 + dom=omj-omjm1 - cumul = 0 + cumul=0 # Calculate the likelihood at each point - pj = self.gaussian_model(fexp, fmod, omj)*factor*priorj - pjm1 = self.gaussian_model( - fexp, fmod, omjm1)*factorm1*priorjm1 + pj=self.gaussian_model(fexp,fmod,omj)*factor*priorj + pjm1=self.gaussian_model(fexp,fmod,omjm1)*factorm1*priorjm1 + + cumul = cumul+(pj+pjm1)/2.0/dom; - cumul = cumul+(pj+pjm1)/2.0/dom + #print(j, fexp, fmod, sigma0, omj,(pj+pjm1)/2.0/dom,cumul) - self.prob_grid[fm, f] = -1.0 * numpy.log(cumul) + self.prob_grid[fm,f]=-1.0 * numpy.log(cumul) + #print "OMEGA: ",fexp,fmod,omega0,cumul,cumul2,cumul2/cumul return self.prob_grid diff --git a/pyext/src/model.py b/pyext/src/model.py index cf09930..6f03379 100644 --- a/pyext/src/model.py +++ b/pyext/src/model.py @@ -1,33 +1,38 @@ """ - Classes that store the representation and the sampled parameters for each - state + Classes that store the representation and the sampled parameters for each state Models are defined for each system state. They contain the sampled parameters for each state. """ - from __future__ import print_function +#import IMP +#import IMP.hdx +import system +import analysis +import hxio import sampling import numpy import scipy +#import scipy.special +from numpy import linalg import sys from copy import deepcopy +import os +import os.path +import time from random import randint class ResidueGridModel(object): ''' - Models the system as individual residues with protection factors along a - grid - Defined by a grid_size and the parameters of the datasets included in - the state + Models the system as individual residues with protection factors along a grid + Defined by a grid_size and the parameters of the datasets included in the state Calculates the grid of observable protection factors and converts between grid values and protection factors. @param state - the state that this model applies to @param grid_size - the size of the sampling grid - @param protection_factors - Boolean. Set to true to calculate protection - factors. False to calculate rates only. + @param protection_factors - Boolean. Set to true to calculate protection factors. False to calculate rates only. ''' def __init__(self, state, grid_size): self.state = deepcopy(state) @@ -46,19 +51,13 @@ def generate_model(self, random=True, value=1, initialize=False): sequence = self.state.get_sequence() model = numpy.zeros(self.length) - if (not random - and (type(value) is not int or value > self.grid_size - or value < 1)): - raise Exception( - "ResidueGridModel.generate_model: Value error. Either allow " - "random assignment or set an integer value from 1 " - "to grid_size") + if not random and (type(value) is not int or value > self.grid_size or value < 1): + raise Exception("ResidueGridModel.generate_model: Value error. Either allow random assignment or set an integer value from 1 to grid_size") for i in range(self.length): if sequence[i] != "P": if random: - model[i] = numpy.random.randint( - 1, high=self.grid_size, size=1)[0] + model[i] = numpy.random.randint(1, high=self.grid_size, size=1)[0] else: model[i] = int(value) @@ -80,34 +79,36 @@ def convert_model_to_exchange_rates(self, model): pass def convert_model_to_protection_factors(self, model): - # For a vector of grid values (the model), return a vector - # of protection values + # For a vector of grid values (the model), return a vector of protection values pf_grids = self.pf_grids - - self.model_protection_factors = [ - pf_grids[i][int(model[i]-1)] for i in range(self.length)] + #print(model) + #for i in range(self.length): + # grid = pf_grids[i] + # pf = grid[int(model[i])] + #self.model_protection_factors=[] + #for i in range(self.length): + # print(i, pf_grids[i], model[i], len(pf_grids), self.length) + # self.model_protection_factors.append(pf_grids[i][int(model[i])]) + + self.model_protection_factors = [pf_grids[i][int(model[i]-1)] for i in range(self.length)] return self.model_protection_factors - def calculate_protection_factor_grids(self, threshold=0.01): + def calculate_protection_factor_grids(self, threshold = 0.01): ''' - We theorize that the protection factor will be constant over multiple - experiments. Therefore, the grid of observable protection factors for - a given residue is dependent on all datasets. - - @param threshold - Means that the grid will start at a protection - factor value where (1-threshold) * 100% - will be observed at the first timepoint and end at a protection - factor value where threshold * 100% of exchange will be - observed for this site at the longest timepoint. + We theorize that the protection factor will be constant over multiple experiments. + Therefore, the grid of observable protection factors for a given residue is dependent + on all datasets. + + @param threshold - Means that the grid will start at a protection factor value where (1-threshold) * 100% + will be observed at the first timepoint and end at a protection factor value where threshold * 100% + of exchange will be observed for this site at the longest timepoint. ''' + bounds = [] - # For each dataset, get the observable protection factors for - # each residue + # For each dataset, get the observable protection factors for each residue observable_pfs = [] if not self.state.has_data: - raise Exception( - "Cannot calculate protection factor grid because state " - + self.state.name + " has no associated datasets") + raise Exception("Cannot calculate protection factor grid because state " + self.state.name + " has no associated datasets") for d in self.state.get_datasets(): d.calculate_observable_rate_bounds(threshold) observable_pfs.append(d.calculate_observable_protection_factors()) @@ -115,7 +116,10 @@ def calculate_protection_factor_grids(self, threshold=0.01): # The observable range is the largest range over all datasets pf_grids = [] for n in range(self.length): - pf_grid = numpy.linspace(0, 10, self.grid_size) + pf_ranges = [pf[n] for pf in observable_pfs] + #print("II", pf_ranges, max([pf[1] for pf in pf_ranges]), min([pf[0] for pf in pf_ranges]), self.grid_size) + #pf_grid = numpy.linspace(min([pf[1] for pf in pf_ranges]), max([pf[0] for pf in pf_ranges]), self.grid_size ) + pf_grid = numpy.linspace( 0, 10, self.grid_size ) pf_grids.append(pf_grid) self.pf_grids = pf_grids @@ -125,26 +129,25 @@ def calculate_protection_factor_grids(self, threshold=0.01): def get_conversion_grid(self): ''' - Returns a list (of length = # of residues) of lists - (of length = grid_size) with the protection factor values for - each residue. + Returns a list (of length = # of residues) of lists (of length = grid_size) + with the protection factor values for each residue. ''' return self.pf_grids def output_conversion_model(self, outfile): model_string = "" for i in self.model: - model_string += " "+str(i) - model_string = str(self.model) + the_string+=" "+str(i) + the_string = str(self.model) # Conversion_string is a list of pf_grid_value_string = str(pf_grid) return pf_grid_value_string, model_string[0:-1] def change_residue(self, residue, value): + #print(value, len(self.model), len(self.model_protection_factors), len(self.pf_grids), len(self.pf_grids[residue-1])) self.model[residue-1] = value - self.model_protection_factors[residue-1] \ - = self.pf_grids[residue-1][int(value-1)] + self.model_protection_factors[residue-1] = self.pf_grids[residue-1][int(value-1)] def get_model_residue(self, residue): return self.model[residue-1] @@ -155,6 +158,7 @@ def get_model(self): def get_masked_model(self, observed_residues): # get model only returning observed values (0 for everything else) mod = [] + #print(observed_residues) for i in range(len(self.model)): if i+1 in observed_residues: mod.append(int(self.model[i])) @@ -163,114 +167,101 @@ def get_masked_model(self, observed_residues): return mod + class MultiExponentialModel(object): - """ A multi-exponential model for exchange that approximates each - independent sector as a sum of exponentials, corresponding to the - number of exchanging amides in that sector. - The exchange rate constant, kex, for each residue is sampled along - a grid of log(kex) values. - - The model can be initialized to a certain bin in the exponential grid, - or to random values. - The Bayesian likelihood of the model is calculated against a set - of hdx.system_setup.Fragment objects which contain experimental data. + """ A multi-exponential model for exchange that approximates each independent sector + as a sum of exponentials, corresponding to the number of exchanging amides in that sector. + The exchange rate constant, kex, for each residue is sampled along a grid of log(kex) values. + + The model can be initialized to a certain bin in the exponential grid, or to random values. + The Bayesian likelihood of the model is calculated against a set of hdx.system_setup.Fragment objects + which contain experimental data. """ - def __init__(self, model, state, frags=None, sigma=10.0, num_exp_bins=10, - init="enumerate", error_model="gaussian", - output_directory=None, noclobber=True, first_tp=10, - last_tp=3600, marginalize=False): + def __init__(self, model, state, frags=None, sigma=10.0, num_exp_bins=10, init="enumerate", + error_model="gaussian", output_directory=None, noclobber=True, first_tp=10, last_tp=3600, marginalize=False): if frags is None: - self.frags = state.frags + self.frags=state.frags else: - self.frags = frags - # print "using frags" - tpoints = [] - times = [] + self.frags=frags + #print "using frags" + tpoints=[] + times=[] # get first and last time points for f in state.frags: tpoints += [tp for tp in f.timepoints] - times += [t.time for t in tpoints] + times+=[t.time for t in tpoints] first_tp = numpy.min(times) last_tp = numpy.max(times) self.first_tp = first_tp self.last_tp = last_tp - self.state = state - self.state.exp_model = self - self.exp_grid = self.calc_exp_grid(num_exp_bins, first_tp, last_tp) - self.exp_seq = [0]*len(model.seq) - self.name = state.state_name + self.state=state + self.state.exp_model=self + self.exp_grid=self.calc_exp_grid(num_exp_bins, first_tp, last_tp) + self.exp_seq=[0]*len(model.seq) + self.name=state.state_name state.add_exp_model(self) - self.offset = model.offset - self.sigma = sigma - self.oldsigma = -1 - self.scores = [] - self.seq = model.seq - self.flag_outliers = True - self.output_directory = output_directory - self.noclobber = noclobber - self.marginalize = marginalize + self.offset=model.offset + self.sigma=sigma + self.oldsigma=-1 + self.scores=[] + self.seq=model.seq + self.flag_outliers=True + self.output_directory=output_directory + self.noclobber=noclobber + self.marginalize=marginalize # These grids contain the values at which the integral of sigma is # computed - - # stdev between 0.01 and 0.5 - self.sigma0_grid = numpy.linspace(0.1, 10, 100) - # stdev between 0.01 and 0.5 - self.omega_grid = numpy.linspace(1, 100, 100) - # fmod between 0.00 and 1.00 - self.fmod_grid = numpy.linspace(0, 100, 101) - # fmod between 0.00 and 1.00 - self.fexp_grid = self.fmod_grid - - # Initialize each log(k) value to one of the grid points. - # If there is no coverage, or there is a proline, make it zero - # so we do not use this position in calculations. Otherwise, use - # a random integer to refer to the exp_grid value to use for log(k). - # Or use a user defined value or just the first value. - - if init == "enumerate": + self.sigma0_grid=numpy.linspace(0.1,10,100) #stdev between 0.01 and 0.5 + self.omega_grid=numpy.linspace(1,100,100) #stdev between 0.01 and 0.5 + self.fmod_grid=numpy.linspace(0,100,101) #fmod between 0.00 and 1.00 + self.fexp_grid=self.fmod_grid #fmod between 0.00 and 1.00 + + #Initialize each log(k) value to one of the grid points. + #If there is no coverage, or there is a proline, make it zero + #so we do not use this position in calculations. Otherwise, use + #a random integer to refer to the exp_grid value to use for log(k). + #Or use a user defined value or just the first value. + + if init=="enumerate": self.guess_init_exp_sequence(5) else: - for n in range(0, model.num_res-1): - if (model.seq[n] == 'P' or model.seq[n] == 'p' - or state.coverage[n] == 0): - self.exp_seq[n] = int(0) - elif init == "random": - self.exp_seq[n] = numpy.random.randint( - len(self.exp_grid)-1)+1 - elif (isinstance(init, int) and init > 0 - and init < len(self.exp_grid)-1): - self.exp_seq[n] = init - elif isinstance(init, int) and len(init) == len(seq): - self.exp_seq = init + for n in range(0,model.num_res-1): + if model.seq[n]=='P' or model.seq[n]=='p' or state.coverage[n]==0: + self.exp_seq[n]=int(0) + elif init=="random": + self.exp_seq[n]=numpy.random.randint(len(self.exp_grid)-1)+1 + elif isinstance(init,int) and init > 0 and init < len(self.exp_grid)-1: + self.exp_seq[n]=init + elif isinstance(init,int) and len(init)==len(seq): + self.exp_seq=init else: - self.exp_seq[n] = 1 + self.exp_seq[n]=1 - state.exp_model = self - self.exp_models = [] + state.exp_model=self + self.exp_models=[] if marginalize: self.calculate_marginal_probabilities() - self.score = self.calculate_bayesian_score( - self.frags, sigma, error_model=error_model) + self.score=self.calculate_bayesian_score(self.frags, sigma, error_model=error_model) - def set_sigma(self, sigma): - self.sigma = sigma + def set_sigma(self,sigma): + self.sigma=sigma - def calculate_bayesian_score(self, frags, sig=None, - error_model="gaussian", report_stats=False): - joint_likelihood = 0 + def calculate_bayesian_score(self, frags, sig=None, error_model="gaussian", report_stats = False): + total_score=0 + joint_likelihood=0 if sig is None: sig = self.sigma # Get the sigma for this timepoint j = self.get_nearest_value(sig, self.sigma0_grid) - # print("SIG:: CBS:", sig, j, self.sigma0_grid[j]) + #print("SIG:: CBS:", sig, j, self.sigma0_grid[j]) for f in frags: @@ -280,82 +271,90 @@ def calculate_bayesian_score(self, frags, sig=None, for tp in f.timepoints: tp.clear_model_deuteration_values() - tp_avg, tp_sd = tp.get_avg_sd() - tp_likelihood = 0 + tp_avg, tp_sd=tp.get_avg_sd() + tp_likelihood=0 - model_deut = 0 + model_deut=0 - sigma = tp.sigma + sigma=tp.sigma # Calculate model deuteration values - for r in range(f.start_res+self.offset+2, - f.end_res+self.offset+1): + # for r in range(f.start_res-self.offset,f.end_res-self.offset-1): + for r in range(f.start_res+self.offset+2,f.end_res+self.offset+1): if self.seq[r] != 'P' and self.seq[r] != 'p': - kr = 10**self.exp_grid[self.exp_seq[r]-1] - model_deut += (1-numpy.exp(-kr*tp.time)) + kr=10**self.exp_grid[self.exp_seq[r]-1] + model_deut+=(1-numpy.exp(-kr*tp.time)) - model_deut = 100.*model_deut/f.num_observable_amides + model_deut=100.*model_deut/f.num_observable_amides tp.add_model_deuteration_value(model_deut) + #print(f.seq, tp.time, sigma, tp_avg, tp_sd, kr, model_deut, self.exp_grid) + + #fm=self.get_nearest_value(model_deut, self.fmod_grid) + # Calculate (or lookup) replicate likelihood for n in range(len(tp.replicates)): frag_replicates += 1 if self.marginalize: - fe = tp.replicates[n].get_nearest_gridpoint( - self.fexp_grid) - replicate_score = self.prob_grid[j, fm, fe] + #print("Marginalized step") + + #delfmod=tp.replicates[n].deut-model_deut + #sig=100.*tp.get_sigma()/f.num_observable_amides + + #fe=self.get_nearest_value(tp.replicates[n].deut, self.fexp_grid) + fe = tp.replicates[n].get_nearest_gridpoint(self.fexp_grid) + + replicate_score=self.prob_grid[j,fm,fe] else: - if error_model == "gaussian": - replicate_score = self.gaussian_model( - tp.replicates[n].deut, model_deut, sigma) - elif error_model == "truncated_gaussian": - replicate_score = self.gaussian_model( - tp.replicates[n].deut, model_deut, sigma) - replicate_score += self.truncated_gaussian_factor( - tp.replicates[n].deut, -10, 120, sigma) - elif error_model == "lognormal": - replicate_score = self.lognormal_model( - tp.replicates[n].deut, model_deut, sigma) + if error_model=="gaussian": + replicate_score = self.gaussian_model(tp.replicates[n].deut, model_deut, sigma) + elif error_model=="truncated_gaussian": + replicate_score = self.gaussian_model(tp.replicates[n].deut, model_deut, sigma) + replicate_score += self.truncated_gaussian_factor(tp.replicates[n].deut,-10,120,sigma) + elif error_model=="lognormal": + replicate_score = self.lognormal_model(tp.replicates[n].deut, model_deut, sigma) else: print("Unrecognized Error Model:", error_model) sys.exit(1) # Hack to deal with log-based overflows from likelihood=0 if replicate_score != 0: - tp_likelihood = \ - tp_likelihood + -1*numpy.log(replicate_score) + tp_likelihood = tp_likelihood + -1*numpy.log(replicate_score) #numpy.log(replicate_score) #+ -1.0*numpy.log(self.get_sigma_score(self.sigma,tp.sigma)) else: tp_likelihood = tp_likelihood + 10000 - frag_likelihood = frag_likelihood+tp_likelihood + + + + #print(f.seq, tp.time, model_deut, [r.deut for r in tp.replicates], sigma, replicate_score, tp_likelihood) + + #avg_diff=(model_deut-numpy.average([r.deut for r in tp.replicates])) + frag_likelihood=frag_likelihood+tp_likelihood f.likelihood = frag_likelihood joint_likelihood = joint_likelihood + frag_likelihood - # print(f.seq, f.start_res, frag_likelihood, joint_likelihood) + #print(f.seq, f.start_res, frag_likelihood, joint_likelihood) return joint_likelihood def guess_init_exp_sequence(self, size_enum): - ''' Guess initial state by finding most likely exp values for each - fragment + ''' Guess initial state by finding most likely exp values for each fragment Returns a sequence of exp values ''' frags = deepcopy(self.frags) - exp_seq = [0]*len(self.exp_seq) - exp_grid = self.calc_exp_grid(size_enum, self.first_tp, self.last_tp) + exp_seq=[0]*len(self.exp_seq) + exp_grid=self.calc_exp_grid(size_enum, self.first_tp, self.last_tp) print("\n##### Initial Guess: Enumerate #####") print("# Enumerating all fragments with gridsize of", size_enum) print("# State:", self.state.state_name, "\n") print("> Fragment start_res end_res exp_grid") for f in frags: - f.exp_vals = sampling.enumerate_fragment( - f, exp_grid, 2.0, num_models=1) + f.exp_vals = sampling.enumerate_fragment(f, exp_grid, 2.0, num_models = 1) print(f.seq, f.exp_vals) - print("Calculating overlaps and casting to gridsize", - len(self.exp_grid)) + print("Calculating overlaps and casting to gridsize", len(self.exp_grid)) sectors = self.state.get_sectors(frags) # sort sectors by size sectors.sort(key=lambda x: x.num_amides, reverse=True) @@ -363,14 +362,11 @@ def guess_init_exp_sequence(self, size_enum): for s in sectors: sector_exp_vals = numpy.zeros(size_enum) sector_frags = deepcopy(s.fragments) - # if there is only one fragment, just place the residues. - # We are done. - if (len(sector_frags) == 1 - and sum(sector_frags[0].exp_vals) == s.num_amides): + # if there is only one fragment, just place the residues. We are done. + if len(sector_frags) == 1 and sum(sector_frags[0].exp_vals) == s.num_amides: method = "Single fragment" sector_exp_vals = sector_frags[0].exp_vals - sector_frags[0].exp_vals = \ - deepcopy(sector_frags[0].exp_vals) - sector_exp_vals + sector_frags[0].exp_vals = deepcopy(sector_frags[0].exp_vals) - sector_exp_vals else: # Populate sector with overlapping exp values # Is total overlap greater or less than s.num_amides? @@ -378,53 +374,53 @@ def guess_init_exp_sequence(self, size_enum): # Overlap is the minimum value for each field # How to deal with the rest? High to low or low to high? # Three possibilities: - # * Equal - assign overlap to each amide. - # Location doesn't matter - # * sum(overlap) > len(seq) - take (random?) subset of - # overlap and assign to sequence - # * sum(overlap) < len(seq) - Assign overlap to sequence. - # Find closest non-overlapping values - - all_exp_vals = [] + # * Equal - assign overlap to each amide. Location doesn't matter + # * sum(overlap) > len(seq) - take (random?) subset of overlap and assign to sequence + # * sum(overlap) < len(seq) - Assign overlap to sequence. Find closest non-overlapping values + + all_exp_vals = [] #numpy.zeros((len(frags),len(self.exp_grid))) for f in sector_frags: all_exp_vals.append(f.exp_vals) overlap = numpy.amin(numpy.array(all_exp_vals), axis=0) - # ------------------ + #------------------ # For equal number of overlap values and amides: - # ------------------ - if numpy.sum(overlap) == s.num_amides: + #------------------ + if numpy.sum(overlap)==s.num_amides: method = "Equal Overlap" sector_exp_vals = overlap + #print("Overlap = num_amides", s.seq, overlap, s.start_res) for f in sector_frags: f.exp_vals = deepcopy(f.exp_vals) - overlap - # ------------------ + #print(exp_seq[s.start_res:s.end_res+1]) + #------------------ # For more overlap values than amides: # - randomly pick amides to assign - # ------------------ + #------------------ elif numpy.sum(overlap) > s.num_amides: method = "Overlap > num_amides" overlap_sub = numpy.zeros(len(overlap)).astype(int) - i = 1 + i=1 while i <= s.num_amides: - picked_val = randint(0, len(overlap)-1) + picked_val = randint(0,len(overlap)-1) + #print(i, picked_val, overlap, overlap[picked_val], overlap_sub) if overlap[picked_val] != 0: overlap_sub[picked_val] += 1 overlap[picked_val] -= 1 - i += 1 + i+=1 + #print("Overlap > num_amides", s.seq, overlap, overlap_sub, s.start_res) sector_exp_vals = overlap_sub for f in sector_frags: f.exp_vals = deepcopy(f.exp_vals) - overlap_sub - # ------------------ + #------------------ # For fewer overlap values than amides: # - assign overlap - # - start from smallest fragment and assign values - # until you hit n_remaining - # ------------------ + # - start from smallest fragment and assign values until you hit n_remaining + #------------------ else: method = "Overlap < num_amides" overlap_sub = numpy.zeros(len(overlap)).astype(int) @@ -432,222 +428,228 @@ def guess_init_exp_sequence(self, size_enum): for f in sector_frags: f.exp_vals = deepcopy(f.exp_vals) - overlap - # sort frags from smallest to largest + #sort frags from smallest to largest sector_frags.sort(key=lambda x: x.num_observable_amides) - i = 1 - # randomly pick value from smallest fragment + #print("***Overlap < num_amides", s.seq, s.num_amides - sum(overlap), overlap, sector_frags[0].exp_vals, s.start_res) + i=1 + #randomly pick value from smallest fragment while i <= n_remaining: - picked_val = randint(0, len(overlap)-1) + picked_val = randint(0,len(overlap)-1) if sector_frags[0].exp_vals[picked_val] != 0: overlap_sub[picked_val] += 1 - sector_frags[0].exp_vals[picked_val] = \ - deepcopy( - sector_frags[0].exp_vals[picked_val])-1 + sector_frags[0].exp_vals[picked_val] = deepcopy(sector_frags[0].exp_vals[picked_val])-1 + #print(i, n_remaining, picked_val, overlap_sub, sector_frags[0].exp_vals) i = i+1 sector_exp_vals = overlap + overlap_sub - # For each other frag, delete exp_vals closest to - # overlap_sub + # For each other frag, delete exp_vals closest to overlap_sub for f in sector_frags[1:-1]: - nonzero_elements = f.exp_vals.nonzero()[0] + nonzero_elements=f.exp_vals.nonzero()[0] for i in range(len(exp_grid)): if overlap_sub[i] != 0: - # if we are at i=0, delete first nonzero - if i == 0: + if i==0: # if we are at i=0, delete first nonzero first_nze = nonzero_elements[0] - f.exp_vals[first_nze] = \ - deepcopy(f.exp_vals[first_nze]) - 1 + f.exp_vals[first_nze] = deepcopy(f.exp_vals[first_nze]) - 1 - # if we are at i=N, delete last nonzero - elif i == len(exp_grid)-1: + elif i==len(exp_grid)-1: # if we are at i=N, delete last nonzero last_nze = nonzero_elements[-1] - f.exp_vals[last_nze] = \ - deepcopy(f.exp_vals[last_nze]) - 1 - - # Otherwise a bit more complicated - else: - # decide whether to go + or - - direction = int(2*(randint(0, 1)) - 1) - # try 1st direction - if f.exp_vals[i+direction] != 0: - f.exp_vals[i+direction] = \ - deepcopy( - f.exp_vals[i+direction]) - 1 - # if not, try second direction - elif f.exp_vals[i-direction] != 0: - f.exp_vals[i-direction] = \ - deepcopy( - f.exp_vals[i-direction]) - 1 - # cop out for this...just use the - # first one. - else: - f.exp_vals[nonzero_elements[0]] = \ - deepcopy( - f.exp_vals[ - nonzero_elements[0]]) - 1 + f.exp_vals[last_nze] = deepcopy(f.exp_vals[last_nze]) - 1 + + else: #Otherwise a bit more complicated + direction = int(2*(randint(0,1)) - 1) #decide whether to go + or - + #print(i, direction, f.exp_vals) + if f.exp_vals[i+direction] != 0: # try 1st direction + f.exp_vals[i+direction] = deepcopy(f.exp_vals[i+direction]) - 1 + elif f.exp_vals[i-direction] != 0: #if not, try second direction + f.exp_vals[i-direction] = deepcopy(f.exp_vals[i-direction]) - 1 + else: #cop out for this...just use the first one. + f.exp_vals[nonzero_elements[0]] = deepcopy(f.exp_vals[nonzero_elements[0]]) - 1 if numpy.sum(sector_exp_vals) != s.num_amides: - print("ERROR: sector ", s.seq, - "does not have right number of amides,", s.num_amides) - print("Overlap not correctly calculated. Sector exp_vals:", - sector_exp_vals) + print("ERROR: sector ", s.seq, "does not have right number of amides,", s.num_amides) + print("Overlap not correctly calculated. Sector exp_vals:", sector_exp_vals) print("Method:", method) sys.exit(1) + # cast exp_grid to self.exp_grid for each sector - long_sector_exp_vals = numpy.zeros(len(self.exp_grid)).astype(int) + long_sector_exp_vals=numpy.zeros(len(self.exp_grid)).astype(int) # first and last are identical for each. - long_sector_exp_vals[0] = sector_exp_vals[0] - long_sector_exp_vals[-1] = sector_exp_vals[-1] + long_sector_exp_vals[0]=sector_exp_vals[0] + long_sector_exp_vals[-1]=sector_exp_vals[-1] - # Cast middle values + #Cast middle values - # something is wrong with these indicies...need to fix. - # long_sector -1 is getting overwritten - factor = float(len(sector_exp_vals)) \ - / (float(len(long_sector_exp_vals))) + #something is wrong with these indicies...need to fix. long_sector -1 is getting overwritten + factor = float(len(sector_exp_vals))/(float(len(long_sector_exp_vals))) - for i in range(1, len(sector_exp_vals)-1): - for j in range(1, len(long_sector_exp_vals)-1): + for i in range(1,len(sector_exp_vals)-1): + for j in range(1,len(long_sector_exp_vals)-1): if j * factor <= i and (j+1)*factor >= i: - ran_int = randint(0, 1) + ran_int = randint(0,1) long_sector_exp_vals[j+ran_int] += sector_exp_vals[i] break + #print(j*factor, i,long_sector_exp_vals[j], sector_exp_vals[i]) + + #print(s.seq, sector_exp_vals, method) + #print(s.seq, long_sector_exp_vals) if numpy.sum(long_sector_exp_vals) != s.num_amides: - print("ERROR: sector ", s.seq, - "does not have right number of amides,", s.num_amides) - print(sector_exp_vals, "not correctly cast to", - long_sector_exp_vals) + print("ERROR: sector ", s.seq, "does not have right number of amides,", s.num_amides) + print(sector_exp_vals, "not correctly cast to", long_sector_exp_vals) print("Method:", method) sys.exit(1) for n in range(s.start_res, s.end_res+1): - if self.seq[n] == 'P' or self.seq == 'p': - exp_seq[n] = 0 + if self.seq[n]=='P' or self.seq=='p': + exp_seq[n]=0 else: for i in range(len(self.exp_grid)): if long_sector_exp_vals[i] != 0: exp_seq[n] = i+1 - long_sector_exp_vals[i] -= 1 + long_sector_exp_vals[i] = long_sector_exp_vals[i] - 1 break + #print("SSSS", s.seq, sector_exp_vals, exp_seq[s.start_res:s.end_res+1]) + + #for f in frags: + #print(f.seq, f.exp_vals) + + #for n in range(len(self.exp_seq)): + #print(n,self.state.seq[n], exp_seq[n]) + #print(exp_seq) self.exp_seq = exp_seq return exp_seq + + #print("Fragment Amides:",s.num_amides,"Overlap:",overlap, "Diff:" , s.num_amides - numpy.sum(overlap), s.seq) + + ''' + + def calc_overlap_score(self, seq): + coverage = self.state.get_coverage(self.frags) + seclist = self.state.get_sectors().sort + sorted_sec = + + if coverage == 0 or self.state.seq[n] == 'P': + seq[n] = 0 + else + for f in self.frags: + if n > f.start_res + 2 and n < f.end_res: + + ''' + + def get_sigma(self, inval, sig): return inval*sig+0.05 - def gaussian_model(self, exp, model, sig): - return numpy.exp(-((model-exp)**2)/(2*sig**2)) \ - / (sig*numpy.sqrt(2*numpy.pi)) + def gaussian_model(self,exp,model,sig): + #print("hello") + return numpy.exp(-((model-exp)**2)/(2*sig**2))/(sig*numpy.sqrt(2*numpy.pi)) - def lognormal_model(self, exp, model, sig): - return numpy.exp(-((numpy.log(model)-numpy.log(exp))**2)/(2*sig**2)) \ - / (sig*numpy.sqrt(numpy.pi)*model) + def lognormal_model(self,exp,model,sig): + return numpy.exp(-( (numpy.log(model)-numpy.log(exp) )**2)/(2*sig**2))/(sig*numpy.sqrt(numpy.pi)*model) - def truncated_gaussian_factor(self, exp, a, b, sig): - return 1 / (0.5 * ( - scipy.special.erf((b-exp)/sig * numpy.sqrt(3.1415)) - - scipy.special.erf((a-exp)/sig * numpy.sqrt(3.1415)))) + def truncated_gaussian_factor(self,exp,a,b,sig): + return 1/ ( 0.5 * ( scipy.special.erf( (b-exp)/sig * numpy.sqrt(3.1415) ) - scipy.special.erf( (a-exp)/sig * numpy.sqrt(3.1415) ) ) ) def set_exp_grid(self, grid): - self.exp_grid = grid + self.exp_grid=grid + + def calc_exp_grid(self, num_bins, first_time_point=10, last_time_point=3600): - def calc_exp_grid(self, num_bins, first_time_point=10, - last_time_point=3600): if num_bins < 2: - raise Exception("get_exp_bins: num_bins must be greater than 1") - # slow_exp_val should satisfy - # 0.99 = exp(-10**slowest_rate_bin*first_time_point) - # fast_exp_val should satisfy - # 0.01 = exp(-10**fastest_rate_bin*last_time_point) + raise Exception("get_exp_bins: num_bins must be greater than 1" ) + # slow_exp_val should satisfy 0.99 = exp(-10**slowest_rate_bin*first_time_point) + # fast_exp_val should satisfy 0.01 = exp(-10**fastest_rate_bin*last_time_point) if first_time_point == 0.0: first_time_point = 1.0 - slowest_rate_bin = numpy.log10(-numpy.log(0.01)/first_time_point) - fastest_rate_bin = numpy.log10(-numpy.log(0.99)/last_time_point) + slowest_rate_bin=numpy.log10(-numpy.log(0.01)/first_time_point) + fastest_rate_bin=numpy.log10(-numpy.log(0.99)/last_time_point) - print("RATE BINS", slowest_rate_bin, fastest_rate_bin, - first_time_point, last_time_point) + print("RATE BINS", slowest_rate_bin, fastest_rate_bin, first_time_point, last_time_point) return numpy.linspace(fastest_rate_bin, slowest_rate_bin, num_bins) def get_sector_averaged_protection_values(self, seq, sectors): - sector_averaged_protection_values = [0.]*len(seq) + sector_averaged_protection_values=[0.]*len(seq) for s in sectors: - sum_exp = 0 + sum_exp=0 for n in range(s.start_res, s.end_res+1): - sum_exp = sum_exp+self.exp_grid[int(seq[n])-1] - frac_exp = sum_exp/s.num_amides + #print n, seq[n], sum_exp, self.exp_seq[n], self.exp_grid[int(seq[n])-1] + sum_exp=sum_exp+self.exp_grid[int(seq[n])-1] + frac_exp=sum_exp/s.num_amides for n in range(s.start_res, s.end_res+1): - if self.seq[n] != 'p' and self.seq[n] != 'P': - sector_averaged_protection_values[n] = frac_exp + if self.seq[n]!='p' and self.seq[n]!='P': + sector_averaged_protection_values[n]=frac_exp return sector_averaged_protection_values - def get_indiv_sector_averaged_protection_values(self, seq, sectors): - sector_averaged_protection_values = [0.]*len(sectors) + def get_indiv_sector_averaged_protection_values(self,seq, sectors): + sector_averaged_protection_values=[0.]*len(sectors) for snum in range(len(sectors)): s = sectors[snum] - sum_exp = 0 + sum_exp=0 for n in range(s.start_res, s.end_res+1): - sum_exp = sum_exp+self.exp_grid[int(seq[n])-1] - frac_exp = sum_exp/s.num_amides - sector_averaged_protection_values[snum] = frac_exp + #print n, seq[n], sum_exp, self.exp_seq[n], self.exp_grid[int(seq[n])-1] + sum_exp=sum_exp+self.exp_grid[int(seq[n])-1] + frac_exp=sum_exp/s.num_amides + sector_averaged_protection_values[snum]=frac_exp return sector_averaged_protection_values - def get_protection_values(self, seq): - protection_values = [0.]*len(seq) + def get_protection_values(self,seq): + protection_values=[0.]*len(seq) for n in range(len(seq)): - protection_values[n] = self.exp_grid[self.exp_seq[n]-1] + protection_values[n]=self.exp_grid[self.exp_seq[n]-1] return protection_values - def get_protection_bins(self, seq): - protection_values = [0.]*len(seq) + def get_protection_bins(self,seq): + protection_values=[0.]*len(seq) for n in range(len(seq)): - protection_values[n] = self.exp_seq[n] + protection_values[n]=self.exp_seq[n] return protection_values def get_frag_deut_comparison(self, frags, exp_seq="none"): - if exp_seq == "none": - exp_seq = self.exp_seq + if exp_seq=="none": + exp_seq=self.exp_seq for f in frags: for t in f.timepoints: - model_deut = 0 - grid = [] - for r in range(f.start_res+2-self.offset, - f.end_res-self.offset): - kr = 10**self.exp_grid[exp_seq[r]-1] - model_deut += (1-numpy.exp(-kr*t.time))*self.sat + model_deut=0 + grid=[] + #print len(grid), f.start_res, f.end_res + for r in range(f.start_res+2-self.offset, f.end_res-self.offset): + kr=10**self.exp_grid[exp_seq[r]-1] + model_deut=model_deut+(1-numpy.exp(-kr*t.time))*self.sat grid.append(exp_seq[r]-1) - # print(f.seq, t.time, grid, model_deut) + #print(f.seq, t.time, grid, model_deut) return model_deut def calc_exp_chiscores_from_file(self, frags, infile, outdir): self.import_model_deuteration_from_file(frags, infile, outdir) for f in frags: - chi = 0 - totd = 0 + chi=0 + totd=0 for t in f.timepoints: t.avg_sd_model() for d in t.deut: - chi = chi+(d-t.model_avg)**2/self.sigma - totd = totd+1 - # print f.seq, chi, t.time, t.model_avg, t.deut - f.chi = chi/totd + chi=chi+(d-t.model_avg)**2/self.sigma + totd=totd+1 + #print f.seq, chi, t.time, t.model_avg, t.deut + f.chi=chi/totd + #print f.seq, chi, [t.model_avg for t in f.timepoints], [t.model_avg for t in f.timepoints] def calc_exp_chiscores_from_gridvals(self, frags, infile, outdir): self.import_model_deuteration_from_gridvals(frags, infile, outdir) for f in frags: - chi = 0 - totd = 0 + chi=0 + totd=0 for t in f.timepoints: t.avg_sd_model() for d in t.deut: - chi = chi+(d-t.model_avg)**2/self.sigma - totd = totd+1 - # print f.seq, chi, t.time, t.model_avg, t.deut - f.chi = chi/totd + chi=chi+(d-t.model_avg)**2/self.sigma + totd=totd+1 + #print f.seq, chi, t.time, t.model_avg, t.deut + f.chi=chi/totd print(f.seq, chi) def get_model_distributions(self, frags): @@ -656,104 +658,112 @@ def get_model_distributions(self, frags): print(t.model) def get_model_deuteration(self, time, frag, exp_seq): - model_deut = 0 + model_deut=0 for r in range(frag.start_res+1, frag.end_res): - # print f.seq, t.time, r, self.seq[r] - if self.seq[r] != 'P' and self.seq[r] != 'p': - kr = 10**self.exp_grid[int(float(exp_seq[r])-1)] - model_deut = model_deut+(1-numpy.exp(-kr*time)) + #print f.seq, t.time, r, self.seq[r] + if self.seq[r] !='P' and self.seq[r] != 'p': + kr=10**self.exp_grid[int(float(exp_seq[r])-1)] + model_deut=model_deut+(1-numpy.exp(-kr*time)) return model_deut/frag.num_observable_amides*100 - def import_model_deuteration_from_file(self, frags, infile, firstline=1, - lastline=-1, append=False): - data = open(infile, "r") - if not append: - self.exp_models = [] + + def import_model_deuteration_from_file(self, frags, infile, firstline=1, lastline=-1, append=False): + data=open(infile,"r") + if append==False: + self.exp_models=[] for line in data: - # print line - exp = map(int, line.split(' ')) + #print line + exp=list(map(int, line.split(' '))) self.exp_models.append(exp) for f in frags: - # print f.seq + #print f.seq for t in f.timepoints: - model_deut = 0 - if not append: + model_deut=0 + if append==False: t.clear_model_deuteration_values() - append = True + append=True for r in range(f.start_res+1, f.end_res): - # print f.seq, t.time, r, self.seq[r] - if self.seq[r] != 'P' and self.seq[r] != 'p': - kr = 10**self.exp_grid[int(float(exp[r])-1)] - model_deut = model_deut+(1-numpy.exp(-kr*t.time)) - model_deut = model_deut/f.num_observable_amides*100 + #print f.seq, t.time, r, self.seq[r] + if self.seq[r] !='P' and self.seq[r] != 'p': + kr=10**self.exp_grid[int(float(exp[r])-1)] + model_deut=model_deut+(1-numpy.exp(-kr*t.time)) + #print f.seq, self.seq[r], r, f.start_res, r-f.start_res-1, (1-numpy.exp(-kr*t.time))*self.sat + model_deut=model_deut/f.num_observable_amides*100 t.add_model_deuteration_value(model_deut) + #print(f.seq, t.time, len(t.models), model_deut, t.get_model_avg(), t.get_model_sd(), [r.deut for r in t.replicates]) + #print f.seq, t.time, t.get_model_avg(), t.get_model_sd(), t.get_avg_sd() def import_model_deuteration_from_gridvals(self, frags, models): for f in frags: - timepoints = f.timepoints + timepoints=f.timepoints for t in timepoints: for m in models: - model_deut = 0 - for r in range(f.start_res+self.offset+2, - f.end_res+self.offset+1): - if self.seq[r] != 'P' and self.seq[r] != 'p': - kr = 10**self.exp_grid[int(float(m[r])-1)] - model_deut += \ - (1-numpy.exp(-kr*t.time))*t.replicates[0].sat - model_deut = model_deut/f.num_observable_amides*100 + model_deut=0 + for r in range(f.start_res+self.offset+2, f.end_res+self.offset+1): + #print r, f.seq, f.start_res+1, f.end_res, self.seq[r] + if self.seq[r] !='P' and self.seq[r] != 'p': + #print(r, self.seq[r], m[r], "|", m) + kr=10**self.exp_grid[int(float(m[r])-1)] + model_deut+=(1-numpy.exp(-kr*t.time))*t.replicates[0].sat + #print(">>", r, m[r], int(float(m[r])), self.exp_grid[int(float(m[r])-1)]) + #print f.seq, self.seq[r], r, f.start_res, r-f.start_res-1, (1-numpy.exp(-kr*t.time))*self.sat + #print(model_deut, t.time, t.sigma, f.seq, f.start_res, f.end_res, m[f.start_res+2+self.offset:f.end_res+1+self.offset], m, self.exp_grid) + model_deut=model_deut/f.num_observable_amides*100 t.add_model_deuteration_value(model_deut) + #t.avg_sd_model() + #print t.time, t.model_avg, t.model_sd, t.deut def import_models_from_gridvals(self, models, append=False): - if append == 'False': - self.exp_models = [] + if append=='False': + self.exp_models=[] for i in models: self.exp_models.append(i) + def import_models_from_file(self, modelfile, append='False'): - data = open(modelfile, "r") - if append == 'False': - self.exp_models = [] + data=open(modelfile,"r") + if append=='False': + self.exp_models=[] for line in data: - self.exp_models.append(map(int, line.split(' '))) + self.exp_models.append(list(map(int, line.split(' ')))) def calc_model_scores(self, frags=None, sig=1.0, error_model="gaussian"): - if self.exp_models == []: + if self.exp_models==[]: print("No models to calculate score") return if frags is None: - frags = self.frags + frags=self.frags - self.model_scores = [] + self.model_scores=[] for m in self.exp_models: - self.exp_seq = deepcopy(m) - score = self.calculate_bayesian_score(frags, sig, error_model) + self.exp_seq=deepcopy(m) + score=self.calculate_bayesian_score(frags,sig,error_model) self.model_scores.append(score) return self.model_scores def import_scores(self, scorefile): try: - data = open(scorefile, "r") + data=open(scorefile,"r") for line in data: self.scores.append(line.split().strip()) - except: # noqa: E722 + except: for i in scorefile: self.scores.append(i) - def get_sigma_score(self, sigma0, sig): - return (2*sigma0 / (numpy.sqrt(numpy.pi) * sig**2)) \ - * numpy.exp(-sigma0**2 / sig**2) + def get_sigma_score(self,sigma0,sig): + return (2*sigma0/ (numpy.sqrt(numpy.pi) * sig**2 ))*numpy.exp(-sigma0**2 /sig**2) def get_model_average(self): - # print(len(self.exp_models)) + #print(len(self.exp_models)) if len(self.exp_models) > 0: - modelavg = numpy.average(numpy.array(self.exp_models), axis=0) + #print(len(self.exp_models), len(self.exp_models[0]), numpy.array(self.exp_models).shape, self.exp_models[0]) + modelavg=numpy.average(numpy.array(self.exp_models), axis=0) return modelavg else: print("No models to get the average of =(") def unimodal_prior(self, omj, sigma0=0.1): - return (2.0 * omj / (numpy.sqrt(numpy.pi) * omj**2)) \ - * numpy.exp(-sigma0**2 / omj**2) + return (2.0 * omj / ( numpy.sqrt(numpy.pi) * omj**2 )) * numpy.exp(-sigma0**2 / omj**2) def get_nearest_value(self, value, array): ''' @@ -761,106 +771,97 @@ def get_nearest_value(self, value, array): ''' return (numpy.abs(array-value)).argmin(), - def calculate_marginal_probabilities(self, - error_model="truncated_gaussian"): - # Calculates marginal probabilities given grids of omega, sigma0, - # fmod, and fexp + def calculate_marginal_probabilities(self, error_model="truncated_gaussian"): + # Calculates marginal probabilities given grids of omega, sigma0, fmod, and fexp - self.prob_grid = numpy.zeros( - (len(self.sigma0_grid), len(self.fmod_grid), len(self.fexp_grid))) + self.prob_grid=numpy.zeros( (len(self.sigma0_grid), len(self.fmod_grid), len(self.fexp_grid)) ) for fm in range(len(self.fmod_grid)): - fmod = self.fmod_grid[fm] + fmod=self.fmod_grid[fm] for f in range(len(self.fexp_grid)): - fexp = self.fexp_grid[f] + fexp=self.fexp_grid[f] for s in range(len(self.sigma0_grid)): - sigma0 = self.sigma0_grid[s] - cumul = 0 - - for j in range(1, len(self.omega_grid)): - # We're going to calculate the likelihoods - # at each point - omj = self.omega_grid[j] - omjm1 = self.omega_grid[j-1] + sigma0=self.sigma0_grid[s] + cumul=0 + + for j in range(1,len(self.omega_grid)): + # We're going to calculate the likelihoods at each point + omj=self.omega_grid[j] + omjm1=self.omega_grid[j-1] priorj = self.unimodal_prior(omj, sigma0) priorjm1 = self.unimodal_prior(omjm1, sigma0) - dom = omj-omjm1 - - if error_model == "truncated_gaussian": - # If we are using the truncated gaussian, - # calculate those pre-factors - # use the truncated gaussian factors at - # -0.1 and 1.2 - factor = self.truncated_gaussian_factor( - fexp, -10, 120, omj) - factorm1 = self.truncated_gaussian_factor( - fexp, -10, 120, omjm1) + dom=omj-omjm1 + + + if error_model=="truncated_gaussian": + # If we are using the truncated gaussian, calculate those pre-factors + # use the truncated gaussian factors at -0.1 and 1.2 + factor=self.truncated_gaussian_factor(fexp,-10,120,omj) + factorm1=self.truncated_gaussian_factor(fexp,-10,120,omjm1) else: - factor = 1.0 - factorm1 = 1.0 + factor=1.0 + factorm1=1.0 - # Calculate the likelihood at each point - pj = self.gaussian_model( - fexp, fmod, omj)*factor*priorj - pjm1 = self.gaussian_model( - fexp, fmod, omjm1)*factorm1*priorjm1 - cumul = cumul+(pj+pjm1)/2.0/dom - self.prob_grid[s, fm, f] = -1.0 * numpy.log(cumul) + # Calculate the likelihood at each point + pj=self.gaussian_model(fexp,fmod,omj)*factor*priorj + pjm1=self.gaussian_model(fexp,fmod,omjm1)*factorm1*priorjm1 + + cumul = cumul+(pj+pjm1)/2.0/dom; + #print(pj, pjm1, omj, omjm1, sigma0, "| ", factor, factorm1, "| ", priorj, priorjm1, cumul) + #print(fexp, fmod, fexp-fmod, sigma0, "| ", cumul, -1.0 * numpy.log(cumul)) + #print j,fexp,fmod,omj,(pj+pjm1)/2.0/dom,cumul, cumul2, cumul2/cumul + self.prob_grid[s,fm,f]=-1.0 * numpy.log(cumul) return self.prob_grid - def get_marginalized_sigma_probabilities(self, sigma0, - error_model="gaussian"): + def get_marginalized_sigma_probabilities(self,sigma0,error_model="gaussian"): ''' - Precomputes marginal probabilities for sigma levels using a - cauchy distribution - @param delfmod - difference between forward model and noise model in - delta pct-D units to normalize for fragment length (num amides) - @param sigma0 - Estimated standard deviation of the instrument - in pct-D units + Precomputes marginal probabilities for sigma levels using a cauchy distribution + @param delfmod - difference between forward model and noise model in delta pct-D + units to normalize for fragment length (num amides) + @param sigma0 - Estimated standard deviation of the instrument in pct-D units ''' - self.prob_grid = numpy.zeros((len(self.fmod_grid), - len(self.fexp_grid))) + self.prob_grid=numpy.zeros( (len(self.fmod_grid), len(self.fexp_grid)) ) + for fm in range(len(self.fmod_grid)): - fmod = self.fmod_grid[fm] + fmod=self.fmod_grid[fm] # Unimodal Distributions for sigma # for f in range(len(self.fexp_grid)): - fexp = self.fexp_grid[f] - if error_model == "truncated_gaussian": - # If we are using the truncated gaussian, calculate - # those pre-factors + fexp=self.fexp_grid[f] + if error_model=="truncated_gaussian": + # If we are using the truncated gaussian, calculate those pre-factors # use the truncated gaussian factors at -0.1 and 1.2 - factor = self.truncated_gaussian_factor( - fexp, -10, 120, omj) - factorm1 = self.truncated_gaussian_factor( - fexp, -10, 120, omjm1) + factor=self.truncated_gaussian_factor(fexp,-10,120,omj) + factorm1=self.truncated_gaussian_factor(fexp,-10,120,omjm1) + o=1 else: - factor = 1.0 - factorm1 = 1.0 + factor=1.0 + factorm1=1.0 for j in range(len(self.omega_grid)): # We're going to calculate the likelihoods at each point - omj = self.omega_grid[j] - omjm1 = self.omega_grid[j-1] + omj=self.omega_grid[j] + omjm1=self.omega_grid[j-1] priorj = self.unimodal_prior(omj, sigma0) priorjm1 = self.unimodal_prior(omjm1, sigma0) - dom = omj-omjm1 + dom=omj-omjm1 - cumul = 0 + cumul=0 # Calculate the likelihood at each point - pj = self.gaussian_model( - fexp, fmod, omj)*factor*priorj - pjm1 = self.gaussian_model( - fexp, fmod, omjm1)*factorm1*priorjm1 + pj=self.gaussian_model(fexp,fmod,omj)*factor*priorj + pjm1=self.gaussian_model(fexp,fmod,omjm1)*factorm1*priorjm1 + + cumul = cumul+(pj+pjm1)/2.0/dom; - cumul = cumul+(pj+pjm1)/2.0/dom + #print(j, fexp, fmod, sigma0, omj,(pj+pjm1)/2.0/dom,cumul) - self.prob_grid[fm, f] = -1.0 * numpy.log(cumul) + self.prob_grid[fm,f]=-1.0 * numpy.log(cumul) + #print "OMEGA: ",fexp,fmod,omega0,cumul,cumul2,cumul2/cumul return self.prob_grid diff --git a/pyext/src/plots.py b/pyext/src/plots.py index 559034a..6db5409 100644 --- a/pyext/src/plots.py +++ b/pyext/src/plots.py @@ -3,46 +3,47 @@ from __future__ import print_function import os +import system_setup +import hdx_models import numpy import math import analysis -from pylab import arange -import matplotlib +from pylab import * import matplotlib.pyplot as plt +import scipy from scipy.stats import gaussian_kde +#from numpy.random import normal +from scipy.integrate import simps + cdict1 = {'red': ((0.0, 0.0, 0.0), - (0.45, 0.0, 0.0), - (1.0, 1.0, 1.0)), - 'green': ((0.0, 0.0, 0.0), - (0.25, 0.0, 0.0), - (0.50, 0.0, 0.0), - (0.75, 0.0, 0.0), - (1.0, 0.0, 0.0)), - 'blue': ((0.0, 1.0, 1.0), - (0.55, 0.0, 0.0), - (1.0, 0.0, 0.0))} + (0.45, 0.0, 0.0), + (1.0, 1.0, 1.0)), + 'green': ((0.0, 0.0, 0.0), + (0.25, 0.0, 0.0), + (0.50, 0.0, 0.0), + (0.75, 0.0, 0.0), + (1.0, 0.0, 0.0)), + 'blue': ((0.0, 1.0, 1.0), + (0.55, 0.0, 0.0), + (1.0, 0.0, 0.0))} def roundup(i, scale): return int(math.ceil(i/float(scale)))*scale - def find_minmax(lists): - flatlist = [item for sublist in lists for item in sublist] - minmax = (math.floor(min(numpy.array(flatlist))), - math.ceil(max(numpy.array(flatlist)))) + flatlist=[item for sublist in lists for item in sublist] + minmax=(math.floor(min(np.array(flatlist))), math.ceil(max(np.array(flatlist))) ) return minmax - def calculate_histogram(list, bins): # Given a list of values, calculate a histogram return numpy.histogram(numpy.array(list), bins=bins) -def plot_apo_lig_dhdx(model, show_plot=True, save_plot=False, - outfile="dhdx.png", outdir=None, noclobber=True): +def plot_apo_lig_dhdx(model, show_plot=True, save_plot=False, outfile="dhdx.png", outdir=None, noclobber=True): """ Takes a sampled model and plots the delta HDX (ligand - Apo) as horizontal line plots for each liganded state Also, outputs raw data for each state to individual dat files @@ -59,69 +60,60 @@ def plot_apo_lig_dhdx(model, show_plot=True, save_plot=False, elif noclobber: raise Exception("Output directory ", outdir, " already exists") - fig = plt.figure(figsize=(8, 2*(len(model.states)-1)), dpi=100) - ax = fig.add_axes([0.1, 0.5, 0.8, 0.4]) + fig=plt.figure(figsize=(8,2*(len(model.states)-1)), dpi=100) + ax=fig.add_axes([0.1,0.5,0.8,0.4]) # calculate bounds # Create color bar - my_cmap = matplotlib.colors.LinearSegmentedColormap( - 'my_colormap', cdict1, 256) + my_cmap=matplotlib.colors.LinearSegmentedColormap('my_colormap',cdict1,256) # Normalize colormap to delta log(k) = 5 as maxes cNorm = matplotlib.colors.Normalize(vmin=-5., vmax=5.) scalarMap = matplotlib.cm.ScalarMappable(norm=cNorm, cmap=my_cmap) # Apo (reference) state is always state 0 - apo_state = model.states[0] + apo_state=model.states[0] # Get plot sizes nres = len(apo_state.seq) # Allow these possible xtics. Want about 4-6 per plot - p_xtics = [5, 10, 20, 50, 100, 200] + p_xtics = [5,10,20,50,100,200] p = list((numpy.array(p_xtics)-nres/4.)**2) xtics = p_xtics[p.index(min(p))] print(xtics) maxx = roundup(nres, xtics) - # Calculate avg and SD log(HDX rate) for apo state; returns vectors of - # length nres - avg_apo, sd_apo = analysis.get_average_sector_values( - apo_state.exp_model.exp_models, apo_state) + # Calculate avg and SD log(HDX rate) for apo state; returns vectors of length nres + avg_apo, sd_apo=analysis.get_average_sector_values(apo_state.exp_model.exp_models, apo_state) - # This is the ligand number being processed. It is the y-value of that - # state's bar in the plot. - nlig = -1 + nlig=-1 # This is the ligand number being processed. It is the y-value of that state's bar in the plot. for s in model.states[1:]: # Find the high and low grid values and calculate the # tolerance for sectors that are too fast/slow to observe - grid_val_hi = s.exp_model.exp_grid[-1] - grid_val_lo = s.exp_model.exp_grid[0] + grid_val_hi=s.exp_model.exp_grid[-1] + grid_val_lo=s.exp_model.exp_grid[0] tol = 0.05*(grid_val_hi-grid_val_lo) # Initialize dat file sname = s.state_name outdat = outdir + sname + "_dhdx.dat" - f = open(outdat, 'w') - f.write("Res# Res state_name dhdx dhdx_z avg_lig avg_apo " - "sd_lig sd_apo flag\n") + f = open(outdat,'w') + f.write("Res# Res state_name dhdx dhdx_z avg_lig avg_apo sd_lig sd_apo flag\n") - nlig = nlig+1 # Increment ligand number - xin = [] # Stores residue numbers (x-values of chart) - yin = (nlig, 1.0) # y value is nlig, bar width is 1.0 - color = [] # stores color value - associated with xin + nlig=nlig+1 # Increment ligand number + xin=[] # Stores residue numbers (x-values of chart) + yin=(nlig,1.0) # y value is nlig, bar width is 1.0 + color=[] # stores color value - associated with xin # Calculate avg and SD log(HDX rate) for liganded state - avg_lig, sd_lig = analysis.get_average_sector_values( - s.exp_model.exp_models, s) + avg_lig, sd_lig=analysis.get_average_sector_values(s.exp_model.exp_models, s) # Calculate Z-score - z = analysis.calculate_zscore( - apo_state.exp_model.exp_models, s.exp_model.exp_models, - apo_state, s) + z=analysis.calculate_zscore(apo_state.exp_model.exp_models, s.exp_model.exp_models, apo_state, s) # Add state name to plot plt.text(nres+1, nlig+0.45, s.state_name, fontsize=14) @@ -129,7 +121,7 @@ def plot_apo_lig_dhdx(model, show_plot=True, save_plot=False, for n in range(len(model.seq)): # Write values to file dhdx = avg_lig[n] - avg_apo[n] - flag = "" + flag="" if avg_apo[n] < grid_val_lo + tol and avg_apo[n] != 0: flag += "*LOW_APO_VAL" if avg_lig[n] < grid_val_lo + tol and avg_lig[n] != 0: @@ -139,55 +131,54 @@ def plot_apo_lig_dhdx(model, show_plot=True, save_plot=False, if avg_lig[n] > grid_val_hi - tol and avg_lig[n] != 0: flag += "*HI_LIG_VAL" - if avg_lig[n] == 0.0 and avg_apo[n] == 0.0: + if avg_lig[n]==0.0 and avg_apo[n]==0.0: flag = "" + #print("NNN", avg_apo[n], avg_lig[n], grid_val_lo, "|", grid_val_lo + tol, grid_val_hi -tol, flag) - f.write("%i %s %s %f %f %f %f %f %f %s\n" - % (n, model.seq[n], s.state_name, dhdx, z[n], avg_lig[n], - avg_apo[n], sd_lig[n], sd_apo[n], flag)) + + f.write("%i %s %s %f %f %f %f %f %f %s\n" % + (n,model.seq[n], s.state_name, dhdx, z[n], avg_lig[n], avg_apo[n], sd_lig[n], sd_apo[n], flag)) # Calculate residue color - rgba = scalarMap.to_rgba(avg_lig[n]-avg_apo[n]) - lst = list(rgba) + rgba=scalarMap.to_rgba(avg_lig[n]-avg_apo[n]) #my_cmap((avg_lig[n]-avg_apo[n])/3.0) + lst=list(rgba) # calculate saturation (|Z|>3.0) is full saturation, linear scale - if abs(z[n]) > 3.0: - lst[3] = 1.0 + if abs(z[n])>3.0: + lst[3]=1.0 else: - lst[3] = abs(z[n])/3.0 + lst[3]=abs(z[n])/3.0 # Append the xin and color tuples - xin.append((n+model.offset, 1)) - rgba = tuple(lst) + xin.append((n+model.offset,1)) + rgba=tuple(lst) color.append(rgba) # Plot bar and close outfile - ax.broken_barh(xin, yin, color=color, lw=0) + ax.broken_barh(xin,yin, color=color, lw=0) f.close() ax.grid(True) - ax.get_xaxis().set_ticks(range(0, maxx, xtics)) + ax.get_xaxis().set_ticks(range(0,maxx,xtics)) ax.get_yaxis().set_ticks([]) - ax.set_xlim([0, maxx]) - ax.set_ylim([0, nlig+1]) + ax.set_xlim([0,maxx]) + ax.set_ylim([0,nlig+1]) ax.set_xlabel('Residue Number', fontsize=10) ax.set_ylabel('Ligand States', fontsize=10) ax.set_title('Delta HDX | Target: ' + model.target_name) # Add axes to bottom [left, bottom, width, height] for colormap - cax = fig.add_axes([0.3, 0.2, 0.4, 0.02]) - cb = matplotlib.colorbar.ColorbarBase( - cax, cmap=my_cmap, norm=cNorm, spacing='proportional', - orientation='horizontal') + cax = fig.add_axes([0.3,0.2,0.4,0.02]) + cb = matplotlib.colorbar.ColorbarBase(cax, cmap=my_cmap, norm=cNorm, spacing='proportional', orientation='horizontal') cb.set_label('DHDX log(k)', fontsize=10) - if show_plot: + if show_plot==True: fig.show() - if save_plot: + if save_plot==True: fig.savefig(outdir + outfile, bbox_inches=0) -def plot_fragment_chi_values(state, sig="model", outfile=None, - show_plot=False, outdir="./output/"): +def plot_fragment_chi_values(state, sig="model", outfile=None, show_plot=False, outdir="./output/"): + if outdir is None: outdir = "./output/fragment_chi_plots/" else: @@ -195,32 +186,34 @@ def plot_fragment_chi_values(state, sig="model", outfile=None, if not os.path.exists(outdir): os.makedirs(outdir) - frags = state.frags - maxchi = 50 + frags=state.frags + maxchi=50 + #for f in frags: + # if f.get_chi_value(sig) > maxchi: + # maxchi = f.chi - nres = len(state.seq) - fig, ax = plt.subplots(figsize=(12, 6)) + nres=len(state.seq) + fig, ax=plt.subplots(figsize=(12,6)) max_overlap = 50 color_map = plt.get_cmap('gnuplot') cNorm = matplotlib.colors.Normalize(vmin=0, vmax=maxchi) scalarMap = matplotlib.cm.ScalarMappable(norm=cNorm, cmap=color_map) sorted_frags = sorted(frags, key=lambda x: x.start_res) - y_val = 1 - maxi = 0 - end_res = [0]*max_overlap + y_val=1 + maxi=0 + end_res=[0]*max_overlap for f in sorted_frags: - i = 1 + i=1 while i < max_overlap: if f.start_res > end_res[i]: - y_val = i - end_res[i] = f.end_res + y_val=i + end_res[i]=f.end_res break - i = i+1 - if i > maxi: - maxi = i - colorVal = scalarMap.to_rgba(f.chi) - ax.hlines(y_val, int(f.start_res)-0.5, int(f.end_res)+0.5, - color=colorVal, lw=25) + i=i+1 + if i> maxi: + maxi=i + colorVal=scalarMap.to_rgba(f.chi) + ax.hlines(y_val, int(f.start_res)-0.5, int(f.end_res)+0.5, color=colorVal, lw=25) if nres < 100: ax.text(int(f.start_res)+1, y_val-0.3, f.seq) ax.text(int(f.start_res)-0.7, y_val, f.start_res) @@ -230,113 +223,103 @@ def plot_fragment_chi_values(state, sig="model", outfile=None, ax.text(1, maxi+0.5, "Max chi value = " + str(maxchi)) ax.set_title("Individual Fragment Fits to Model - " + state.state_name) ax.set_xlabel('Residue Number') - ax.set_ylim([0, maxi+1]) - ax.set_xlim([0, nres+1]) + ax.set_ylim([0,maxi+1]) + ax.set_xlim([0,nres+1]) - cax = fig.add_axes([0.95, 0.2, 0.02, 0.6]) - cb = matplotlib.colorbar.ColorbarBase( - cax, cmap=color_map, norm=cNorm, spacing='proportional') + + cax = fig.add_axes([0.95,0.2,0.02,0.6]) + cb = matplotlib.colorbar.ColorbarBase(cax, cmap=color_map, norm=cNorm, spacing='proportional') cb.set_label('Fragment Chi') if outfile is None: outfile = state.state_name+"_fragment_chi_fits.png" - if outfile is None: + if outfile==None: plt.show() - elif not show_plot: + elif show_plot==False: plt.savefig(outdir+outfile, bbox_inches=0, format="png") else: plt.savefig(outdir+outfile, bbox_inches=0, format="png") plt.show() +def plot_fragment_avg_model_fits(state, sig, frags="all", outfile=None, show_plot=False): + if frags=="all": + frags=state.frags -def plot_fragment_avg_model_fits(state, sig, frags="all", outfile=None, - show_plot=False): - if frags == "all": - frags = state.frags - - fig = plt.figure() - ax = plt.gca() + fig=plt.figure() + ax=plt.gca() for f in frags: - x = [] - yavg = [] - yerror = [] - # print f.seq, x + + x=[] + yavg=[] + yerror=[] + #print f.seq, x for t in f.timepoints: - # print t.model + #print t.model if t.get_model_avg() is not None and t.get_model_sd() is not None: x.append(t.time) - xt = [int(t.time)]*len(t.replicates) - yt = [float(r.deut) for r in t.replicates] + xt=[int(t.time)]*len(t.replicates) + yt=[float(r.deut) for r in t.replicates] plt.scatter(xt, yt) - chi = f.get_chi_value(sig) + chi=f.get_chi_value(sig) for t in f.timepoints: if t.get_model_avg() is not None and t.get_model_sd() is not None: yavg.append(t.model_avg) yerror.append(t.model_sd) - # print f.seq, t.time, t.model_avg, t.model_sd, len(t.models) - # plt.show() + #print f.seq, t.time, t.model_avg, t.model_sd, len(t.models) + #plt.show() plt.errorbar(x, yavg, yerr=yerror) ax.set_xscale('log') - # ax.set_xlim=(1,3600) - plt.axis = (1, 3600, 0, 100) - plt.text(1, 1, chi) - fig.title = (str(f.seq) + "_" + str(chi)) - if outfile is None: + #ax.set_xlim=(1,3600) + plt.axis=(1,3600,0,100) + plt.text(1,1,chi) + fig.title=(str(f.seq)+"_"+str(chi)) + if outfile==None: plt.show() - elif not show_plot: + elif show_plot==False: plt.savefig(outfile, bbox_inches=0) else: plt.show() plt.savefig(outfile, bbox_inches=0) - def plot_model_scores(kinetic_model, show_plot=True, outfile=None): '''plots model vs. score for a sorted list''' try: scores = kinetic_model.model_scores - except: # noqa: E722 + except: scores = numpy.sort(kinetic_model.calc_model_scores()) fig = plt.figure() - x = range(len(scores)) + x=range(len(scores)) plt.xlabel("Model Rank") plt.ylabel("Model Score") - plt.scatter(x, scores) + plt.scatter(x,scores) - fig.title = (str(kinetic_model.state.state_name) + "_top_models") + fig.title=(str(kinetic_model.state.state_name)+"_top_models") if outfile is not None: plt.savefig(outfile, bbox_inches=0) - if show_plot: + if show_plot==True: plt.show() -def plot_2state_fragment_avg_model_fits( - state1, state2, sig, num_best_models=100, outdir=None, - write_file=True, show_plot=False): - '''For two states (apo and 1 ligand, e.g.), plot the fits to model - for all fragments. - Output an individual plot for each fragment fit +def plot_2state_fragment_avg_model_fits(state1, state2, sig, num_best_models=100, outdir=None, write_file=True, show_plot=False): + ''' For two states (apo and 1 ligand, e.g.), plot the fits to model for all fragments. + Output an individual plot for each fragment fit ''' for s in [state1, state2]: - bsm, scores = analysis.get_best_scoring_models( - s.modelfile, s.scorefile, num_best_models=num_best_models, - prefix=s.state_name, write_file=write_file) + bsm, scores=analysis.get_best_scoring_models(s.modelfile, s.scorefile, num_best_models=num_best_models, prefix=s.state_name, write_file=write_file) s.exp_model.import_model_deuteration_from_gridvals(s.frags, bsm) s.exp_model.import_models_from_gridvals(bsm) - # takes a model and score file and writes a new model file with the best - # X scoring models. - # This new file can then be imported into an HDXModel class for analysis + #takes a model and score file and writes a new model file with the best X scoring models. + #This new file can then be imported into an HDXModel class for analysis if outdir is None: - outdir = "./output/fragment_fit-to-data_" + str(state1.state_name) \ - + "-" + str(state2.state_name) + "/" + outdir = "./output/fragment_fit-to-data_"+ str(state1.state_name) + "-" + str(state2.state_name) +"/" else: - outdir = outdir + "fragment_fit-to-data_" + str(state1.state_name) \ - + "-" + str(state2.state_name) + "/" + outdir = outdir + "fragment_fit-to-data_"+ str(state1.state_name) + "-" + str(state2.state_name) +"/" if not os.path.exists(outdir): os.makedirs(outdir) @@ -344,66 +327,182 @@ def plot_2state_fragment_avg_model_fits( # Get the state2 frag. If this is not there, it will return "" s2frag = state2.get_frag(s1frag.seq, s1frag.start_res) - outfile = str(s1frag.seq) + "_model_fits.png" + outfile=str(s1frag.seq) + "_model_fits.png" - fig = plt.figure() - ax = plt.gca() - x = [] - yavg = [] - yerror = [] - # print f.seq, x + fig=plt.figure() + ax=plt.gca() + x=[] + yavg=[] + yerror=[] + #print f.seq, x + s1avg=numpy.average(state1.exp_model.get_model_average()[s1frag.start_res+1:s1frag.end_res+1]) for t in s1frag.timepoints: - # print t.model + #print t.model if t.get_model_avg() is not None and t.get_model_sd() is not None: x.append(t.time) - xt = [int(t.time)]*len(t.replicates) - yt = [float(r.deut) for r in t.replicates] + xt=[int(t.time)]*len(t.replicates) + yt=[float(r.deut) for r in t.replicates] plt.scatter(xt, yt, c='b') - chi = s1frag.get_chi_value(sig) + chi=s1frag.get_chi_value(sig) for t in s1frag.timepoints: if t.get_model_avg() is not None and t.get_model_sd() is not None: yavg.append(t.model_avg) yerror.append(t.model_sd) + #print(s1frag.seq, t.time, t.model_avg, t.model_sd, len(t.models)) plt.errorbar(x, yavg, yerr=yerror, c='b') if s2frag != "": - x2 = [] - yavg2 = [] - yerror2 = [] - # print(s2frag, len(s2frag.timepoints)) + s2avg=numpy.average(state2.exp_model.get_model_average()[s2frag.start_res+1:s2frag.end_res]) + x2=[] + yavg2=[] + yerror2=[] + #print(s2frag, len(s2frag.timepoints)) for t in s2frag.timepoints: # Plot experimental data - # print(t.models) - if (t.get_model_avg() is not None - and t.get_model_sd() is not None): + #print(t.models) + if t.get_model_avg() is not None and t.get_model_sd() is not None: x2.append(t.time) - xt = [int(t.time)]*len(t.replicates) - yt = [float(r.deut) for r in t.replicates] + xt=[int(t.time)]*len(t.replicates) + yt=[float(r.deut) for r in t.replicates] plt.scatter(xt, yt, c='r') - chi2 = s2frag.get_chi_value(sig) + chi2=s2frag.get_chi_value(sig) for t in s2frag.timepoints: - # Plot model average and SD errorbars - if (t.get_model_avg() is not None - and t.get_model_sd() is not None): + #Plot model average and SD errorbars + if t.get_model_avg() is not None and t.get_model_sd() is not None: yavg2.append(t.model_avg) yerror2.append(t.model_sd) + #print f.seq, t.time, t.model_avg, t.model_sd, len(t.models) plt.errorbar(x2, yavg2, yerr=yerror2, c='r') avg = sum(yavg2)/len(yavg2) - fig.title = (str(s1frag.seq)+"_"+str(chi)+"_"+str(chi2)) + fig.title=(str(s1frag.seq)+"_"+str(chi)+"_"+str(chi2)) else: - fig.title = (str(s1frag.seq)+"_"+str(chi)) + fig.title=(str(s1frag.seq)+"_"+str(chi)) ax.set_xscale('log') ax.set_xlabel('Time') ax.set_ylabel('%D Incorporation') - plt.axis = (1, 3600, 0, 100) - plt.text(2, avg+10, s1frag.seq) + #ax.set_xlim=(1,3600) + plt.axis=(1,3600,0,100) + plt.text(2,avg+10,s1frag.seq) + #plt.text(2,avg+5,"Apo chi:"+str(chi)) + #plt.text(2,avg,"Lilly Diff: "+str(lilly_diff)) + #plt.text(2,avg-5,"Sali Diff: "+str(sali_diff)) + #plt.text(1,avg-5,"Sali Exp Diff: "+str(exp_diff)) + + if show_plot==False: + plt.savefig(outdir+outfile, bbox_inches=0, format="png") + else: + plt.show() + plt.savefig(outdir+outfile, bbox_inches=0, format="png") + plt.close() + fig.clear() - if not show_plot: + +def plot_2state_fragment_avg_model_fits_num_model_per_state(state1, state2, sig, num_best_models=[100, 100], + outdir=None, write_file=True, show_plot=False): + ''' + For two states (apo and 1 ligand, e.g.), plot the fits to model for all fragments. + num_best_models is an array with 2 elements contaning the number of best models for each state + Output an individual plot for each fragment fit + ''' + + bsm, scores = analysis.get_best_scoring_models(state1.modelfile, state1.scorefile, + num_best_models=num_best_models[0], + prefix=state1.state_name, write_file=write_file) + state1.exp_model.import_model_deuteration_from_gridvals(state1.frags, bsm) + state1.exp_model.import_models_from_gridvals(bsm) + + bsm, scores = analysis.get_best_scoring_models(state2.modelfile, state2.scorefile, + num_best_models=num_best_models[1], + prefix=state2.state_name, write_file=write_file) + state2.exp_model.import_model_deuteration_from_gridvals(state2.frags, bsm) + state2.exp_model.import_models_from_gridvals(bsm) + + #takes a model and score file and writes a new model file with the best X scoring models. + #This new file can then be imported into an HDXModel class for analysis + + if outdir is None: + outdir = "./output/fragment_fit-to-data_"+ str(state1.state_name) + "-" + str(state2.state_name) +"/" + else: + outdir = outdir + "fragment_fit-to-data_"+ str(state1.state_name) + "-" + str(state2.state_name) +"/" + if not os.path.exists(outdir): + os.makedirs(outdir) + + for s1frag in state1.frags: + # Get the state2 frag. If this is not there, it will return "" + s2frag = state2.get_frag(s1frag.seq, s1frag.start_res) + + outfile=str(s1frag.seq) + "_model_fits.png" + + fig=plt.figure() + ax=plt.gca() + x=[] + yavg=[] + yerror=[] + #print f.seq, x + s1avg=numpy.average(state1.exp_model.get_model_average()[s1frag.start_res+1:s1frag.end_res+1]) + + for t in s1frag.timepoints: + #print t.model + if t.get_model_avg() is not None and t.get_model_sd() is not None: + x.append(t.time) + xt=[int(t.time)]*len(t.replicates) + yt=[float(r.deut) for r in t.replicates] + plt.scatter(xt, yt, c='b') + chi=s1frag.get_chi_value(sig) + + for t in s1frag.timepoints: + if t.get_model_avg() is not None and t.get_model_sd() is not None: + yavg.append(t.model_avg) + yerror.append(t.model_sd) + #print(s1frag.seq, t.time, t.model_avg, t.model_sd, len(t.models)) + plt.errorbar(x, yavg, yerr=yerror, c='b') + + if s2frag != "": + s2avg=numpy.average(state2.exp_model.get_model_average()[s2frag.start_res+1:s2frag.end_res]) + x2=[] + yavg2=[] + yerror2=[] + #print(s2frag, len(s2frag.timepoints)) + for t in s2frag.timepoints: + # Plot experimental data + #print(t.models) + if t.get_model_avg() is not None and t.get_model_sd() is not None: + x2.append(t.time) + xt=[int(t.time)]*len(t.replicates) + yt=[float(r.deut) for r in t.replicates] + plt.scatter(xt, yt, c='r') + chi2=s2frag.get_chi_value(sig) + + for t in s2frag.timepoints: + #Plot model average and SD errorbars + if t.get_model_avg() is not None and t.get_model_sd() is not None: + yavg2.append(t.model_avg) + yerror2.append(t.model_sd) + #print f.seq, t.time, t.model_avg, t.model_sd, len(t.models) + plt.errorbar(x2, yavg2, yerr=yerror2, c='r') + avg = sum(yavg2)/len(yavg2) + + fig.title=(str(s1frag.seq)+"_"+str(chi)+"_"+str(chi2)) + else: + fig.title=(str(s1frag.seq)+"_"+str(chi)) + + ax.set_xscale('log') + ax.set_xlabel('Time') + ax.set_ylabel('%D Incorporation') + #ax.set_xlim=(1,3600) + plt.axis=(1,3600,0,100) + plt.text(2,avg+10,s1frag.seq) + #plt.text(2,avg+5,"Apo chi:"+str(chi)) + #plt.text(2,avg,"Lilly Diff: "+str(lilly_diff)) + #plt.text(2,avg-5,"Sali Diff: "+str(sali_diff)) + #plt.text(1,avg-5,"Sali Exp Diff: "+str(exp_diff)) + + if show_plot==False: plt.savefig(outdir+outfile, bbox_inches=0, format="png") else: plt.show() @@ -412,24 +511,23 @@ def plot_2state_fragment_avg_model_fits( fig.clear() -def get_cdf(ax, data, pos, bp=False): +def get_cdf(ax,data,pos, bp=False): ''' create violin plots on an axis ''' dist = max(pos)-min(pos) - w = min(0.15*max(dist, 1.0), 0.5) - for d, p in zip(data, pos): - # calculates the kernel density - k = gaussian_kde(d, bw_method="silverman") - m = k.dataset.min() # lower bound of violin - M = k.dataset.max() # upper bound of violin - x = arange(m, M, (M-m)/100.) # support for violin - v = k.evaluate(x) # violin profile (density curve) - v = v/v.max()*w # scaling the violin to the available space - ax.fill_betweenx(x, p, v+p, facecolor='y', alpha=0.3) - ax.fill_betweenx(x, p, -v+p, facecolor='y', alpha=0.3) + w = min(0.15*max(dist,1.0),0.5) + for d,p in zip(data,pos): + k = gaussian_kde(d, bw_method="silverman") #calculates the kernel density + m = k.dataset.min() #lower bound of violin + M = k.dataset.max() #upper bound of violin + x = arange(m,M,(M-m)/100.) # support for violin + v = k.evaluate(x) #violin profile (density curve) + v = v/v.max()*w #scaling the violin to the available space + ax.fill_betweenx(x,p,v+p,facecolor='y',alpha=0.3) + ax.fill_betweenx(x,p,-v+p,facecolor='y',alpha=0.3) if bp: - ax.boxplot(data, notch=1, positions=pos, vert=1) + ax.boxplot(data,notch=1,positions=pos,vert=1) def calculate_shannon_bits(hist): @@ -437,7 +535,7 @@ def calculate_shannon_bits(hist): # a numpy array of probabilities, calculate the shannon information # gain over a uniform distribution - if sum(hist) == 0: + if sum(hist)==0: return 0 nbins = len(hist) @@ -450,7 +548,6 @@ def calculate_shannon_bits(hist): return base_info - hist_info - def import_output_file(model_file): ''' Import an output file @@ -462,16 +559,15 @@ def import_output_file(model_file): ''' return 0 - -def plot_residue_rate_distributions( - model_files, rate_bins=None, resrange=None, plot_prior=True): +def plot_residue_rate_distributions(model_files, rate_bins = None, resrange=None, plot_prior=True): # Input is a standard output model file and the rate bins # Should note the sectors as well, along with overlap. # Outputs a plot with a sing + import csv colors = ["red", "blue", "yellow", "green"] - resnum_label_skip = 10 + resnum_label_skip=10 # Get data and place into list of lists. @@ -486,8 +582,7 @@ def plot_residue_rate_distributions( nmod = len(d_list[0]) maxbin = int(numpy.max(d_list[0])) - # How to calculate xlim? Keep the uniform prior the same proportion of - # the window + # How to calculate xlim? Keep the uniform prior the same proportion of the window # So proportional to 1/bins. Say 3 or 4 times this value? xlim = 1.0/maxbin * 3 @@ -498,38 +593,35 @@ def plot_residue_rate_distributions( else: resrange = range(resrange[0], resrange[1]+1) - bins = range(1, maxbin+1) + bins = range(1,maxbin+1) # What is the optimal figsize? - fig, ax = plt.subplots(1, len(resrange), sharey='row', figsize=(20, 2)) + fig, ax = plt.subplots(1, len(resrange), sharey='row', figsize=(20,2)) - # print(bins, maxbin) + #print(bins, maxbin) data = [] if rate_bins is not None: if len(rate_bins) < maxbin: - raise Exception( - "Number of inputted rate_bins is less than the maximum bin " - "in output file.") + raise Exception("Number of inputted rate_bins is less than the maximum bin in output file.") x = bins else: + x = bins # Calculate the histograms for n in resrange: - d_hists = [] + d_hists=[] for d in d_list: - nums = d[:, n] - h = calculate_histogram( - list(nums), numpy.array(range(maxbin+1))+0.5) + nums = d[:,n] + h = calculate_histogram(list(nums), numpy.array(range(maxbin+1))+0.5) hist = 1.0 * h[0] / nmod d_hists.append(hist) data.append(d_hists) - # data is a list of lists. Outer index is resnum, inner index is - # the dataset. + # data is a list of lists. Outer index is resnum, inner index is the dataset. # Figure out some way to determine the best ytick rate - ytick_rate = 2 + ytick_rate=2 # Outer loop over all residues for nd in resrange: @@ -543,57 +635,50 @@ def plot_residue_rate_distributions( ax[n].set_xticks([]) - if nd % resnum_label_skip == 0: + if nd%resnum_label_skip == 0: ax[n].set_title(str(nd), fontsize=12) ax[n].set_xticks([0]) ax[n].xaxis.set_ticks_position("top") - ax[n].plot([0, 0], [x[0]-0.2, x[-1]+0.2], - color="black", lw=0.5) + ax[n].plot([0,0],[x[0]-0.2,x[-1]+0.2], color="black", lw=0.5) # Calculate bits of information bits = calculate_shannon_bits(arr) - # print(n, i, bits) + #print(n, i, bits) - ax[n].set_xlim((-xlim, xlim)) - ax[n].set_ylim((-numpy.log(len(x))+1, x[-1]+0.2)) + ax[n].set_xlim((-xlim,xlim)) + ax[n].set_ylim((-numpy.log(len(x))+1,x[-1]+0.2)) if plot_prior: - ax[n].fill_betweenx( - x, 0, 1.0*numpy.ones(len(x))/len(x), facecolor='grey', - alpha=0.5, lw=0) - ax[n].fill_betweenx( - x, 0, -1.0*numpy.ones(len(x))/len(x), facecolor='grey', - alpha=0.5, lw=0) + ax[n].fill_betweenx(x,0,1.0*numpy.ones(len(x))/len(x),facecolor='grey',alpha=0.5, lw=0) + ax[n].fill_betweenx(x,0,-1.0*numpy.ones(len(x))/len(x),facecolor='grey',alpha=0.5, lw=0) if sum(arr) != 0: # Fill in the prior probability (uniform for now) - ax[n].fill_betweenx( - x, 0, arr, facecolor=colors[i], alpha=0.5, lw=0) - ax[n].fill_betweenx( - x, 0, -arr, facecolor=colors[i], alpha=0.5, lw=0) + ax[n].fill_betweenx(x,0,arr,facecolor=colors[i],alpha=0.5, lw=0) + ax[n].fill_betweenx(x,0,-arr,facecolor=colors[i],alpha=0.5, lw=0) # Add in lower bar for information content - ax[n].barh(bottom=-1*bits+1, - width=2*xlim/len(x_lists), height=bits, - left=-1*xlim+i*2*xlim/len(x_lists), - color=colors[i], alpha=0.7, lw=0) - for yval in range(1, int(numpy.max(d))+1, ytick_rate): - ax[n].axhline(y=yval, ls='-', lw=0.5) + ax[n].barh(bottom=-1*bits+1,width=2*xlim/len(x_lists), height=bits,left=-1*xlim+i*2*xlim/len(x_lists),color=colors[i], alpha=0.7, lw=0) + for yval in range(1, int(numpy.max(d))+1,ytick_rate): + ax[n].axhline(y=yval,ls='-', lw=0.5) + + #ax[n].set_xticks([str(nd)]) + ax[n].tick_params(axis='x', which='major', labelsize=0, color="grey") - # ax[n].set_xticks([str(nd)]) - ax[n].tick_params(axis='x', which='major', labelsize=0, - color="grey") ax[n].set_frame_on(False) + #ax[n].spines['top'].set_visible(False) + #ax[n].spines['right'].set_visible(False) + #ax[n].spines['bottom'].set_visible(False) + #ax[n].spines['left'].set_visible(False) ax[0].set_ylabel("HX Rate Bin") ax[0].set_yticks(x[::ytick_rate]) - # ax[0].tick_params(axis='y', which='major', labelsize=8) + #ax[0].tick_params(axis='y', which='major', labelsize=8) plt.savefig("test_violins.png", dpi=300, format="png") plt.show() - def plot_residue_protection_factors(parse_output, rate_bins=None, resrange=None, plot_prior=True, resnum_skip=10, num_best_models=100, @@ -603,10 +688,11 @@ def plot_residue_protection_factors(parse_output, rate_bins=None, # true_vals is a list of residue numbers and protection factors. # These will be sorted and plotted as a red line. + import csv colors = ["red", "blue", "yellow", "green"] - resnum_label_skip = resnum_skip + resnum_label_skip=resnum_skip # Get data and place into list of lists. @@ -615,10 +701,10 @@ def plot_residue_protection_factors(parse_output, rate_bins=None, data_list = [] for po in parse_output: - data_list.append( - po.get_best_scoring_models(num_best_models, return_pf=True)) + #print(po, len(po.get_best_scoring_models(num_best_models, return_pf=True))) + data_list.append(po.get_best_scoring_models(num_best_models, return_pf=True)) - # print(len(data_list), len(data_list[0])) + #print(len(data_list), len(data_list[0])) nres = len(data_list[0][1][1][0]) nmod = len(data_list[0]) @@ -634,11 +720,11 @@ def plot_residue_protection_factors(parse_output, rate_bins=None, pfs = [] for i in d: pfs.append(numpy.array(i[1][0])) - # print(i[1][0]) + #print(i[1][0]) pf_list.append(numpy.array(pfs)) + #print(len(data_list),len(pf_list), len(pf_list[0]), len(pf_list[0][0])) - # How to calculate xlim? Keep the uniform prior the same proportion - # of the window + # How to calculate xlim? Keep the uniform prior the same proportion of the window # So proportional to 1/bins. Say 3 or 4 times this value? xlim = 1.0/maxbin * 3 @@ -649,45 +735,40 @@ def plot_residue_protection_factors(parse_output, rate_bins=None, else: resrange = range(resrange[0], resrange[1]+1) - bins = range(1, maxbin+1) + bins = range(1,maxbin+1) # What is the optimal figsize? - fig, ax = plt.subplots(1, len(resrange), sharey='row', figsize=(20, 2)) + fig, ax = plt.subplots(1, len(resrange), sharey='row', figsize=(20,2)) - # print(bins, maxbin) + #print(bins, maxbin) data = [] if rate_bins is not None: if len(rate_bins) < maxbin: - raise Exception( - "Number of inputted rate_bins is less than the maximum bin " - "in output file.") + raise Exception("Number of inputted rate_bins is less than the maximum bin in output file.") x = bins else: - x = numpy.linspace(minbin, maxbin, parse_output[0].grid_size-1) + x = numpy.linspace(minbin,maxbin,parse_output[0].grid_size-1) # Calculate the histograms for n in resrange: - d_hists = [] + d_hists=[] for d in pf_list: - # print(type(d), d) - nums = d[:, n-1] - # print(n, nums, nums[0]) + #print(type(d), d) + nums = d[:,n-1] + #print(n, nums, nums[0]) if math.isnan(nums[0]): d_hists.append(numpy.zeros(parse_output[0].grid_size-1)) else: - h = calculate_histogram( - list(nums), - numpy.linspace(minbin, maxbin, parse_output[0].grid_size)) - # print(n-1, nums, h) + h = calculate_histogram(list(nums), numpy.linspace(minbin,maxbin,parse_output[0].grid_size)) + #print(n-1, nums, h) hist = 1.0 * h[0] / nmod d_hists.append(hist) data.append(d_hists) - # data is a list of lists. Outer index is resnum, inner index is - # the dataset. + # data is a list of lists. Outer index is resnum, inner index is the dataset. # Figure out some way to determine the best ytick rate - ytick_rate = 10 + ytick_rate=10 # Outer loop over all residues for nd in resrange: @@ -701,60 +782,58 @@ def plot_residue_protection_factors(parse_output, rate_bins=None, ax[n].set_xticks([]) - if nd % resnum_label_skip == 0: + if nd%resnum_label_skip == 0: ax[n].set_title(str(nd), fontsize=12) ax[n].set_xticks([0]) ax[n].xaxis.set_ticks_position("top") - ax[n].plot([0, 0], [x[0]-0.2, x[-1]+0.2], - color="black", lw=0.5) + ax[n].plot([0,0],[x[0]-0.2,x[-1]+0.2], color="black", lw=0.5) # Calculate bits of information bits = calculate_shannon_bits(arr) - # print(n, i, bits) + #print(n, i, bits) - ax[n].set_xlim((-xlim, xlim)) - ax[n].set_ylim((-numpy.log(len(x))+1, x[-1]+0.2)) + ax[n].set_xlim((-xlim,xlim)) + ax[n].set_ylim((-numpy.log(len(x))+1,x[-1]+0.2)) if plot_prior: # Fill in the prior probability (uniform for now) - ax[n].fill_betweenx( - x, 0, 1.0*numpy.ones(len(x))/len(x), facecolor='grey', - alpha=0.5, lw=0) - ax[n].fill_betweenx( - x, 0, -1.0*numpy.ones(len(x))/len(x), facecolor='grey', - alpha=0.5, lw=0) + ax[n].fill_betweenx(x,0,1.0*numpy.ones(len(x))/len(x),facecolor='grey',alpha=0.5, lw=0) + ax[n].fill_betweenx(x,0,-1.0*numpy.ones(len(x))/len(x),facecolor='grey',alpha=0.5, lw=0) + if not math.isnan(arr[0]): - # print(nd, x, arr, len(arr), len(x), numpy) - ax[n].fill_betweenx( - x, 0, arr, facecolor=colors[i], alpha=0.5, lw=0) - ax[n].fill_betweenx( - x, 0, -arr, facecolor=colors[i], alpha=0.5, lw=0) + #print(nd, x, arr, len(arr), len(x), numpy) + ax[n].fill_betweenx(x,0,arr,facecolor=colors[i],alpha=0.5, lw=0) + ax[n].fill_betweenx(x,0,-arr,facecolor=colors[i],alpha=0.5, lw=0) # Add in lower bar for information content - ax[n].barh(bottom=-1*bits+minbin, - width=2*xlim/len(x_lists), height=bits, - left=-1*xlim+i*2*xlim/len(x_lists), - color=colors[i], alpha=0.7, lw=0) + ax[n].barh(bottom=-1*bits+minbin,width=2*xlim/len(x_lists), height=bits,left=-1*xlim+i*2*xlim/len(x_lists),color=colors[i], alpha=0.7, lw=0) for yval in range(minbin, maxbin, ytick_rate): - ax[n].axhline(y=yval, ls='-', lw=0.5) + ax[n].axhline(y=yval,ls='-', lw=0.5) else: - ax[n].fill_betweenx( - x, 0, numpy.zeros(parse_output[0].grid_size-0), - facecolor=colors[i], alpha=0.5, lw=0) - ax[n].tick_params(axis='x', which='major', labelsize=0, - color="grey") + ax[n].fill_betweenx(x,0,numpy.zeros(parse_output[0].grid_size-0),facecolor=colors[i],alpha=0.5, lw=0) + #ax[n].fill_betweenx(x,0,-arr,facecolor=colors[i],alpha=0.5, lw=0) + + #ax[n].set_xticks([str(nd)]) + ax[n].tick_params(axis='x', which='major', labelsize=0, color="grey") + ax[n].set_frame_on(False) + #ax[n].spines['top'].set_visible(False) + #ax[n].spines['right'].set_visible(False) + #ax[n].spines['bottom'].set_visible(False) + #ax[n].spines['left'].set_visible(False) if true_vals: if nd in true_vals.keys: print("TVal", nd, true_vals[nd]) - ax[n].fill_between(-1, 1, 10**true_vals[nd], facecolor='red', - alpha=0.7, lw=0) + ax[n].fill_between(-1,1,10**true_vals[nd], facecolor='red', alpha=0.7, lw=0) ax[0].set_ylabel("Log(Protection Factor)") ax[0].set_yticks(x[::ytick_rate]) - # ax[0].tick_params(axis='y', which='major', labelsize=8) + #ax[0].tick_params(axis='y', which='major', labelsize=8) plt.savefig("test_violins_pf.png", dpi=300, format="png") plt.show() + + + #fig, ax = plt.subplots(nrows=1, ncols=nres, figsize=(20,5)) diff --git a/v2/bin/hdxworkbench.py b/v2/bin/hdxworkbench.py index f15182d..93ed99d 100755 --- a/v2/bin/hdxworkbench.py +++ b/v2/bin/hdxworkbench.py @@ -8,6 +8,31 @@ import model import hxio import analysis +from operator import itemgetter + + +def get_best_scoring_models(o, minsize=100): + new_pof = analysis.deepcopy(o.pof1) + pof_all = analysis.concatenate_pofs(new_pof, o.pof2) + smt = pof_all.get_models(return_pf=False) + smt.sort(key=itemgetter(0)) + sum2=0.0 + sum1=0.0 + finali=minsize #HB store the final i; default is the minimum number of points to use, obviously! + minstderr2=1.0E+34 #HB 23-May-2018. Can never be negative, so a dumb value to initialize for debugging purposes. + print('Processing {} models'.format(len(smt))) + for i in range(0,len(smt),1): + sum1+=smt[i][0] + sum2+=smt[i][0]*smt[i][0] + avg=sum1/(i+1) + var=(sum2/(i+1)) - (avg * avg) + stderr2=var/(i+1) #square of stderr, monotonic with stderr and faster to calculate + if (i >= minsize) and (stderr2 < minstderr2): + # print('Adjusting stderr: {}'.format(stderr2)) + finali=i+1 #HB this is the i'th entry in smt, which is the current minimum in the stderr + minstderr2=stderr2 #HB 23-May-2018 + return (finali, minstderr2) + if __name__ == '__main__': parser = argparse.ArgumentParser( @@ -17,6 +42,8 @@ parser.add_argument('--mol_name', help='Molecule name', required=True) parser.add_argument('-o','--outputdir', help='Output directory.', required=True) + parser.add_argument('--control', help='Control sample name', required=True) + parser.add_argument('--ligand', help='Ligand sample name', required=True) parser.add_argument('--init', help='How to initialize - either "random" or "enumerate". ' 'Enumerate is slower but sampling will converge faster. ' 'Default: enumerate', @@ -112,7 +139,7 @@ sampler = sampling.MCSampler(sys, sigma_sample_level="timepoint") # First, run a short minimization step - sampler.run(50, 0.0001, write=True) + sampler.run(100, 0.0001, write=True) # Slowly cool system sampler.run(args.annealing_steps, 3) @@ -128,31 +155,42 @@ sampler.run(args.nsteps, 1, write=True) files = [f for dr, ds, files in os.walk(args.outputdir) for f in files if f.endswith('.dat')] - if len(files) >= 2: - os.chdir(args.outputdir) - oa = analysis.OutputAnalysis([files[0]]) - oa1 = analysis.OutputAnalysis([files[1]]) - - conv = oa.get_convergence(20) - conv1 = oa1.get_convergence(20) - - distmat = conv.get_distance_matrix(num_models=20) - distmat1 = conv1.get_distance_matrix(num_models=20) - - cutoff_list = conv.get_cutoffs_list(1.0) - cutoff_list1 = conv1.get_cutoffs_list(1.0) - - pvals, cvs, percents = conv.get_clusters(cutoff_list) - pvals1, cvs1, percents1 = conv1.get_clusters(cutoff_list1) - - sampling_precision,pval_converged,cramersv_converged,percent_converged = conv.get_sampling_precision(cutoff_list, pvals, cvs, percents) - sampling_precision1,pval_converged1,cramersv_converged1,percent_converged1 = conv1.get_sampling_precision(cutoff_list1, pvals1, cvs1, percents1) - - pofs = conv.cluster_at_threshold_and_return_pofs(sampling_precision) - pofs1 = conv1.cluster_at_threshold_and_return_pofs(sampling_precision1) - - dhdx = analysis.DeltaHDX(pofs[0], pofs1[0]) - diff, Z, mean1, mean2, sd1, sd2 = dhdx.calculate_dhdx() - dhdx.write_dhdx_file() - - + control_files = [] + ligand_files = [] + for f in files: + if args.control in f: + control_files.append(os.path.join(args.outputdir, f)) + elif args.ligand in f: + ligand_files.append(os.path.join(args.outputdir, f)) + print('Using control files: {}'.format(control_files)) + print('Using ligand files: {}'.format(ligand_files)) + oa = analysis.OutputAnalysis(control_files) + oa1 = analysis.OutputAnalysis(ligand_files) + + num_models, minstderr = get_best_scoring_models(oa, 100) + num_models1, minstderr1 = get_best_scoring_models(oa1, 100) + print('Control best num of models {} with stderr: {}'.format(num_models, minstderr)) + print('Ligand best num of models {} with stderr: {}'.format(num_models1, minstderr1)) + + conv = oa.get_convergence(num_models) + conv1 = oa1.get_convergence(num_models1) + + distmat = conv.get_distance_matrix(num_models=num_models) + distmat1 = conv1.get_distance_matrix(num_models=num_models1) + + cutoff_list = conv.get_cutoffs_list(1.0) + cutoff_list1 = conv1.get_cutoffs_list(1.0) + + pvals, cvs, percents = conv.get_clusters(cutoff_list) + pvals1, cvs1, percents1 = conv1.get_clusters(cutoff_list1) + + sampling_precision,pval_converged,cramersv_converged,percent_converged = conv.get_sampling_precision(cutoff_list, pvals, cvs, percents) + sampling_precision1,pval_converged1,cramersv_converged1,percent_converged1 = conv1.get_sampling_precision(cutoff_list1, pvals1, cvs1, percents1) + + pofs = conv.cluster_at_threshold_and_return_pofs(sampling_precision) + pofs1 = conv1.cluster_at_threshold_and_return_pofs(sampling_precision1) + + dhdx = analysis.DeltaHDX(pofs[0], pofs1[0]) + diff, Z, mean1, mean2, sd1, sd2 = dhdx.calculate_dhdx() + dhdx.write_dhdx_file(prefix='{}/'.format(args.outputdir)) + \ No newline at end of file diff --git a/v2/examples/HDXWorkbench_example/modeling_new.py b/v2/examples/HDXWorkbench_example/modeling_new.py index 9cae753..6978a29 100644 --- a/v2/examples/HDXWorkbench_example/modeling_new.py +++ b/v2/examples/HDXWorkbench_example/modeling_new.py @@ -24,7 +24,7 @@ ### output directory for this simulation. -outputdir = "./output_v2/" +outputdir = "./output_v2_test/" ### ########################################## ### Experimental Input Parameters diff --git a/v2/pyext/src/analysis.py b/v2/pyext/src/analysis.py index 3b1dcc2..5908889 100644 --- a/v2/pyext/src/analysis.py +++ b/v2/pyext/src/analysis.py @@ -2,6 +2,7 @@ Analysis functions for HDX simulations """ from __future__ import print_function +from __future__ import division import hxio #from scipy.stats import cumfreq #from scipy.stats import chi2_contingency @@ -59,8 +60,8 @@ def get_models(self, num_gsm="all"): mod2 = [g[1] for g in m2] else: mod1 = [g[1] for g in m1[0:num_gsm]] - mod2 = [g[1] for g in m2[0:num_gsm]] - return mod1, mod2 + mod2 = [g[1] for g in m2[0:num_gsm]] + return mod1, mod2 def total_score_pvalue_and_cohensd(self, num_gsm="all"): @@ -79,7 +80,7 @@ def total_score_pvalue_and_cohensd(self, num_gsm="all"): return pvalue, cohens_d def residue_pvalue_and_cohensd(self, num_gsm="all"): - + bsm1, bsm2 = self.get_models(num_gsm) output=[] @@ -132,12 +133,12 @@ def precision_cluster(self, threshold): neighbors.append([count]) # model is a neighbor of itself for i in range(num_models-1): - for j in range(i+1,num_models): + for j in range(i+1,num_models): if distmat[i][j]<=threshold: neighbors[i].append(j) neighbors[j].append(i) - #print(i,j,distmat[i][j],len(neighbors[i]),len(neighbors[j])) + #print(i,j,distmat[i][j],len(neighbors[i]),len(neighbors[j])) # 2). Get the weightiest cluster, and iterate @@ -159,16 +160,16 @@ def precision_cluster(self, threshold): for eachu in unclustered: # if multiple clusters have same maxweight this tie is broken arbitrarily! if len(neighbors[eachu])>max_neighbors: max_neighbors=len(neighbors[eachu]) - currcenter=eachu - + currcenter=eachu + #form a new cluster with u and its neighbors cluster_centers.append(currcenter) - cluster_members.append([n for n in neighbors[currcenter]]) + cluster_members.append([n for n in neighbors[currcenter]]) - #update neighbors + #update neighbors for n in neighbors[currcenter]: #removes the neighbor from the pool - unclustered.remove(n) #first occurence of n is removed. + unclustered.remove(n) #first occurence of n is removed. boolUnclustered[n]=False # clustered for n in neighbors[currcenter]: @@ -176,7 +177,7 @@ def precision_cluster(self, threshold): if not boolUnclustered[unn]: continue neighbors[unn].remove(n) - + return cluster_centers, cluster_members @@ -202,7 +203,7 @@ def get_clusters(self, cutoffs_list): pvals.append(pval) cvs.append(cramersv) percents.append(percent_explained) - + f1.write(str(c)+", "+str(pval)+", "+str(cramersv)+", "+str(percent_explained)+"\n") return pvals, cvs, percents @@ -238,7 +239,7 @@ def get_contingency_table(self, num_clusters,cluster_members,all_models,run1_mod reduced_ctable=[] retained_clusters=[] - + for i in range(num_clusters): if full_ctable[i][0]<=10.0 or full_ctable[i][1]<=10.0: #if full_ctable[i][0]<=0.10*numModelsRun1 and full_ctable[i][1] <= 0.10*numModelsRun2: @@ -252,14 +253,14 @@ def test_sampling_convergence(self, contingency_table, total_num_models): if len(contingency_table)==0: return 0.0,1.0 - + ct = numpy.transpose(contingency_table) [chisquare,pvalue,dof,expected]=scipy.stats.chi2_contingency(ct) if dof==0.0: cramersv=0.0 else: cramersv=math.sqrt(chisquare/float(total_num_models)) - + return(pvalue,cramersv) def percent_ensemble_explained(self, ctable,total_num_models): @@ -419,15 +420,15 @@ def parse_header(self): ''' f = open(self.output_file, "r") for line in f.readlines(): - + # > means model data (so header is over.) if line[0]==">": break - + # #-symbol means datasets elif line[0:2]=="# ": self.datafiles.append( (line[2:].split("|")[0].strip(), float(line[2:].split("|")[2].strip())) ) - + # @-symbol means sectors elif line[0:2]=="@ ": for s_string in line[2:].strip().split("|"): @@ -498,7 +499,7 @@ def get_all_models(self, return_pf=False): f = open(self.output_file, "r") models = [] # Cycle over all lines - for line in f.readlines(): + for line in f.readlines(): if line[0]==">": score = float(line.split("|")[1].strip()) @@ -509,7 +510,7 @@ def get_all_models(self, return_pf=False): if return_pf: ml1 = model_list model_list = self.models_to_protection_factors(model_list) - models.append((score, model_list)) + models.append((score, model_list)) self.models = models @@ -704,12 +705,12 @@ def split_into_two_POFs(self): random.shuffle(n_output_files) pof1 = ParseOutputFile(self.output_files[n_output_files[0]]) - for i in n_output_files[1:int(len(n_output_files)/2)]: + for i in n_output_files[1:len(n_output_files)//2]: new_pof = ParseOutputFile(self.output_files[n_output_files[i]]) concatenate_pofs(pof1, new_pof) - pof2 = ParseOutputFile(self.output_files[n_output_files[int(len(n_output_files)/2)]]) - for i in n_output_files[int(len(n_output_files)/2)+1:]: + pof2 = ParseOutputFile(self.output_files[n_output_files[len(n_output_files)//2]]) + for i in n_output_files[len(n_output_files)//2+1:]: new_pof = ParseOutputFile(self.output_files[n_output_files[i]]) concatenate_pofs(pof2, new_pof) @@ -773,19 +774,20 @@ def get_best_scoring_models(modelfile, scorefile, num_best_models=100, prefix=No top_models=[] top_score_indices=[] top_scores=[] - infile=open(scorefile, "r") - for line in infile: - scores.append(float(line.split()[0].strip())) - infile=open(modelfile, "r") - - for line in infile: - models.append(line) + with open(scorefile, "r") as infile: + for line in infile: + scores.append(float(line.split()[0].strip())) + with open(modelfile, "r") as infile: + for line in infile: + models.append(line) for i in range(num_best_models): top_score_indices.append(scores.index(min(scores))) top_scores.append(min(scores)) #print(scores, min(scores), scores.index(min(scores)), top_score_indices) scores[scores.index(min(scores))]=max(scores)+1 + if i > len(scores): # To get as many models as in the scorefile + break if write_file: output_file=open(outfile, "w") return top_models, top_scores @@ -830,7 +832,7 @@ def get_residue_rate_probabilities(modelfile, scorefile, sectors, seq, grid, num if hasattr(grid, '__iter__'): grid=len(grid) - + best_models, best_scores=get_best_scoring_models(modelfile, scorefile, num_models, write_file=False) best_models=numpy.array(best_models) @@ -882,9 +884,9 @@ def get_cdf(exp_models): A=numpy.array(exp_models) y=numpy.linspace(1./len(exp_models),1,len(exp_models)) #print(len(exp_models[0])) - for i in range(len(exp_models[0])): - counts, edges = numpy.histogram(A[:,i], len(A), range=(-6,0), density=False) - #print i,A[:,i],counts,numpy.cumsum(counts*1.0/len(A)) + for i in range(len(exp_models[0])): + counts, edges = numpy.histogram(A[:,i], len(A), range=(-6,0), density=False) + #print i,A[:,i],counts,numpy.cumsum(counts*1.0/len(A)) exp_model_edf[:,i]=numpy.cumsum(counts*1.0/len(A)) return exp_model_edf @@ -896,7 +898,7 @@ def get_chisq(exp_models1, exp_models2, nbins): #print(len(exp_models1[0])) for i in range(269,len(exp_models1[0])): meanA = numpy.mean(A[:,i]) - ssd = numpy.std(A[:,i])**2 + numpy.std(B[:,i])**2 + ssd = numpy.std(A[:,i])**2 + numpy.std(B[:,i])**2 sstdev = numpy.sqrt( ssd / 5000 ) meanB = numpy.mean(B[:,i]) t = 1.96 diff --git a/v2/pyext/src/plots.py b/v2/pyext/src/plots.py index e41bf6b..3e76322 100644 --- a/v2/pyext/src/plots.py +++ b/v2/pyext/src/plots.py @@ -2,19 +2,20 @@ """ from __future__ import print_function + +import math import os -import system + +# from pylab import * +import matplotlib import numpy -import math + import analysis -#from pylab import * -import matplotlib + matplotlib.use('Agg') import matplotlib.pyplot as plt -import scipy #from scipy.stats import gaussian_kde #from numpy.random import normal -from scipy.integrate import simps import tools @@ -462,7 +463,6 @@ def plot_residue_rate_distributions(model_files, rate_bins = None, resrange=None # Input is a standard output model file and the rate bins # Should note the sectors as well, along with overlap. # Outputs a plot with a sing - import csv colors = ["red", "blue", "yellow", "green"] @@ -585,7 +585,6 @@ def plot_residue_protection_factors(parse_output, rate_bins=None, # Input is a standard output model file and the rate bins # Should note the sectors as well, along with overlap. # Outputs a plot with a sing - import csv colors = ["red", "blue", "yellow", "green"]