started alignment

Former-commit-id: 08032a3da12d833fbbb72d701a942f5655322b7c

Author: bbengfort

Makefile

@@ -26,3 +26,8 @@ test:
 # Targets for installation
 install:
 	$(PYTHON_BIN)/python setup.py install
+
+# Targets for running fixture
+runfixture:
+	rm -rf fixtures/output
+	-spark-submit $(CURDIR)/apps/brisera_align.py $(CURDIR)/fixtures/s_suis.ser $(CURDIR)/fixtures/100k.ser $(CURDIR)/fixtures/output

apps/brisera_align.py

@@ -0,0 +1,67 @@
+"""
+Spark application that takes a preprocessed query file and DNA sequence
+and uses seed-and-reduce (e.g. BLAST) to find alignments in a distributed
+fashion (similar to CloudBurst's implementation on Hadoop).
+
+Spark implements RDDs - an efficient way of cacheing repeated data in
+memory, saving some of the overhead of disk IO that is requied in Hadoop,
+as such, Spark is a much faster processing platform for distributed
+alignments than Hadoop and Cloudburst.
+"""
+
+##########################################################################
+## Imports
+##########################################################################
+
+import sys
+import brisera
+
+from brisera.utils import timeit
+from brisera.config import settings
+from pyspark import SparkConf, SparkContext
+
+@timeit
+def run_brisera_alignment(sc, refpath, qrypath, outpath):
+    """
+    Runs the complete alignment
+    """
+    # Execute the alignments
+    alignments, adelta = brisera.align_all(sc, refpath, qrypath)
+
+    # Filter best alignments
+    if settings.filter_align:
+        alignments, fdelta = brisera.filter_alignments(sc, alignments)
+    else:
+        fdelta = 0
+
+    # Write alignments to disk
+    alignments.saveAsTextFile(outpath)
+
+    return alignments, adelta, fdelta
+
+if __name__ == "__main__":
+
+    if len(sys.argv) != 4:
+        sys.stderr.write("Usage: convert_fasta.py refpath qrypath outpath\n")
+        sys.stderr.write("    all other settings are stored in brisera.yaml\n")
+        sys.exit(-1)
+
+    conf = SparkConf().setAppName("Brisera Alignment")
+    sc   = SparkContext(conf=conf)
+
+    refpath = sys.argv[1]
+    qrypath = sys.argv[2]
+    outpath = sys.argv[3]
+
+    if settings.redundancy < 1:
+        raise brisera.ImproperlyConfigured("Minimum redundancy is 1")
+
+    if settings.max_read_len > settings.overlap:
+        raise brisera.ImproperlyConfigured("Increase overlap for %i length reads"
+            " and reconvert FASTA file.", settings.max_read_len)
+
+    result, delta = 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

brisera/__init__.py

@@ -4,3 +4,7 @@
 """
 
 __version__ = "1.0"
+
+from brisera.exceptions import *
+from brisera.align import align_all
+from brisera.filter import filter_alignments

brisera/align.py

@@ -0,0 +1,169 @@
+"""
+Implements MerReduce alignment (seed-and-reduce) in a distributed fashion
+"""
+
+##########################################################################
+## Imports
+##########################################################################
+
+from brisera.utils import *
+from brisera.records import *
+from brisera.config import settings
+from brisera.exceptions import *
+
+N = '.'
+
+##########################################################################
+## MerAlignment
+##########################################################################
+
+class MerAlignment(object):
+    """
+    Performs both map and reduce alignment akin to CloudBurst
+    """
+
+    def __init__(self, **kwargs):
+        setting = lambda name: kwargs.get(name, getattr(settings, name))
+
+        self.min_read_len = setting('min_read_len')
+        self.max_read_len = setting('max_read_len')
+        self.seed_len     = setting('seed_len')
+        self.flank_len    = setting('flank_len')
+        self.k            = setting('k')
+        self.redundancy   = setting('redundancy')
+
+    def parse_record(self, data):
+        key, val = data     # Expand the key, value pair
+        record     = deserialize_record(val)
+        sequence   = record[0]
+        offset     = record[1]
+        is_last     = record[2]
+        seqlen     = len(sequence)
+
+        return key, sequence, offset, is_last, seqlen
+
+    def map_reference(self, arg):
+        """
+        Input (id, (sequence, offset, tag))
+        Yields tuples:
+            (mer, (id, pos, tag, left, right, r))
+        for each seed in the sequences that are passed in
+        if tag = 0, also output the reverse complement sequences
+        """
+
+        key, sequence, offset, is_last, seqlen = self.parse_record(arg)
+
+        start = 0
+        if offset != 0:
+            # Not the first chunk, shift for room on left flank
+            start  = settings.overlap + 1 - self.flank_len - self.seed_len
+            offset += start
+
+        # stop so the last mer will fit
+        end = seqlen - self.seed_len + 1
+
+        if not is_last:
+            # If not the last chunk, stop so the right flank fits as well
+            end -= self.flank_len
+
+        # Emit the mers starting at every position
+        for idx in xrange(start, end):
+
+            seed   = sequence[start:self.seed_len]
+            if N in seed:
+                continue
+
+            offset += 1
+            start  += 1
+
+            leftstart = start - self.flank_len
+            if leftstart < 0:
+                leftstart = 0
+            leftlen = start-leftstart
+
+            rightstart = start + self.seed_len
+            rightend   = rightstart + self.flank_len
+            if rightend > seqlen:
+                rightend = seqlen
+            rightlen = rightend-rightstart
+
+            seed = sequence[start:start+self.seed_len]
+            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))
+            else:
+                yield (seed, (key, offset, is_last, leftstart, leftlen, rightstart, rightlen, 0))
+
+    def map_queries(self, arg):
+        """
+        Input (id, (sequence, offset, tag))
+        Yields tuples:
+            (mer, (id, pos, tag, left, right, r))
+        for each seed in the sequences that are passed in
+        if tag = 0, also output the reverse complement sequences
+        """
+        key, sequence, offset, is_last, seqlen = self.parse_record(arg)
+
+        if seqlen < self.min_read_len:
+            raise ReadLengthException("read length %i < minimum read length %i", seqlen, self.min_read_len)
+        elif seqlen > self.max_read_len:
+            raise ReadLengthException("read length %i > maximum read length %i", seqlen, self.max_read_len)
+
+        numN = sum(1 for char in sequence if char == N)
+
+        for rc in xrange(2):
+
+            if numN > self.k:
+                break
+
+            if rc == 1:
+                # Reverse complement the sequence
+                sequence = revc(sequence)
+                is_rc = True
+            else:
+                is_rc = False
+
+            # emit non-overlapping mers
+            for idx in xrange(0, seqlen, self.seed_len):
+                seed   = sequence[idx:idx+self.seed_len]
+                if N in seed:
+                    continue
+
+                rightstart = idx+self.seed_len
+                rightlen = seqlen - rightstart
+
+                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))
+                else:
+                    yield (seed, (key, idx, is_last, 0, idx, rightstart, rightlen, 0))
+
+##########################################################################
+## Spark Functionality
+##########################################################################
+
+@timeit
+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()
+
+    # Perform mapping
+    reference = reference.flatMap(alignment.map_reference)
+    return reference
+
+if __name__ == '__main__':
+    from brisera.convert import *
+    # path = fixture('s_suis.fa', 'cloudburst')
+    path = fixture('100k.fa', 'cloudburst')
+    chunker = FastaChunker(path)
+    aligner = MerAlignment()
+    for chunk in chunker.convert():
+        # for record in aligner.map_reference(chunk):
+        for record in aligner.map_queries(chunk):
+            print record
+        # break

brisera/config.py

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

brisera/convert.py

@@ -2,10 +2,19 @@
 Handles the conversion of a FASTA sequence into a sequence format
 """
 
+##########################################################################
+## Imports
+##########################################################################
+
 import cPickle
 
 from brisera.utils import fasta
 from brisera.config import settings
+from brisera.records import serialize_record
+
+##########################################################################
+## Chunker for RDDs
+##########################################################################
 
 class FastaChunker(object):
 
@@ -62,7 +71,7 @@ def convert(self):
 
         for idx, seq in self:
             for record in self.chunk(seq):
-                yield (idx, cPickle.dumps(record, cPickle.HIGHEST_PROTOCOL))
+                yield (idx, serialize_record(record))
 
     def __iter__(self):
         """

brisera/exceptions.py

@@ -0,0 +1,21 @@
+"""
+Class hierarchy for exceptions in Brisera
+"""
+
+class BriseraException(Exception):
+    """
+    Top level class for Brisera exceptions
+    """
+    pass
+
+class ImproperlyConfigured(BriseraException):
+    """
+    Something is wrong with a setting or configuration
+    """
+    pass
+
+class ReadLengthException(BriseraException):
+    """
+    The read is not in bounds of the minimum and maximum read lengths
+    """
+    pass

brisera/filter.py

@@ -0,0 +1,18 @@
+"""
+Filters the best alignments from the computed alignments.
+"""
+
+##########################################################################
+## Imports
+##########################################################################
+
+from brisera.utils import *
+from brisera.config import settings
+
+##########################################################################
+## Helper functions
+##########################################################################
+
+@timeit
+def filter_alignments(sc, alignments):
+    return alignments

brisera/records.py

@@ -2,6 +2,16 @@
 Utilities to help create and serialize records in binary format
 """
 
+##########################################################################
+## Imports
+##########################################################################
+
+import cPickle
+
+##########################################################################
+## Module Constants
+##########################################################################
+
 DNA_BYTES = {
     'A': 0x00,
     'C': 0x01,
@@ -14,6 +24,10 @@
 
 BYTES_DNA = dict((v, k) for (k,v) in DNA_BYTES.items())
 
+##########################################################################
+## Helper functions
+##########################################################################
+
 def dna_from_seq(dna, pos, length):
     if length == 0:
         return ""
@@ -46,6 +60,13 @@ def dna_from_seq(dna, pos, length):
 
     return string
 
+def repseed(seq, start, slen):
+    first = seq[start]
+    for idx in xrange(slen):
+        if seq[idx+start] != first:
+            return False
+    return True
+
 def record_from_bytes(raw):
 
     last_chunk = raw[0] == 1
@@ -59,6 +80,19 @@ def record_from_bytes(raw):
 
     return sequence, offset, last_chunk
 
+def serialize_record(record):
+    """
+    Convert a tuple into a binary string for use with SequenceFiles
+    """
+    return cPickle.dumps(record, 0)
+
+def deserialize_record(record):
+    """
+    Read a binary record object and return the tuple
+    """
+    record = record.encode('utf-8')
+    return cPickle.loads(record)
+
 if __name__ == '__main__':
     value = bytearray(b'\x01\x00\x00\x00\x00!\x14$AD@\x10B\x04DD"A@$$\x04"')
     print record_from_bytes(value)