richard
/
geinf
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Compare commits

base: richard:bbc06411e9e5c72ab9a761203aabe483819a3d36
Branches Tags
richard:master
..
compare: richard:592b16e7f30906df01c36f4a62afcd485c186bbf
Branches Tags
richard:master

No commits in common. "bbc06411e9e5c72ab9a761203aabe483819a3d36" and "592b16e7f30906df01c36f4a62afcd485c186bbf" have entirely different histories.
bbc06411e9 ... 592b16e7f3

2 changed files with 8 additions and 19 deletions
Show Stats Expand all files Collapse all files

14
sfs_tools.py
View File

@ -208,7 +208,6 @@ def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed =
# else: # else:
# sfs_val.append(ksi) # sfs_val.append(ksi)
if transformed: if transformed:
ylab = r'$ \phi_i $'
if folded: if folded:
val = ((k*(2*n - k)) / (2*n))*(ksi) val = ((k*(2*n - k)) / (2*n))*(ksi)
else: else:
@ -216,9 +215,6 @@ def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed =
else: else:
val = ksi val = ksi
sfs_val.append(val) sfs_val.append(val)
if not transformed and not normalized:
ylab = r'$ \eta_i $'
#terminal case, same for folded or unfolded #terminal case, same for folded or unfolded
if transformed: if transformed:
@ -227,21 +223,19 @@ def barplot_sfs(sfs, xlab, ylab, folded=True, title = "Barplot", transformed =
last_bin = list(sfs.values())[n-1] last_bin = list(sfs.values())[n-1]
sfs_val[-1] = last_bin sfs_val[-1] = last_bin
if normalized: if normalized:
#ylab = "Fraction of SNPs " ylab = "Fraction of SNPs"
ylab = r'$ \phi_i $'
sum_val = sum(sfs_val) sum_val = sum(sfs_val)
for k, sfs_bin in enumerate(sfs_val): for k, sfs_bin in enumerate(sfs_val):
sfs_val[k] = sfs_bin / sum_val sfs_val[k] = sfs_bin / sum_val
#print(sum(sfs_val)) #print(sum(sfs_val))
#build the plot #build the plot
title = title+" (n="+str(len(sfs_val))+") [folded="+str(folded)+"]"+" [transformed="+str(transformed)+"]" title = title+" (n="+str(len(sfs_val)+1)+") [folded="+str(folded)+"]"+" [transformed="+str(transformed)+"]"
print("SFS =", sfs) print("SFS =", sfs)
if folded: if folded:
xlab = "Minor allele frequency" xlab = "Minor allele frequency"
if transformed: if transformed:
print("Transformed SFS ( n =",len(sfs_val), ") :", sfs_val) print("Transformed SFS ( n =",len(sfs_val)+1, ") :", sfs_val)
#plt.axhline(y=1/n, color='r', linestyle='-') plt.axhline(y=1/n, color='r', linestyle='-')
else: else:
if normalized: if normalized:
# then plot a theoritical distribution as 1/i # then plot a theoritical distribution as 1/i

13
swp2.py
View File

@ -27,19 +27,17 @@ def return_x_y_from_stwp_theta_file(stwp_theta_file, breaks, mu, tgen):
#### END of parsing #### END of parsing
# quit this file if the number of dimensions is incorrect # quit this file if the number of dimensions is incorrect
if dim < breaks+1: if dim < breaks+1:
return 0,0,0 return 0,0,0
# get n, the last bin of the last group # get n, the last bin of the last group
#n = int(groups[-1].split(",")[-1]) n = int(groups[-1].split(",")[-1])
# revert the list of groups as the most recent times correspond
# to the closest and last leafs of the coal. tree.
groups = groups[::-1]
# initiate the dict of times # initiate the dict of times
t = {} t = {}
# list of thetas # list of thetas
theta_L = [] theta_L = []
sum_t = 0 sum_t = 0
#groups.sort(reverse=True)
for group_nb, group in enumerate(groups): for group_nb, group in enumerate(groups):
print(group_nb, group, theta_site[group_nb], len(theta_site)) print(group_nb, group)
# store all the thetas one by one, with one theta per group # store all the thetas one by one, with one theta per group
theta_L.append(float(theta_site[group_nb])) theta_L.append(float(theta_site[group_nb]))
# if the group is of size 1 # if the group is of size 1
@ -71,14 +69,11 @@ def return_x_y_from_stwp_theta_file(stwp_theta_file, breaks, mu, tgen):
# build the time x axis # build the time x axis
x = [0] x = [0]
# need to fix: t.values are sorted? and in reverse order because the groups and the times are sorted in opposite order
for time in range(0, len(t.values())-1): for time in range(0, len(t.values())-1):
x.append(list(t.values())[time]) x.append(list(t.values())[time])
x.append(list(t.values())[time]) x.append(list(t.values())[time])
x.append(list(t.values())[len(t.values())-1]) x.append(list(t.values())[len(t.values())-1])
#x = x[::-1]
#y = y[::-1]
#x.sort(reverse=True) #x.sort(reverse=True)
#y.sort(reverse=True) #y.sort(reverse=True)