#!/usr/bin/env python # coding: utf-8 # In[1]: get_ipython().run_line_magic('load_ext', 'autoreload') get_ipython().run_line_magic('autoreload', '2') get_ipython().run_line_magic('matplotlib', 'inline') # In[2]: get_ipython().run_line_magic('run', 'nb_init.py') # In[3]: all_kwargs = {'seq_kwargs': {'num_cells': 30, ## The number of apoptotice cells 'width_apopto':60., ## The width of the apoptotic region 'amp': 0., ## The importance of the ventral bias 'seed': 10, ## Seed for the random process 'num_steps': 10, ## The number of step to perform apoptosis 'ventral_bias': True, 'random': True, 'gamma': 1, 'theta_sorted': False}, 'apopto_kwargs': {'vol_reduction':0.5, 'contractility': 1., 'radial_tension': 0.}, 'post_kwargs': {'max_ci':3., 'rate_ci':1.4, 'span_ci':2} } # In[4]: eptm = lj.Epithelium(graphXMLfile='../saved_graphs/xml/before_apoptosis.xml', identifier='joint', copy=True) # In[5]: eptm.isotropic_relax() # In[6]: ax = lj.plot_avg_rho(eptm, bin_width=20) # In[7]: seq_kwargs = all_kwargs['seq_kwargs'] apopto_kwargs = all_kwargs['apopto_kwargs'] post_kwargs = all_kwargs['post_kwargs'] apopto_cells, fold_cells, apopto_sequence = lj.get_apoptotic_cells(eptm, **seq_kwargs) print('%i apoptotic cell(s)' % len(apopto_cells)) print('surrounded by %i fold cells' % len(fold_cells)) # In[8]: eptm.set_local_mask(None) lj.local_slice(eptm, theta_amp=2*np.pi, zed_amp=seq_kwargs['width_apopto']) eptm.update_rhotheta() d_theta = 0. z_angle = np.pi / 6 pseudo_x = eptm.ixs.copy() pseudo_y = eptm.ixs.copy() pseudo_x.a = eptm.zeds.a * np.cos(z_angle) - eptm.rhos.a * np.sin( eptm.thetas.a + d_theta) * np.sin(z_angle) pseudo_y.a = eptm.rhos.a * np.cos(eptm.thetas.a + d_theta) is_apopto = eptm.is_cell_vert.copy() is_apopto.a[:] = 0 color_dead = eptm.zeds.copy() color_dead.a[:] = 0. for cell in apopto_cells: color_dead[cell] = 1. is_apopto[cell] = 1 for jv in cell.out_neighbours(): is_apopto[jv] = 1 # In[12]: ax = lj.plot_eptm_generic(eptm, pseudo_x, pseudo_y, local=True, cell_kwargs={'cell_colors':color_dead, 'cmap': 'Reds', 'alpha':0.8}, edge_kwargs={'c':'g', 'lw':1, 'alpha':0.4}) plt.savefig('../doc/imgs/apoptosis_repartition.svg') # In[10]: mean_rho = eptm.rhos.a.mean() thetas = np.linspace(0, 2 * np.pi, 60) fig, ax = plt.subplots(figsize=(8, 4)) num_rows = 12 delta = 60. for n, cells in enumerate(apopto_sequence[1:]): if len(cells) == 0: break dx = + (n % num_rows) * delta dy = - (n // num_rows) * delta ax.plot(mean_rho * np.cos(thetas) + dx, mean_rho * np.sin(thetas) + dy, 'k-', alpha=0.4) ax.text(dx, dy, str(n + 1), fontsize=10) for cell in cells: ax.plot(eptm.wys[cell] + dx, eptm.ixs[cell] + dy, 'ro', alpha=0.5) ax.set_aspect('equal') ax.set_xticks([]) ax.set_yticks([]) dx = + ((n + 1) % num_rows) * delta dy = - ((n + 1) // num_rows) * delta ax.plot(mean_rho * np.cos(thetas) + dx, mean_rho * np.sin(thetas) + dy, 'k-', alpha=0.4) ax.text(dx, dy, 'all', fontsize=10) for cell in apopto_cells: ax.plot(eptm.wys[cell] + dx, eptm.ixs[cell] + dy, 'ko', alpha=0.7) ax.set_title('Sequence of apoptoses around the joint') plt.show() # In[11]: lj.gradual_apoptosis(eptm, seq_kwargs, apopto_kwargs, post_kwargs) # In[13]: eptm.graph # In[15]: (eptm.is_local_vert.a * (eptm.is_cell_vert.a)).sum() # In[28]: eptm.zeds.fa.size # In[16]: fig, ax = plt.subplots() ax.plot(eptm.zeds.fa, eptm.ixs.fa, 'o', alpha=0.5) # In[47]: for je in eptm.graph.edges(): print(je.source(), je.target(), eptm.zeds[je.source()], eptm.zeds[je.target()], eptm.is_cell_vert[je.source()], eptm.is_cell_vert[je.target()]) # In[22]: eptm.update_geometry() # In[24]: np.all(np.isfinite(eptm.ixs.fa)) # In[25]: eptm.stamp # In[26]: import hdfgraph # In[28]: latest = hdfgraph.graph_from_hdf(eptm.paths['hdf'], stamp=6746) # In[29]: latest # In[44]: eptm.graph # In[33]: vertex_df, edge_df = hdfgraph.frames_from_hdf(eptm.paths['hdf'], stamp=6746) # In[41]: vertex_df.drop_duplicates().xs(0, level='vertex_index') # In[31]: eptm.graph.set_vertex_filter(None) eptm.graph.set_edge_filter(None) # In[43]: eptm.graph.set_vertex_filter(eptm.is_local_vert) eptm.graph.set_edge_filter(eptm.is_local_edge) # In[13]: ax = lj.plot_avg_rho(eptm, bin_width=20) # In[3]: eptm = lj.Epithelium(graphXMLfile='../saved_graphs/joint_2014-03-04T14_02_13/xml/post_optimize.xml', identifier='joint') # In[29]: eptm.identifier # In[28]: # In[4]: lj.draw(eptm, output3d=os.path.join(eptm.paths['png'], 'post_apopto_3d.png'), output2d=os.path.join(eptm.paths['png'], 'post_apopto_2d.png')) # In[5]: import IPython.core.display as disp # In[6]: disp.Image(os.path.join(eptm.paths['png'], 'post_apopto_3d.png')) # In[7]: disp.Image(os.path.join(eptm.paths['png'], 'post_apopto_2d.png')) # In[23]: ax = lj.plot_avg_rho(eptm, bin_width=20) # In[20]: anisotropies, alignments = eptm.cells.get_anisotropies() # In[30]: rel_c = anisotropies #alignments lj.local_slice(eptm, zed_amp=15, theta_amp=np.pi) ax_zs, ax_xy = lj.plot_2pannels(eptm, cell_kwargs={'c_text':False, 'cell_colors':rel_c, 'alpha':0.4}, edge_kwargs={'j_text':False, 'c':'k', 'lw':0.5, 'alpha':0.8}) # In[49]: lj.local_slice(eptm, zed_amp=15, theta_amp=np.pi/2., theta_c=0) ax_zs = lj.plot_eptm_generic(eptm, eptm.zeds, eptm.ixs, cell_kwargs={'c_text':False, 'cell_colors':eptm.zeds, 'alpha':0.4}, edge_kwargs={'j_text':False, 'c':'k', 'lw':0.5, 'alpha':0.8}) # In[24]: fig, axes = plt.subplots(2, 1, sharex=True) eptm.graph.set_vertex_filter(eptm.is_cell_vert) axes[0].plot(eptm.zeds.fa, anisotropies.fa, 'o', alpha=0.4) axes[0].set_ylim(0.8, 2) axes[1].plot(eptm.zeds.fa, alignments.fa, 'o', alpha=0.4) axes[1].set_ylim(-0.1, 1.1) eptm.graph.set_vertex_filter(None) # In[34]: disp.Image('../saved_graphs/joint_2014-03-04T14_02_13/png/eptm3d_4649.png') # In[35]: disp.Image('../saved_graphs/joint_2014-03-04T14_02_13/png/eptm2d_4649.png') # In[27]: import pandas as pd store_path = os.path.join(eptm.paths['hdf']) with pd.get_store(store_path) as store: print(store) all_areas = store.select('vertices', columns=['areas', 'zeds', 'is_cell_vert', 'is_alive']) # In[1]: from mpl_toolkits.mplot3d import Axes3D all_areas = all_areas[all_areas.is_cell_vert > 0] all_areas = all_areas[all_areas.is_alive > 0] stamp = all_areas.index.get_level_values('stamp').values colors = plt.cm.jet(stamp) fig, ax = plt.subplots(subplot_kw={'projection':'3d'}) ax.scatter(stamp, all_areas.zeds.values, zs=all_areas.areas.values, c=stamp, cmap='jet', alpha=0.1) # In[49]: all_areas.index.get_level_values('stamp').unique().size # In[9]: lj.draw(eptm, output3d=os.path.join(eptm.paths['png'], 'after_apopto.png')) # In[ ]: avg_dists = [] for cell in eptm.cells: neighbs = eptm.cells.get_neighbor_cells(cell) xyz = np.array([[eptm.ixs[n] - eptm.ixs[cell], eptm.wys[n] - eptm.wys[cell], eptm.zeds[n] - eptm.zeds[cell]] for n in neighbs]) if not len(xyz): continue avg_dist = np.sqrt((xyz**2).sum(axis=1)).mean() avg_dists.append(avg_dist) avg_dists = np.array(avg_dists) num_sigma = 2*np.pi*eptm.rhos.a.mean() / avg_dists.mean() print('Number of cells around the leg: %.1f' % num_sigma) # In[ ]: width_apopto = fixed_params['width_apopto'] p0 = 1. amp = fixed_params['amp'] thetas = np.linspace(-np.pi, np.pi, 512) zeds = np.linspace(-20, 20, 512) fig, ax = plt.subplots()#subplot_kw={'projection':'3d'}) thetas, zeds = np.meshgrid(thetas, zeds) probas = p_apopto(zeds, thetas, width_apopto=width_apopto, p0=p0, amp=0.4) im = ax.imshow(probas.T, vmin=0, vmax=1, cmap='gray', aspect=20/180., extent=(-20, 20, 0, 360)) _ = plt.colorbar(im) #plt.savefig('../doc/imgs/apoptosis_proba_width_{}_p0_{}_amp_{}.svg'.format(width_apopto, p0, amp)) # In[ ]: fig, axes = plt.subplots(1, 4, figsize=(16, 4)) radii = np.ones(8) width = np.pi/4 avg_hists *= norm for i in range(4): num_depths = avg_hists[i, :] ax = axes[i] # bars = ax.bar(np.arange(0, 2*np.pi, np.pi/4), radii, # width=width, bottom=0.0, ) # for nd, bar in zip(num_depths, bars): # bar.set_facecolor([0, 0, 0]) # bar.set_alpha(nd/avg_hists.max()) labels = ['%1.1f' % nd for nd in avg_hists[i, :]] sizes = np.ones(8) colors = [(1 - nd/avg_hists.max(), 1 - nd/avg_hists.max(), 1 - nd/avg_hists.max()) for nd in avg_hists[i, :]] ax.pie(sizes, colors=colors, labels=labels, startangle=-90) ax.axis('equal') fig.savefig('../doc/imgs/repartition_pie_plot.svg') plt.draw() # In[ ]: all_hists = [] n_cells = [] for seed in range(10): print('Random seed: %i' % seed) apopto_cells = get_apoptotic_cells(eptm, seed=seed, z0=0., width_apopto=width_apopto, p0=p0, amp=amp) apopto_sequence = get_sequence(apopto_cells, num_steps=10) n_seq = len(apopto_sequence) n_cells.append(len(apopto_cells)) apopto_4steps = [np.concatenate( apopto_sequence[n_seq//4 * i: n_seq//4 * i + n_seq//4] ) for i in range(4)] thetas_4steps = [np.unique([eptm.thetas[cell] for cell in seq]) for seq in apopto_4steps] hists_4steps = [np.histogram(thetas, bins=8, range=(-np.pi, np.pi))[0] for thetas in thetas_4steps] all_hists.append(hists_4steps) all_hists = np.array(all_hists) n_cells= np.array(n_cells) print('Average number of apoptic cells: %.1f ± %.2f' % (n_cells.mean(), n_cells.std())) avg_hists = all_hists.mean(axis=0) norm = n_cells.mean() / avg_hists.sum() print('repartition: ') print(avg_hists * norm) print('Summed:') print(avg_hists.sum(axis=0) * norm) in_vivo_hist = np.array([[0.6, 0.6, 0.2, 0.0, 0.3, 0.6, 0.6, 1.0], [1.5, 0.8, 0.3, 0.5, 0.3, 0.7, 2.7, 2.9], [1.2, 2.4, 0.6, 0.8, 1.3, 1.8, 3.5, 1.0], [0.0, 0.0, 0.0, 1.3, 1.4, 0.6, 0.4, 0.0]]) # In[ ]: eptm_unbiased = eptm # In[16]: def average_rho(eptm, bin_width=10): eptm.update_rhotheta() zeds = eptm.zeds.a rhos = eptm.rhos.a rhos = rhos[np.argsort(zeds)] zeds = np.sort(zeds) rhos_cliped = rhos[: -(rhos.size % bin_width)] rhos_cliped = rhos_cliped.reshape((rhos_cliped.size // bin_width, bin_width)) rhos_avg = rhos_cliped.mean(axis=1) rhos_max = rhos_cliped.max(axis=1) rhos_min = rhos_cliped.min(axis=1) zeds_cliped = zeds[: -(zeds.size % bin_width)] zeds_cliped = zeds_cliped.reshape((zeds_cliped.size // bin_width, bin_width)) zeds_avg = zeds_cliped.mean(axis=1) return zeds_avg, rhos_avg, rhos_max, rhos_min def plot_avg_rho(eptm, bin_width, ax=None, retall=False, ls='r-'): if ax is None: fig, ax = plt.subplots(figsize=(12,4)) else: fig = ax.get_figure() zeds_avg, rhos_avg, rhos_max, rhos_min = average_rho(eptm, bin_width) ax.fill_between(zeds_avg, rhos_max, rhos_min, facecolor='0.5', edgecolor='0.9') ax.plot(zeds_avg, rhos_avg, ls, lw=2, alpha=0.7) ax.set_aspect('equal') max_zed = ax.get_ylim()[1] ax.set_ylim(0, max_zed) if not retall: return ax return ax, (zeds_avg, rhos_avg, rhos_max, rhos_min) # In[17]: out = plot_avg_rho(eptm, 20) # In[17]: thetas = np.linspace(0, 2*np.pi, 180)# + np.pi/2 i = eptm.stamp + 1 for n, theta in enumerate(thetas): #output_nt= '../saved_graphs/png/after_apopto_nt_3d/angle_%03i.png' %n output3d= os.path.join(eptm.paths['png'], 'eptm3d_%03i.png' % i) lj.draw(eptm, d_theta=theta, output3d=output3d) i += 1 # In[17]: def average_aniso(eptm, bin_width=10): anisotropies, alignments = eptm.cells.get_anisotropies() eptm.graph.set_vertex_filter(eptm.is_cell_vert) zeds = eptm.zeds.fa anisotropies = anisotropies.fa eptm.graph.set_vertex_filter(None) anisotropies = anisotropies[np.argsort(zeds)] zeds = np.sort(zeds) anisotropies_cliped = anisotropies[: -(anisotropies.size % bin_width)].clip(0, 4) anisotropies_cliped = anisotropies_cliped.reshape((anisotropies_cliped.size // bin_width, bin_width)) anisotropies_avg = anisotropies_cliped.mean(axis=1) anisotropies_max = anisotropies_cliped.max(axis=1) anisotropies_min = anisotropies_cliped.min(axis=1) zeds_cliped = zeds[: -(zeds.size % bin_width)] zeds_cliped = zeds_cliped.reshape((zeds_cliped.size // bin_width, bin_width)) zeds_avg = zeds_cliped.mean(axis=1) return zeds_avg, anisotropies_avg, anisotropies_max, anisotropies_min def plot_avg_aniso(eptm, bin_width, ax=None, retall=False, ls='r-'): if ax is None: fig, ax = plt.subplots(figsize=(12,4)) else: fig = ax.get_figure() zeds_avg, anisotropies_avg, anisotropies_max, anisotropies_min = average_aniso(eptm, bin_width) ax.fill_between(zeds_avg, anisotropies_max, anisotropies_min, facecolor='0.5', edgecolor='0.9') ax.plot(zeds_avg, anisotropies_avg, ls, lw=2, alpha=0.7) max_zed = ax.get_ylim()[1] ax.set_ylim(0, max_zed) if not retall: return ax return ax, (zeds_avg, anisotropies_avg, anisotropies_max, anisotropies_min) # In[19]: ax = plot_avg_aniso(eptm, bin_width=10, ax=None, retall=False, ls='r-') # In[20]: anisotropies, alignments = eptm.cells.get_anisotropies() eptm.graph.set_vertex_filter(eptm.is_cell_vert) fig, ax = plt.subplots() ax.plot(eptm.zeds.fa, anisotropies.fa, 'ro', alpha=0.2) ax.plot(eptm.zeds.fa, alignments.fa, 'bo', alpha=0.2) ax.set_ylim(-2, 4) eptm.graph.set_vertex_filter(None) # In[ ]: