1. Mark Howison
  2. gabi

Commits

Mark Howison  committed 1fc39c7

new option to sample to load an existing assembly as the initial state

  • Participants
  • Parent commits 74b922e
  • Branches chloroplast, devel 4
    1. hiv
    2. issue-6
    3. issue-7
    4. issue7a

Comments (0)

Files changed (2)

File gabi/graph.py

View file
                 break
         self.write_dot('random.dot', active=True)
 
+    def load_assembly(self, edges):
+        # activate each node in the state
+        map(self.activate, edges)
+        # split any junctions
+        for n in self.A.nodes_iter():
+            if self.A.out_degree(n) > 1:
+                for n1 in self.A.successors(n)[:-1]:
+                    self.A.remove_edge(n, n1)
+            if self.A.in_degree(n) > 1:
+                for n1 in self.A.predecessors(n)[:-1]:
+                    self.A.remove_edge(n1, n)
+        self.write_dot('random.dot', active=True)
+

File gabi/sample.py

View file
 
 from time import time
 
-import trace
-from sampler import Sampler
-from graph import Graph
+from gabi import trace
+from gabi.graph import Graph
+from gabi.sampler import Sampler
 
 from biolite import diagnostics
 from biolite.pipeline import IlluminaPipeline
 pipe.add_arg('--seed', type=long, default=-1, help="""
     Deterministic (instead of random) sampling, starting from a seed value.""")
 
+pipe.add_path_arg('--assembly', default=None, help="""
+    Use an assembly as the initial state (provided as a text file of comma
+    separated edge IDs).""")
+
 # Required arguments.
 
 pipe.add_arg('--genome_size', '-s', type=int, metavar='SIZE', help="""
 @pipe.stage
 def sample(
     data, graphml, genome_size, contigs, scale, nodata, flatprior, seed,
-    insert_mean, insert_stdev, samples, perturbations, outdir):
+    assembly, insert_mean, insert_stdev, samples, perturbations, outdir):
     """
     Perform MCMC sampling: iteratively perturb the graph, re-calculate
     likelihood with LAP, and accept new states that pass the odds ratio.
     diagnostics.log_path(trace_path, 'trace')
 
     if last is None:
-        # Random initial state.
-        graph.randomize_state()
+        if assembly:
+            edges = map(int, open(assembly).read().strip().split(','))
+            graph.load_assembly(edges)
+        else:
+            # Random initial state.
+            graph.randomize_state()
         # Start with an initial assembly from the unperturbed graph.
         current = sampler.sample(graph.contigs(), 'initial-contigs.fa')
         save_state = graph.get_state()