Browse Source

correction of transformed sfs plot

tforest 1 year ago
parent
commit
f3fc152df3
1 changed files with 17 additions and 11 deletions
  1. 17 11
      sfs_tools.py

+ 17 - 11
sfs_tools.py View File

192
     return SFS_values, count_pluriall
192
     return SFS_values, count_pluriall
193
 
193
 
194
 
194
 
195
-def barplot_sfs(sfs,  xlab, ylab, folded=True, title = "Barplot", transformed = False):
195
+def barplot_sfs(sfs,  xlab, ylab, folded=True, title = "Barplot", transformed = False, normalized = False):
196
     sfs_val = []
196
     sfs_val = []
197
     n = len(sfs.values())
197
     n = len(sfs.values())
198
-    print("n =", n)
198
+    sum_sites = sum(list(sfs.values()))
199
     for k, ksi in sfs.items():
199
     for k, ksi in sfs.items():
200
         #ksi = list(sfs.values())[k-1]
200
         #ksi = list(sfs.values())[k-1]
201
         # k+1 because k starts from 0
201
         # k+1 because k starts from 0
209
         #         sfs_val.append(ksi)
209
         #         sfs_val.append(ksi)
210
         if transformed:
210
         if transformed:
211
             if folded:
211
             if folded:
212
-                #sfs_val.append(ksi * k * (2*n - k))
213
-                sfs_val.append(((k*(2*n - k)) / (2*n))*ksi) 
212
+                val = ((k*(2*n - k)) / (2*n))*(ksi)
214
             else:
213
             else:
215
-                sfs_val.append(ksi * k)
214
+                val = ksi * k
216
         else:
215
         else:
217
-             sfs_val.append(ksi)
216
+            val = ksi
217
+        sfs_val.append(val)
218
             
218
             
219
     #terminal case, same for folded or unfolded
219
     #terminal case, same for folded or unfolded
220
     if transformed:
220
     if transformed:
221
-        sfs_val[-1] = list(sfs.values())[n-1] * n
221
+        last_bin = list(sfs.values())[n-1] * n/2
222
     else:
222
     else:
223
-        sfs_val[-1] = list(sfs.values())[n-1]
223
+        last_bin = list(sfs.values())[n-1]
224
+    sfs_val[-1] = last_bin
225
+    if normalized:
226
+        ylab = "Fraction of SNPs"
227
+        sum_val = sum(sfs_val)
228
+        for k, sfs_bin in enumerate(sfs_val):
229
+            sfs_val[k] = sfs_bin / sum_val
224
     #build the plot
230
     #build the plot
225
-    title = title+" [folded="+str(folded)+"]"+" [transformed="+str(transformed)+"]"
231
+    title = title+" (n="+str(len(sfs_val)+1)+") [folded="+str(folded)+"]"+" [transformed="+str(transformed)+"]"
226
     print("SFS =", sfs)
232
     print("SFS =", sfs)
227
     if transformed:
233
     if transformed:
228
-        print("Transformed SFS ( n =",len(sfs_val), ") :", sfs_val)
229
-    customgraphics.barplot(x = sfs.keys(), y= sfs_val, xlab = xlab, ylab = ylab, title = title)
234
+        print("Transformed SFS ( n =",len(sfs_val)+1, ") :", sfs_val)
235
+    customgraphics.barplot(x = [x for x in list(sfs.keys())], y= sfs_val, xlab = xlab, ylab = ylab, title = title)
230
     plt.show()
236
     plt.show()
231
 
237
 
232
 if __name__ == "__main__":
238
 if __name__ == "__main__":