-
Notifications
You must be signed in to change notification settings - Fork 1
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
added visualize_MCchain() function to visualization module
changed init.ipy to import visualization module added script to process output of MC simulations scripts for senior thesis / presentation figures script for group meeting plots
- Loading branch information
Deepti Kannan
committed
Jun 3, 2019
1 parent
4fb8776
commit d2a4dd3
Showing
7 changed files
with
1,931 additions
and
30 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,115 @@ | ||
| """Script to analyze output of MC simulations of nucleosome chains.""" | ||
| import re | ||
| from pathlib import Path | ||
| import pickle | ||
|
|
||
| import matplotlib.cm as cm | ||
| import numpy as np | ||
| import seaborn as sns | ||
| import pandas as pd | ||
| from matplotlib import pyplot as plt | ||
| import scipy | ||
| from scipy import stats | ||
| from scipy.optimize import curve_fit | ||
|
|
||
| from nuc_chain import geometry as ncg | ||
| from nuc_chain import linkers as ncl | ||
| from MultiPoint import propagator | ||
| from nuc_chain import fluctuations as wlc | ||
| from nuc_chain.linkers import convert | ||
| from nuc_chain import visualization as vis | ||
| from mpl_toolkits.axes_grid1 import make_axes_locatable | ||
|
|
||
| #These follow Andy's plotting preferences for | ||
| params = {'axes.edgecolor': 'black', 'axes.grid': True, 'axes.titlesize': 20.0, | ||
| 'axes.linewidth': 0.75, 'backend': 'pdf','axes.labelsize': 18,'legend.fontsize': 18, | ||
| 'xtick.labelsize': 18,'ytick.labelsize': 18,'text.usetex': False,'figure.figsize': [7, 5], | ||
| 'mathtext.fontset': 'stixsans', 'savefig.format': 'pdf', 'xtick.bottom':True, 'xtick.major.pad': 5, 'xtick.major.size': 5, 'xtick.major.width': 0.5, | ||
| 'ytick.left':True, 'ytick.right':False, 'ytick.major.pad': 5, 'ytick.major.size': 5, 'ytick.major.width': 0.5, 'ytick.minor.right':False, 'lines.linewidth':2} | ||
|
|
||
| plt.rcParams.update(params) | ||
|
|
||
| teal_flucts = '#387780' | ||
| red_geom = '#E83151' | ||
| dull_purple = '#755F80' | ||
| rich_purple = '#e830e8' | ||
|
|
||
| dir_default='csvs/MC/41bp_notrans' | ||
|
|
||
| def extract_MC(dir=dir_default, timestep=55, v=1): | ||
| """Extract the positions and orientation traids of each nucleosome in a MC simulation.""" | ||
| dir = Path(dir) | ||
| if dir.is_dir() is False: | ||
| raise ValueError(f'The director {dir} does not exist') | ||
| entry_pos = np.loadtxt(dir/f'r{timestep}v{v}') | ||
| entry_us = np.loadtxt(dir/f'u{timestep}v{v}') | ||
| entry_t3 = entry_us[:, 0:3] #U | ||
| entry_t1 = entry_us[:, 3:] #V | ||
| entry_t2 = np.cross(entry_t3, entry_t1, axis=1) #UxV | ||
| num_nucs = entry_pos.shape[0] | ||
| entry_rots = [] | ||
| for i in range(num_nucs): | ||
| #t1, t2, t3 as columns | ||
| rot = np.eye(3) | ||
| rot[:, 0] = entry_t1[i, :] #V | ||
| rot[:, 1] = entry_t2[i, :] #UxV | ||
| rot[:, 2] = entry_t3[i, :] #U | ||
| entry_rots.append(rot) | ||
| entry_rots = np.array(entry_rots) | ||
| return entry_rots, entry_pos | ||
| #^above files saved in csvs/MLC in npy format | ||
|
|
||
| def check_equilibration(maxtimestep=55): | ||
| """Plot the end-to-end distance of the simulated structure vs. time to see if it's equilibrated""" | ||
| MC_sqrtR2 = np.zeros((maxtimestep+1,)) | ||
|
|
||
| for i in range(maxtimestep): | ||
| rots, pos = extract_MC(timestep=i, v=1) | ||
| pos = pos - pos[0, :] | ||
| r2 = np.array([pos[i, :].T@pos[i, :] for i in range(pos.shape[0])]) | ||
| MC_sqrtR2[i] = np.sqrt(r2[-1]) | ||
|
|
||
| fig, ax = plt.subplots() | ||
| ax.plot(np.arange(0, maxtimestep+1), MC_sqrtR2) | ||
| plt.xlabel('Time step') | ||
| plt.ylabel('End-to-end distance (nm)') | ||
|
|
||
|
|
||
| def plot_R2_MC_vs_theory(links=np.tile(36, 6999), unwrap=0, num_sims=10, maxtimestep=55, plot_sims=False, simdir='csvs/MC/36bp_notrans'): | ||
| """Plot the end-to-end distance predicted by the theory against the R^2 of the simulated structure.""" | ||
| link = links[0] #linker length (36bp) | ||
| R2, Rmax, kuhn = wlc.R2_kinked_WLC_no_translation(links, plotfig=False, unwraps=unwrap) | ||
| r2 = np.concatenate(([0], R2)) | ||
| rmax = np.concatenate(([0], Rmax)) | ||
| fig, ax = plt.subplots() | ||
|
|
||
| MCr_allsim = np.ones((rmax.size, num_sims)) | ||
| #plot results from the 55th time step of 10 different simulations | ||
| for i in range(1, num_sims): | ||
| rots, pos = extract_MC(dir=simdir, timestep=maxtimestep, v=i) | ||
| #shift MC positions so beginning of polymer coincides with the origin | ||
| pos = pos - pos[0, :] | ||
| #calculate end-to-end distance | ||
| MCr2 = np.array([pos[i, :].T@pos[i, :] for i in range(pos.shape[0])]) | ||
| MCr_allsim[:, i-1] = np.sqrt(MCr2) | ||
| #plot each of the individual simulated end-to-end distances | ||
| if plot_sims: | ||
| if i==1: | ||
| mylabel = 'Simulation' | ||
| else: | ||
| mylabel = None | ||
| plt.loglog(rmax, np.sqrt(MCr2), color=dull_purple, alpha=0.4, label=mylabel) | ||
|
|
||
| plt.loglog(rmax, np.mean(MCr_allsim, axis=1), color=dull_purple, lw=3, label='Simulation') | ||
| plt.loglog(rmax, np.sqrt(r2), color='k', lw=3, label='Theory') | ||
| plt.legend() | ||
| plt.xlabel(r'Total linker length (nm)') | ||
| plt.ylabel(r'$\sqrt{\langle R^2 \rangle}$') | ||
| plt.title('Homogeneous Chain 36bp linkers') | ||
| plt.subplots_adjust(left=0.15, bottom=0.16, top=0.93, right=0.98) | ||
| #plt.savefig('plots/end-to-end-distance-theory-vs-simulation_36bp_0unwraps.png') | ||
|
|
||
|
|
||
|
|
||
|
|
||
|
|
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,81 @@ | ||
| #plotting parameters | ||
| import matplotlib.cm as cm | ||
| import numpy as np | ||
| import seaborn as sns | ||
| import pandas as pd | ||
| from matplotlib import pyplot as plt | ||
| import scipy | ||
| from scipy import stats | ||
| from scipy.optimize import curve_fit | ||
| from nuc_chain import geometry as ncg | ||
| from nuc_chain import linkers as ncl | ||
| from MultiPoint import propagator | ||
| from nuc_chain import fluctuations as wlc | ||
| from nuc_chain.linkers import convert | ||
| from mpl_toolkits.axes_grid1 import make_axes_locatable | ||
|
|
||
| #These follow Andy's plotting preferences for | ||
| params = {'axes.edgecolor': 'black', 'axes.grid': False, 'axes.facecolor': 'white', 'axes.titlesize': 20.0, | ||
| 'axes.linewidth': 0.75, 'backend': 'pdf','axes.labelsize': 18,'legend.fontsize': 18, | ||
| 'xtick.labelsize': 18,'ytick.labelsize': 18,'text.usetex': False,'figure.figsize': [7, 5], | ||
| 'mathtext.fontset': 'stixsans', 'savefig.format': 'pdf', 'xtick.bottom':True, 'xtick.major.pad': 5, 'xtick.major.size': 5, 'xtick.major.width': 0.5, | ||
| 'ytick.left':True, 'ytick.right':False, 'ytick.major.pad': 5, 'ytick.major.size': 5, 'ytick.major.width': 0.5, 'ytick.minor.right':False, 'lines.linewidth':2} | ||
|
|
||
| plt.rcParams.update(params) | ||
|
|
||
| #Mouse data | ||
| def plot_looping_within_contact_radius(a=10, lp=50, loglog=True): | ||
| """Plot probability that 2 ends will form a loop within contact radius a, in nm, | ||
| as a function of dimensionless chain length N=Rmax/2lp, where lp is in nm. | ||
| Plots kinked model vs. bare WLC looping probabilities for both Rlinks and Llinks.""" | ||
|
|
||
| #convert a and lp to basepairs | ||
| a_in_bp = a / ncg.dna_params['lpb'] | ||
| lp_in_bp = lp / ncg.dna_params['lpb'] | ||
|
|
||
| #load in linker lengths used to simulate nucleosome positions in mice (mu=45bp) | ||
| links = np.load('csvs/Bprops/0unwraps/heterogenous/Sarah/mice2/linker_lengths_101nucs.npy') | ||
| #Rlinks -- linkers right of TSS, starting with 60bp between -1 and +1 nuc | ||
| Rlinks = links[50:] | ||
| #Llinks -- linkers left of TSS, starting with 15bp between -2 and -1 nuc | ||
| Llinks = links[0:50] | ||
| #reverse so linker from -1 to -2 nuc is first | ||
| Llinks_rev = Llinks[::-1] | ||
| total_links = np.concatenate((Llinks_rev, Rlinks)) | ||
|
|
||
| #cumulative chain length including burried basepairs | ||
| unwrap = 0 | ||
| #plot as positive distance from TSS in bp | ||
| ldna_Rlinks = convert.genomic_length_from_links_unwraps(Rlinks, unwraps=unwrap) #max WLC chain length in bp | ||
| #plot as negative distance from TSS in bp | ||
| ldna_Llinks = -1*convert.genomic_length_from_links_unwraps(Llinks_rev, unwraps=unwrap) #max WLC chain length in bp | ||
|
|
||
| #load in calculated looping probabilities | ||
| loops = pd.read_csv('../../data_for_Sarah/mice2_looping_probs_a=10nm_102nucs.csv') | ||
|
|
||
| #50th nucleosome or row 49, corresponds to -2nuc | ||
| loops_to_plot = loops['49'] | ||
| #only interested in plotting looping with nucleation site (-2 nuc) | ||
| Lloops = np.concatenate((loops_to_plot[0:49], [loops_to_plot[50]])) | ||
| Lloops = Lloops[::-1] | ||
| #linkers left of -2 nuc | ||
| #ldna_leftnuc = ldna_Llinks[1:] | ||
|
|
||
| #linkers right of -2 nuc | ||
| Rloops = loops_to_plot[51:] | ||
| #ldna_rightnuc = np.concatenate(([ldna_Llinks[0]], ldna_Rlinks)) | ||
|
|
||
|
|
||
| fig, ax = plt.subplots(figsize=(6.25, 4.89)) | ||
| colors = sns.color_palette("BrBG", 9) | ||
| ax.plot(ldna_Rlinks, Rloops, '-o', lw=2, markersize=4, color=colors[1], label='Right of TSS') | ||
| ax.plot(ldna_Llinks, Lloops, '-o', lw=2, markersize=4, color=colors[-2], label='Left of TSS') | ||
| plt.xlabel(r'Distance from TSS (bp)') | ||
| plt.ylabel(f'Looping probability, a={a}nm') | ||
| plt.subplots_adjust(left=0.16, bottom=0.15, top=0.98, right=0.98) | ||
| plt.yscale('log') | ||
| plt.legend(loc=0) | ||
| #plt.ylim([10**-4.5, 10**-1.8]) | ||
| plt.savefig(f'plots/lp{lp}nm_a{a}nm_left_right_contact_probability_mice2.png') | ||
|
|
||
|
|
Oops, something went wrong.