Namita

fittingParameters.py

@@ -30,7 +30,7 @@ def __init__(self, concentrations=None, f_abs_green_max=None,
         self.frac_bound_lowerbound = 0.999
         self.fdr_cutoff = 0.05
 
-        self.fitParameters = pd.DataFrame(columns=['fmax', 'dG', 'fmin', 'toff'],
+        self.fitParameters = pd.DataFrame(columns=['fmax', 'dG', 'fmin', 'toff', 'ton'],
                                           index=['lowerbound', 'initial', 'upperbound'])
         if concentrations is not None:
             self.concentrations = concentrations
@@ -72,7 +72,23 @@ def __init__(self, concentrations=None, f_abs_green_max=None,
                 self.fitParameters[currParam]['lowerbound']= 0
                 self.fitParameters[currParam]['initial']= np.nan # this will be defined per cluster
                 self.fitParameters[currParam]['upperbound'] = self.find_fmax_lowerbound(f_abs_green_nonbinders, self.fdr_cutoff) # this is futher reducedor relaxed if the first point in the binding curve was fit
+
+            if fittype == 'onrate':
+                currParam = 'ton'
+                self.fitParameters[currParam]['lowerbound'] = 1E-1
+                self.fitParameters[currParam]['initial']= 1E4 
+                self.fitParameters[currParam]['upperbound'] = 1E6
+
+                currParam = 'fmin'  
+                self.fitParameters[currParam]['lowerbound'] = 0
+                self.fitParameters[currParam]['initial'] = np.nan # this will be defined per cluster
+                self.fitParameters[currParam]['upperbound'] = self.find_fmax_lowerbound(f_abs_green_nonbinders, self.fdr_cutoff) # this is futher reducedor relaxed if the first point in the binding curve was fit
 
+                currParam = 'fmax'
+                self.fitParameters[currParam]['lowerbound']= 0
+                self.fitParameters[currParam]['initial']= np.nan # this will be defined per cluster
+                self.fitParameters[currParam]['upperbound'] = self.find_fmax_lowerbound(f_abs_green_nonbinders, self.fdr_cutoff) # this is futher reducedor relaxed if the first point in the binding curve was fit
+
 
 
             # estimate conversion of f_abs_red to f_abs_green

libs/CurveFitFun.m

@@ -28,5 +28,13 @@
          function q = findFDR(score, null_scores)
              q = sum(abs(null_scores) > abs(score))/length(null_scores);
          end
+
+         % function to find on rate curve
+         function fracbound = findOnRate(x, time)
+             fmax = x(1);
+             ton  = x(2);
+             fmin = x(3);
+             fracbound = fmin + fmax*(1 - exp(-time/ton));
+         end
      end
 end

libs/IMlibs.py

@@ -340,18 +340,21 @@ def fitSetKds(fitParametersFilenameParts, bindingSeriesFilenameParts, initialFit
     fitParameters = joinTogetherFitParts(fitParametersFilenameParts)
     return fitParameters
 
-def fitSetKoff(fitParametersFilenameParts, bindingSeriesFilenameParts, initialFitParameters, scale_factor):
+def fitSetKoff(fitParametersFilenameParts, bindingSeriesFilenameParts, initialFitParameters, scale_factor, fittype=None):
+    if fittype is None: fittype = 'offrate'
+    if fittype == 'offrate': parameter = 'toff'
+    if fittype == 'onrate': parameter = 'ton'
     workerPool = multiprocessing.Pool(processes=len(bindingSeriesFilenameParts)) #create a multiprocessing pool that uses at most the specified number of cores
     for i, bindingSeriesFilename in bindingSeriesFilenameParts.items():
         result = workerPool.apply_async(findKoff, args=(bindingSeriesFilename, fitParametersFilenameParts[i],
                                                      initialFitParameters['fmax']['lowerbound'], initialFitParameters['fmax']['upperbound'], initialFitParameters['fmax']['initial'],
-                                                     initialFitParameters['toff']['lowerbound'], initialFitParameters['toff']['upperbound'], initialFitParameters['toff']['initial'],
+                                                     initialFitParameters[parameter]['lowerbound'], initialFitParameters[parameter]['upperbound'], initialFitParameters[parameter]['initial'],
                                                      initialFitParameters['fmin']['lowerbound'], initialFitParameters['fmin']['upperbound'], initialFitParameters['fmin']['initial'],
-                                                     scale_factor,)
+                                                     scale_factor, fittype)
                                )
     workerPool.close()
     workerPool.join()
-    fitParameters = joinTogetherFitParts(fitParametersFilenameParts, parameter='toff')
+    fitParameters = joinTogetherFitParts(fitParametersFilenameParts, parameter=parameter)
     return fitParameters
 
 def findKds(bindingSeriesFilename, outputFilename, fmax_min, fmax_max, fmax_initial, kd_min, kd_max, kd_initial, fmin_min, fmin_max, fmin_initial, scale_factor):
@@ -368,14 +371,14 @@ def findKds(bindingSeriesFilename, outputFilename, fmax_min, fmax_max, fmax_init
     except Exception,e:
         return(matlabFunctionCallString,'Python excpetion generated in findKds: ' + e.message + e.stack)
 
-def findKoff(bindingSeriesFilename, outputFilename, fmax_min, fmax_max, fmax_initial, kd_min, kd_max, kd_initial, fmin_min, fmin_max, fmin_initial, scale_factor):
-    matlabFunctionCallString = "fitOffRateCurve('%s', [%4.2f, %4.2f, %4.2f], [%4.2f, %4.2f, %4.2f], '%s', [%4.2f, %4.2f, %4.2f], %4.2f );"%(bindingSeriesFilename,
+def findKoff(bindingSeriesFilename, outputFilename, fmax_min, fmax_max, fmax_initial, kd_min, kd_max, kd_initial, fmin_min, fmin_max, fmin_initial, scale_factor, fittype):
+    matlabFunctionCallString = "fitOffRateCurve('%s', [%4.2f, %4.2f, %4.2f], [%4.2f, %4.2f, %4.2f], '%s', [%4.2f, %4.2f, %4.2f], %4.2f, '%s' );"%(bindingSeriesFilename,
 
                                                                                                                 fmax_min, kd_min, fmin_min,
                                                                                                                 fmax_max, kd_max, fmin_max,
                                                                                                                 outputFilename,
                                                                                                                 fmax_initial, kd_initial, fmin_initial,
-                                                                                                                scale_factor)
+                                                                                                                scale_factor, fittype)
     try:
         logString = spawnMatlabJob(matlabFunctionCallString)
         return (matlabFunctionCallString, logString)
@@ -563,11 +566,15 @@ def loadOffRatesCurveFromCPsignal(filename, timeStampDict, numCores=None):
         xvalues[tiles==tile] = timeDelta
     return binding_series, np.array(table['all_cluster_signal']), xvalues, tiles  
 
-def loadFittedCPsignal(filename):
+def loadFittedCPsignal(filename, index_by_cluster=None):
+    if index_by_cluster is None: index_by_cluster = False
     f = open(filename); header = np.array(f.readline().split()); f.close()
     index_start = np.where(header=='variant_number')[0][0]
     index_end = len(header)
-    table =  pd.read_table(filename, usecols=tuple(range(index_start,index_end)))
+    if index_by_cluster:
+        table =  pd.read_table(filename, usecols=tuple([0]+range(index_start,index_end)), index_col=0)
+    else:
+        table =  pd.read_table(filename, usecols=tuple(range(index_start,index_end)))
     binding_series, all_cluster_image = loadBindingCurveFromCPsignal(filename)
     for col in range(binding_series.shape[1]):
         table[col] = binding_series[:, col]
@@ -628,7 +635,7 @@ def findVariantTable(table, parameter=None, numCores=None, variants=None):
 def filterFitParameters(sub_table):
     not_nan_binding_points = [0,1] # if first two binding points are NaN, filter
     nan_filter = np.all(np.isfinite(sub_table[not_nan_binding_points]), 1)
-    barcode_filter = np.array(sub_table['fraction_consensus'] > 50)
+    barcode_filter = np.array(sub_table['fraction_consensus'] >= 67)
     fit_filter = np.array(sub_table['rsq'] > 0)
 
     # use barcode_filter, nan_filter to reduce sub_table