diff --git a/environment.yml b/environment.yml index 54c92ee..54b2d4a 100644 --- a/environment.yml +++ b/environment.yml @@ -1,4 +1,4 @@ -name: crossover +name: crossover2 channels: - defaults dependencies: diff --git a/extra_panels/new_model.svg b/extra_panels/new_model.svg index 14f1405..abfe6c5 100644 --- a/extra_panels/new_model.svg +++ b/extra_panels/new_model.svg @@ -1,6 +1,4 @@ - - @@ -227,9 +225,9 @@ borderopacity="1.0" inkscape:pageopacity="0.0" inkscape:pageshadow="2" - inkscape:zoom="2.9685662" - inkscape:cx="183.2647" - inkscape:cy="925.85453" + inkscape:zoom="1.7208721" + inkscape:cx="629.94698" + inkscape:cy="350.38816" inkscape:document-units="mm" inkscape:current-layer="g2765" showgrid="true" @@ -242,7 +240,8 @@ inkscape:snap-bbox="true" inkscape:snap-nodes="false" inkscape:snap-text-baseline="false" - inkscape:bbox-nodes="true"> + inkscape:bbox-nodes="true" + inkscape:document-rotation="0"> @@ -273,7 +272,7 @@ style="fill:#008000;fill-opacity:1;stroke:#008000;stroke-width:0.16277526;stroke-linecap:butt;stroke-linejoin:round;stroke-miterlimit:4;stroke-dasharray:none;stroke-dashoffset:0" /> SC + style="font-size:3.52778px;stroke-width:0.264583;stroke-linecap:butt;stroke-linejoin:round">SC + id="text5565" /> + ns1:pdfconverter="inkscape" + ns1:text="\\begin{tabular}{|c|c|c|}\n\\hline\n$D$ & $1.1 \\, \\mu \\textrm{m}^2 \\textrm{s}^{-1}$ & HEI10 diffusion coefficient \\\\\n$\\alpha$ & $2.1 \\, \\mu \\textrm{m}\\, \\textrm{s}^{-1}$ & RI HEI10 absorption rate \\\\\n$\\beta_C$ & $0.5\\,\\textrm{s}^{-1}$ & RI HEI10 escape rate \\\\\n$\\gamma$ & $1.25$ & Escape rate Hill coefficient \\\\\n$K_C$ & $1 \\, \\textrm{arb. units}$ & Hill function threshold \\\\\n\\hline\n$c_0$ & $1.2 \\,\\textrm{arb. units} \\, \\mu \\textrm{m}^{-1}$ & HEI10 initial loading on bivalent\n\\\\\n$C_0$ & $6.8\\, \\textrm{arb. units}$ & HEI10 initial RI loading \\\\\n$\\sigma$ & $2.2 \\, \\textrm{arb. units}$ & Noise in HEI10 RI loading\\\\\n$x_e$ & $0.1$ & Telomere length (fraction of bivalent)\\\\\n$f_e$ & $2.0$ & Telomere weighting for RI loading\\\\\n$\\rho$ & $ 0.5 \\, \\mu \\textrm{m}^{-1} $ & RI density\\\\\n\\hline\n$T$ & $36000\\, \\textrm{s}$ & Simulation duration\\\\\n\\hline\n\\end{tabular}\n\n\n" + ns1:preamble="/home/foz/.config/inkscape/extensions/textext/default_packages.tex" + ns1:scale="1.11319614602" + ns1:alignment="middle center" + ns1:inkscapeversion="0.92.4" + ns1:jacobian_sqrt="0.392711" + id="g3545"> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + id="id-a9f19014-bbea-4aaf-aed1-d4097e09f6f6" + transform="translate(-148.512, -124.801)"> + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.0010625,0.001 H 298.94425" + transform="matrix(1 0 0 -1 148.712 125.001)" + id="id-63c53655-29c0-4266-b72c-29ee6bda4f42" /> - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.0010625,-2.5e-4 V 11.956781" + transform="matrix(1 0 0 -1 148.712 137.156)" + id="id-f5ad667f-c703-49db-82c8-fdad3923788f" /> + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.001875,-2.5e-4 V 11.956781" + transform="matrix(1 0 0 -1 175.83 137.156)" + id="id-cd9305f9-0a51-4314-9f47-dfead5795554" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 0.00175,-2.5e-4 V 11.956781" + transform="matrix(1 0 0 -1 272.092 137.156)" + id="id-f88bf9a0-3734-44ae-bf52-cdc51718f092" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 2.5e-4,-2.5e-4 V 11.956781" + transform="matrix(1 0 0 -1 447.656 137.156)" + id="id-5f921b3d-59b1-47a5-8a1d-0adb62e1e66f" /> - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.0010625,0.001625 V 11.95475" + transform="matrix(1 0 0 -1 148.712 149.111)" + id="id-4716ef7f-45c1-4a50-b57a-b872c7ef571a" /> + + + + - - - - - - - - - - - - - - - - - - - - - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.001875,0.001625 V 11.95475" + transform="matrix(1 0 0 -1 175.83 149.111)" + id="id-29fec620-8969-4478-812e-d182cfa06ddd" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 0.00175,0.001625 V 11.95475" + transform="matrix(1 0 0 -1 272.092 149.111)" + id="id-aa7f9e60-b897-4eb5-bdb1-eb3b92cc2bed" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 2.5e-4,0.001625 V 11.95475" + transform="matrix(1 0 0 -1 447.656 149.111)" + id="id-377ace4e-f750-4c6a-945e-2640ebddef0c" /> - - - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.0010625,-4.0625e-4 V 11.956625" + transform="matrix(1 0 0 -1 148.712 161.066)" + id="id-91696451-cf47-4849-8f55-408ade6d280a" /> + + + + + + + + + - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.001875,-4.0625e-4 V 11.956625" + transform="matrix(1 0 0 -1 175.83 161.066)" + id="id-5ab67f53-d0de-4649-b396-d87822f66e44" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + + + + + - - + + + + + + + + + + + + + + + + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + + + + + + + + + + - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.001875,-5.625e-4 V 11.956469" + transform="matrix(1 0 0 -1 175.83 184.976)" + id="id-187cc60d-acfd-4b25-9a0a-662f56cebcd0" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - + + + + + + + + + + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 2.5e-4,0.001875 V 11.955" + transform="matrix(1 0 0 -1 447.656 197.33)" + id="id-3ec1581a-ee63-4976-8f99-7cd05f943cf0" /> - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.0010625,-1.5625e-4 V 11.956875" + transform="matrix(1 0 0 -1 148.712 209.285)" + id="id-ebd56ee9-d2d1-4c65-a253-420b79241f99" /> + + + + + + + + + - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.001875,-1.5625e-4 V 11.956875" + transform="matrix(1 0 0 -1 175.83 209.285)" + id="id-8bc63801-1b11-4547-a8ce-ca13107f1c6a" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - + + + + + + - - - - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.001875,0.00171875 V 11.954844" + transform="matrix(1 0 0 -1 175.83 221.24)" + id="id-0808e5d3-cb35-4c3a-84a6-1c066a3d665a" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - + + + + + + + + + + + - - - + + + + + + + + + + + + + + + - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - + + + + + + + + + + + - - - - - - - + + + + + + + + + + + + + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 2.5e-4,-0.00134375 V 11.955688" + transform="matrix(1 0 0 -1 447.656 245.151)" + id="id-59ab4699-43aa-4fcb-82e2-e8ed37b368d1" /> - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.0010625,5.3125e-4 V 11.953656" + transform="matrix(1 0 0 -1 148.712 257.106)" + id="id-ba12c382-f34b-4361-9c26-6b63b595bc13" /> + + + + - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 0.00175,5.3125e-4 V 11.953656" + transform="matrix(1 0 0 -1 272.092 257.106)" + id="id-fd20a08d-39e4-4792-8515-821fe9896fc0" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.001875,-9.375e-4 V 11.956094" + transform="matrix(1 0 0 -1 175.83 269.46)" + id="id-bc872b53-1903-4152-8f5c-8068e3530c25" /> + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 0.00175,-9.375e-4 V 11.956094" + transform="matrix(1 0 0 -1 272.092 269.46)" + id="id-ea785c00-cb48-456f-8c89-d1e5c58b0de8" /> + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M 2.5e-4,-9.375e-4 V 11.956094" + transform="matrix(1 0 0 -1 447.656 269.46)" + id="id-adcc6b2a-fbe4-4e92-ab43-346781761f77" /> + style="fill:none;stroke:#000000;stroke-width:0.398;stroke-linecap:butt;stroke-linejoin:miter;stroke-miterlimit:10;stroke-opacity:1" + d="M -0.0010625,-0.00115625 H 298.94425" + transform="matrix(1 0 0 -1 148.712 269.659)" + id="id-2dcd8e24-2501-45a8-88cb-50d64d48f5a5" /> diff --git a/image_analysis/data_preprocess.py b/image_analysis/data_preprocess.py index 8bd1940..f589ed4 100644 --- a/image_analysis/data_preprocess.py +++ b/image_analysis/data_preprocess.py @@ -20,7 +20,9 @@ ##### Edit this path to point to downloaded dataset data_path = '/media/foz/a92bd7ce-a444-4635-8659-b03fda836a5e/JIC/' -os.makedirs('data_output', exist_ok=True) +data_output_path = '../output/data_output/' + +os.makedirs(data_output_path, exist_ok=True) # Combine multiple dictionaries @@ -251,8 +253,9 @@ def process_data(base_dir, csv_fn, output_fn, category): pickle.dump(data, f, protocol=3) -data_output_path = 'data_output/' -process_data(data_path +'WT Col-0/', '200406.csv', data_output_path+'test.pkl', 'wt') -process_data(data_path+'ox/HEI10 overexpressor/', 'OX.csv', data_output_path+'test_ox.pkl', 'ox') -process_data(data_path+'underexpressor/HEI10 underexpressor/', 'UX.csv', data_output_path+'test_ux.pkl', 'ux') + + +process_data(data_path +'WT Col-0/', '../input_data/200406.csv', data_output_path+'test.pkl', 'wt') +process_data(data_path+'ox/HEI10 overexpressor/', '../input_data/OX.csv', data_output_path+'test_ox.pkl', 'ox') +process_data(data_path+'underexpressor/HEI10 underexpressor/', '../input_data/UX.csv', data_output_path+'test_ux.pkl', 'ux') diff --git a/plotting/CO_length.py b/plotting/CO_length.py new file mode 100644 index 0000000..354d193 --- /dev/null +++ b/plotting/CO_length.py @@ -0,0 +1,74 @@ + + +import svgutils.transform as sg +import base64 +import sys + +def generate_figure(input_julia_path, input_data_path, output_path): + fig0 = sg.SVGFigure( "210mm", "297mm") + + + MPL_SCALE = 0.4*2 + JL_SCALE = 0.08 + + def get_file(fn, scale=MPL_SCALE, pos=None): + fig = sg.fromfile(fn) + plot = fig.getroot() + plot.scale_xy(scale, scale) + if pos is not None: + plot.moveto(*pos) + return plot + + + YS = 160*2 + XS = 185*2 + + + + def make_row(data): + panels = [] + for i, (f, l, c) in enumerate(data): + panels.append(get_file(f, pos=(i*XS, 0))) + panels.append(sg.TextElement(i*XS-10, 0, l, size=20, weight="bold")) + panels.append(sg.TextElement(i*XS+0.5*XS, 0, c, size=24, anchor="middle")) + + return sg.GroupElement(panels) + + + def make_single(data, label=True, label_offset_x=-10): + panels = [] + if label: + f, l = data + panels.append(get_file(f, pos=(0, 0))) + if l: + panels.append(sg.TextElement(label_offset_x, 10, l, size=20, weight="bold")) + else: + f=data + panels.append(get_file(f, pos=(0, 0))) + return sg.GroupElement(panels) + + + + # load matpotlib-generated figures + + + row = [ ( input_data_path+'/violin_number_length.svg', 'a', 'Experimental data'), + ( input_julia_path+'/new_end_nco_vs_length.svg', 'b', 'Simulation output') ] + + g = make_row(row) + g.moveto(0,70) + + + + gpage = sg.GroupElement([g])#, gh]) + gpage.moveto(30,30) + + + fig0.append([gpage]) + + + + # save generated SVG files + fig0.save(output_path+"/CO_length.svg") + +generate_figure(*sys.argv[1:]) diff --git a/plotting/SC_length_vs_CO_number.py b/plotting/SC_length_vs_CO_number.py index d172d12..e202ea8 100644 --- a/plotting/SC_length_vs_CO_number.py +++ b/plotting/SC_length_vs_CO_number.py @@ -24,9 +24,11 @@ utils = importr('utils') -utils.chooseCRANmirror(ind=1) -utils.install_packages('dunn.test') -utils.install_packages('kruskal') +## Uncomment if necessary + +#utils.chooseCRANmirror(ind=1) +#utils.install_packages('dunn.test') +#utils.install_packages('kruskal') import rpy2.robjects as robjects @@ -43,14 +45,13 @@ 'axes.labelsize':28, 'savefig.edgecolor': 'none', 'savefig.facecolor': 'none', - + 'svg.fonttype' : 'none', }) LARGE_FS=32 from matplotlib import lines -LARGE_FS = 32 import cycler @@ -95,7 +96,7 @@ def lin_fit(x, y, min_x=None, of=None, r2=False): else: return X, yy, est2.rsquared -## Analysis of each dataset +## Analysis of each dataset (similar code as for data_analysis_plots) def analyse_data(A, data): all_peak_data = data['all_peak_data'] @@ -241,7 +242,6 @@ def plot_corr(A, image_output_path, data_fn): print(A.keys()) - # Restrict data to SCs from late cells with good quality traces A_late = A.loc[(A.all_stage=='late') & (A.all_quality==1)] @@ -271,17 +271,13 @@ def plot_corr(A, image_output_path, data_fn): print('Data from {} SC'.format(len(A_late))) - # Absolute SC length per crossover number plt.figure() data = [] - NN = 3#np.max(A_late['num_sig_peaks']) + NN = 3 - dl = ' '.join([ str(len(A_late['SC length'][A_late['num_sig_peaks']==i])) for i in range(1, np.max(A_late['num_sig_peaks'])+1) ]) - - print('NN', NN) @@ -294,12 +290,6 @@ def plot_corr(A, image_output_path, data_fn): plt.violinplot(data, positions=range(1,NN+1), showmeans=True, vert=False) -# X, yy, r2 = lin_fit(total_SC_length, total_CO_number, r2=True) -# plt.plot(X, yy, 'r-') - -# plt.title(r'$r^2 = '+ f'{r2:.2f}'+r'$') - -# plt.text(0.1, 0.8, f'r^2 = {r2:.2f}', transform=plt.gca().transAxes) for i in range(3): @@ -310,161 +300,27 @@ def plot_corr(A, image_output_path, data_fn): plt.xticks([20,40,60,80]) plt.yticks([1,2,3]) plt.xlim(15,85) -# plt.subplots_adjust(left=0.2, right=0.9, top=0.9, bottom=0.2) - plt.savefig('violin_number_length.svg') - - plt.title(dl) - - plt.savefig('violin_number_length.png') - - # Relative SC length vc crossover number - - rel_lengths = [] - CO_number = [] - - data = [[] for i in range(5)] - - for i in range(0, len(A), 5): - if A.Stage[i]=='late' and A.all_good[i]=='y': - lengths = A['SC length'][i:i+5] - lengths = lengths/np.sum(lengths) - num_peaks = A['num_sig_peaks'][i:i+5] - for l, n in zip(lengths, num_peaks): - if 1<=n<6: - data[n-1].append(l) - - plt.figure() - plt.violinplot(data, positions=range(1,6), showmeans=True) - plt.xlabel('CO Number') - plt.ylabel('SC rel length') - plt.subplots_adjust(left=0.2, right=0.9, top=0.9, bottom=0.2) - plt.savefig('violin_number_rel_length.svg') - plt.savefig('violin_number_rel_length.png') - - - - - # Total SC length vs total crossover number for each cell? - - total_SC_length = [] - total_CO_number = [] - - for i in range(0, len(A), 5): - if A.Stage[i]=='late' and A.all_good[i]=='y': - total_SC_length.append(np.sum(A['SC length'][i:i+5])) - total_CO_number.append(np.sum(A['num_sig_peaks'][i:i+5])) - - plt.figure() - plt.scatter(total_SC_length, total_CO_number) - - X, yy, r2 = lin_fit(total_SC_length, total_CO_number, r2=True) - plt.plot(X, yy) - - plt.title(f'R2 = {r2:.2f}') - plt.xlabel('total SC length') - plt.ylabel('total CO number') - plt.subplots_adjust(left=0.2, right=0.9, top=0.9, bottom=0.2) - plt.savefig('total_SC_CO.svg') - plt.savefig('total_SC_CO.png') - - - - rel_lengths = [] - abs_lengths = [] - CO_number = [] - - - for i in range(0, len(A), 5): - if A.Stage[i]=='late' and A.all_good[i]=='y': - lengths = A['SC length'][i:i+5] - abs_lengths += list(lengths) - lengths = lengths/np.sum(lengths) - num_peaks = A['num_sig_peaks'][i:i+5] - rel_lengths+= list(lengths) - CO_number += list(num_peaks) - - - print('Max CO number', np.max(CO_number)) - - plt.figure() - plt.scatter(rel_lengths, CO_number) - - X, yy, r2 = lin_fit(rel_lengths, CO_number, r2=True) - plt.plot(X, yy) - plt.title(f'R2 = {r2:.2f}') - plt.xlabel('rel SC length') - plt.ylabel('CO number') - plt.subplots_adjust(left=0.2, right=0.9, top=0.9, bottom=0.2) - plt.savefig('rel_SC_CO.svg') - plt.savefig('rel_SC_CO.png') - - plt.figure() - plt.scatter(abs_lengths, CO_number) - - X, yy, r2 = lin_fit(abs_lengths, CO_number, r2=True) - plt.plot(X, yy) - plt.title(f'R2 = {r2:.2f}') - plt.xlabel('SC length') - plt.ylabel('CO number') - plt.subplots_adjust(left=0.2, right=0.9, top=0.9, bottom=0.2) - plt.savefig('abs_SC_CO.svg') - plt.savefig('abs_SC_CO.png') - - - - - # Single / double / triple total HEI10 per SC - ncells = 0 - norm_hei10 = [] - norm_peaks = dict((i,[]) for i in range(1,4)) - for i in range(0, len(A), 5): - if A.Stage[i]=='late' and A.all_good[i]=='y': - ncells +=1 - hei10 = [] - hei10_all = [] - for j in range(i, i+5): - hei10.append(np.sum(A['new_peak_hei10'][j])) - hei10_all.append(A['new_peak_hei10'][j]) - hei10 = np.array(hei10) - sc_hei10 = hei10/np.sum(hei10) - tot_hei10 = np.sum(hei10) - lengths = A['SC length'][i:i+5] - lengths = lengths/np.sum(lengths) - num_peaks = A['num_sig_peaks'][i:i+5] - for j, h in zip(num_peaks, hei10_all): - #print(j, h) - if j in norm_peaks: - norm_peaks[j].append(np.sum(h)/tot_hei10) - - label_map = { 1:'single', 2:'double', 3:'triple'} - - plt.figure() - for i in range(1,4): - plt.hist(norm_peaks[i], histtype='step', label=label_map[i], density=True, lw=2) - - - plt.xlabel('Relative HEI10 focus intensity') - plt.ylabel('Frequency density') - plt.legend() - plt.subplots_adjust(left=0.2, right=0.9, top=0.9, bottom=0.2) - plt.savefig('single_double_triple_SC_tot_new_intensities.svg') - - - - - plt.show() + plt.savefig(output_path+'/violin_number_length.svg') + + with open(output_path+'/violin_number_length_data.csv', 'w') as f: + f.write('num_sig_peaks, SC_length\n') + for i in range(1, NN+1): + d = data[i-1] + for v in np.array(d): + f.write(str(i) + ', ' + str(v) + '\n') + - +import sys -data_output_path = 'data_output/' -data_output_path2 = 'data_output/' -image_output_path='data_output/' +input_path = '../input_data/' +data_input_path = sys.argv[1] +output_path = sys.argv[2] for image_output_base, csv_fn, data_fn in [ - ( image_output_path, '200406.csv', data_output_path+'test.pkl') + ( output_path, input_path+'/200406.csv', data_input_path+'/test.pkl') ]: A = pd.read_csv(csv_fn) plot_corr(A, image_output_base, data_fn) diff --git a/plotting/all_plots.sh b/plotting/all_plots.sh index 8a50d50..8ecaf6f 100644 --- a/plotting/all_plots.sh +++ b/plotting/all_plots.sh @@ -1,8 +1,10 @@ +set -e + base="../output/" julia_dir=${base}"julia_plots/" kymo_dir=${base}"kymo/" -data_dir="data_output/" +data_dir=${base}"/data_output/" output_dir=${base}"figures/" mkdir -p ${julia_dir} @@ -13,7 +15,9 @@ mkdir -p ${output_dir} python julia_plots.py ${base} ${julia_dir} -python data_analysis_plots.py ${base} +python data_analysis_plots.py ${base}"/data_output/" + +python SC_length_vs_CO_number.py ${data_dir} ${data_dir} for s in "" no_end_ exp_ escape_; do diff --git a/plotting/data_analysis_plots.py b/plotting/data_analysis_plots.py index 10b76d5..064055f 100644 --- a/plotting/data_analysis_plots.py +++ b/plotting/data_analysis_plots.py @@ -26,6 +26,8 @@ utils = importr('utils') +## Uncomment these to install R packages + #utils.chooseCRANmirror(ind=1) #utils.install_packages('dunn.test') #utils.install_packages('kruskal') @@ -37,17 +39,18 @@ stats = importr('stats') # Update matplotlib parameters -mpl.rcParams.update({ #'figure.figsize': (6.0,4.0), - 'figure.facecolor': 'none', #(1,1,1,0), # play nicely with white background in the Qt and notebook -# 'axes.facecolor': 'none', #(1,1,1,0), # play nicely with white background in the Qt and notebook +mpl.rcParams.update({ + 'figure.facecolor': 'none', +# 'axes.facecolor': 'none', 'figure.edgecolor': 'none', - 'font.size': 20, # 12pt labels get cutoff on 6x4 logplots, so use 10pt. - 'figure.dpi': 72, # 72 dpi matches SVG/qtconsole - 'figure.subplot.bottom' : .15, # 10pt still needs a little more room on the xlabel + 'font.size': 20, + 'figure.dpi': 72, + 'figure.subplot.bottom' : .15, 'axes.labelsize':28, 'savefig.edgecolor': 'none', 'savefig.facecolor': 'none', - + 'svg.fonttype' : 'none', + }) @@ -69,7 +72,7 @@ def paler(col_str): # Criterion used for identifying peak as a CO - normalized (with mean and s.d.) -# hei10 levels being above 0.4 time maximum peak level +# hei10 levels being above 0.4 times maximum peak level def pc(peaks, alpha=alpha_max, criterion='max'): pos, v, l = peaks idx = (v>=alpha*np.max(v)) @@ -105,7 +108,6 @@ def analyse_data(A, data): o_hei10_traces = data['o_hei10_traces'] orig_trace_lengths = data['orig_trace_lengths'] - #print(len(A), len(all_peak_data), len(o_hei10_traces), len(orig_trace_lengths)) stages = [] cell_id = [] @@ -117,7 +119,7 @@ def analyse_data(A, data): orig_peak_hei10 = [] - ### For each cell in turn, process the + ### For each cell in turn for i in range(0, len(A), 5): stage = A.iloc[(i//5)*5].Stage @@ -143,8 +145,6 @@ def analyse_data(A, data): orig_peak_hei10.append(h[pp[0]]) h = h-np.median(h) - #h=(h-median_all) #/std_all - #h = (h-np.mean(h))/np.std(h) new_peak_hei10.append(h[pp[0]]) @@ -173,7 +173,7 @@ def analyse_data(A, data): # Sort the 5 chromosomes in each cell by length - # Fo each cell, sort the chromosomes in order + # For each cell, sort the chromosomes in order def analyse_SC_lengths(df): ch_sort_idx = [] @@ -273,7 +273,7 @@ def hist_peak_no(df): # Histograms of the positions of all of the hei10 peaks # (positions normalized to length of SC - # Data can be made symmetric about 0.5 + # Data can be made symmetric about 0.5 (not used in paper) # Also make stacked histograms (showing the distributions of the 1st / 2nd / 3rd / etc crossovers def hist_peak_pos(df, make_symmetric=False, stacked=0, summary_file=None): all = [] @@ -326,24 +326,20 @@ def all_summary(df, summary_file): print('Num late SC ', np.sum(A.all_stage=='late'), file=summary_file) - # Number of cells this data is drawn from for reporting summary + # Number of cells this data (i.e. everything in A_late) is drawn from for reporting summary nc = 0 for i in range(0, len(A.index), 5): if A.all_stage[i] == 'late' and any(A.quality[j]==1 for j in range(i, i+5)): nc += 1 - #print(A.iloc[i], file=summary_file) print("USED", i, A.all_stage[i], 'All quality:', A.all_quality[i], 'All good:', A.all_good[i], file=summary_file) else: print("NOT USED", i, A.all_stage[i], 'All quality: ', A.all_quality[i], 'All good: ', A.all_good[i], file=summary_file) -# print(A.iloc[i], file=summary_file) -# print(A.quality[i:i+5], file=summary_file) print('SCs from ' + str(nc) + ' cells', file=summary_file) - # Restrict data to SCs from late cells with good quality traces A_late = A.loc[(A.all_stage=='late') & (A.quality==1)] @@ -468,8 +464,6 @@ def all_summary(df, summary_file): plt.plot(data_c, data_x, 'rx') plt.plot(X, yy, 'b') plt.title('N='+str(len(data_x)) + f' $R^2$ = {r2:.2f}') - # plt.text(0.1, 0.9, , transform = plt.gca().transAxes) - #plt.xlabel('Bivalent fraction closest to CO') plt.xlabel('Midpoint relative position') plt.ylabel('Left CO relative intensity', size=22) plt.savefig(image_output_base+'new_centre.svg') @@ -542,12 +536,10 @@ def all_summary(df, summary_file): plt.legend() plt.savefig(image_output_base+'single_double_triple_peak_new_intensities.svg') - + + # Write normalized peaks into text file - with open(image_output_base+'single_double_triple.json', 'w') as of123: - json.dump(norm_peaks, of123) - @@ -595,21 +587,12 @@ def all_summary(df, summary_file): p_adj = pvals_adj[idx] c = compare[idx] -# _, p_val, _ = ttest_ind(norm_peaks[i0+1], norm_peaks[i1+1]) -# u_val, p_val = mannwhitneyu(norm_peaks[i0+1], norm_peaks[i1+1], alternative='two-sided') - -# print('ttest_ind', i0+1, i1+1, ttest_ind(norm_peaks[i0+1], norm_peaks[i1+1]), 'n=', len(norm_peaks[i0+1]), len(norm_peaks[i1+1]), file=summary_file) -# print('mann_whitney', i0+1, i1+1, mannwhitneyu(norm_peaks[i0+1], norm_peaks[i1+1]), 'n=', len(norm_peaks[i0+1]), len(norm_peaks[i1+1]), file=summary_file) - fig.text(0.87, 0.5*(y_pos[i0]+y_pos[i1]), f'p={p_adj:.1e}', fontsize=15, rotation=90, ha='center', va='center') - #fig.text(0.88, 0.5*(y_pos[i0]+y_pos[i1]), f'Z={Z_val:.1f} ', fontsize=8, rotation=90, ha='center', va='center') fig.lines.append(l) for i0, i1, idx in [[0,2, 1]]: l = lines.Line2D([0.81, 0.91, 0.91, 0.81], [y_pos[i0]+0.02, y_pos[i0]+0.02, y_pos[i1]-0.02, y_pos[i1]-0.02], transform = fig.transFigure, lw=2, color='k') - #_, p_val, _ = ttest_ind(norm_peaks[i0+1], norm_peaks[i1+1]) -# u_val, p_val = mannwhitneyu(norm_peaks[i0+1], norm_peaks[i1+1], alternative='two-sided') Z_val = Z[idx] p_val = pvals[idx] @@ -618,19 +601,20 @@ def all_summary(df, summary_file): fig.text(0.95, 0.5*(y_pos[i0]+y_pos[i1]), f'p={p_adj:.1e}', fontsize=15, rotation=90, ha='center', va='center') -# fig.text(0.96, 0.5*(y_pos[i0]+y_pos[i1]), f'Z={Z_val:.2f}', fontsize=8, rotation=90, ha='center', va='center') fig.lines.append(l) - # print('ttest_ind', i0+1, i1+1, ttest_ind(norm_peaks[i0+1], norm_peaks[i1+1]), 'n=', len(norm_peaks[i0+1]), len(norm_peaks[i1+1]), file=summary_file) -# print('mann_whitney', i0+1, i1+1, mannwhitneyu(norm_peaks[i0+1], norm_peaks[i1+1]), 'n=', len(norm_peaks[i0+1]), len(norm_peaks[i1+1]), file=summary_file) ax= fig.add_axes([0.15, 0.15, 0.7, 0.8], frame_on=False) plt.tick_params(labelcolor="none", bottom=False, left=False) - #ax.set_xlabel('Relative HEI10 focus intensity', size=20) - #ax.set_ylabel('Frequency density', size=20) plt.savefig(image_output_base+'single_double_triple_peak_new_intensities_stacked.svg') + with open(image_output_base+'single_double_triple_peak_new_intensities_stacked_data.csv', 'w') as f: + f.write('Number_of_foci_on_SC,Relative_HEI10_focus_intensity\n') + for i, d in norm_peaks.items(): + for dd in d: + f.write(str(i)+ ',' + str(dd)+'\n') + print('single double triple kruskal', kruskal(norm_peaks[1], norm_peaks[2], norm_peaks[3]), file=summary_file) print('single double triple anova', f_oneway(norm_peaks[1], norm_peaks[2], norm_peaks[3]), file=summary_file) @@ -641,7 +625,7 @@ def all_summary(df, summary_file): for prefix in ['new']: for num_co, title, out_fn in [ ((-1, 10000), 'All late SC', '_rel_length.svg'), ((2,2), 'All late SC with double CO', '_rel_length2.svg' ) ]: - + # Normalize sum of hei10 peaks for each cell plt.figure() norm_lengths = [] @@ -665,23 +649,23 @@ def all_summary(df, summary_file): print(title+' Relative bivalent length vs Relative total focus HEI10 per bivalent', file=summary_file) xx, yy, r2 = lin_fit(norm_lengths, norm_hei10, min_x=0, of=summary_file, r2=True) plt.plot(xx, yy,'r-') -# plt.title(title) plt.xlim(xmin=0) if num_co!=2: plt.ylim(ymin=0) plt.text(0.1, 0.9, f'$R^2$ = {r2:.2f}', transform = plt.gca().transAxes) -# plt.xlabel('Relative bivalent length') -# plt.ylabel('Relative total focus HEI10 per bivalent') -# plt.tight_layout() plt.savefig(image_output_base+prefix+out_fn) + with open(image_output_base+prefix+out_fn[:-4]+'_data.csv', 'w') as f: + f.write('Relative_bivalent_SC_length,Relative_total_focal_HEI10_per_bivalent\n') + for x,y in zip(norm_lengths, norm_hei10): + f.write(str(x) + ',' + str(y) + '\n') + source_data = {} # Fig 1b left panel plt.figure() - - + used_sc = {} for stage, criterion, col in [ ('late', lambda s: s=='late', 'r'), ('mid', lambda s: s=='mid' or s=='Mid', 'b'), @@ -704,11 +688,8 @@ def all_summary(df, summary_file): for j in range(i, i+5): h = o_hei10_traces[j] - #h = np.maximum(h - np.median(h), 0) - h = h - np.median(h) #median_all + h = h - np.median(h) - #h = h-median_all - #h = np.maximum(h - median_all, 0) hei10.append(np.sum(h)) used_sc[stage] +=1 @@ -722,23 +703,27 @@ def all_summary(df, summary_file): print(stage, np.mean(norm_lengths), np.std(norm_lengths), file=summary_file) plt.plot(norm_lengths, norm_hei10, col+'o', markersize=2) + source_data[stage] = (norm_lengths, norm_hei10) + print(stage + ' Relative bivalent length vs Relative total focus HEI10 per bivalent', file=summary_file) xx, yy, r2 = lin_fit(norm_lengths, norm_hei10, min_x=0, of=summary_file, r2=True) -# plt.plot(xx,yy,col+'-', label=stage + f' $R^2$ = {r2:.2f}' ) plt.plot(xx,yy,col+'-', label=stage + f' {r2:.2f}' ) - # plt.title('Total HEI10 above median') plt.xlim(xmin=0) plt.ylim(ymin=0) -# plt.xlabel('Relative bivalent SC length') -# plt.ylabel('Relative total HEI10 per SC') plt.legend(fontsize=18) - # plt.tight_layout() plt.savefig(image_output_base+'rel_length_tot_compare.svg') + with open(image_output_base+'rel_length_tot_compare_data.csv', 'w') as f: + f.write('stage, Relative_bivalent_SC_length,Relative_total_HEI10_per_bivalent\n') + for stage, (lengths, hei10) in source_data.items(): + for p in zip(lengths, hei10): + f.write(stage + ',' + str(p[0]) + ',' + str(p[1]) + '\n') + + print('rel_length_tot_compare SC numbers', used_sc, file=summary_file) ### Calulate late total foci number per cell and total SC length per cell. @@ -746,23 +731,19 @@ def all_summary(df, summary_file): cell_idx = [] cell_total_SC = [] -# cell_num_peaks = [] sc_length_all = [] for i in range(0, len(A), 5): if A.Stage[i]=='late' and A.all_good[i]=='y': total_SC = 0 -# total_peaks = 0 sc_length = [] for j in range(i, i+5): total_SC += A['SC length'][j] -# total_peaks += A['num_sig_peaks'][j] sc_length.append(A['SC length'][j]) sc_length.sort() cell_idx.append(i//5) cell_total_SC.append(total_SC) -# cell_num_peaks.append(total_peaks) sc_length_all.append(sc_length) @@ -773,7 +754,6 @@ def all_summary(df, summary_file): for i in range(0, len(A), 5): if True: #A.Stage[i]=='late': total_SC = 0 -# total_peaks = 0 for j in range(i, i+5): total_SC += A['SC length'][j] cell_total_SC_all.append(total_SC) @@ -844,20 +824,22 @@ def plot_traces(data_fn, csv_fn, idx_list, image_output_path): base = sys.argv[1] -data_output_path = base+'/data_output/' -image_output_path='data_output/' + +input_path = '../input_data/' +data_output_path = base +image_output_path = base for image_output_base, csv_fn, data_fn, max_n, stacked_bins, summary_fn, summary_fn2 in [ - ( image_output_path, '200406.csv', data_output_path+'test.pkl', 4, np.linspace(0, 0.3, 30), data_output_path+'summary.txt', data_output_path+'cell_summary.txt'), - ( image_output_path+'ox_', 'OX.csv', data_output_path+'test_ox.pkl', 4, np.linspace(0, 0.12, 11), data_output_path+'summary_ox.txt', data_output_path+'cell_summary_ox.txt'), - ( image_output_path+'ux_', 'UX.csv', data_output_path+'test_ux.pkl', 3, np.linspace(0, 0.3, 10), data_output_path+'summary_ux.txt', data_output_path+'cell_summary_ux.txt'), + ( image_output_path, input_path+'200406.csv', data_output_path+'test.pkl', 4, np.linspace(0, 0.3, 30), data_output_path+'summary.txt', data_output_path+'cell_summary.txt'), + ( image_output_path+'ox_', input_path+'OX.csv', data_output_path+'test_ox.pkl', 4, np.linspace(0, 0.12, 11), data_output_path+'summary_ox.txt', data_output_path+'cell_summary_ox.txt'), + ( image_output_path+'ux_', input_path+'UX.csv', data_output_path+'test_ux.pkl', 3, np.linspace(0, 0.3, 10), data_output_path+'summary_ux.txt', data_output_path+'cell_summary_ux.txt'), ]: A = pd.read_csv(csv_fn) plot_data(A, image_output_base, data_fn, max_n, stacked_bins, summary_fn, summary_fn2) print(image_output_path) # Plots of HEI10 traces along individual chromosomes (used for panels in Fig1A, 3A and 4A) -plot_traces(data_output_path+'test.pkl', '200406.csv', [675, 664, 492 ], image_output_path) -plot_traces(data_output_path+'test_ox.pkl', 'OX.csv', [ 260], image_output_path) -plot_traces(data_output_path+'test_ux.pkl', 'UX.csv', [ 115], image_output_path) +plot_traces(data_output_path+'test.pkl', input_path+'200406.csv', [675, 664, 492 ], image_output_path) +plot_traces(data_output_path+'test_ox.pkl', input_path+'OX.csv', [ 260], image_output_path) +plot_traces(data_output_path+'test_ux.pkl', input_path+'UX.csv', [ 115], image_output_path) diff --git a/plotting/f_plot.py b/plotting/f_plot.py index 87ef23d..e4cf094 100644 --- a/plotting/f_plot.py +++ b/plotting/f_plot.py @@ -27,6 +27,6 @@ plt.xlim(0,1) plt.xticks([0,1]) plt.yticks([0,1,2]) -plt.savefig('f_plot.svg') +plt.savefig('../output/julia_plots/f_plot.svg') diff --git a/plotting/fig1_wide.py b/plotting/fig1_wide.py index 3e45f9a..49af041 100644 --- a/plotting/fig1_wide.py +++ b/plotting/fig1_wide.py @@ -32,12 +32,12 @@ def get_mpl_base(fn, scale=MPL_SCALE, pos=None): return plot - chris_fig = get_file('extra_panels/HEI10 figure 1.14.svg', scale=1.0) + chris_fig = get_file('../extra_panels/HEI10 figure 1.14.svg', scale=1.0) - fig_panels = [ get_mpl_base('data_output/rel_length_tot_compare.svg', pos=(14,123)), - get_mpl_base('data_output/new_rel_length.svg', pos=(65,123)), - get_mpl_base('data_output/new_rel_length2.svg', pos=(111,123)), - get_mpl_base('data_output/single_double_triple_peak_new_intensities_stacked.svg', pos=(155,123))] + fig_panels = [ get_mpl_base(input_data_path+'/rel_length_tot_compare.svg', pos=(14,123)), + get_mpl_base(input_data_path+'/new_rel_length.svg', pos=(65,123)), + get_mpl_base(input_data_path+'/new_rel_length2.svg', pos=(111,123)), + get_mpl_base(input_data_path+'/single_double_triple_peak_new_intensities_stacked.svg', pos=(155,123))] gpage = sg.GroupElement([chris_fig]+fig_panels) diff --git a/plotting/fig2_wide.py b/plotting/fig2_wide.py index 5d017e8..564f04b 100644 --- a/plotting/fig2_wide.py +++ b/plotting/fig2_wide.py @@ -5,6 +5,7 @@ def generate_figure(input_julia_path, input_data_path, output_path): fig0 = sg.SVGFigure( "210mm", "297mm") + fig0.root.set("viewBox", "0 0 210 297") MPL_SCALE = 0.4 @@ -140,7 +141,7 @@ def make_single(data, label=True, label_offset_x=-10): gall.scale_xy(0.9, 0.9) gall.moveto(30,280) - model = get_file("extra_panels/new_model.svg", scale=4) + model = get_file("../extra_panels/new_model.svg", scale=4) model.moveto(-10, -60) @@ -150,7 +151,8 @@ def make_single(data, label=True, label_offset_x=-10): gmodel = sg.GroupElement([model, alabel, blabel]) gpage = sg.GroupElement([gmodel, gall]) - + gpage.scale_xy(0.26,0.26) + fig0.append([gpage]) diff --git a/plotting/fig3_wide.py b/plotting/fig3_wide.py index fecf605..e324c1c 100644 --- a/plotting/fig3_wide.py +++ b/plotting/fig3_wide.py @@ -9,6 +9,7 @@ def generate_figure(input_julia_path, input_data_path, output_path): fig0 = sg.SVGFigure( "210mm", "297mm") + fig0.root.set("viewBox", "0 0 210 297") MPL_SCALE = 0.4 @@ -68,10 +69,10 @@ def make_single(data, label=True, label_offset_x=-10): # load matpotlib-generated figures - image_data = base64.b64encode(open("extra_panels/HEI10 fig 3A.png", "rb").read()) + image_data = base64.b64encode(open("../extra_panels/HEI10 fig 3A.png", "rb").read()) - image_panel = sg.ImageElement(open("extra_panels/HEI10 fig 3A.png", "rb"), 1584, 444) + image_panel = sg.ImageElement(open("../extra_panels/HEI10 fig 3A.png", "rb"), 1584, 444) image_panel.moveto(10,10) image_panel.scale_xy(0.48, 0.48) @@ -114,8 +115,9 @@ def make_single(data, label=True, label_offset_x=-10): gall = sg.GroupElement([g, g2, gh, gh2]) gall.moveto(30,270) - + gpage = sg.GroupElement([g_image, gall]) + gpage.scale_xy(0.26, 0.26) fig0.append([gpage]) diff --git a/plotting/fig4_wide.py b/plotting/fig4_wide.py index 0302078..18c3bea 100644 --- a/plotting/fig4_wide.py +++ b/plotting/fig4_wide.py @@ -9,6 +9,7 @@ def generate_figure(input_julia_path, input_data_path, output_path): fig0 = sg.SVGFigure( "210mm", "297mm") + fig0.root.set("viewBox", "0 0 210 297") MPL_SCALE = 0.4 @@ -70,7 +71,7 @@ def make_single(data, label=True, label_offset_x=-10): # load matpotlib-generated figures - image_panel = sg.ImageElement(open("extra_panels/HEI10 fig 4A.png", "rb"), 1584, 444) + image_panel = sg.ImageElement(open("../extra_panels/HEI10 fig 4A.png", "rb"), 1584, 444) image_panel.moveto(10,10) image_panel.scale_xy(0.48, 0.48) @@ -117,7 +118,8 @@ def make_single(data, label=True, label_offset_x=-10): gall.moveto(30,270) gpage = sg.GroupElement([g_image, gall]) - + gpage.scale_xy(0.26, 0.26) + fig0.append([gpage]) diff --git a/plotting/julia_plots.py b/plotting/julia_plots.py index 73448a8..d620f60 100644 --- a/plotting/julia_plots.py +++ b/plotting/julia_plots.py @@ -1,10 +1,16 @@ import numpy as np +import os import sys import matplotlib.pyplot as plt import matplotlib as mpl import csv import random +import pickle + +import statsmodels.api as sm +import pandas as pd + mpl.rcParams.update({ 'figure.facecolor': 'none', @@ -15,20 +21,19 @@ 'axes.labelsize':28, 'savefig.edgecolor': 'none', 'savefig.facecolor': 'none', + 'svg.fonttype' : 'none', }) LARGE_FS=32 -import pickle - -import statsmodels.api as sm - +# Peak criterion def pc(v, alpha=0.4): idx = (v>=alpha*np.max(v)) return idx +# Code to read output data file def read_file(fn): all_data = [] @@ -65,11 +70,10 @@ def read_file(fn): return head_dict, all_data -import os +# Make histograms, length and foci intensity arrays for the data from a single file def summary(data, L, K, tp=-1, max_k=4): - print('tp', tp) head, all_data = data all_foci = [] @@ -81,42 +85,32 @@ def summary(data, L, K, tp=-1, max_k=4): for L, pos, hei10 in all_data: if len(pos) and len(hei10) and tp=max_k: @@ -124,11 +118,12 @@ def summary(data, L, K, tp=-1, max_k=4): hist_nn[max_k], _ = np.histogram(f, bins=np.linspace(0,1,11)) - hist, _ = np.histogram(foci, bins = np.linspace(0,1,11)) + hist, _ = np.histogram(foci, bins = np.linspace(0,1,11)) # hist_nn[k] = Histogram of all relative foci positions return hist_n, hist_nn, hist, len(foci), len(all_foci), L_array, foci_intensities, orig_foci_intensities, all_foci +### Linear regression def lin_fit(x, y, r2=False): X = np.array(x) @@ -146,6 +141,7 @@ def lin_fit(x, y, r2=False): return X, yy, est2.rsquared +# Plot for fig 2k - for SCs with two COs, rel intensity of left peak vs position of centre def plot_centre(all_data, tp=-1): c = [] @@ -154,7 +150,6 @@ def plot_centre(all_data, tp=-1): for L, pos, hei10 in all_data: if len(pos) and len(hei10) and tp2*np.mean(hei10[tp]): idx = pc(hei10[tp]) idx = np.where(idx)[0] @@ -173,8 +168,6 @@ def plot_centre(all_data, tp=-1): c = np.array(c) v = np.array(v) - with open('centre.pkl', 'wb') as f: - pickle.dump([c,v,s], f) plt.plot(c[s], v[s], 'rx') @@ -187,42 +180,22 @@ def plot_centre(all_data, tp=-1): -def plot_diff(all_data, order, tp=-1): - - - d = [] - - for L, pos, hei10 in all_data: - - if len(pos) and len(hei10) and tp2*np.mean(hei10[tp]): - idx = pc(hei10[tp]) - idx = np.where(idx)[0] - - if len(idx)==order: - p = pos[idx]/L - d += list(np.diff(p)) - plt.figure() - plt.hist(d) - plt.title("{:.4g} {:.4g}".format(np.mean(d), np.std(d))) - plt.autoscale(enable=True, axis='x', tight=True) - +# Plot spacing between consecutive foci def plot_diff_all(all_data, rel=False, tp=-1): - d = [] for L, pos, hei10 in all_data: if len(pos) and len(hei10) and tp2*np.mean(hei10[tp]): idx = pc(hei10[tp]) idx = np.where(idx)[0] - if rel: - p = pos[idx]/L - else: - p = pos[idx] - d += list(np.diff(p)) + if len(idx)=len(survey_base) and f[:len(survey_base)] == survey_base: - print(f) params = f[len(survey_base):-4] all_data_ext[params] = read_file(survey_dir+'/'+f) @@ -260,15 +233,14 @@ def get_class(params): new_data = {} for k in all_data_ext: - p = get_class(k)[0:8] + get_class(k)[9:] + p = get_class(k)[0:8] + get_class(k)[9:] # Omit the starting index from file if p in new_data: new_data[p][1] += all_data_ext[k][1] else: h, d = all_data_ext[k] new_data[p] = [h, d] - print(list(new_data)) - + return new_data @@ -278,55 +250,35 @@ def to_val(v): except ValueError: return float(v) -def line_hist(ax, x, data, label=None): - data = [0] + list(np.repeat(data, 2)) + [0] - xx = np.repeat(x, 2) - print(data, len(data)) - print(xx, len(xx)) - ax.plot(xx, data, label=label) - # plt.ylim(0, np.max(data)*1.05) - -import pandas as pd + +## Process all output files for WT/OX/UC etc def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, max_n=4): # Load preprocessed pkl file new_data = load_data(*data_path) k = next(iter(new_data)) - print(k, new_data[k][0]) - print('loaded data') - print(list(new_data.values())[0]) - + # Organize data into groups new_data_keys = [ (idx, to_val(h['L']), to_val(h['density'])) for idx, (h, v) in enumerate(new_data.values()) ] idx, L, density = map(list, zip(*new_data_keys)) df = pd.DataFrame(data={'idx':idx, 'L':L, 'density':density}) df2 = df.sort_values(['L', 'density']) - - print(df2) df2['density_idx'] = list(range(3))*5 - print(df2) - - - grouped = df2.groupby(['L']) - for name, group in grouped: - print(name) - print(group) - print('***') new_data_list = list(new_data.values()) - - all_hist_n = np.zeros((19,)) all_hist_nn = np.zeros((4, 10)) all_hist = np.zeros(10,) + + all_lengths = [] all_intensities = [] @@ -337,27 +289,23 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, nf_tot = [] + all_nco = [] + + all_data_1 = [] all_co = [] for group_idx, (name, group) in enumerate(grouped): - print('name', name) s = group.iloc[-1,:] -# print('s', s) - print("s['idx']", s['idx']) h, v = new_data_list[int(s['idx'])] all_data_1 += v - print('h, ', h) - L = to_val(h['L']) density = to_val(h['density']) K = to_val(h['K']) - print('L, density, K', L, density, K) - hist_n, hist_nn, hist, nf, nch, L_array, intensities, orig_foci_intensities, co = summary((h,v), L, K, tp=-1) all_intensities.append(orig_foci_intensities) @@ -369,15 +317,14 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, all_co += co - - print(hist, len(hist)) - + all_lengths += list(L_array) + all_nco += [len(z) for z in co] + dsb_density = [] nf_all = [] for i in range(3): s = group.iloc[i,:] - print(s) j = int(s['idx']) dsb_density.append(s['density']) _, _, _, nf, nch, _, _, _, _ = summary(new_data_list[j], L, K, tp=-1) @@ -387,6 +334,38 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, + + plt.figure() + + length_by_nco = [ np.array(all_lengths)[np.array(all_nco)==n] for n in range(1, np.max(all_nco)+1) ] + + N_data = [str(len(p)) for p in length_by_nco] + + + NN = 3 + + length_by_nco = [ np.array(all_lengths)[np.array(all_nco)==n] for n in range(1, NN+1) ] + + + plt.violinplot(length_by_nco, positions=np.arange(1,4), showmeans=True, vert=False) + + for i in range(3): + plt.annotate('n='+str(len(length_by_nco[i])), (20, i+1+0.2), fontsize=18) + + plt.xlabel('SC length ($\mu$m)') + plt.ylabel('CO number') + plt.xticks([20, 40,60,80]) + plt.yticks([1,2,3]) + plt.xlim(15,85) + + plt.savefig(output_prefix+"nco_vs_length.svg") + + with open(output_prefix+'nco_vs_length_data.csv', 'w') as f: + f.write('CO_number,SC_length\n') + for i in range(3): + for v in length_by_nco[i]: + f.write(str(i+1) + ','+str(v) + '\n') + mean_co = np.mean([len(c) for c in co]) if centre_plot: @@ -401,11 +380,11 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, - with open('H2A.csv') as f: + with open('../input_data/H2A.csv') as f: reader = csv.reader(f) H2A = list(reader) - with open('MHL1_diakinesis.csv') as f: + with open('../input_data/MHL1_diakinesis.csv') as f: reader = csv.reader(f) MHL = list(reader) @@ -413,7 +392,6 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, H2A = np.array([float(p[1]) for p in H2A[:3]]) MHL = np.array([float(p[1]) for p in MHL[:3]]) - print(H2A, MHL) dsb_data = H2A/H2A[0] CO_data = MHL @@ -455,10 +433,6 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, plt.savefig(output_prefix+'julia_all_hist.svg') - print(all_hist_n) - print('mean peaks all', np.sum(all_hist_n*np.arange(19))/np.sum(all_hist_n), np.sum(all_hist_n*np.arange(19))/np.sum(all_hist_n)*5) - - print('MHL', MHL) fig, ax = plt.subplots(2,2) for i in range(4): @@ -530,7 +504,6 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, plt.close() - print(H2A, MHL) dsb_density = np.array(dsb_density) dsb_density = dsb_density/dsb_density[-1] @@ -540,16 +513,12 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, CO_data = CO_data/CO_data[0] - - print('nf tot', nf_tot) plt.figure() nf2 = np.mean(np.array(nf_tot), axis=0) CO_sim = nf2/1000.0 CO_sim = CO_sim/CO_sim[-1] - - print(dsb_density, CO_sim) plt.plot(dsb_density, CO_sim) plt.plot(dsb_data, CO_data, 'rx') @@ -576,10 +545,6 @@ def count_CO(data_path): grouped = df2.groupby(['L']) - for name, group in grouped: - print(name) - print(group) - print('***') new_data_list = list(new_data.values()) @@ -590,19 +555,14 @@ def count_CO(data_path): for group_idx, (name, group) in enumerate(grouped): print('name', name) s = group.iloc[-1,:] -# print('s', s) - print("s['idx']", s['idx']) h, v = new_data_list[int(s['idx'])] - print('h, ', h) L = to_val(h['L']) density = to_val(h['density']) K = to_val(h['K']) - print('L, density, K', L, density, K) - hist_n, hist_nn, hist, nf, nch, L_array, intensities, orig_foci_intensities, co = summary((h,v), L, K, tp=-1) co_tot += nf @@ -612,9 +572,9 @@ def count_CO(data_path): def main(sim_data_path, output_path): - make_plots((sim_data_path+'/survey_escape/', 'at_'), output_path+'/escape_', intensity_bins=np.linspace(0,0.3,30)) make_plots((sim_data_path+'/survey_julia_new_ends/', 'at_'), output_path+'/new_end_', intensity_bins=np.linspace(0,0.3,30)) - quit() + make_plots((sim_data_path+'/survey_escape/', 'at_'), output_path+'/escape_', intensity_bins=np.linspace(0,0.3,30)) + make_plots((sim_data_path+'/survey_julia_ox/', 'at_'), output_path+'/ox_new_end_', centre_plot=False, intensity_bins=np.linspace(0,0.12,11), max_n=7) make_plots((sim_data_path+'/survey_julia_ux/', 'at_'), output_path+'/ux_new_end_', intensity_bins=np.linspace(0,0.3,10), max_n=3) make_plots((sim_data_path+'/survey_julia_no_ends/', 'at_'), output_path+'/no_end_', intensity_bins=np.linspace(0,0.3,30)) diff --git a/plotting/kymo.py b/plotting/kymo.py index 267962c..c7ad57e 100644 --- a/plotting/kymo.py +++ b/plotting/kymo.py @@ -7,7 +7,7 @@ import matplotlib as mpl -### Script to make matplotlib kymograph plot from +### Script to make matplotlib kymograph plot mpl.rcParams.update({ #'figure.figsize': (6.0,4.0), diff --git a/plotting/kymo_tc.py b/plotting/kymo_tc.py index e216061..1bbef87 100644 --- a/plotting/kymo_tc.py +++ b/plotting/kymo_tc.py @@ -16,6 +16,7 @@ 'axes.labelsize':28, 'savefig.edgecolor': 'none', 'savefig.facecolor': 'none', + 'svg.fonttype' : 'none', }) diff --git a/simulation/run_simulations.py b/simulation/run_simulations.py index 3aa6700..4e723cb 100644 --- a/simulation/run_simulations.py +++ b/simulation/run_simulations.py @@ -70,7 +70,7 @@ julia_path='/home/foz/JuliaPro-1.4.0-1/Julia/bin/julia' # Path to Julia executable julia_options = [] # Possible to precompile Julia packages to speed-up startup -#julia_options = ["--sysimage", "nsim.so" ] # Cached +julia_options = ["--sysimage", "nsim.so" ] # Cached import numpy as np diff --git a/source_data_scripts/copy_files.sh b/source_data_scripts/copy_files.sh new file mode 100644 index 0000000..517ff1e --- /dev/null +++ b/source_data_scripts/copy_files.sh @@ -0,0 +1,23 @@ + + +cp ../output/kymo/kymo1.dat ../source_data/fig2j_kymo1.csv +cp ../output/kymo/kymo2.dat ../source_data/fig2j_kymo2.csv +cp ../output/kymo/kymo3.dat ../source_data/fig2j_kymo3.csv + +cp ../output/kymo/no_end_kymo1.dat ../source_data/figS3h_kymo1.csv +cp ../output/kymo/no_end_kymo2.dat ../source_data/figS3h_kymo2.csv +cp ../output/kymo/no_end_kymo3.dat ../source_data/figS3h_kymo3.csv + +cp ../output/kymo/escape_kymo1.dat ../source_data/figS6h_kymo1.csv +cp ../output/kymo/escape_kymo2.dat ../source_data/figS6h_kymo2.csv +cp ../output/kymo/escape_kymo3.dat ../source_data/figS6h_kymo3.csv + +cp ../output/data_output/rel_length_tot_compare_data.csv ../source_data/fig1b_left.csv +cp ../output/data_output/new_rel_length_data.csv ../source_data/fig1b_left_mid.csv +cp ../output/data_output/new_rel_length2_data.csv ../source_data/fig1b_right_mid.csv +cp ../output/data_output/single_double_triple_peak_new_intensities_stacked_data.csv ../source_data/fig1b_right.csv + +cp ../output/data_output/violin_number_length_data.csv ../source_data/figS5a.csv +cp ../output/julia_plots/new_end_nco_vs_length_data.csv ../source_data/figS5b.csv + + diff --git a/source_data_scripts/data_make_table.py b/source_data_scripts/data_make_table.py index 54e7f62..3aecb11 100644 --- a/source_data_scripts/data_make_table.py +++ b/source_data_scripts/data_make_table.py @@ -152,7 +152,7 @@ def analyse_SC_lengths(df): -def process_data(csv_fn, pkl_fn, category): +def process_data(csv_fn, pkl_fn, output_file, output_trace_file, specific_traces, trace_output_paths): # Read CSV file with data specification A = pd.read_csv(csv_fn) @@ -160,7 +160,7 @@ def process_data(csv_fn, pkl_fn, category): # Remove some unused columns and rename one field A = A.drop(A.columns[[6,7,8,9]], axis=1) - A = A.drop(columns=['foci count', 'foci per chromosome', 'comments']) + A = A.drop(columns=[x for x in ('foci count', 'foci per chromosome', 'comments') if x in A]) A = A.rename(columns={'Plant ':'Plant'}) @@ -183,17 +183,28 @@ def process_data(csv_fn, pkl_fn, category): print(A.iloc[0]) - A.to_csv('test.csv') + A.to_csv(output_file) # Write original hei10 traces to file - with open('traces.txt', 'w') as f: + with open(output_trace_file, 'w') as f: for i, v in data['o_hei10_traces'].items(): f.write(', '.join(map(str, v)) + '\n') + + for i, fn in zip(specific_traces, trace_output_paths): + with open(fn, 'w') as f: + t = data['o_hei10_traces'][i] + for v in t: + f.write(f'{v}' + '\n') + -data_output_path = 'data_output/' -process_data('200406.csv', data_output_path+'test.pkl', 'wt') + +data_output_path = '../output/data_output/' + +process_data('../input_data/200406.csv', data_output_path+'test.pkl', '../source_data/fig2_cytology.csv', '../source_data/fig2_cytology_raw_traces.csv', [675, 664, 492], [f'../source_data/fig1a_HEI10_trace_{s}.csv' for s in ['upper', 'mid', 'lower']]) +process_data('../input_data/OX.csv', data_output_path+'test_ox.pkl', '../source_data/fig3_cytology.csv', '../source_data/fig3_cytology_raw_traces.csv', [260], ['../source_data/fig3a_HEI10_trace.csv']) +process_data('../input_data/UX.csv', data_output_path+'test_ux.pkl', '../source_data/fig4_cytology.csv', '../source_data/fig4_cytology_raw_traces.csv', [115], ['../source_data/fig4a_HEI10_trace.csv']) diff --git a/source_data_scripts/data_table_plots.py b/source_data_scripts/data_table_plots.py index b4216b2..793717a 100644 --- a/source_data_scripts/data_table_plots.py +++ b/source_data_scripts/data_table_plots.py @@ -733,48 +733,23 @@ def all_summary(df, summary_file): print('Mean SC lengths', np.mean(sc_length_all, axis=0), np.std(sc_length_all, axis=0, ddof=1), len(sc_length_all), file=summary_file) -## Plot hei10 (mean) intensity traces along a specific SC -def plot_traces(data_fn, csv_fn, idx_list, image_output_path): - A = pd.read_csv(csv_fn) + - with open(data_fn, 'rb') as f: - data = pickle.load(f) - all_peak_data, orig_trace_lengths = data['all_peak_data'], data['orig_trace_lengths'] - o_hei10_traces = data['o_hei10_traces'] # Don't have dapi data for UX. - for i in idx_list: - plt.figure() - v = np.array(o_hei10_traces[i]) - n = len(v) - plt.plot(np.linspace(0,1,n), v/65535) - plt.savefig(image_output_path+f'trace_{i}.png') - plt.close() - plt.figure() - plt.plot(np.linspace(0,1,n), v/65535) - plt.plot(np.linspace(0,1,n), np.median(v)*np.ones_like(v)/65535) - plt.savefig(image_output_path+f'trace_median_{i}.png') - plt.close() - print(A.iloc[i]) - - with open(image_output_path+f'trace_{i}.txt', 'w') as f: - for u in v: - f.write(f'{u:.02f}\n') - +input_data_path = '../input_data/' +source_data_path = '../source_data/' +test_output_path = '../source_data_test/' - -data_output_path = 'data_output/' -data_output_path2 = 'data_output_test/' -image_output_path='data_output_test/' +os.makedirs(test_output_path, exist_ok=True) for image_output_base, csv_fn, max_n, stacked_bins, summary_fn, summary_fn2 in [ - ( image_output_path, 'test.csv',4, np.linspace(0, 0.3, 30), data_output_path2+'summary.txt', data_output_path2+'cell_summary.txt'), - ]: + ( test_output_path, source_data_path+'fig2_cytology.csv', 4, np.linspace(0, 0.3, 30), test_output_path+'summary.txt', test_output_path+'cell_summary.txt'), + ( test_output_path+'ox_', source_data_path+'fig3_cytology.csv', 4, np.linspace(0, 0.12, 11), test_output_path+'summary_ox.txt', test_output_path+'cell_summary_ox.txt'), + ( test_output_path+'ux_', source_data_path+'fig4_cytology.csv', 3, np.linspace(0, 0.3, 10), test_output_path+'summary_ux.txt', test_output_path+'cell_summary_ux.txt'), + + +]: A = pd.read_csv(csv_fn) plot_data(A, image_output_base, max_n, stacked_bins, summary_fn, summary_fn2) - print(image_output_path) -# Plots of HEI10 traces along individual chromosomes (used for panels in Fig1A, 3A and 4A) -plot_traces(data_output_path+'test.pkl', '200406.csv', [675, 664, 492 ], image_output_path) -plot_traces(data_output_path+'test_ox.pkl', 'OX.csv', [ 260], image_output_path) -plot_traces(data_output_path+'test_ux.pkl', 'UX.csv', [ 115], image_output_path) diff --git a/source_data_scripts/julia_convert_to_source_data.py b/source_data_scripts/julia_convert_to_source_data.py index 784bbfb..37cb356 100644 --- a/source_data_scripts/julia_convert_to_source_data.py +++ b/source_data_scripts/julia_convert_to_source_data.py @@ -106,7 +106,7 @@ def load_data(survey_dir, survey_base): import pandas as pd -def make_data_files(data_path, output_prefix): +def make_data_files(data_path, output_prefix, only_wt_density=False): # Load preprocessed pkl file new_data = load_data(*data_path) @@ -159,6 +159,8 @@ def make_data_files(data_path, output_prefix): L = to_val(h['L']) density = to_val(h['density']) + if density != 0.5 and only_wt_density: + continue K = to_val(h['K']) print('L, density, K', L, density, K) @@ -198,7 +200,8 @@ def make_data_files(data_path, output_prefix): dataset = pd.DataFrame(dataset) - print(np.sum(dataset['Number COs']==0)) + if not only_wt_density: + print('Zero crossover SCs', np.sum(dataset['Number COs']==0)) dataset.to_csv(output_prefix+'density_{}.csv'.format(name), index=False) @@ -213,14 +216,15 @@ def main(sim_data_path, output_path): #make_plots((sim_data_path+'/survey_escape/', 'at_'), output_path+'/escape_', intensity_bins=np.linspace(0,0.3,30)) #return - make_data_files((sim_data_path+'/survey_julia_new_ends/', 'at_'), output_path+'/wt_') - make_data_files((sim_data_path+'/survey_julia_ox/', 'at_'), output_path+'/ox_') - make_data_files((sim_data_path+'/survey_julia_ux/', 'at_'), output_path+'/ux_') - make_data_files((sim_data_path+'/survey_julia_no_ends/', 'at_'), output_path+'/no_end_') - make_data_files((sim_data_path+'/survey_exp/', 'at_'), output_path+'/exp_') + make_data_files((sim_data_path+'/survey_julia_new_ends/', 'at_'), output_path+'/fig2_simulations_') + make_data_files((sim_data_path+'/survey_julia_ox/', 'at_'), output_path+'/fig3_simulations_') + make_data_files((sim_data_path+'/survey_julia_ux/', 'at_'), output_path+'/fig4_simulations_') + make_data_files((sim_data_path+'/survey_julia_no_ends/', 'at_'), output_path+'/figS3_simulations_') - make_data_files((sim_data_path+'/survey_female/', 'at_'), output_path+'/female_') + make_data_files((sim_data_path+'/survey_female/', 'at_'), output_path+'/figS4b_simulations_', only_wt_density=True) + + make_data_files((sim_data_path+'/survey_escape/', 'at_'), output_path+'/figS6_simulations_') -main(sys.argv[1], sys.argv[2]) +main('../output', '../source_data') diff --git a/source_data_scripts/julia_plots_from_table.py b/source_data_scripts/julia_plots_from_table.py index da1c609..6f527ae 100644 --- a/source_data_scripts/julia_plots_from_table.py +++ b/source_data_scripts/julia_plots_from_table.py @@ -141,9 +141,9 @@ def line_hist(ax, x, data, label=None): intensity_4 NaN """ -def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, max_n=4, max_k=4): +def make_plots(input_csv_path, output_prefix, centre_plot=True, intensity_bins=None, max_n=4, max_k=4, do_density_plot=True): # Load preprocessed pkl file - new_data = pd.read_csv('test/wt_density_0.5.csv') + new_data = pd.read_csv(input_csv_path+'_density_0.5.csv') all_hist_n = np.zeros((19,)) @@ -175,7 +175,7 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, s = new_data.iloc[i] print(s) - nco = s['Number COs'] + nco = int(s['Number COs']) L = s['L'] all_nco.append(nco) @@ -278,11 +278,11 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, plt.savefig(output_prefix+'diff.svg') - with open('H2A.csv') as f: + with open('../input_data/H2A.csv') as f: reader = csv.reader(f) H2A = list(reader) - with open('MHL1_diakinesis.csv') as f: + with open('../input_data/MHL1_diakinesis.csv') as f: reader = csv.reader(f) MHL = list(reader) @@ -437,8 +437,6 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, plt.plot(c[s], v[s], 'rx') - with open('centre2.pkl', 'wb') as f: - pickle.dump([c,v,s], f) X, yy, r = lin_fit(c[s], v[s], r2=True) @@ -453,51 +451,61 @@ def make_plots(data_path, output_prefix, centre_plot=True, intensity_bins=None, print(H2A, MHL) - dsb_density = [] - CO_sim = [] - - for d in ['0.25', '0.375', '0.5']: - data = pd.read_csv(f'test/wt_density_{d}.csv') - nf = np.sum(data['Number COs']) - nch = len(data) - dsb_density.append(float(d)) - CO_sim.append(nf/nch) + if do_density_plot: + dsb_density = [] + CO_sim = [] + for d in ['0.25', '0.375', '0.5']: + data = pd.read_csv(input_csv_path+f'_density_{d}.csv') + nf = np.sum(data['Number COs']) + nch = len(data) + dsb_density.append(float(d)) + CO_sim.append(nf/nch) - - dsb_density = np.array(dsb_density) - dsb_density = dsb_density/dsb_density[-1] - - dsb_data = H2A/H2A[0] - CO_data = MHL[:3]/5 - CO_data = CO_data/CO_data[0] + dsb_density = np.array(dsb_density) + dsb_density = dsb_density/dsb_density[-1] - print('nf tot', nf_tot) - - plt.figure() - CO_sim = np.array(CO_sim) - CO_sim = CO_sim/CO_sim[-1] + dsb_data = H2A/H2A[0] + CO_data = MHL[:3]/5 + CO_data = CO_data/CO_data[0] - print(dsb_density, CO_sim) - - plt.plot(dsb_density, CO_sim) - plt.plot(dsb_data, CO_data, 'rx') - plt.xlabel('Relative RI density') - plt.ylabel('Relative number of \nCOs per bivalent') + print('nf tot', nf_tot) + + plt.figure() + CO_sim = np.array(CO_sim) + CO_sim = CO_sim/CO_sim[-1] - plt.savefig(output_prefix+'julia_hs.svg') - + print(dsb_density, CO_sim) + plt.plot(dsb_density, CO_sim) + plt.plot(dsb_data, CO_data, 'rx') + plt.xlabel('Relative RI density') + plt.ylabel('Relative number of \nCOs per bivalent') -def main(sim_data_path, output_path): - make_plots((sim_data_path+'/survey_julia_new_ends/', 'at_'), output_path+'/wt_', intensity_bins=np.linspace(0,0.3,30)) + + plt.savefig(output_prefix+'julia_hs.svg') -main(sys.argv[1], sys.argv[2]) + + +source_data_path = '../source_data/' +test_output_path = '../source_data_test_julia/' + +os.makedirs(test_output_path, exist_ok=True) + +def main(): + make_plots(source_data_path+'fig2_simulations', test_output_path+'/new_end_', intensity_bins=np.linspace(0,0.3,30)) + make_plots(source_data_path+'figS6_simulations', test_output_path+'/escape_', intensity_bins=np.linspace(0,0.3,30)) + make_plots(source_data_path+'fig3_simulations', test_output_path+'/ox_new_end_', intensity_bins=np.linspace(0,0.12,11), max_n=7, centre_plot=False) + make_plots(source_data_path+'fig4_simulations', test_output_path+'/ux_new_end_', intensity_bins=np.linspace(0,0.3,30), max_n=3) + make_plots(source_data_path+'figS3_simulations', test_output_path+'/no_end_', intensity_bins=np.linspace(0,0.3,30)) + make_plots(source_data_path+'figS4b_simulations', test_output_path+'/female_', intensity_bins=np.linspace(0,0.3,30), do_density_plot=False) + +main()