Klauda Lab Simulation Wiki

This page is read only. You can view the source, but not change it. Ask your administrator if you think this is wrong.
<code python>
# determines the frequency of peptide 
# uses local_defect output from blob script and truDiff coordinates output 
# Written by Ky Wildermuth, 11/16. Feel free to edit. Feel free to ask questions: kywildermuth@gmail.com


import pandas as pd
import numpy as np
import glob
import re 


#### USER INPUT PARAMETERS #####
# define molecular interaction cutoff in Angstroms (choose large max_cutoff (ex 100) and small min_cutoff (ex -100) if don't want)
max_cutoff = 20
min_cutoff = 0
# define the frame at which the analysis ends (typically before binding is initiated)
prebind_limit = 1800 
# define the frame at which analysis begins (can set to 0 but allows for more control of range)
start_limit = 0
# define the number of dcd 
num_dcd = 1
#define the number of frames per dcd
num_frame = 4000
# define the defect area (d) and movement (m) threshold ranges and intervals. d is measured in Ang^2 and m in Ang. 
d_thresh_arr = range(0,301)
m_thresh_arr = np.arange(0,5,0.1)

######################


m_thresh_length = len(m_thresh_arr)
# compute and define the theoretical total number of data points per coor file (will be used to check if files have same number of data entries)
theo_length = num_dcd*num_frame

# load the data files
dist_path = "binding_distances.dat"
freq_path = "local_defect_top.dat"

# define number of rows, numrows
numrows = prebind_limit - start_limit

# open data files for line check
dist = open(dist_path)
freq = open(freq_path)


files = [freq, dist]

# define a natural sort function type for numerical file combination
def natural_key(string_):
	return [int(s) if s.isdigit() else s for s in re.split(r'(\d+)', string_)]

#check to make sure data files have the same number of lines
def file_len(fname):
	with fname as f:
	    for i, l in enumerate(f):
	        pass
	return i + 1

for fname in files:
	length = file_len(fname)
	print fname,",", length

yes = set(['yes','y', 'ye', ''])
no = set(['no','n'])
print 'Confirm the coordinate files all have a length of', theo_length, '[y/n]'
choice = raw_input().lower()
if choice in yes:
   print "confirmed"
elif choice in no:
    print('Check input files and try again')
    exit(1)

# loop through movement threshold array
m_count = start_limit
for m_thresh in m_thresh_arr:
	m_count += 1
	print('%d of %d' % (m_count, m_thresh_length))
	stat_sum_df_path = "%d_primary_local.dat" % m_count

	# load the coordinate data files to be used (phosphate and protein coordinates)  - only load up to defined prebind_limit
	dist_dat = pd.read_csv(dist_path, delim_whitespace = True, skiprows=start_limit, nrows=numrows, header = True, names = ["com", "custom", "min"])


	# extract the coordinates from each file to be used (z coord)
	c1 = dist_dat['min']

	# concatenate the coordinates to be used to one data frame and label the columns
	defect_cond_dat = list(c1)
	df1 = pd.DataFrame(data = defect_cond_dat, columns = ['dist'])

	# define a new column, dist. Value is +1 if distance between peptide and phosphate layer is less than or equal to cutoff distance, otherwise dist is 0. 
	df1['sample_space'] = 0
	df1.loc[(df1['dist'] <= max_cutoff) & (df1['dist'] >= min_cutoff), 'sample_space'] = 1

	# define a new column, dir. value is +1 if peptide moves towards primary leaflet, -1 if moves away. 
	# maybe want to make it 0 if dist is less than 0 (because then in the membrane, and net movement isn't very interesting)
	df1['dir'] = -1
	df1.loc[df1['dist'].shift(-1) < df1['dist'], 'dir'] = 1

	# define a new column, dDiff. Value is magnitude of net z-motion from current frame n, to next frame n+1, 0 if cross PBC
	# note that a negative dDiff value means that the peptide moves toward the membrane
	df1['dDiff'] = df1['dist'].shift(-1) - df1['dist'] 

	# define new column, M. value is +1 if and only if peptide moves toward primary leaflet by an interval more than m_thresh AND starts within cutoff range. 
	# Value is -1 if doesnt meet threshold but still within cutoff range
	# move towards membrane if dDiff is negative, thus if dDiff*(-1) > threshold, this means the peptide move towards membrane by at least the threshold amount
	df1['M'] = 0
	df1.loc[(df1['dDiff']*(-1) > m_thresh), 'M'] = 1
	df1.loc[(df1['dDiff']*(-1) < m_thresh), 'M'] = -1
	df1.loc[df1['sample_space'] != 1, 'M'] = 0

	# if using frequency data, load the freq data files to be used
	#defect = pd.read_csv(freq_path, delim_whitespace = True, header = True, names = ["dyn_num", "defect_area", "comp_val"])

	# if using local defect data, load the data file to be used
	defect = pd.read_csv(freq_path, delim_whitespace = True, header = False, names = ['dyn_num', 'defect_area', 'protein_area', 'comp_val'])
	#print(defect[:10])
	df1['max_defect'] = 0
	
	# define var frame_num, starting at start_limit frame
	frame_num = start_limit + 1
	for index, row in df1.iterrows():
		#print(frame_num)
		# make a set of all of the defects for the given frame number
		frame_defects = defect[defect['dyn_num'] == frame_num]
		# define the max defect area from above set
		# print a list of max defects and compare timeseries to local defect
		max_defect_area = frame_defects['defect_area'].max()
		#print(max_defect_area)
		# insert max defect area of binding surface in max_defect column
		df1.set_value(index, 'max_defect', max_defect_area)
		frame_num = frame_num + 1


	# build columns of new dataframe to write statistical results to
	# can clean this up with a loop
	thresh_d_s = pd.Series(666, index=d_thresh_arr)
	thresh_m_s = pd.Series(m_thresh, index=d_thresh_arr)
	M_s = pd.Series(666, index=d_thresh_arr)
	D_s = pd.Series(666, index=d_thresh_arr)
	MD_s = pd.Series(666, index=d_thresh_arr)
	sample_space_s = pd.Series(666, index=d_thresh_arr)


	# loop over threshold values and conduct D stats
	for d_thresh in d_thresh_arr:
		D = 'D_%s' % d_thresh
		MD = 'MD_m%s_d%s' % (m_thresh, d_thresh)

		# define new variable D. Value is 1 if max defect is larger than current defect threshold and event is in sample space
		df1[D] = 0
		df1.loc[(df1['max_defect'] >= d_thresh) & (df1['sample_space'] == 1), D] = 1

		# define new variable MD. Value is 1 if and only if M is true and D is true
		df1[MD] = 0
		df1.loc[(df1['M'] == 1) & (df1[D] == 1) , MD] = 1		

		#df1.fillna(0, inplace=True)

		# M_one_count is number of events that net movement of peptide was towards membrane AND meets distance threshold
		try:
			M_count = len(df1.groupby(['M']).groups[1])
		except KeyError:
			M_count = 0			
		# D_count is number of events that max defect area is larger than or equal to defect area threshold, restricting to those events in sample space
		try:	
			D_count = len(df1.groupby([D]).groups[1])
		except KeyError:
			D_count = 0	
		# MD_count is number of events that both M and D_ss are 1
		try:	
			MD_count = len(df1.groupby([MD]).groups[1])
		except KeyError:
			MD_count = 0				
		#sample_space is number of events that do not cross PBC boundary and/or begin within cutoff distance
		try:
			sample_space_count = len(df1.groupby(['sample_space']).groups[1])
		except KeyError:
			sample_space_count = 0	
					
		thresh_d_s.set_value(d_thresh, d_thresh)
		M_s.set_value(d_thresh, M_count)
		D_s.set_value(d_thresh, D_count)
		MD_s.set_value(d_thresh, MD_count)
		sample_space_s.set_value(d_thresh, sample_space_count)
	

	# concatenate the coordinates to be used to one data frame and label the columns
	stat_sum = list(zip(thresh_d_s, thresh_m_s, sample_space_s, M_s, D_s, MD_s))
	df2 = pd.DataFrame(data = stat_sum, columns = ['d_thresh','m_thresh','sample_space', 'M', 'D','MD'])
	df2.to_csv(stat_sum_df_path, sep = ' ')

	# remove header if it is not the first file
	if m_count != 1:
		with open(stat_sum_df_path, 'r') as fin:
			data = fin.read().splitlines(True)
		with open(stat_sum_df_path, 'w') as fout:
			fout.writelines(data[1:])	

# combine all primary-surface-only files into one 
read_files_1 = glob.glob("*primary_local.dat")
read_files_1
read_files = sorted(read_files_1, key=natural_key)
read_files
with open("freq_stat_total.dat", "wb") as outfile:
	for f in read_files:
		with open(f, "rb") as infile:
			outfile.write(infile.read())


</code>