removing pruning threshold

README

@@ -1,4 +1,4 @@
-uDance version 1.3.0
+uDance version 1.3.1
 
 AUTHORS:
 Metin Balaban
@@ -9,6 +9,8 @@ Unix/Linux
 Anaconda
 
 Changelog:
+1.3.1:
+-Getting rid of pruning thresholds (automated finding)
 
 1.3.0:
 -Pruning Large partitions

config.yaml

@@ -39,8 +39,6 @@ prep_config:
   fraglength: 75
   # large clusters are pruned. ASTRAL's limit is around 9000 for 400+ genes.
   pruneafter: 9000
-  # [0-1]. lower threshold = stricter pruning
-  prune_thr: 0.96
 
 infer_config:
   # [raxml-ng,iqtree,raxml-8]

prune_similar.py

@@ -1,8 +1,7 @@
 import json
 from multiprocessing import cpu_count
 from optparse import OptionParser
-from os.path import join, isfile
-import shutil
+from os.path import join
 
 from uDance.subsample_partition import subsample_partition
 
@@ -15,8 +14,6 @@
     parser.add_option("-T", "--threads", type=int, dest="num_thread", default=0,
                       help="number of cores used in placement. "
                            "0 to use all cores in the running machine", metavar="NUMBER")
-    parser.add_option("-c", "--cutoff", type=float, dest="cutoff_threshold", metavar='NUMBER', default=0.96,
-                      help="threshold number of dissimilar genes.")
     parser.add_option("-S", "--size", type=int, dest="minimum_size", metavar='NUMBER', default=9000,
                       help="partition size requirement for pruning.")
 
@@ -37,12 +34,10 @@
             if numspecies < options.minimum_size:
                 continue
         print(numspecies)
-        res = subsample_partition(partition_output_dir, options.cutoff_threshold)
+        res = subsample_partition(partition_output_dir, options.minimum_size)
         if res:
             dupmapstrs.append(res)
 
-    # if isfile(join(options.output_fp, "rm_map.txt")):
-    #     shutil.copyfile(join(options.output_fp, "rm_map.txt"), join(options.output_fp, "dedupe_map.txt"))
     if len(dupmapstrs) > 0:
         with open(join(options.output_fp, "dedupe_map.txt"), "a") as f:
             for st in dupmapstrs:

uDance/subsample_partition.py

@@ -8,7 +8,7 @@
 import treeswift as ts
 
 
-def subsample_partition(partition_output_dir, cutoff):
+def subsample_partition(partition_output_dir, limit):
     with open(join(partition_output_dir, "species.txt")) as f:
         species = set(map(lambda x: x.strip(), f.readlines()))
     t = ts.read_tree_newick(join(partition_output_dir, "astral_constraint.nwk"))
@@ -65,11 +65,25 @@ def subsample_partition(partition_output_dir, cutoff):
     x = counts_ij / np.minimum(counts_i[..., np.newaxis], counts_i[np.newaxis, ...])
     x.dump(join(partition_output_dir, "adj_mat.pkl"), protocol=4)
     #print(x)
-    y = (x >= cutoff)
-    print("redo %.3f." % (time.time() - start))
-    start = time.time()
-    n, components = connected_components(y)
-    print(n, components)
+
+    clow = 1
+    chigh = 100
+
+    while chigh > clow or chigh == 2:
+        cutoff = (chigh + clow + 1)//2
+        print(clow, cutoff, chigh)
+        y = (x >= (cutoff/100))
+        n, components = connected_components(y)
+        if limit < n:
+            chigh = cutoff-1
+        elif limit > n:
+            clow = cutoff
+        else:
+            clow = cutoff
+            chigh = cutoff # making sure clow=chigh always holds at the exit. we will print that.
+            break
+
+    print("Partition " + partition_output_dir + " is pruned to %d taxa at the automatic cutoff %f." % (n, clow/100))
 
     print("components %.3f." % (time.time() - start))
 

udance.smk

@@ -198,7 +198,6 @@ checkpoint decompose:
         sub=config["prep_config"]["sublength"],
         frag=config["prep_config"]["fraglength"],
         pra=config["prep_config"]["pruneafter"],
-        prt=config["prep_config"]["prune_thr"],
         char=config["chartype"]
 
     resources: cpus=config["resources"]["cores"],
@@ -212,7 +211,7 @@ checkpoint decompose:
             else
                 python run_udance.py decompose -p -s {input.ind} -o {outdir}/udance -t {params.size} -j {input.j} -m {params.method} -T {resources.cpus} -l {params.sub} -f {params.frag}
             fi
-            python prune_similar.py -T {resources.cpus} -o {outdir}/udance -S {params.pra} -c {params.prt}
+            python prune_similar.py -T {resources.cpus} -o {outdir}/udance -S {params.pra}
             if [  -f {outdir}/udance/dedupe_map.txt ]; then 
                 cat {outdir}/udance/dedupe_map.txt > {outdir}/dedupe_map.txt 
             fi