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

sfs plotting

master
tforest 2023-11-10 10:28:31 +01:00
parent 592b16e7f3
commit ecc6fd23fc
1 changed files with 10 additions and 4 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

14
sfs_tools.py
View File

@ -208,6 +208,7 @@ 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:
@ -215,6 +216,9 @@ 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:
@ -223,19 +227,21 @@ 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)+1)+") [folded="+str(folded)+"]"+" [transformed="+str(transformed)+"]" title = title+" (n="+str(len(sfs_val))+") [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)+1, ") :", sfs_val) print("Transformed SFS ( n =",len(sfs_val), ") :", 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