version utilisée pour la thèse / PCIA
parent
1848102140
commit
33a120ed7e
|
|
@ -114,7 +114,7 @@ Github: http://github.com/JohnLonginotto/pybam
|
|||
my_bam = pybam.read('/my/data.bam',decompressor='internal')
|
||||
|
||||
[ Parse Words (hah) ]'''
|
||||
wat += '\n'+''.join([('\n===============================================================================================\n\n ' if code is 'file_alignments_read' or code is 'sam' else ' ')+(code+' ').ljust(25,'-')+description+'\n' for code,description in sorted(parse_codes.items())]) + '\n'
|
||||
wat += '\n'+''.join([('\n===============================================================================================\n\n ' if code == 'file_alignments_read' or code == 'sam' else ' ')+(code+' ').ljust(25,'-')+description+'\n' for code,description in sorted(parse_codes.items())]) + '\n'
|
||||
|
||||
class read():
|
||||
'''
|
||||
|
|
@ -155,7 +155,7 @@ class read():
|
|||
|
||||
if fields is not False:
|
||||
print(fields)
|
||||
if type(fields) is not list or len(fields) is 0:
|
||||
if type(fields) is not list or len(fields) == 0:
|
||||
raise PybamError('\n\nFields for the static parser must be provided as a non-empty list. You gave a ' + str(type(fields)) + '\n')
|
||||
else:
|
||||
for field in fields:
|
||||
|
|
@ -167,7 +167,7 @@ class read():
|
|||
|
||||
if decompressor:
|
||||
if type(decompressor) is str:
|
||||
if decompressor is not 'internal' and '{}' not in decompressor: raise PybamError('\n\nWhen a custom decompressor is used and the input file is a string, the decompressor string must contain at least one occurence of "{}" to be substituted with a filepath by pybam.\n')
|
||||
if decompressor != 'internal' and '{}' not in decompressor: raise PybamError('\n\nWhen a custom decompressor is used and the input file is a string, the decompressor string must contain at least one occurence of "{}" to be substituted with a filepath by pybam.\n')
|
||||
else: raise PybamError('\n\nUser-supplied decompressor must be a string that when run on the command line decompresses a named file (or stdin), to stdout:\ne.g. "lzma --decompress --stdout {}" if pybam is provided a path as input file, where {} is substituted for that path.\nor just "lzma --decompress --stdout" if pybam is provided a file object instead of a file path, as data from that file object will be piped via stdin to the decompression program.\n')
|
||||
|
||||
## First we make a generator that will return chunks of uncompressed data, regardless of how we choose to decompress:
|
||||
|
|
@ -203,7 +203,7 @@ class read():
|
|||
|
||||
elif magic == b"\x1f\x8b\x08\x04": # The user has passed us compressed gzip/bgzip data, which is typical for a BAM file
|
||||
# use custom decompressor if provided:
|
||||
if decompressor is not False and decompressor is not 'internal':
|
||||
if decompressor != False and decompressor != 'internal':
|
||||
if type(f) is str: self._subprocess = subprocess.Popen( decompressor.replace('{}',f), shell=True, stdout=subprocess.PIPE, stderr=DEVNULL)
|
||||
else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + decompressor, stdin=self._file, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL)
|
||||
self.file_decompressor = decompressor
|
||||
|
|
@ -230,7 +230,7 @@ class read():
|
|||
use = 'gzip'
|
||||
except OSError:
|
||||
use = 'internal'
|
||||
if use is not 'internal' and decompressor is not 'internal':
|
||||
if use != 'internal' and decompressor != 'internal':
|
||||
if type(f) is str: self._subprocess = subprocess.Popen([ use , '--decompress','--stdout', f ], stdout=subprocess.PIPE, stderr=DEVNULL)
|
||||
else: self._subprocess = subprocess.Popen('{ printf "'+magic+'"; cat; } | ' + use + ' --decompress --stdout', stdin=f, shell=True, stdout=subprocess.PIPE, stderr=DEVNULL)
|
||||
time.sleep(1)
|
||||
|
|
@ -307,7 +307,7 @@ class read():
|
|||
yield b''.join(internal_cache)
|
||||
return
|
||||
|
||||
elif decompressor is not False and decompressor is not 'internal':
|
||||
elif decompressor != False and decompressor != 'internal':
|
||||
# It wouldn't be safe to just print to the shell four random bytes from the beginning of a file, so instead it's
|
||||
# written to a temp file and cat'd. The idea here being that we trust the decompressor string as it was written by
|
||||
# someone with access to python, so it has system access anyway. The file/data, however, should not be trusted.
|
||||
|
|
|
|||
25
sfs_tools.py
25
sfs_tools.py
|
|
@ -241,7 +241,9 @@ def sfs_from_parsed_vcf(n, vcf_dict, folded = True, diploid = True, phased = Fal
|
|||
return SFS_values, count_pluriall
|
||||
|
||||
|
||||
def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed = False, normalized = False, ploidy = 2, output = None):
|
||||
|
||||
def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed = False,
|
||||
normalized = False, ploidy = 2, output = None, step = 8):
|
||||
sfs_val = []
|
||||
n = len(sfs.values())
|
||||
sum_sites = sum(list(sfs.values()))
|
||||
|
|
@ -300,8 +302,8 @@ def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed =
|
|||
|
||||
if transformed:
|
||||
print("Transformed SFS ( n =",n_title, ") :", sfs_val)
|
||||
#plt.axhline(y=1/n, color='r', linestyle='-')
|
||||
plt.bar([x+0.2 for x in list(sfs.keys())], [1/n]*n, fill=False, hatch="///", linestyle='-', width = 0.4, label= "H0 Theoric constant")
|
||||
plt.axhline(y=1/n, color='r', linestyle='-', label= "H0 Theoric constant")
|
||||
#plt.bar([x+0.2 for x in list(sfs.keys())], [1/n]*n, fill=False, hatch="///", linestyle='-', width = 0.4, label= "H0 Theoric constant")
|
||||
|
||||
else:
|
||||
if normalized:
|
||||
|
|
@ -311,12 +313,27 @@ def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed =
|
|||
print(expected_y)
|
||||
plt.bar([x+0.2 for x in list(sfs.keys())], expected_y, fill=False, hatch="///", linestyle='-', width = 0.4, label= "H0 Theoric constant")
|
||||
print(sum(expected_y))
|
||||
|
||||
# graphics parameters
|
||||
if output is not None:
|
||||
# if write in a file, don't open the window dynamically
|
||||
plot = False
|
||||
else:
|
||||
plot = True
|
||||
customgraphics.barplot(x = [x-0.2 for x in X_axis], width=0.4, y= sfs_val, xlab = xlab, ylab = ylab, title = title, label = "H1 Observed spectrum", xticks =list(sfs.keys()), plot = plot )
|
||||
# customgraphics.barplot(x = [x-0.2 for x in X_axis], width=0.4, y= sfs_val, xlab = xlab, ylab = ylab, title = title, label = "H1 Observed spectrum", xticks =list(sfs.keys()), plot = plot )
|
||||
customgraphics.barplot(x = X_axis, width=0.4, y= sfs_val, xlab = xlab, ylab = ylab, title = title, label = "H1 Observed spectrum", xticks =list(sfs.keys()), plot = plot )
|
||||
|
||||
|
||||
|
||||
|
||||
# graphics parameters
|
||||
if n > 12:
|
||||
# plt.locator_params(axis='x', nbins=10)
|
||||
step = n//step
|
||||
x_ticks = [k for k in range(1,n,step)]+[n]
|
||||
print(x_ticks)
|
||||
plt.xticks(x_ticks)
|
||||
|
||||
if output:
|
||||
plt.savefig(f"{output}/{original_title}_SFS.pdf")
|
||||
else:
|
||||
|
|
|
|||
180
swp2.py
180
swp2.py
|
|
@ -78,18 +78,20 @@ def parse_stwp_theta_file(stwp_theta_file, breaks, mu, tgen, relative_theta_scal
|
|||
j = int(group.split(",")[-1])
|
||||
t[i] = 0
|
||||
#t =
|
||||
#T += y[i]*2 / (x[i]*(x[i]-1))
|
||||
|
||||
if len(group.split(',')) == 1:
|
||||
k = i
|
||||
if relative_theta_scale:
|
||||
t[i] += ((theta_L[group_nb] ) / (k*(k-1)))
|
||||
t[i] += ((theta_L[group_nb]*2 ) / (k*(k-1)))
|
||||
else:
|
||||
t[i] += ((theta_L[group_nb] ) / (k*(k-1)) * tgen) / mu
|
||||
t[i] += ((theta_L[group_nb]*2 ) / (k*(k-1)) ) / mu
|
||||
else:
|
||||
for k in range(j, i-1, -1 ):
|
||||
if relative_theta_scale:
|
||||
t[i] += ((theta_L[group_nb] ) / (k*(k-1)))
|
||||
t[i] += ((theta_L[group_nb]*2 ) / (k*(k-1)))
|
||||
else:
|
||||
t[i] += ((theta_L[group_nb] ) / (k*(k-1)) * tgen) / mu
|
||||
t[i] += ((theta_L[group_nb]*2 ) / (k*(k-1)) ) / mu
|
||||
# we add the cumulative times at the end
|
||||
t[i] += sum_t
|
||||
sum_t = t[i]
|
||||
|
|
@ -126,7 +128,7 @@ def plot_straight_x_y(x,y):
|
|||
|
||||
def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title",
|
||||
theta_scale = True, ax = None, input = None, output = None):
|
||||
my_dpi = 500
|
||||
my_dpi = 800
|
||||
L = full_dict["L"]
|
||||
if ax is None:
|
||||
# intialize figure
|
||||
|
|
@ -143,14 +145,14 @@ def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title",
|
|||
plot_handles = []
|
||||
best_plot = full_dict['all_epochs']['best']
|
||||
p0, = ax1.plot(best_plot[0], best_plot[1], linestyle = "-",
|
||||
alpha=1, lw=2, label = str(best_plot[2])+' brks | Lik='+best_plot[3])
|
||||
alpha=1, lw=2, label = str(best_plot[2])+' brks | Lik='+str(round(float(best_plot[3]), 4)))
|
||||
plot_handles.append(p0)
|
||||
#ax1.grid(True,which="both", linestyle='--', alpha = 0.3)
|
||||
for k, plot_Lk in enumerate(full_dict['all_epochs']['plots']):
|
||||
plot_Lk = str(full_dict['all_epochs']['plots'][k][3])
|
||||
# plt.rcParams['font.size'] = fnt_size
|
||||
p, = ax1.plot(full_dict['all_epochs']['plots'][k][0], full_dict['all_epochs']['plots'][k][1], linestyle = "-",
|
||||
alpha=1/(k+1), lw=1.5, label = str(full_dict['all_epochs']['plots'][k][2])+' brks | Lik='+plot_Lk)
|
||||
alpha=1/(k+1), lw=1.5, label = str(full_dict['all_epochs']['plots'][k][2])+' brks | Lik='+str(round(float(plot_Lk), 4)))
|
||||
plot_handles.append(p)
|
||||
if theta_scale:
|
||||
ax1.set_xlabel("Coal. time", fontsize=fnt_size)
|
||||
|
|
@ -172,11 +174,12 @@ def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title",
|
|||
# ax1.set_xticklabels([f'{k}\n{k/(mu)}\n{k/(mu)*tgen}' for k in x_ticks], fontsize = fnt_size*0.8)
|
||||
# plt.rcParams['font.size'] = fnt_size
|
||||
# print(fnt_size, "rcParam font.size=", plt.rcParams['font.size'])
|
||||
ax1.legend(handles = plot_handles, loc='best', fontsize = fnt_size*0.5)
|
||||
ax1.legend(handles = plot_handles, fontsize = fnt_size*0.5, bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0.)
|
||||
ax1.set_title(title)
|
||||
breaks = len(full_dict['all_epochs']['plots'])
|
||||
if ax is None:
|
||||
plt.savefig(title+'_best_'+str(breaks+1)+'_epochs.pdf')
|
||||
fig.tight_layout()
|
||||
plt.savefig(title+'_best_'+str(breaks+1)+'_epochs.pdf', bbox_inches='tight')
|
||||
# plot likelihood against nb of breakpoints
|
||||
if ax is None:
|
||||
fig, ax2 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi)
|
||||
|
|
@ -192,12 +195,13 @@ def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title",
|
|||
ax2.scatter(full_dict['Ln_Brks'][0], full_dict['Ln_Brks'][1], s=50, c=colors, marker='o', zorder=2)
|
||||
ax2.axhline(y=full_dict['best_Ln'], linestyle = "-.", color = "red", label = "$-\log\mathcal{L}$ = "+str(round(full_dict['best_Ln'], 2)))
|
||||
ax2.set_yscale('log')
|
||||
ax2.set_xlabel("# breakpoints", fontsize=fnt_size)
|
||||
ax2.set_xlabel("# breakpoints", fontsize=fnt_size * 0.8)
|
||||
ax2.set_ylabel("$-\log\mathcal{L}$", fontsize=fnt_size)
|
||||
ax2.legend(loc='best', fontsize = fnt_size*0.5)
|
||||
ax2.set_title(title+" Likelihood gain from # breakpoints")
|
||||
if ax is None:
|
||||
plt.savefig(title+'_Breakpts_Likelihood.pdf')
|
||||
fig.tight_layout()
|
||||
plt.savefig(title+'_Breakpts_Likelihood.pdf', bbox_inches='tight')
|
||||
# AIC
|
||||
if ax is None:
|
||||
fig, ax3 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi)
|
||||
|
|
@ -212,11 +216,12 @@ def plot_all_epochs_thetafolder(full_dict, mu, tgen, title = "Title",
|
|||
ax3.axhline(y=full_dict['best_AIC'], linestyle = "-.", color = "red",
|
||||
label = "Min. AIC = "+str(round(full_dict['best_AIC'], 2)))
|
||||
ax3.set_yscale('log')
|
||||
ax3.set_xlabel("# breakpoints", fontsize=fnt_size)
|
||||
ax3.set_ylabel("AIC")
|
||||
ax3.set_xlabel("# breakpoints", fontsize=fnt_size*0.8)
|
||||
ax3.set_ylabel("AIC", fontsize=fnt_size*0.8)
|
||||
ax3.legend(loc='best', fontsize = fnt_size*0.5)
|
||||
ax3.set_title(title+" AIC")
|
||||
if ax is None:
|
||||
fig.tight_layout()
|
||||
plt.savefig(title+'_Breakpts_Likelihood_AIC.pdf')
|
||||
else:
|
||||
# return plots
|
||||
|
|
@ -307,7 +312,7 @@ def save_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc
|
|||
for i in range(2, my_n):
|
||||
an +=1.0/i
|
||||
|
||||
print("an=", an, "theta_w", S/an, "theta_w_p_site", (S/an)/L)
|
||||
#print("an=", an, "theta_w", S/an, "theta_w_p_site", (S/an)/L)
|
||||
# compute Ln
|
||||
Ln = log_facto(S+S0) - log_facto(S0) + np.log(float(S0)/(S+S0)) * S0
|
||||
for xi in range(0, len(SFS_stored)):
|
||||
|
|
@ -319,10 +324,13 @@ def save_all_epochs_thetafolder(folder_path, mu, tgen, title = "Title", theta_sc
|
|||
AIC = []
|
||||
for brk in np.array(brkpt_lik)[:, 0]:
|
||||
brk = int(brk)
|
||||
AIC.append((2*brk+1)+2*np.array(brkpt_lik)[brk, 1].astype(float))
|
||||
AIC.append((2*(brk+1))+(2*np.array(brkpt_lik)[brk, 1].astype(float)))
|
||||
AIC_Brks = [list(np.array(brkpt_lik)[:, 0]), AIC]
|
||||
# AIC = 2*k - 2ln(L) ; where k is the number of parameters, here brks+1
|
||||
AIC_ln = 2*(len(brkpt_lik)+1) - 2*Ln
|
||||
# AIC_new = []
|
||||
# for brk in range(len(brkpt_lik)):
|
||||
# AIC_new.append(2*brk+1 - 2*Ln)
|
||||
best_AIC = AIC_ln
|
||||
selected_brks_nb = AIC.index(min(AIC))
|
||||
# to return : plots ; Ln_Brks ; AIC_Brks ; best_Ln ; best_AIC
|
||||
|
|
@ -463,7 +471,7 @@ def save_k_theta(folder_path, mu, tgen, title = "Title", theta_scale = True,
|
|||
x[i] = int(x[i])
|
||||
# compute the times as: theta_k / (k*(k-1))
|
||||
for i in range(0, len(x)):
|
||||
T += y[i]*2 / (x[i]*(x[i]-1))
|
||||
T += y[i]*4 / (x[i]*(x[i]-1))
|
||||
x_2.append(T)
|
||||
# Save plotting (fig 2)
|
||||
# x_2 = [0]+x_2
|
||||
|
|
@ -501,8 +509,10 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
nb_epochs = len(plot_lines)
|
||||
# fig 2 & 3
|
||||
if ax is None:
|
||||
my_dpi = 500
|
||||
my_dpi = 800
|
||||
fnt_size = 18
|
||||
# make sure that everything is closed before...
|
||||
plt.close()
|
||||
fig2, ax2 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi)
|
||||
fig3, ax3 = plt.subplots(figsize=(5000/my_dpi, 2800/my_dpi), dpi=my_dpi)
|
||||
else:
|
||||
|
|
@ -519,12 +529,46 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
swp2_lines[0][k] = swp2_lines[0][k]/tgen
|
||||
for k in range(len(swp2_lines[1])):
|
||||
swp2_lines[1][k] = swp2_lines[1][k]
|
||||
# x2_plot, y2_plot = plot_straight_x_y(swp2_lines[0],swp2_lines[1])
|
||||
x2_plot, y2_plot = swp2_lines[0], swp2_lines[1]
|
||||
p2, = ax2.plot(x2_plot, y2_plot, linestyle="-", alpha=0.75, lw=2, label = 'swp2', color="black")
|
||||
#x2_plot, y2_plot = swp2_lines[0], swp2_lines[1]
|
||||
x2_plot, y2_plot = plot_straight_x_y(swp2_lines[0],swp2_lines[1])
|
||||
|
||||
title_swp2 = "swp2"
|
||||
|
||||
swp2_x = swp2_lines[0]
|
||||
swp2_y = swp2_lines[1]
|
||||
# SIM ONLY
|
||||
# sim1
|
||||
# title_swp2 = "sim1"
|
||||
# x2_plot = [0, 100, 100, 4600, 4600, 10000]
|
||||
# y2_plot = [5000, 5000, 25000, 25000, 125000, 125000]
|
||||
|
||||
# sim2
|
||||
# title_swp2 = "sim2"
|
||||
# x2_plot = [0, 100, 100, 280, 280, 10000]
|
||||
# y2_plot = [5000, 5000, 1000, 1000, 200, 200]
|
||||
|
||||
# # sim3
|
||||
# title_swp2 = "sim3"
|
||||
# x2_plot = [0, 100, 100, 280, 280, 10000]
|
||||
# y2_plot = [5000, 5000, 1000, 1000, 5000, 5000]
|
||||
|
||||
|
||||
# # sim4
|
||||
# title_swp2 = "sim4"
|
||||
# x2_plot = [0, 100, 100, 4600, 4600, 10000]
|
||||
# y2_plot = [5000, 5000, 25000, 25000, 5000, 5000]
|
||||
|
||||
|
||||
|
||||
# # to comment if not using sim!!
|
||||
# swp2_lines[0] = x2_plot
|
||||
# swp2_lines[1] = y2_plot
|
||||
## END OF SIM ONLY
|
||||
|
||||
p2, = ax2.plot(x2_plot, y2_plot, linestyle="-", alpha=0.75, lw=2, label = title_swp2, color="black")
|
||||
lines_fig2.append(p2)
|
||||
# Plotting (fig 3) which is the same but log scale for x
|
||||
p3, = ax3.plot(x2_plot, y2_plot, linestyle="-", alpha=0.75, lw=2, label = 'swp2', color="black")
|
||||
p3, = ax3.plot(x2_plot, y2_plot, linestyle="-", alpha=0.75, lw=2, label = title_swp2, color="black")
|
||||
lines_fig3.append(p3)
|
||||
min_x = 1
|
||||
min_y = 1
|
||||
|
|
@ -535,6 +579,14 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
x2_plot, y2_plot = plot_straight_x_y(x,y)
|
||||
if subset is not None:
|
||||
if breaks in subset:
|
||||
with open(title+"_theo.csv", "w") as filin:
|
||||
filin.write("xtheo,ytheo\n")
|
||||
for i in range(len(swp2_x)):
|
||||
filin.write(str(swp2_x[i]) + "," + str(swp2_y[i]) + "\n")
|
||||
with open(title+"_obs.csv", "w") as filin:
|
||||
filin.write("xobs,yobs\n")
|
||||
for i in range(len(x)):
|
||||
filin.write(str(x2_plot[i]) + "," + str(y2_plot[i]) + "\n")
|
||||
masking_alpha = 0.75
|
||||
autoscale = True
|
||||
min_x = min(min_x, min(x2_plot))
|
||||
|
|
@ -554,20 +606,31 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
Ne_max_below_limit = y[min(x.index(t_max_below_limit)+1, len(y)-1)]
|
||||
Ne_min_below_limit = y[x.index(t_min_below_limit)]
|
||||
if recent_change:
|
||||
print(f"\n{breaks} breaks ; This is below the recent limit of {recent_limit_years} years:\n",
|
||||
f"t_min (most recent time point under the limit) : {t_min_below_limit/mu*tgen:.1f} t_max (most ancient time point under the limit) : {t_max_below_limit/mu*tgen:.1f}",
|
||||
f"\nNe_min (effective size at t_min) : {Ne_min_below_limit/(4*mu):.1f} Ne_max (effective size at t_max) : {Ne_max_below_limit/(4*mu):.1f}",
|
||||
f"\nNe_min/Ne_max = {(Ne_min_below_limit/(4*mu)) / (Ne_max_below_limit/(4*mu)):.1f}",
|
||||
f"\nEvolution: {((Ne_min_below_limit/(4*mu)) - (Ne_max_below_limit/(4*mu)))/((Ne_max_below_limit/(4*mu)))*100:.1f}%")
|
||||
else:
|
||||
print(f"Recent event under {recent_limit_years} years: NA")
|
||||
# need to compute the last change and when it occured
|
||||
tmin = x[1]
|
||||
tmin_plus_1 = x[2]
|
||||
Ne_min = y[1]
|
||||
Ne_min_plus_1 = y[2]
|
||||
print(f"Last was {tmin/mu*tgen:.1f} years ago. And was of {((Ne_min/(4*mu)) - (Ne_min_plus_1/(4*mu)))/(Ne_min_plus_1/(4*mu))*100:.1f}%")
|
||||
# print(t_max_below_limit, y[min(x.index(t_max_below_limit)+1, len(y)-1)])
|
||||
recent_events_x = x[0:min(x.index(t_max_below_limit)+1, len(y)-1)+1]
|
||||
recent_events_y = y[0:min(x.index(t_max_below_limit)+1, len(y)-1)+1]
|
||||
size_changes = []
|
||||
for i in range(len(recent_events_y)-1):
|
||||
size = recent_events_y[i]
|
||||
# size of t+1
|
||||
size_1 = recent_events_y[i+1]
|
||||
size_changes.append(((size) - (size_1))/((size_1))*100)
|
||||
|
||||
print(f"\n{breaks} breaks ; This is below the recent limit of {recent_limit_years} years:\n",
|
||||
f"t_min (most recent time point under the limit) : {t_min_below_limit*tgen:.1f} t_max (most ancient time point under the limit) : {t_max_below_limit*tgen:.1f}",
|
||||
f"\nNe_min (effective size at t_min) : {Ne_min_below_limit:.1f} Ne_max (effective size at t_max) : {Ne_max_below_limit:.1f}",
|
||||
f"\nNe_min/Ne_max = {(Ne_min_below_limit) / (Ne_max_below_limit):.1f}",
|
||||
f"\nEvolution: {((Ne_min_below_limit) - (Ne_max_below_limit))/((Ne_max_below_limit))*100:.1f}%",
|
||||
f"\nDetail: {len(size_changes)} events since the last {recent_limit_years} years (each event size change in %) : {[round(k,2) for k in size_changes]}")
|
||||
|
||||
else:
|
||||
print(f"{breaks} breaks ; Last recent event under {recent_limit_years} years: NA")
|
||||
# get the last change and when it occured in all cases (under the recent limit or not)
|
||||
tmin = x[1]
|
||||
tmin_plus_1 = x[2]
|
||||
Ne_min = y[1]
|
||||
Ne_min_plus_1 = y[2]
|
||||
print(f"{breaks} breaks ; Last event was {tmin*tgen:.1f} years ago. And size change was of {((Ne_min) - (Ne_min_plus_1))/(Ne_min_plus_1)*100:.1f}%")
|
||||
else:
|
||||
masking_alpha = 0
|
||||
autoscale = False
|
||||
|
|
@ -594,12 +657,15 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
if ax is None:
|
||||
# if not ax, then use the plt syntax, not ax...
|
||||
plt.xlabel(xlabel, fontsize=fnt_size)
|
||||
plt.ylabel(ylabel, fontsize=fnt_size)
|
||||
plt.ylabel(ylabel, fontsize=fnt_size*0.8)
|
||||
plt.gca().set_xlim(0, recent_limit * 3)
|
||||
if recent_change:
|
||||
plt.ylim(Ne_min_below_limit/3, Ne_max_below_limit *3)
|
||||
else:
|
||||
plt.ylim(y2_plot[0]/3, y2_plot[0])
|
||||
plt.ylim(y[0]/3, y[0]+y[0]*0.5)
|
||||
plt.gca().set_xlim(0, x[0] * 3)
|
||||
|
||||
|
||||
# plt.ylim(0, max(max_y+(max_y*0.05), max(swp2_lines[1])+(max(swp2_lines[1])*0.05)))
|
||||
#plt.xlim(0, recent_limit * 3)
|
||||
#xlim_val = plt.gca().get_xlim()
|
||||
|
|
@ -607,7 +673,7 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
# plt.xlim(min(min_x,min(swp2_lines[0])), max(max(swp2_lines[0]), max_x))
|
||||
# x_ticks = list(plt.gca().get_xticks())
|
||||
# plt.gca().set_xticks(x_ticks)
|
||||
# plt.xticks(x_ticks)
|
||||
plt.xticks(x_ticks)
|
||||
# plt.gca().set_xlim(xlim_val)
|
||||
# plt.gca().set_xticklabels([f'{k:.0e}\n{k/(mu):.0e}\n{k/(mu)*tgen:.0e}' for k in x_ticks], fontsize = fnt_size*0.5)
|
||||
plt.gca().set_xticklabels([f'{k:.1f}\n{k*tgen:.1f}' for k in x_ticks], fontsize = fnt_size*0.5)
|
||||
|
|
@ -623,20 +689,21 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
# plt.title(title, fontsize=fnt_size)
|
||||
# plt.legend(handles=lines_fig2, loc='best', fontsize = fnt_size*0.5)
|
||||
# # plt.text(-0.13, -0.135, 'Coal. time\nGen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes)
|
||||
plt.text(-0.13, -0.135, 'Gen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes)
|
||||
plt.text(-0.15, -0.155, 'Gen.\nYears', ha='left', va='bottom', transform=ax3.transAxes)
|
||||
|
||||
plt.subplots_adjust(bottom=0.2) # Adjust the value as needed
|
||||
plt.tight_layout()
|
||||
plt.savefig(title+'_plotB_'+str(nb_epochs)+'_epochs.pdf')
|
||||
# close fig2 to save memory
|
||||
plt.close(fig2)
|
||||
else:
|
||||
# when ax subplotting is used
|
||||
ax2.set_xlabel(xlabel, fontsize=fnt_size)
|
||||
ax2.set_ylabel(ylabel, fontsize=fnt_size)
|
||||
ax2.set_ylabel(ylabel, fontsize=fnt_size*0.8)
|
||||
ax2.set_title(title, fontsize=fnt_size)
|
||||
ax2.legend(handles=lines_fig2, loc='best', fontsize = fnt_size*0.5)
|
||||
ax3.set_xlabel(xlabel, fontsize=fnt_size)
|
||||
ax3.set_ylabel(ylabel, fontsize=fnt_size)
|
||||
ax3.set_ylabel(ylabel, fontsize=fnt_size*0.8)
|
||||
ax3.set_title(title, fontsize=fnt_size)
|
||||
ax3.legend(handles=lines_fig3, loc='best', fontsize = fnt_size*0.5)
|
||||
ax3.set_xscale('log')
|
||||
|
|
@ -644,10 +711,13 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
# Scale the x-axis
|
||||
# x_ticks = list(ax3.get_xticks())
|
||||
# ax3.set_xticks(x_ticks)
|
||||
# x_ticks = [i for i in range(0.1,max(max_x, max(swp2_lines[0]))), ]
|
||||
# x_ticks = [i for i in range(0.1,max(max_x, max(swp2_lines[0]))) ]
|
||||
|
||||
# ax3.set_xticks(x_ticks)
|
||||
ax3.set_xlim(0.1, max(max_x, max(swp2_lines[0])))
|
||||
ax3.set_xlim(1, max(max_x, max(swp2_lines[0])))
|
||||
x_ticks = ax3.get_xticks()
|
||||
ax3.set_xticks(x_ticks)
|
||||
ax3.set_xlim(1, max(max_x, max(swp2_lines[0])))
|
||||
# ax3.set_xlim(min(min(x_ticks), min(swp2_lines[0])), max(max_x, max(swp2_lines[0])))
|
||||
# ax3.set_xlim(1, max(max_x, max(swp2_lines[0])))
|
||||
# ax3.set_xticklabels([f'{k:.0e}\n{k/(mu):.0e}\n{k/(mu)*tgen:.0e}' for k in x_ticks], fontsize = fnt_size*0.5)
|
||||
|
|
@ -664,11 +734,14 @@ def plot_scaled_theta(plot_lines, prop, title, mu, tgen, swp2_lines = None, ax =
|
|||
# ax3.set_yticklabels([f'{k/(4*mu):.0e}' for k in y_ticks], fontsize = fnt_size*0.5)
|
||||
# plt.text(-0.13, -0.135, 'Coal. time\nGen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes)
|
||||
# plt.text(-0.13, -0.135, 'Gen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes)
|
||||
plt.text(-0.13, -0.085, 'Gen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes)
|
||||
# plt.text(-0.13, -0.085, 'Gen. time\nYears', ha='left', va='bottom', transform=ax3.transAxes)
|
||||
|
||||
plt.subplots_adjust(bottom=0.2) # Adjust the value as needed
|
||||
plt.subplots_adjust(top=0.9) # Adjust the value as needed
|
||||
|
||||
if ax is None:
|
||||
# nb of plot_lines represent the number of epochs stored (len(plot_lines) = #breaks+1)
|
||||
# plt.tight_layout()
|
||||
plt.savefig(title+'_plotC_'+str(nb_epochs)+'_epochs_log.pdf')
|
||||
# close fig3 to save memory
|
||||
plt.close(fig3)
|
||||
|
|
@ -700,7 +773,9 @@ def plot_raw_stairs(plot_lines, prop, title, ax = None, n_ticks = 10, rescale =
|
|||
x,y = plot
|
||||
x_plot, y_plot = plot_straight_x_y(x,y)
|
||||
p, = ax1.plot(x_plot, y_plot, 'o', linestyle="-", alpha=0.75, lw=2, label = str(breaks)+' brks')
|
||||
print("breaks=", breaks, "theta0", y[0])
|
||||
print("breaks=", breaks)
|
||||
# for xi in range(len(x)):
|
||||
# print("times", xi, "theta", y[xi])
|
||||
# add plot to the list of all plots to superimpose
|
||||
plots.append(p)
|
||||
x_ticks = x
|
||||
|
|
@ -767,18 +842,31 @@ def combined_plot(folder_path, mu, tgen, breaks, title = "Title", theta_scale =
|
|||
|
||||
# Start of Parsing real swp2 output
|
||||
folder_splitted = folder_path.split("/")
|
||||
swp2_blueprint = "/".join(folder_splitted[:-2])+".blueprint"
|
||||
# parsing blueprint and get some parameters from the run:
|
||||
with open(swp2_blueprint) as filin:
|
||||
for line in filin:
|
||||
if line.startswith("year_per_generation"):
|
||||
blueprint_tgen = float(line.strip().split(":")[-1].strip())
|
||||
|
||||
#print("BLUEPRINT",swp2_blueprint, "tgen", blueprint_tgen )
|
||||
swp2_summary = "/".join(folder_splitted[:-2])+'/'+folder_splitted[-3]+".final.summary"
|
||||
swp2_vals = parse_stairwayplot_output_summary(stwplt_out = swp2_summary)
|
||||
swp2_x, swp2_y = swp2_vals[0], swp2_vals[1]
|
||||
if blueprint_tgen != tgen:
|
||||
# need to rescale swp2 output to the new tgen
|
||||
for i in range(len(swp2_x)):
|
||||
swp2_x[i] = swp2_x[i] / blueprint_tgen * tgen
|
||||
remove_back_and_forth_points(swp2_x, swp2_y)
|
||||
# End of Parsing real swp2 output
|
||||
plot_raw_stairs(plot_lines = loaded_data['raw_stairs'],
|
||||
prop = loaded_data['prop'], title = title, ax = None, max_breaks = breaks)
|
||||
plot_scaled_theta(plot_lines = loaded_data['scaled_stairs'], mu = mu, tgen = tgen, subset=[loaded_data['best_epoch_by_AIC']]+selected_breaks,
|
||||
#plot_scaled_theta(plot_lines = loaded_data['scaled_stairs'], mu = mu, tgen = tgen, subset=[loaded_data['best_epoch_by_AIC']]+selected_breaks,
|
||||
plot_scaled_theta(plot_lines = loaded_data['scaled_stairs'], mu = mu, tgen = tgen, subset=selected_breaks,
|
||||
# plot_scaled_theta(plot_lines = loaded_data['scaled_stairs'], subset=list(range(0,3))+[loaded_data['best_epoch_by_AIC']]+selected_breaks,
|
||||
prop = loaded_data['prop'], title = title, swp2_lines = [swp2_x, swp2_y], ax = None)
|
||||
plot_all_epochs_thetafolder(loaded_data, mu, tgen, title, theta_scale, ax = None)
|
||||
|
||||
plt.close()
|
||||
# plt.close(fig1)
|
||||
# plt.close(fig2)
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue