alignment with seed reduce

Former-commit-id: ef438849c031ae4282b9286c1d9eea8d3f682563

Author: bbengfort

apps/brisera_align.py

@@ -14,12 +14,21 @@
 ##########################################################################
 
 import sys
+import csv
 import brisera
+import StringIO
 
 from brisera.utils import timeit
 from brisera.config import settings
 from pyspark import SparkConf, SparkContext
 
+def write_records(records):
+    output = StringIO.StringIO()
+    writer = csv.writer(output)
+    for record in records:
+        writer.writerow(record)
+    return [output.getvalue()]
+
 @timeit
 def run_brisera_alignment(sc, refpath, qrypath, outpath):
     """
@@ -35,7 +44,7 @@ def run_brisera_alignment(sc, refpath, qrypath, outpath):
         fdelta = 0
 
     # Write alignments to disk
-    alignments.saveAsTextFile(outpath)
+    alignments.mapPartitions(write_records).saveAsTextFile(outpath)
 
     return alignments, adelta, fdelta
 
@@ -64,4 +73,4 @@ def run_brisera_alignment(sc, refpath, qrypath, outpath):
     alignments, adelta, fdelta = result
     print "Total execution time: %0.3f seconds" % delta
     print "    Alignment time: %0.3f seconds" % adelta
-    print "    Filtering time: %0.3f seconds" % fdelta
+    # print "    Filtering time: %0.3f seconds" % fdelta

brisera/align.py

@@ -6,6 +6,7 @@
 ## Imports
 ##########################################################################
 
+from brisera.lv import *
 from brisera.utils import *
 from brisera.records import *
 from brisera.config import settings
@@ -91,9 +92,9 @@ def map_reference(self, arg):
             if self.redundancy > 1 and repseed(sequence, start, self.seed_len):
                 for rdx in xrange(self.redundancy):
                     r   = rdx % self.redundancy
-                    yield (seed, (key, offset, is_last, leftstart, leftlen, rightstart, rightlen, r))
+                    yield (seed, (key, offset, is_last, sequence[leftstart:leftlen], sequence[rightstart:rightlen], r))
             else:
-                yield (seed, (key, offset, is_last, leftstart, leftlen, rightstart, rightlen, 0))
+                yield (seed, (key, offset, is_last, sequence[leftstart:leftlen], sequence[rightstart:rightlen], 0))
 
     def map_queries(self, arg):
         """
@@ -135,9 +136,43 @@ def map_queries(self, arg):
 
                 if self.redundancy > 1 and repseed(sequence, idx, self.seed_len):
                     r = key % self.redundancy
-                    yield (seed, (key, idx, is_last, 0, idx, rightstart, rightlen, r))
+                    yield (seed, (key, idx, is_last, sequence[0:idx], sequence[rightstart:rightlen], r))
                 else:
-                    yield (seed, (key, idx, is_last, 0, idx, rightstart, rightlen, 0))
+                    yield (seed, (key, idx, is_last, sequence[0:idx], sequence[rightstart:rightlen], 0))
+
+    def extend(self, arg):
+        """
+        Performs the extend portion of the BLAST algorithm.
+        """
+        key, (query, reference) = arg
+        refstart    = reference[1]
+        refend      = refstart + self.seed_len
+        differences = 0
+
+        try:
+            # Align left flanks
+            if len(query[3]) != 0:
+                a = LandauVishkin().kdifference(reference[3], query[3], self.k)
+                if a[0] == -1:
+                    return NO_ALIGNMENT
+                if not is_bazea_yates_seed(a, len(query[3]), self.seed_len):
+                    return NO_ALIGNMENT
+
+                refstart -= a[0]
+                differences = a[1]
+
+            # Align right flanks
+            if len(query[4]) != 0:
+                b = LandauVishkin.kdifference(reference[4], query[4], self.k-differences)
+                if a[0] == -1:
+                    return NO_ALIGNMENT
+
+                refend += b[0]
+                differences += b[1]
+
+            return reference[0], refstart, refend, differences
+        except:
+            return NO_ALIGNMENT
 
 ##########################################################################
 ## Spark Functionality
@@ -148,13 +183,21 @@ def align_all(sc, refpath, qrypath):
     """
     Returns an RDD of alignments (no writes to disk)
     """
-    reference = sc.sequenceFile(refpath)
-    queries   = sc.sequenceFile(qrypath)
-    alignment = MerAlignment()
+    reference  = sc.sequenceFile(refpath)
+    queries    = sc.sequenceFile(qrypath)
+    alignment  = MerAlignment()
 
     # Perform mapping
-    reference = reference.flatMap(alignment.map_reference)
-    return reference
+    reference  = reference.flatMap(alignment.map_reference)
+    queries    = queries.flatMap(alignment.map_queries)
+
+    # Find shared seeds
+    shared     = queries.join(reference)
+
+    # Compute alignments with Landau-Viskin
+    alignments = shared.map(alignment.extend).filter(lambda a: a != NO_ALIGNMENT )
+
+    return alignments
 
 if __name__ == '__main__':
     from brisera.convert import *

brisera/config.py

@@ -44,7 +44,7 @@ class BriseraConfiguration(confire.Configuration):
     max_read_len = 36
     k            = 3
     allow_diff   = False
-    filter_align = True
+    filter_align = False
     block_size   = 128
     redundancy   = 16
 

brisera/filter.py

@@ -1,5 +1,7 @@
 """
 Filters the best alignments from the computed alignments.
+
+Was unable to implement due to time constraints.
 """
 
 ##########################################################################

conf/brisera-example.yaml

@@ -14,9 +14,9 @@ min_read_len: 36 # the minimum length of the reads
 max_read_len: 36 # the maximum length of the reads
 
 
-k: 3               # number of mismatches to allow (bigger means more time)
-allow_diff: false  # False is mismatches only, True indels as well
-fliter_align: true # False uses all alignments, True only report unambiguous best alignment
+k: 3                # number of mismatches to allow (bigger means more time)
+allow_diff: false   # False is mismatches only, True indels as well
+fliter_align: false # False uses all alignments, True only report unambiguous best alignment
 
 
 block_size: 128    # number of qry and ref tuples to consider at a time in the reduce phase