Browse Source

Modify the sfs plot

tforest 1 year ago
parent
commit
20d66bf1cc
2 changed files with 13 additions and 1 deletions
  1. 1 1
      __init__.py
  2. 12 0
      sfs_tools.py

+ 1 - 1
__init__.py View File

1
-from frst import sfs_tools, customgraphics, vcf_utils, sfs_tools, stats_sfs, dependences
1
+from frst import sfs_tools, customgraphics, vcf_utils, sfs_tools, stats_sfs, dependences, swp2
2
 
2
 

+ 12 - 0
sfs_tools.py View File

227
         sum_val = sum(sfs_val)
227
         sum_val = sum(sfs_val)
228
         for k, sfs_bin in enumerate(sfs_val):
228
         for k, sfs_bin in enumerate(sfs_val):
229
             sfs_val[k] = sfs_bin / sum_val
229
             sfs_val[k] = sfs_bin / sum_val
230
+        #print(sum(sfs_val))
230
     #build the plot
231
     #build the plot
231
     title = title+" (n="+str(len(sfs_val)+1)+") [folded="+str(folded)+"]"+" [transformed="+str(transformed)+"]"
232
     title = title+" (n="+str(len(sfs_val)+1)+") [folded="+str(folded)+"]"+" [transformed="+str(transformed)+"]"
232
     print("SFS =", sfs)
233
     print("SFS =", sfs)
234
+    if folded:
235
+        xlab = "Minor allele frequency"
233
     if transformed:
236
     if transformed:
234
         print("Transformed SFS ( n =",len(sfs_val)+1, ") :", sfs_val)
237
         print("Transformed SFS ( n =",len(sfs_val)+1, ") :", sfs_val)
238
+        plt.axhline(y=1/n, color='r', linestyle='-')
239
+    else:
240
+        if normalized:
241
+            # then plot a theoritical distribution as 1/i
242
+            expected_y = [1/(2*x+1) for x in list(sfs.keys())]
243
+            print(sum(expected_y))
244
+            plt.plot([x for x in list(sfs.keys())], expected_y, color='r', linestyle='-')
245
+            #print(expected_y)
246
+            
235
     customgraphics.barplot(x = [x for x in list(sfs.keys())], y= sfs_val, xlab = xlab, ylab = ylab, title = title)
247
     customgraphics.barplot(x = [x for x in list(sfs.keys())], y= sfs_val, xlab = xlab, ylab = ylab, title = title)
236
     plt.show()
248
     plt.show()
237
 
249