%load_ext autoreload %autoreload 2 %matplotlib inline %run nb_init.py import IPython.core.display as disp import hdfgraph import pandas as pd !ls -lth /media/data/Simulations/ graph_dir = "/media/data/Simulations/" def get_lists(graph_dir, cond): dir_list = os.listdir(graph_dir) dir_list.sort() dir_list = [d for d in dir_list if cond(d)] png_dir_list = [os.path.join(graph_dir, d, 'png') for d in dir_list] svg_list = [os.path.join(graph_dir, d, 'svg', 'avg_rho_%s.svg' % d) for d in dir_list] json_list = [os.path.join(graph_dir, d, 'params_%s.json' %d) for d in dir_list] xml_list = [os.path.join(graph_dir, d, 'xml', 'after_apopto.xml') for d in dir_list] return (dir_list, png_dir_list, svg_list, json_list, xml_list) cond = lambda d: d.startswith('100') or d.startswith('000') dir_list, png_dir_list, svg_list, json_list, xml_list = get_lists(graph_dir, cond) conditions = pd.DataFrame(index=np.arange(len(json_list)), columns=['max_ci', 'radial_tension', 'num_cells', 'width_apopto', 'json']) for i, json_fname in enumerate(json_list): #continue with open(json_fname) as json_file: kwargs = json.load(json_file) print('index: {} => \n'.format(i)) conditions.loc[i, 'json'] = json_fname print('\t max_ci: {}'.format( kwargs['post_kwargs']['max_ci'])) conditions.loc[i, 'max_ci'] = kwargs['post_kwargs']['max_ci'] print('\t radial_tension : {}'.format( kwargs['apopto_kwargs']['radial_tension'])) conditions.loc[i, 'radial_tension'] = kwargs['apopto_kwargs']['radial_tension'] print('\t num cells : {}'.format( kwargs['seq_kwargs']['num_cells'])) conditions.loc[i, 'num_cells'] = kwargs['seq_kwargs']['num_cells'] print('\t width : {}'.format( kwargs['seq_kwargs']['width_apopto'])) conditions.loc[i, 'width_apopto'] = kwargs['seq_kwargs']['width_apopto'] print('\n') conditions xml_paths = {'before': '../leg_joint/data/graphs/before_apoptosis.xml', '5_random': xml_list[0], '10_random': xml_list[1], '20_random': xml_list[2], 'rt_1_ct2': xml_list[3], 'rt_0_ct2': xml_list[4], 'rt_05_ct2': xml_list[5], 'rt_2_ct2': '/media/data/Simulations/0006_2014-05-24T09_04_38/xml/after_apopto.xml', ### backup run 'ectopic': '/media/data/Simulations/joint_2014-06-06T08_48_11/xml/after_apopto.xml', 'rt_1_ct1': xml_list[8], 'rt_1_ct15': xml_list[13], 'rt_1_ct3': xml_list[15] } xml_paths2 = {'before': '../saved_graphs/xml/before_apoptosis.xml', 'intermediate': '../saved_graphs/xml/single_cell_before_elimination.xml', 'single_after': '../saved_graphs/xml/single_cell_after_elimination.xml', 'intermediate_nf': '../saved_graphs/xml/single_cell_before_elimination_noforce.xml', 'single_after_nf': '../saved_graphs/xml/single_cell_after_elimination_noforce.xml'} epithelia = {name: lj.Epithelium(xmlfile, copy=False) for name, xmlfile in [('before', xml_paths['before']),]} eptm_b = epithelia['before'] lj.local_slice(eptm_b, theta_amp=np.pi/8, theta_c=0, zed_amp=2, zed_c=0) local_cells = eptm_b.is_cell_vert.copy() local_cells.a = local_cells.a * eptm_b.is_local_vert.a eptm_b.graph.set_vertex_filter(local_cells) avg_area = eptm_b.cells.areas.fa.mean() avg_rho = eptm_b.rhos.fa.mean() # fig, ax = plt.subplots() # ax.plot(eptm_b.zeds.fa, eptm_b.cells.areas.fa, 'o') eptm_b.graph.set_vertex_filter(None) eptm_b.is_local_vert.a.sum() avg_area, avg_rho def new_2pannels(eptm, rot=45): eptm.update_rhotheta() rot = -(rot + 90) * np.pi / 180 eptm.rotate(rot) lj.local_slice(eptm, theta_amp=np.pi/8, theta_c=0, zed_amp=4, zed_c=0) local_patch = eptm.is_local_vert.copy() lj.local_slice(eptm, theta_amp=np.pi/5, theta_c=0, zed_amp=10, zed_c=0) colors = eptm.rhos.copy() eptm.graph.set_vertex_filter(eptm.is_cell_vert) colors.fa = colors.fa - 5 colors.fa = colors.fa / 25#colors.fa.max() colors.a *= local_patch.a eptm.graph.set_vertex_filter(None) fig, (ax_zs, ax_3d) = plt.subplots(1, 2, figsize=(12, 4), sharey=True) ax_zs = lj.plot_eptm_generic(eptm, eptm.zeds, eptm.proj_sigma(), ax=ax_zs, cell_kwargs={'c_text':False, 'cell_colors':colors, #'c':'k', 'alpha':0.4, 'normalize':False}, edge_kwargs={'j_text':False, 'c':'k', 'lw':0.4, 'alpha':0.5}) eptm.rotate(- rot) rhos, thetas, zeds = eptm.rhos, eptm.thetas, eptm.zeds z_angle=np.pi/24 d_theta = - 6 * np.pi / 7 pseudo_x = eptm.graph.new_vertex_property('float') pseudo_y = eptm.graph.new_vertex_property('float') pseudo_x.a = zeds.a * np.cos(z_angle) - rhos.a * np.cos( thetas.a + d_theta) * np.sin(z_angle) pseudo_y.a = rhos.a * np.sin(thetas.a + d_theta) edge_alpha = eptm.dzeds.copy() v_depth = eptm.zeds.copy() v_depth.a = rhos.a * (1 - np.cos(thetas.a + d_theta)) v_prop = v_depth alpha0 = 0. edge_alpha.a = np.array([(v_prop[s] + v_prop[t]) for s, t in eptm.graph.edges()]) edge_alpha.a -= edge_alpha.a.min() edge_alpha.a += alpha0 * edge_alpha.a.max() edge_alpha.a /= edge_alpha.a.max() ax_3d = lj.plot_edges_generic(eptm, pseudo_x, pseudo_y, ax=ax_3d, efilt=None, edge_alpha=edge_alpha, edge_color=edge_alpha, j_text=False, c='g', lw=0.6) ax_3d.plot([-40, -40], [-15, -20], '-k', lw=2) ax_zs.plot([-20, -20], [-15, -20], '-k', lw=2) for ax in (ax_zs, ax_3d): ax.set_xticks([]) ax.set_yticks([]) ax.set_frame_on(False) ax.set_ylim(-35, 35) plt.savefig(os.path.join(eptm.paths['svg'], 'two_pannels.svg')) import pandas as pd dfs = {} def new_aniso(eptm, name, dfs, rot=45): rot = -(rot + 90) * np.pi / 180 eptm.rotate(rot) sigmas = eptm.wys #ixs #eptm.sigmas #eptm.proj_sigma() coords = [eptm.zeds, sigmas] anisotropies, orientations = eptm.cells.get_anisotropies([eptm.zeds, eptm.sigmas]) areas = eptm.cells.areas rhos = eptm.rhos eptm.graph.set_vertex_filter(eptm.is_cell_vert) n_cells = eptm.is_cell_vert.a.sum() df = pd.DataFrame(np.zeros((n_cells, 5)), columns=['zeds', 'anisotropies', 'orientations', 'areas', 'rhos']) df['zeds'] = eptm.zeds.fa df['anisotropies'] = anisotropies.fa df['rhos'] = rhos.fa / avg_rho df['areas'] = areas.fa / avg_area df['orientations'] = orientations.fa eptm.graph.set_vertex_filter(None) dfs[name] = df lj.local_slice(eptm, theta_amp=np.pi/12, theta_c=0, zed_amp=2, zed_c=0) local_patch = eptm.is_local_vert.copy() lj.local_slice(eptm, theta_amp=np.pi/2, theta_c=0, zed_amp=10, zed_c=0) fig, axes = plt.subplots(1, 4, figsize=(10, 4)) ax_depth, ax_aniso, ax_area, ax_orient = axes ### Depth colors = eptm.rhos.copy() eptm.graph.set_vertex_filter(eptm.is_cell_vert) colors.fa = colors.fa - 5 colors.fa = colors.fa / 25 colors.a = colors.a * local_patch.a eptm.graph.set_vertex_filter(None) ax_depth = lj.plot_eptm_generic(eptm, eptm.zeds, sigmas, ax=ax_depth, cell_kwargs={'c_text': False, 'cell_colors': colors, 'cmap': 'afmhot_r', 'alpha': 0.5,},# 'c':'b', edge_kwargs={'j_text':False, 'c':'k', 'lw':0.4, 'alpha':0.5}) ax_depth.set_title('{}: Depth'.format(name)) ## Anisotropy colors = anisotropies.copy() eptm.graph.set_vertex_filter(eptm.is_cell_vert) colors.fa = colors.fa.clip(0, 4) colors.fa = colors.fa / 4 eptm.graph.set_vertex_filter(None) colors.a = colors.a * local_patch.a ax_aniso = lj.plot_eptm_generic(eptm, eptm.zeds, sigmas, ax=ax_aniso, cell_kwargs={'c_text':False, 'cell_colors':colors, 'cmap':'afmhot_r', 'alpha':1.}, edge_kwargs={'j_text':False, 'c':'k', 'lw':0.4, 'alpha':0.5}) ax_aniso.set_title('Anisotropy') ### Areas colors = eptm.cells.areas.copy() eptm.graph.set_vertex_filter(eptm.is_cell_vert) colors.fa = colors.fa / avg_area colors.fa = colors.fa.clip(0, 1.5) colors.fa = colors.fa / 1.5 eptm.graph.set_vertex_filter(None) colors.a *= local_patch.a ax_area = lj.plot_eptm_generic(eptm, eptm.zeds, sigmas, ax=ax_area, cell_kwargs={'c_text':False, 'cell_colors':colors, 'cmap':'afmhot_r', 'alpha':1.}, edge_kwargs={'j_text':False, 'c':'k', 'lw':0.4, 'alpha':0.5}) ax_area.set_title('Area (µm²)') ### Orientation colors = orientations.copy() eptm.graph.set_vertex_filter(eptm.is_cell_vert) colors.fa = (colors.fa.clip(30, 90) - 30) / 60 eptm.graph.set_vertex_filter(None) colors.a = colors.a * local_patch.a ax_orient = lj.plot_eptm_generic(eptm, eptm.zeds, sigmas, ax=ax_orient, cell_kwargs={'c_text':False, 'cell_colors':colors, 'cmap':'afmhot_r', 'alpha':0.8}, edge_kwargs={'j_text':False, 'c':'k', 'lw':0.4, 'alpha':0.5}) eptm.rotate(- rot) for ax in axes: ax.set_xticks([]) ax.set_yticks([]) ax.set_xlim(-7, 7) ax.set_ylim(-12, 12) ax.set_aspect('equal') ax_orient.set_title('Orientation (°)') saved_name = os.path.join('/home/guillaume/Projets/Droso/Simulations/', name+'_two_pannels.svg') plt.savefig(saved_name) print('saved to {}'.format(saved_name)) for name, path in xml_paths.items(): try: eptm = lj.Epithelium(path, copy=False) except OSError: print('******* Problem with {} *********\n' 'File {} not found \n' '************************\n'.format(name, path)) continue new_aniso(eptm, name, dfs) width = 2. panel = pd.Panel.from_dict({name: df[(- 2 < df['zeds']) & (df['zeds'] < 2)] for name, df in dfs.items()}) panel = panel.swapaxes(0, 2) #panel = panel.swapaxes(1, 2) for name, df in dfs.items(): print(name, df[(- width/2 < df['zeds']) & (df['zeds'] < width/2)].shape) a = np.zeros((3, 4)) b = np.zeros((6, 6)) b[:3, :3] = a panel['anisotropies'] = panel['anisotropies'].clip(0, 7) conditions = {'ncells': ['before', '5_random', '10_random', '20_random', 'rt_1_ct2', 'ectopic'], 'ab_force': ['before', 'rt_0_ct2', 'rt_05_ct2', 'rt_1_ct2', 'rt_2_ct2'], 'contract': ['before', 'rt_1_ct1', 'rt_1_ct15', 'rt_1_ct2', 'rt_1_ct3']} xticks = {'ncells': [0, 5, 10, 20, 30, 'random\npattern'], 'ab_force': ['Ø', 0.0, 0.5, 1.0, 2.0], 'contract': ['Ø', 0.0, 0.5, 1.0, 2.0]} xlabels = {'ncells': 'Number of apoptotic cells', 'ab_force': r'Apico-basal force (units of $\Lambda$)', 'contract': r'Contractility increase (units of $\Gamma$)'} ylabels = {'rhos': 'Normalized radius', 'anisotropies': 'Anisotropy', 'areas': 'Normalized area', 'orientations': 'Orientation (°)'} ylims = {'rhos': (0, 1.4), 'anisotropies': (0, 7.1), 'areas': (0, 1.8), 'orientations': (0, 90)} ncells = {'ncells': conditions['ncells'],} for cond, cols in ncells.items(): for item in ylabels: fig = plt.figure(figsize=(4, 3)) bp = panel[item][cols].boxplot() ax = plt.gca() ax.set_xticklabels(xticks[cond]) ax.set_ylabel(ylabels[item], fontsize=12) ax.set_xlabel(xlabels[cond], fontsize=12) ax.set_ylim(*ylims[item]) fig.tight_layout() plt.savefig('../doc/imgs/{}_bp_{}.svg'.format(item, cond)) def apical_sagital(eptm, ax_a, ax_s, theta_c=0, zed_amp=8, theta_amp=np.pi/12 ): lj.local_slice(eptm, theta_amp=theta_amp, theta_c=theta_c, zed_amp=zed_amp, zed_c=0) colors = eptm.cells.contractilities.copy() eptm.graph.set_vertex_filter(eptm.is_cell_vert) colors.fa = colors.fa - colors.fa.min() colors.fa = colors.fa / 111 colors.fa = colors.fa.clip(0, 1) eptm.graph.set_vertex_filter(None) edge_width = eptm.dzeds.copy() edge_width.a[:] = 0 eptm.graph.set_edge_filter(eptm.is_junction_edge) edge_width.fa = np.array([(eptm.cells.contractilities[eptm.junctions.adjacent_cells[je][0]] + eptm.cells.contractilities[eptm.junctions.adjacent_cells[je][1]])/2 for je in eptm.junctions]) edge_width.fa = (1.4 * edge_width.fa / 400)**2 edge_width.fa = edge_width.fa.clip(0, 2) eptm.graph.set_edge_filter(None) cell_kwargs={'c_text':False, 'cell_colors':colors, 'cmap':'Greens', #'c':'b', 'alpha':0.8, 'normalize':False} edge_kwargs={'j_text':False, 'edge_width':edge_width, 'lw':1., 'c':'g', 'alpha':1.} ax_a = lj.plot_eptm_generic(eptm, eptm.zeds, eptm.wys, ax=ax_a, cell_kwargs=cell_kwargs, edge_kwargs=edge_kwargs) edge_alpha = edge_width.copy() edge_alpha.a[:] = 0. for edge in eptm.junctions.local_junctions(): edge_alpha[edge] = eptm.wys[edge.source()] if ((eptm.wys[edge.source()] > 0) and (eptm.wys[edge.target()] > 0)): eptm.is_local_edge[edge] = 0 edge_alpha[edge] = 0 for cell in eptm.cells.local_cells(): if eptm.thetas[cell] > 0: eptm.is_local_vert[cell] = 0 edge_alpha.a = edge_alpha.a - edge_alpha.a.min() +0.2 edge_alpha.a = edge_alpha.a / edge_alpha.a.max() edge_kwargs['edge_alpha'] = edge_alpha ax_s = lj.plot_eptm_generic(eptm, eptm.zeds, eptm.ixs, ax=ax_s, cell_kwargs=cell_kwargs, edge_kwargs=edge_kwargs) fig, axes = plt.subplots(2, 6, sharex='col', sharey='row', figsize=(12, 4)) eptm_b = epithelia['before'] eptm_i = epithelia['intermediate'] eptm_a = epithelia['single_after'] eptm_inf = epithelia['intermediate_nf'] eptm_anf = epithelia['single_after_nf'] eptms = [eptm_b, eptm_i, eptm_a, eptm_b, eptm_inf, eptm_anf] theta_a_cell = eptm_b.thetas[eptm_a.graph.vertex(3596)] theta_shift = - theta_a_cell for eptm, [ax_s, ax_a] in zip(eptms, axes.T): eptm.rotate(theta_shift) apical_sagital(eptm, ax_a, ax_s, theta_c=0, zed_amp=5, theta_amp=np.pi/12) eptm.rotate(-theta_shift) for ax in axes.ravel(): ax.set_xticks([]) ax.set_yticks([]) ax.set_frame_on(False) plt.savefig('../doc/imgs/single_cell.svg')