diff --git a/Bio/Alignment/AAlign.hs b/Bio/Alignment/AAlign.hs
--- a/Bio/Alignment/AAlign.hs
+++ b/Bio/Alignment/AAlign.hs
@@ -26,17 +26,17 @@
 -- Edit distances
 
 -- | Calculate global edit distance (Needleman-Wunsch alignment score)
-global_score :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> a
+global_score :: (Num a, Ord a) => SubstMx t a -> (a,a) -> Sequence t -> Sequence t -> a
 global_score mx g s1 s2 = uncurry max . last . last 
                           $ columns (score_select minf mx g) (0,fst g) s1 s2
 
 -- | Calculate local edit distance (Smith-Waterman alignment score)
-local_score :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> a
+local_score :: (Num a, Ord a) => SubstMx t a -> (a,a) -> Sequence t -> Sequence t -> a
 local_score mx g s1 s2 = maximum . map (uncurry max) . concat 
                          $ columns (score_select 0 mx g) (0,fst g) s1 s2
 
 -- | Generic scoring and selection function for global and local scoring
-score_select :: (Num a,Ord a) => a -> SubstMx a -> (a,a) -> Selector (a,a)
+score_select :: (Num a,Ord a) => a -> SubstMx t a -> (a,a) -> Selector (a,a)
 score_select minf mx (go,ge) cds = 
     let (reps,ids) = partition (isRepl.snd) cds 
         s = maximum $ minf:[max sub gap +mx (x,y) | ((sub,gap),Repl x y) <- reps]
@@ -51,12 +51,12 @@
 fp (x,ax) (s,e) = (x+s,e:ax)
 
 -- | Calculate global alignment (Needleman-Wunsch)
-global_align :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> (a,EditList)
+global_align :: (Num a, Ord a) => SubstMx t a -> (a,a) -> Sequence t -> Sequence t -> (a,EditList)
 global_align mx g s1 s2 = revsnd . uncurry max' . last . last 
                $ columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2
 
 -- | Calculate local alignmnet (Smith-Waterman)
-local_align :: (Num a, Ord a) => SubstMx a -> (a,a) -> Sequence -> Sequence -> (a,EditList)
+local_align :: (Num a, Ord a) => SubstMx t a -> (a,a) -> Sequence t -> Sequence t -> (a,EditList)
 local_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst)
                          . map (uncurry max') . concat
                          $ columns (align_select 0 mx g) ((0,[]),(fst g,[])) s1 s2
@@ -64,7 +64,7 @@
 revsnd (s,a) = (s,reverse a)
 
 -- | Generic scoring and selection for global and local alignment
-align_select :: (Num a, Ord a) => a -> SubstMx a -> (a,a) -> Selector ((a,EditList),(a,EditList))
+align_select :: (Num a, Ord a) => a -> SubstMx t a -> (a,a) -> Selector ((a,EditList),(a,EditList))
 align_select minf mx (go,ge) cds = 
     let (reps,ids) = partition (isRepl.snd) cds 
         s = maximumBy (compare `on` fst) 
diff --git a/Bio/Alignment/ACE.hs b/Bio/Alignment/ACE.hs
--- a/Bio/Alignment/ACE.hs
+++ b/Bio/Alignment/ACE.hs
@@ -25,6 +25,8 @@
           DS (phred header? left empty by CAP3)
           RD read2 ...
    @
+
+   As far as I know, this is only used for nucleotide sequences.
 -}
 
 {-# LANGUAGE CPP #-}
@@ -47,11 +49,11 @@
 import Control.Monad (liftM) -- ,when?
 import Data.Char (chr)
 
-data Assembly = Asm { contig :: (Sequence,Gaps), fragments :: Alignment }
+data Assembly = Asm { contig :: (Sequence Nuc,Gaps), fragments :: Alignment Nuc}
                 deriving Show
 
 {-# DEPRECATED reads "Stupid name, replaced by 'fragments'." #-}
-reads :: Assembly -> Alignment
+reads :: Assembly -> Alignment Nuc
 reads = fragments -- deprecated, stupid, stupid name.
 
 type Str = ByteString
@@ -121,7 +123,7 @@
                                 _ -> Nothing) <?> "empty line"
 
 -- | parse the contig and quality information (CO, BQ)
-ctg :: AceParser (Sequence,Gaps)
+ctg :: AceParser (Sequence Nuc,Gaps)
 ctg = do
   name <- co
   sd <- sdata
@@ -163,7 +165,7 @@
 bq = parse1 (\t -> case t of BQ -> Just (); _ -> Nothing) <?> "BQ"
 
 -- | Given the CO info, get the AFS'es
-asm :: (Sequence,Gaps) -> AceParser Assembly
+asm :: (Sequence Nuc,Gaps) -> AceParser Assembly
 asm cg = do
   many blank
   afs <- many1 af
@@ -191,14 +193,14 @@
 readInt' :: Str -> Maybe Int
 readInt' = liftM (fromIntegral . fst) . B.readInt
 
-rds :: (Sequence,Gaps) -> [(Str,Dir,Offset)] -> AceParser Assembly
+rds :: (Sequence Nuc,Gaps) -> [(Str,Dir,Offset)] -> AceParser Assembly
 rds cg xs = do
      r <- many1 rseq
      -- todo: check the number and merge with the afs
      let f (_name,d,off) (s,gs) = (off,d,s,gs)
      return $ Asm { contig = cg, fragments = zipWith f xs r }
 
-rseq :: AceParser (Sequence,Gaps)
+rseq :: AceParser (Sequence Nuc,Gaps)
 rseq = do
   (rn,_len,_,_) <- rd
   (s,gaps) <- return . extractGaps =<< sdata
diff --git a/Bio/Alignment/AlignData.hs b/Bio/Alignment/AlignData.hs
--- a/Bio/Alignment/AlignData.hs
+++ b/Bio/Alignment/AlignData.hs
@@ -37,7 +37,7 @@
 
 data Dir = Fwd | Rev deriving (Eq,Show)
 type Gaps = [Offset]
-type Alignment = [(Offset,Dir,Sequence,Gaps)]
+type Alignment a = [(Offset,Dir,Sequence a,Gaps)]
 
 -- | Gaps are coded as '*'s, this function removes them, and returns
 --   the sequence along with the list of gap positions.
@@ -84,11 +84,11 @@
 isRepl _ = False
 
 -- | A substitution matrix gives scores for replacing a character with another.
---   Typically, it will be symmetric.
-type SubstMx a = (Chr,Chr) -> a
+--   Typically, it will be symmetric.  It is type-tagged with the alphabet - Nuc or Amino.
+type SubstMx t a = (Chr,Chr) -> a
 
 -- | Evaluate an Edit based on SubstMx and gap penalty
-eval :: SubstMx a -> a -> Edit -> a
+eval :: SubstMx t a -> a -> Edit -> a
 eval mx g c = case c of Ins _ -> g; Del _ -> g; Repl x y -> mx (x,y)
 
 -- | A Selector consists of a zero element, and a funcition
@@ -99,7 +99,7 @@
 -- | Calculate a set of columns containing scores
 --   This represents the columns of the alignment matrix, but will only require linear space
 --   for score calculation.
-columns :: Selector a -> a -> Sequence -> Sequence -> [[a]]
+columns :: Selector a -> a -> Sequence b -> Sequence b -> [[a]]
 columns f z (Seq _ s1 _) (Seq _ s2 _) = columns' f z s1 s2
 
 columns' :: Selector a -> a -> SeqData -> SeqData -> [[a]]
diff --git a/Bio/Alignment/Bowtie.hs b/Bio/Alignment/Bowtie.hs
new file mode 100644
--- /dev/null
+++ b/Bio/Alignment/Bowtie.hs
@@ -0,0 +1,153 @@
+{-| This module provides a data type to represent an alignment
+produced by the Bowtie short-read alignment tool (see
+<http://bowtie-bio.sourceforge.net/index.shtml>).
+
+The simple accessors recapitulate the details of the Bowtie alignment
+output.  The position of the alignment is given by the \"0-based offset
+into the reference sequence where leftmost character of the alignment
+occurs\".  Thus, for forward-strand alignments this is the 5\' end of
+the query sequence while for reverse-complement alignments this is the
+3\' end of the query sequence.  Similarly, the query sequence and query
+quality are shown in reference forward strand orientation, and thus
+may be reverse complemented.
+
+-}
+
+module Bio.Alignment.Bowtie ( -- * Data type and basic accessors
+                              Align(..), Mismatch(..), length, nmismatch, querySequ, queryQual
+
+                            -- * Sequence positions of alignments
+                            , refCLoc, refCSeqLoc, refSeqLoc, refSeqPos, mismatchSeqPos
+
+                            -- * Parsing Bowtie output
+                            , parse
+
+                            -- * Other utilities
+                            , sameRead
+                            )
+    where 
+
+import Prelude hiding (length)
+import Control.Monad.Error
+import Control.Monad.State
+import qualified Data.ByteString.Lazy as LBSW
+import qualified Data.ByteString.Lazy.Char8 as LBS
+import Data.Char
+import qualified Data.List as List (length)
+import Data.Maybe
+
+import qualified Bio.Location.ContigLocation as CLoc
+import qualified Bio.Location.Location as Loc
+import qualified Bio.Location.Position as Pos
+import Bio.Location.OnSeq
+import qualified Bio.Location.SeqLocation as SeqLoc
+import Bio.Location.Strand
+import Bio.Sequence.SeqData
+
+data Align = Align { name :: !SeqName -- ^ Name of the query sequence
+                   , strand :: !Strand -- ^ Strand of the alignment on the reference sequence
+                   , refname :: !SeqName -- ^ Name of the reference sequence
+                   , leftoffset :: !Offset -- ^ Zero-based offset of the left-most aligned position in the reference
+                   , sequ :: !SeqData -- ^ Query sequence, in the reference forward strand orientation
+                   , qual :: !QualData -- ^ Query quality, in the reference forward strand orientation
+                   , mismatches :: ![Mismatch] -- ^ Mismatches
+                   } deriving (Read, Show, Eq, Ord)
+
+-- | Returns the length of the query sequence
+length :: Align -> Offset
+length = LBS.length . sequ
+
+-- | Returns the number of mismatches in the alignment
+nmismatch :: Align -> Int
+nmismatch = List.length . mismatches
+
+-- | Parses a line of Bowtie output to produce a 'Align'
+parse :: LBS.ByteString -> Either String Align
+parse bstr
+    = case LBS.split '\t' bstr of
+        [nameBStr,strandBStr,refnameBStr,leftoffBStr,sequBStr,qualBStr,_,mismatchesBStr]
+            -> do str <- parseStrand strandBStr
+                  loff <- parseInt64 leftoffBStr
+                  mms <- mapM parseMismatch . LBS.split (',') $ mismatchesBStr
+                  return $ Align (LBS.copy nameBStr) str (LBS.copy refnameBStr) loff (LBS.copy sequBStr) (LBS.copy qualBStr) mms
+        _ -> throwError $ strMsg $ "Malformed Bowtie alignment " ++ (show . LBS.unpack) bstr
+    where parseStrand str | str == (LBS.singleton '+') = return Fwd
+                          | str == (LBS.singleton '-') = return RevCompl
+                          | otherwise = throwError $ "Unknown strand " ++ (show . LBS.unpack $ str)
+
+-- | Query sequence as given in the query file
+querySequ :: Align -> SeqData
+querySequ ba = stranded (strand ba) (sequ ba)
+
+-- | Query quality as given in the query file
+queryQual :: Align -> QualData
+queryQual ba = case strand ba of
+                 Fwd -> qual ba
+                 RevCompl -> LBSW.reverse $ qual ba
+
+-- | Returns the sequence position of the start of the query sequence
+-- alignment.  This will include the strand of the alignment and will
+-- not be the same as the position computed from 'leftoffset' when the
+-- alignment is on the reverse complement strand.
+refSeqPos :: Align -> SeqLoc.SeqPos
+refSeqPos ba = OnSeq (refname ba) $ CLoc.startPos $ refCLoc ba
+
+-- | Returns the sequence location covered by the query in
+-- the alignment.  This will be a sequence location on the reference
+-- sequence and may run on the forward or the reverse complement
+-- strand.
+refCSeqLoc :: Align -> SeqLoc.ContigSeqLoc
+refCSeqLoc ba = OnSeq (refname ba) (refCLoc ba)
+
+-- | As 'refCSeqLoc' but without the reference sequence name.
+refCLoc :: Align -> CLoc.ContigLoc
+refCLoc ba = CLoc.ContigLoc (leftoffset ba) (length ba) (strand ba)
+
+-- | Returns the sequence location covered by the query, as
+-- 'refCSeqLoc', as a 'SeqLoc.SeqLoc' location.
+refSeqLoc :: Align -> SeqLoc.SeqLoc 
+refSeqLoc ba = OnSeq (refname ba) (Loc.Loc [ refCLoc ba ])
+
+-- | Returns true when two alignments were derived from the same
+-- sequencing read.  As Bowtie writes alignments of query sequences in
+-- their order in the query file, all alignments of a given read are
+-- grouped together and the lists of all alignments for each read can
+-- be gathered with
+-- 
+-- > groupBy sameRead
+sameRead :: Align -> Align -> Bool
+sameRead ba1 ba2 = (name ba1) == (name ba2)
+
+-- | Representation of a single mismatch in a bowtie alignment
+data Mismatch = Mismatch { mmoffset :: !Offset -- ^ Offset of the mismatch site from the 5\' end of the query
+                         , refbase :: !Char -- ^ Reference nucleotide
+                         , readbase :: !Char -- ^ Query nucleotide
+                         } deriving (Read, Show, Eq, Ord)
+
+-- | Sequence position of a mismatch on the reference sequence.
+mismatchSeqPos :: Align -> Mismatch -> SeqLoc.SeqPos
+mismatchSeqPos ba mm = OnSeq (refname ba) mmpos
+    where mmpos = fromMaybe badOffset $ (Pos.Pos (mmoffset mm) Fwd) `CLoc.posOutof` (refCLoc ba)
+          badOffset = error $ "Bad mismatch offset: " ++ show (ba, mm)
+
+parseMismatch :: LBS.ByteString -> Either String Mismatch
+parseMismatch = evalStateT parser
+    where parser = liftM3 Mismatch parseOffset parseRefNt parseReadNt
+          parseOffset = parseTo ':' >>= lift . parseInt64
+          parseRefNt = parseTo '>' >>= parseChar
+          parseReadNt = parseToEnd >>= parseChar
+          parseTo ch = StateT $ \str ->
+                       case LBS.span (/= ch) str of
+                         (out, rest) -> do (_, afterCh) <- maybe missingChError Right . LBS.uncons $ rest
+                                           return (out, afterCh)
+              where missingChError = Left ("Failed to find " ++ show ch)
+          parseToEnd = StateT $ \str -> return (str, LBS.empty)
+          parseChar str = lift $ case LBS.uncons str of
+                                   Just (ch, rest) | LBS.null rest -> return ch
+                                   _ -> Left ("Malformed character " ++ (show . LBS.unpack) str) 
+
+parseInt64 :: LBS.ByteString -> Either String Offset
+parseInt64 zstr = case LBS.readInteger zstr of
+                    Just (z, rest) | LBS.null rest -> return $ fromIntegral z
+                    _ -> throwError $ "Malformed integer " ++ show zstr
+
diff --git a/Bio/Alignment/Multiple.hs b/Bio/Alignment/Multiple.hs
--- a/Bio/Alignment/Multiple.hs
+++ b/Bio/Alignment/Multiple.hs
@@ -14,7 +14,7 @@
 -- | Progressive multiple alignment.
 --   Calculate a tree from agglomerative clustering, then align
 --   at each branch going bottom up.  Returns a list of columns (rows?).
-progressive :: (Sequence -> Sequence -> (Double,EditList)) -> [Sequence] -> [String]
+progressive :: (Sequence a -> Sequence a -> (Double,EditList)) -> [Sequence a] -> [String]
 progressive = undefined
 
 -- |  Derive alignments indirectly, i.e. calculate A|C using alignments A|B and B|C.
diff --git a/Bio/Alignment/QAlign.hs b/Bio/Alignment/QAlign.hs
--- a/Bio/Alignment/QAlign.hs
+++ b/Bio/Alignment/QAlign.hs
@@ -49,10 +49,10 @@
 type QSelector a = [(a,Edit,Qual,Qual)] -> a
 
 -- using 21 as default corresponds to blastn default ratio (1,-3)
-columns :: QSelector a -> a -> Sequence -> Sequence -> [[a]]
+columns :: QSelector a -> a -> Sequence t -> Sequence t -> [[a]]
 columns sel z s1 s2 = columns' sel (z,r0,c0) s1' s2'
     where (s1',s2') = (zup s1, zup s2)
-          zup :: Sequence -> [(Chr,Qual)]
+          zup :: Sequence t -> [(Chr,Qual)]
           zup (Seq _ sd Nothing) = map (\c -> (c,22)) $ B.unpack sd
           zup (Seq _ sd (Just qd)) = zip (B.unpack sd) (B.unpack qd)
 
@@ -88,7 +88,7 @@
 minf :: Double
 minf = -100000000
 
-type QualMx a = Qual -> Qual -> SubstMx a
+type QualMx t a = Qual -> Qual -> SubstMx t a
 
 qualMx :: Qual -> Qual -> (Chr,Chr) -> Double
 qualMx q1 q2 (x,y) = if isN x || isN y then 0.0 else
@@ -116,12 +116,12 @@
 -- Edit distances
 
 -- | Calculate global edit distance (Needleman-Wunsch alignment score)
-global_score :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> Double
+global_score :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> Double
 global_score mx g s1 s2 = uncurry max . last . last
                           $ columns (score_select minf mx g) (0,fst g) s1 s2
 
 -- | Calculate local edit distance (Smith-Waterman alignment score)
-local_score :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> Double
+local_score :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> Double
 local_score mx g s1 s2 = maximum . map (uncurry max) . concat
                          $ columns (score_select 0 mx g) (0,fst g) s1 s2
 
@@ -132,7 +132,7 @@
 -- Oh. local score_select will not work, since we cannot replace any matrix entry
 -- with zero in order to initiate a new alignment.  So we need 'minf', except in the
 -- initial row/column.  Damn.
-overlap_score :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> Double
+overlap_score :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> Double
 overlap_score mx g s1 s2 = maximum $ map (uncurry max) $ sel cols
     where cols   = columns overlap_score_select (0,fst g) s1 s2 
           sel cs = map last cs ++ last cs
@@ -141,7 +141,7 @@
           overlap_score_select [_] = (0,minf)
 
 -- | Generic scoring and selection function for global and local scoring
-score_select :: Double -> QualMx Double -> (Double,Double) -> QSelector (Double,Double)
+score_select :: Double -> QualMx t Double -> (Double,Double) -> QSelector (Double,Double)
 score_select minf mx (go,ge) cds =
     let (reps,ids) = partition (isRepl.snd') cds
         s = maximum $ minf:[max sub gap + mx q1 q2 (x,y) | ((sub,gap),Repl x y,q1,q2) <- reps]
@@ -161,13 +161,13 @@
 fp (x,ax) (s,e) = (x+s,e:ax)
 
 -- | Calculate global alignment (Needleman-Wunsch)
-global_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList)
+global_align :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> (Double,EditList)
 global_align mx g s1 s2 = revsnd . uncurry max' . last . last
                $ columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2
 
 -- | Calculate local alignment (Smith-Waterman)
 --   (can we replace uncurry max' with fst - a local alignment must always end on a subst, no?)
-local_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double, EditList)
+local_align :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> (Double, EditList)
 local_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst)
                          . map (uncurry max') . concat
                          $ columns (align_select 0 mx g) ((0,[]),(fst g,[])) s1 s2
@@ -175,7 +175,7 @@
 -- | Calucalte best overlap score, where gaps at the edges are free
 --   The starting point is like for local score (0 cost for initial indels),
 --   the result is the maximum anywhere in the last column or bottom row of the matrix.
-overlap_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList)
+overlap_align :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> (Double,EditList)
 overlap_align mx g s1 s2 = revsnd . maximumBy (compare `on` fst) . map (uncurry max') $ sel cols
     where cols   = columns overlap_align_select ((0,[]),(minf,[])) s1 s2 
           sel cs = map last cs ++ last cs
@@ -184,7 +184,7 @@
           overlap_align_select [_] = ((0,[]),(minf,[]))
 
 -- | Variant that retains indels to retain the entire sequence in the result
-overlap_align' :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,EditList)
+overlap_align' :: QualMx t Double -> (Double,Double) -> Sequence t -> Sequence t -> (Double,EditList)
 overlap_align' mx g s1 s2 = revsnd . maximumBy (compare `on` fst) . map (uncurry max') $ sel cols
     where cols   = columns overlap_align_select ((0,[]),(fst g,[])) s1 s2 
           sel cs = map last cs ++ last cs
@@ -196,7 +196,7 @@
 revsnd (s,a) = (s,reverse a)
 
 -- | Generic scoring and selection for global and local alignment
-align_select :: Double -> QualMx Double -> (Double,Double) -> QSelector ((Double,EditList),(Double,EditList))
+align_select :: Double -> QualMx t Double -> (Double,Double) -> QSelector ((Double,EditList),(Double,EditList))
 align_select minf mx (go,ge) cds =
     let (reps,ids) = partition (isRepl.snd') cds
         s = maximumBy (compare `on` fst)
diff --git a/Bio/Alignment/SAlign.hs b/Bio/Alignment/SAlign.hs
--- a/Bio/Alignment/SAlign.hs
+++ b/Bio/Alignment/SAlign.hs
@@ -22,19 +22,19 @@
 -- Edit distances
 
 -- | Calculate global edit distance (Needleman-Wunsch alignment score)
-global_score :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> a
+global_score :: (Num a, Ord a) => SubstMx t a -> a -> Sequence t -> Sequence t -> a
 global_score mx g s1 s2 = last . last $ columns (g_score mx g) 0 s1 s2
 
 -- | Scoring\/selection function for global alignment
-g_score:: (Num a,Ord a) => SubstMx a -> a -> Selector a
+g_score:: (Num a,Ord a) => SubstMx t a -> a -> Selector a
 g_score mx g cds = maximum [s+eval mx g e | (s,e) <- cds]
 
 -- | Calculate local edit distance (Smith-Waterman alignment score)
-local_score :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> a
+local_score :: (Num a, Ord a) => SubstMx t a -> a -> Sequence t -> Sequence t -> a
 local_score mx g s1 s2 = maximum . concat $ columns (l_score mx g) 0 s1 s2
 
 -- | Scoring\/selection funciton for local alignmnet
-l_score :: (Num a,Ord a) => SubstMx a -> a -> Selector a
+l_score :: (Num a,Ord a) => SubstMx t a -> a -> Selector a
 l_score mx g cds = maximum (0:[s+eval mx g e | (s,e) <- cds])
 
 -- ------------------------------------------------------------
@@ -46,17 +46,17 @@
 -- the score in each cell.  Unreachable cells can then be GC'ed.
 
 -- | Calculate alignments.
-global_align :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> EditList
+global_align :: (Num a, Ord a) => SubstMx t a -> a -> Sequence t -> Sequence t -> EditList
 global_align mx g s1 s2 = reverse . snd . last . last 
                           $ columns (g_align mx g) (0,[]) s1 s2
 
-g_align :: (Num a, Ord a) => SubstMx a -> a -> Selector (a,EditList)
+g_align :: (Num a, Ord a) => SubstMx t a -> a -> Selector (a,EditList)
 g_align mx g cds = maximumBy (compare `on` fst) [(s+eval mx g e,e:a) | ((s,a),e) <- cds]
 
-local_align :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> EditList 
+local_align :: (Num a, Ord a) => SubstMx t a -> a -> Sequence t -> Sequence t -> EditList 
 local_align mx g s1 s2 = reverse . snd . maximumBy (compare `on` fst) . concat 
                          $ columns (l_align mx g) (0,[]) s1 s2
 
-l_align :: (Num a, Ord a) => SubstMx a -> a -> Selector (a,EditList)
+l_align :: (Num a, Ord a) => SubstMx t a -> a -> Selector (a,EditList)
 l_align mx g cds = maximumBy (compare `on` fst) $ (0,[]):[(s+eval mx g e,e:a) | ((s,a),e) <- cds]
 
diff --git a/Bio/GFF3/FeatureHierSequences.hs b/Bio/GFF3/FeatureHierSequences.hs
--- a/Bio/GFF3/FeatureHierSequences.hs
+++ b/Bio/GFF3/FeatureHierSequences.hs
@@ -32,7 +32,7 @@
 parse str = do (feats, seqlines) <- F.parseWithFasta str
                fromLists feats $ mkSeqs seqlines
 
-fromLists :: (Error e, MonadError e m) => [F.Feature] -> [Sequence] -> m FeatureHierSequences
+fromLists :: (Error e, MonadError e m) => [F.Feature] -> [Sequence a] -> m FeatureHierSequences
 fromLists feats seqs = let !seqmap = M.fromList $ map seqAssoc seqs
                            !featSeqids = S.fromList $ map F.seqid feats
                            !missingSeqids = featSeqids `S.difference` M.keysSet seqmap
@@ -44,7 +44,7 @@
 features :: FeatureHierSequences -> S.Set F.Feature
 features = FH.features . featureHier
 
-sequences :: FeatureHierSequences -> [Sequence]
+sequences :: FeatureHierSequences -> [Sequence a]
 sequences = M.foldWithKey mkSeqList [] . sequenceMap
     where mkSeqList seqname sequ = ((Seq seqname sequ Nothing) :)
 
@@ -65,7 +65,7 @@
 getSequence fhs seqid = maybe seqidNotFound return $ M.lookup seqid $ sequenceMap fhs
     where seqidNotFound = throwError $ strMsg $ "SeqID " ++ show (LBS.unpack seqid) ++ " not found"
 
-featureSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m Sequence
+featureSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m (Sequence a)
 featureSequence fhs f = do sequ <- seqData fhs $ F.seqLoc f
                            let seqname = fromMaybe (LBS.pack $ SeqLoc.display $ F.seqLoc f) $ listToMaybe $ F.ids f
                            return $ Seq seqname sequ Nothing
diff --git a/Bio/GFF3/SGD.hs b/Bio/GFF3/SGD.hs
--- a/Bio/GFF3/SGD.hs
+++ b/Bio/GFF3/SGD.hs
@@ -46,7 +46,7 @@
     where isRRNA = any (flip S.member $ S.fromList rRNAFeatureIDs) . F.ids
           rRNAFeatureIDs = map LBS.pack ["RDN25-1", "RDN58-1", "RDN18-1", "RDN5-1"]
 
-geneSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m Sequence
+geneSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m (Sequence a)
 geneSequence fhs g = do sequ <- geneSeqLoc fhs g >>= SeqLoc.seqData (getSequence fhs)
                         seqname <- maybe (throwError $ strMsg $ "No gene ID for " ++ show g) return $ listToMaybe $ F.ids g
                         return $ Seq seqname sequ Nothing
@@ -58,7 +58,7 @@
 geneCDSes fhs g = filter isCDS $ children fhs g
     where isCDS = (== cdsType) . F.ftype          
 
-noncodingSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m Sequence
+noncodingSequence :: (Error e, MonadError e m) => FeatureHierSequences -> F.Feature -> m (Sequence a)
 noncodingSequence fhs nc = do sequ <- noncodingSeqLoc fhs nc >>= SeqLoc.seqData (getSequence fhs)
                               seqname <- maybe (throwError $ strMsg $ "No gene ID for " ++ show nc) return $ listToMaybe $ F.ids nc
                               return $ Seq seqname sequ Nothing
diff --git a/Bio/Location/ContigLocation.hs b/Bio/Location/ContigLocation.hs
--- a/Bio/Location/ContigLocation.hs
+++ b/Bio/Location/ContigLocation.hs
@@ -1,15 +1,29 @@
--- |Data types for working with locations in a sequence.  Zero-based
---  'Int64' indices are used throughout, to facilitate direct use of
---  indexing functions on 'SeqData'.
+{-| Data type for a sequence location consiting of a contiguous range
+of positions on the sequence.
 
-module Bio.Location.ContigLocation ( ContigLoc(..), fromStartEnd, fromPosLen
-                                   , bounds, startPos, endPos
-                                   , slide, extend, posInto, posOutof
-                                   , seqData, seqDataPadded, isWithin, overlaps
-                                   , display
-                                   )
-    where 
+Throughout, /sequence position/ refers to a 'Pos.Pos' which includes a
+strand.  An index into a sequence is referred to as an /offset/, and
+is generally of type 'Offset'.
 
+ -}
+
+module Bio.Location.ContigLocation ( 
+  -- * Sequence locations
+  ContigLoc(..), fromStartEnd, fromPosLen
+
+  -- * Locations and positions
+  , bounds, startPos, endPos, posInto, posOutof, isWithin, overlaps
+
+  -- * Extracting subsequences
+  , seqData, seqDataPadded
+
+  -- * Transforming locations
+  , slide, extend
+
+  -- * Displaying locations
+  , display
+                                   ) where
+
 import Prelude hiding (length)
 
 import Control.Monad.Error
@@ -20,56 +34,69 @@
 import qualified Bio.Location.Position as Pos
 import Bio.Location.Strand
 
--- | Contiguous set of positions in a sequence
-data ContigLoc = ContigLoc { offset5 :: !Offset   -- ^ 5' end of region on target sequence, 0-based index
-                           , length :: !Offset -- ^ length of region on target sequence
-                           , strand :: !Strand -- ^ strand of region
+-- | Contiguous sequence location defined by a span of sequence
+-- positions, lying on a specific strand of the sequence.
+data ContigLoc = ContigLoc { offset5 :: !Offset   -- ^ The offset of the 5\' end of the location, as a 0-based index
+                           , length :: !Offset    -- ^ The length of the location
+                           , strand :: !Strand    -- ^ The strand of the location
                            } deriving (Eq, Ord, Show)
 
 instance Stranded ContigLoc where
     revCompl (ContigLoc seq5 len str) = ContigLoc seq5 len $ revCompl str
 
--- | Create a 'ContigLoc' from 0-based starting and ending positions.
---   When 'start' is less than 'end' the position will be on the 'Fwd'
---   'Strand', otherwise it will be on the 'RevCompl' strand.
+-- | Create a sequence location lying between 0-based starting and
+-- ending offsets.  When @start < end@, the location
+-- be on the forward strand, otherwise it will be on the
+-- reverse complement strand.
+
 fromStartEnd :: Offset -> Offset -> ContigLoc
 fromStartEnd start end
     | start < end = ContigLoc start (1 + end - start) Fwd
     | otherwise   = ContigLoc end   (1 + start - end) RevCompl
 
--- | Create a 'ContigLoc' from a Pos.Pos defining the start
--- ('ContigLoc' 5 prime end) position on the sequence and the length.
+-- | Create a sequence location from the sequence position of the
+-- start of the location and the length of the position.  The strand
+-- of the location, and the direction it extends from the starting
+-- position, are determined by the strand of the starting position.
 fromPosLen :: Pos.Pos -> Offset -> ContigLoc
 fromPosLen (Pos.Pos off5 Fwd)      len = ContigLoc off5               len Fwd
 fromPosLen (Pos.Pos off3 RevCompl) len = ContigLoc (off3 - (len - 1)) len RevCompl
 
--- | The bounds of a 'ContigLoc', a pair of the lowest and highest
---   sequence indices covered by the region, which ignores the strand
---   of the 'ContigLoc'.  The first element of the pair will always be
---   lower than the second.
+-- | The bounds of a sequence location.  This is a pair consisting of
+-- the lowest and highest sequence offsets covered by the region.  The
+-- bounds ignore the strand of the sequence location, and the first
+-- element of the pair will always be lower than the second.
 bounds :: ContigLoc -> (Offset, Offset)
 bounds (ContigLoc seq5 len _) = (seq5, seq5 + len - 1)
 
--- | 0-based starting (5' in the region orientation) position
+-- | Sequence position of the start of the location.  This is the 5'
+-- end on the location strand, which will have a higher offset than
+-- 'endPos' if the location is on the 'RevCompl' strand.
 startPos :: ContigLoc -> Pos.Pos
 startPos (ContigLoc seq5 len str) 
     = case str of
         Fwd      -> Pos.Pos seq5             str
         RevCompl -> Pos.Pos (seq5 + len - 1) str
 
--- | 0-based ending (3' in the region orientation) position
+-- | Sequence position of the end of the location, as described in 'startPos'.
 endPos :: ContigLoc -> Pos.Pos
 endPos (ContigLoc seq5 len str) 
     = case str of
         Fwd      -> Pos.Pos (seq5 + len - 1) str
         RevCompl -> Pos.Pos seq5             str
 
--- | Move a 'ContigLoc' region by a specified offset
+-- | Returns a location resulting from sliding the original location
+-- along the sequence by a specified offset.  A positive offset will
+-- move the location away from the 5\' end of the forward stand of the
+-- sequence regardless of the strand of the location itself.  Thus,
+-- 
+-- > slide (revCompl cloc) off == revCompl (slide cloc off)
 slide :: Offset -> ContigLoc -> ContigLoc
 slide dpos (ContigLoc seq5 len str) = ContigLoc (seq5 + dpos) len str
 
--- | Subsequence 'SeqData' for a 'ContigLoc', provided that the region
---   is entirely within the sequence.
+-- | Extract the nucleotide 'SeqData' for the sequence location.  If
+-- any part of the location lies outside the bounds of the sequence,
+-- an error results.
 seqData :: (Error e, MonadError e m) => SeqData -> ContigLoc -> m SeqData
 seqData sequ (ContigLoc seq5 len str)
     | seq5 < 0 = outOfBounds
@@ -78,7 +105,9 @@
                            | otherwise -> outOfBounds
     where outOfBounds = throwError $ strMsg $ "contig seq loc " ++ show (seq5, seq5 + len - 1) ++ " out of SeqData bounds"
 
--- | Subsequence 'SeqData' for a 'ContigLoc', padded as needed with Ns
+-- | As 'seqData', extract the nucleotide subsequence for the
+-- location.  Any positions in the location lying outside the bounds
+-- of the sequence are returned as @N@ rather than producing an error.
 seqDataPadded :: SeqData -> ContigLoc -> SeqData
 seqDataPadded sequ (ContigLoc seq5 len str) = stranded str fwdseq
     where fwdseq
@@ -90,9 +119,12 @@
               | sublen <= LBS.length subsequ = LBS.take sublen subsequ
               | otherwise = subsequ `LBS.append` LBS.replicate (sublen - LBS.length subsequ) 'N'
 
--- | For a 'Pos' and a 'ContigLoc' on the same sequence, find the
---   corresponding 'Pos' relative to the 'ContigLoc', provided it is
---   within the 'ContigLoc'.
+-- | Given a sequence position and a sequence location relative to the
+-- same sequence, compute a new position representing the original
+-- position relative to the subsequence defined by the location.  If
+-- the sequence position lies outside of the sequence location,
+-- @Nothing@ is returned; thus, the offset of the new position will
+-- always be in the range @[0, length cloc - 1]@.
 posInto :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos
 posInto (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)
     | pos < seq5 || pos >= seq5 + len = Nothing
@@ -100,9 +132,14 @@
                            Fwd      -> Pos.Pos (pos - seq5)              pStrand
                            RevCompl -> Pos.Pos (seq5 + len  - (pos + 1)) (revCompl pStrand)
 
--- | For a 'Pos' specified relative to a 'ContigLoc', find the
---   corresponding 'Pos' relative to the outer sequence, provided that
---   the 'Pos' is within the bounds of the 'ContigLoc'.
+
+-- | Given a sequence location and a sequence position within that
+-- location, compute a new position representing the original position
+-- relative to the outer sequence.  If the sequence position lies
+-- outside the location, @Nothing@ is returned.
+-- 
+-- This function inverts 'posInto' when the sequence position lies
+-- within the position is actually within the location.
 posOutof :: Pos.Pos -> ContigLoc -> Maybe Pos.Pos
 posOutof (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)
     | pos < 0 || pos >= len = Nothing
@@ -110,26 +147,34 @@
                            Fwd      -> Pos.Pos (pos + seq5)             pStrand
                            RevCompl -> Pos.Pos (seq5 + len - (pos + 1)) (revCompl pStrand)
 
--- | 'ContigLoc' extended on the 5' and 3' ends.
+-- | Returns a sequence location produced by extending the original
+-- location on each end, based on a pair of (/5\' extension/, /3\'
+-- extension/).  The 5\' extension is applied to the 5\' end of the
+-- location on the location strand; if the location is on the
+-- 'RevCompl' strand, the 5\' end will have a higher offset than the
+-- 3\' end and this offset will increase by the amount of the 5\'
+-- extension.  Similarly, the 3\' extension is applied to the 3\' end
+-- of the location.
 extend :: (Offset, Offset) -> ContigLoc -> ContigLoc
 extend (ext5, ext3) (ContigLoc seq5 len str)
     = case str of
         Fwd -> ContigLoc (seq5 - ext5) (len + ext5 + ext3) str
         RevCompl -> ContigLoc (seq5 - ext3) (len + ext5 + ext3) str
 
--- | For a 'Pos' and a 'ContigLoc' on the same sequence, is the 'Pos'
---   within the 'ContigLoc'.
+-- | Returns @True@ when  a sequence position lies within a sequence
+-- location on the same sequence, and occupies the same strand.
 isWithin :: Pos.Pos -> ContigLoc -> Bool
 isWithin (Pos.Pos pos pStrand) (ContigLoc seq5 len cStrand)
     = (pos >= seq5) && (pos < seq5 + len) && (cStrand == pStrand)
 
--- | For a pair of 'ContigLoc' regions on the same sequence, indicates
---   if they overlap at all.
+-- | Returns @True@ when two sequence locations overlap at any
+-- position.
 overlaps :: ContigLoc -> ContigLoc -> Bool
 overlaps contig1 contig2
     = case (bounds contig1, bounds contig2) of
         ((low1, high1),(low2, high2)) -> (strand contig1 == strand contig2)
                                          && (low1 <= high2) && (low2 <= high1)
 
+-- | Display a human-friendly, zero-based representation of a sequence location.
 display :: ContigLoc -> String
 display cloc = show (Pos.offset $ startPos cloc) ++ "to" ++ show (Pos.offset $ endPos cloc)
diff --git a/Bio/Location/LocMap.hs b/Bio/Location/LocMap.hs
--- a/Bio/Location/LocMap.hs
+++ b/Bio/Location/LocMap.hs
@@ -1,12 +1,28 @@
--- | A structure to allow fast lookup of objects whose sequence
---   location lines up with a give position.
-module Bio.Location.LocMap (LocMap, mkLocMap, defaultZonesize, fromList
-                           , lookupWithin, lookupOverlaps
-                           , delete, deleteBy, insert
-                           , checkInvariants
-                           )
-    where 
+{-| Efficient lookup of sequence positions and locations in a large
+map of target locations.  For example, target locations might
+represent a collection of genes annotated on a chromosome.  The
+'LocMap' would efficiently find which gene(s) overlapped a sequence
+position on that chromosome.
 
+Target locations are assigned to one or more zones based on
+'Loc.bounds'.  Query locations are then tested only against the target
+locations in the relevant zones.
+-}
+
+module Bio.Location.LocMap (
+  -- * Constructing location lookup maps
+  LocMap, fromList
+
+  -- * Searching for target locations
+  , lookupWithin, lookupOverlaps
+
+  -- * Modifying location lookup maps
+  , delete, deleteBy, insert
+
+  -- * Utilities
+  , checkInvariants
+  ) where 
+
 import qualified Data.IntMap as IM
 import qualified Data.IntSet as IS
 import Data.List hiding (insert, delete, deleteBy)
@@ -20,23 +36,20 @@
 defaultZonesize :: Offset
 defaultZonesize = 65536
 
--- | Collection mapping a collection of 'Loc' locations, possibly
---   overlapping, binned for efficient lookup by position.
+-- | Data structure allowing efficient lookup of target sequence
+-- locations that overlap a query location.  Target locations can be
+-- paired with an arbitrary object.
 data LocMap a = LocMap !Offset !(IM.IntMap (Loc.Loc, a)) !(IM.IntMap IS.IntSet)
 
 instance Monoid (LocMap a) where
-    mempty = mkLocMap defaultZonesize
+    mempty = LocMap defaultZonesize IM.empty IM.empty
     mappend lm0 (LocMap _ keyToLoc1 _) = foldl' (\lm (l,x) -> insert l x lm) lm0 $ IM.elems keyToLoc1
 
--- | Create an empty 'LocMap' with a specified position bin size
-mkLocMap :: Offset -> LocMap a
-mkLocMap zonesize = LocMap zonesize IM.empty IM.empty
-
--- | Create a 'LocMap' from an associated list.
+-- | Create a 'LocMap' from an association list of target locations.
 fromList :: Offset -> [(Loc.Loc, a)] -> LocMap a
-fromList zonesize = foldl' (\lm0 (l,x) -> insert l x lm0) (mkLocMap zonesize)
+fromList zonesize = foldl' (\lm0 (l,x) -> insert l x lm0) mempty
 
--- | Add an object with an associated 'Loc' sequence region
+-- | Insert a new target association into a target location map.
 insert :: Loc.Loc -> a -> LocMap a -> LocMap a
 insert l x (LocMap zonesize keyToLoc zoneKeys) 
     = let !key = case IM.maxViewWithKey keyToLoc of
@@ -50,28 +63,29 @@
                                                 Just currZoneKeySet -> IM.insert z (IS.insert key currZoneKeySet) currZoneKeys
           duplicateError k (l1, _) (l2, _) = error $ "LocMap: key " ++ show k ++ " duplicated: " ++ show (l1, l2)
 
--- | Find the (possibly empty) list of sequence regions and associated
---   objects that contain a 'Pos' position, in the sense of 'withinLoc'
+-- | Find the (possibly empty) list of target locations and associated
+-- objects that contain a sequence position, in the sense of
+-- 'Loc.isWithin'
 lookupWithin :: Pos.Pos -> LocMap a -> [(Loc.Loc, a)]
 lookupWithin pos (LocMap zonesize keyToLoc zoneKeys) 
     = let !zoneKeySet = IM.findWithDefault IS.empty (posZone zonesize pos) zoneKeys
       in filter ((Loc.isWithin pos) . fst) $ keySetLocs keyToLoc zoneKeySet
 
--- | Find the (possibly empty) list of sequence regions and associated
---   objects that overlap a 'Loc' region, in the sense of 'overlapsLoc'
+-- | Find the (possibly empty) list of target locations and associated
+-- objects that overlap a sequence location, in the sense of 'Loc.overlaps'
 lookupOverlaps :: Loc.Loc -> LocMap a -> [(Loc.Loc, a)]
 lookupOverlaps loc (LocMap zonesize keyToLoc zoneKeys)
     = let !zonesKeySet = IS.unions $ map (\z -> IM.findWithDefault IS.empty z zoneKeys) $ locZones zonesize loc
       in filter ((Loc.overlaps loc) . fst) $ keySetLocs keyToLoc zonesKeySet
 
--- | Remove a region / object association from the map, if it is
+-- | Remove a target location and object association from the map, if it is
 -- present.  If it is present multiple times, only the first
 -- occurrence will be deleted.
 delete :: (Eq a) => (Loc.Loc, a) -> LocMap a -> LocMap a
 delete target = deleteBy (== target)
 
--- | Remove the first region / object association satisfying a
--- predicate function.
+-- | Generalized version of 'delete' that removes the first target
+-- location / object association that satisfies a predicate function.
 deleteBy :: ((Loc.Loc, a) -> Bool) -> LocMap a -> LocMap a
 deleteBy isTarget lm0@(LocMap zonesize keyToLoc zoneKeys) 
     = case filter (isTarget . snd) $ IM.toList keyToLoc of
diff --git a/Bio/Location/Location.hs b/Bio/Location/Location.hs
--- a/Bio/Location/Location.hs
+++ b/Bio/Location/Location.hs
@@ -1,21 +1,35 @@
--- |Data types for working with locations in a sequence.  Zero-based
---  'Int64' indices are used throughout, to facilitate direct use of
---  indexing functions on 'SeqData'.
+{-| Data type for a more general sequence location consiting of
+potentially disjoint ranges of positions on the sequence.
 
-module Bio.Location.Location ( Loc(..), bounds, length, startPos, endPos
-                             , extend, posInto, posOutof, isWithin, overlaps
-                             , seqData, seqDataPadded
-                             , display
-                             )
-    where 
+Throughout, /sequence position/ refers to a 'Pos.Pos' which includes a
+strand.  An index into a sequence is referred to as an /offset/, and
+is generally of type 'Offset'.
 
+ -}
+
+module Bio.Location.Location ( 
+  -- * Sequence locations
+  Loc(..)
+
+  -- * Locations and positions
+  , bounds, length, startPos, endPos, posInto, posOutof, isWithin, overlaps
+
+  -- * Extracting subsequences
+  , seqData, seqDataPadded
+
+  -- * Transforming locations
+  , extend
+
+  -- * Displaying locations
+  , display
+  ) where 
+
 import Prelude hiding (length)
 
 import Control.Arrow ((***))
 import Control.Monad.Error
 import qualified Data.ByteString.Lazy.Char8 as LBS
 import Data.List (intercalate)
---import Data.Maybe
 
 import qualified Bio.Location.ContigLocation as CLoc
 import qualified Bio.Location.Position as Pos
@@ -23,49 +37,65 @@
 import Bio.Sequence.SeqData
 
 -- | General (disjoint) sequence region consisting of a concatenated
---   set of contiguous regions
+--   set of contiguous regions (see 'CLoc.ContigLoc').
 newtype Loc = Loc [CLoc.ContigLoc] deriving (Eq, Ord, Show)
 
 instance Stranded Loc where
     revCompl (Loc contigs) = Loc $ reverse $ map revCompl contigs
 
--- | Length of the region
+-- | Returns the length of the region
 length :: Loc -> Offset
 length (Loc contigs) = sum $ map CLoc.length contigs
 
--- | The bounds of a 'Loc', consisting of the lowest & highest
---   sequence indices lying within the region.  The first element of
---   the pair will always be lower than the second.
+-- | The bounds of a sequence location.  This is a pair consisting of
+-- the lowest and highest sequence offsets covered by the region.  The
+-- bounds ignore the strand of the sequence location, and the first
+-- element of the pair will always be lower than the second.  Even if
+-- the positions in the location do not run monotonically through the
+-- location, the overall lowest and highest sequence offsets are returned.
 bounds :: Loc -> (Offset, Offset)
 bounds (Loc []) = error "locBounds on zero-contig Loc"
 bounds (Loc contigs) = (minimum *** maximum) $ unzip $ map CLoc.bounds contigs
 
--- | 0-based starting (5' in the region orientation) offset of the
---   region on its sequence.
+-- | Sequence position of the start of the location.  This is the 5'
+-- end on the location strand, which will have a higher offset than
+-- 'endPos' if the location is on the 'RevCompl' strand.
 startPos :: Loc -> Pos.Pos
 startPos (Loc [])      = error "startPos: zero-contig Loc"
 startPos (Loc contigs) = CLoc.startPos $ head contigs
 
--- | 0-based ending (3' in the region orientation) offset of the
---   region on its sequence.
+-- | Sequence position of the end of the location, as described in 'startPos'.
 endPos :: Loc -> Pos.Pos
 endPos (Loc [])      = error "endPos: zero-contig Loc"
 endPos (Loc contigs) = CLoc.endPos $ last contigs
 
--- | Subsequence 'SeqData' for a 'Loc', provided that the region is
---   entirely within the sequence.
+-- | Extract the nucleotide 'SeqData' for the sequence location.  If
+-- any part of the location lies outside the bounds of the sequence,
+-- an error results.
 seqData :: (Error e, MonadError e m) => SeqData -> Loc -> m SeqData
 seqData sequ (Loc contigs)
     = liftM LBS.concat $ mapM (CLoc.seqData sequ) contigs
 
+-- | As 'seqData', extract the nucleotide subsequence for the
+-- location.  Any positions in the location lying outside the bounds
+-- of the sequence are returned as @N@ rather than producing an error.
 seqDataPadded :: SeqData -> Loc -> SeqData
 seqDataPadded sequ (Loc contigs)
     = LBS.concat $ map (CLoc.seqDataPadded sequ) contigs
 
--- | For a 'Pos' and a 'Loc' region on the same sequence, find the
---   corresponding 'Pos' relative to the region, if the 'Pos' is
---   within the region.  If the 'Loc' region has redundant positions
---   for a given sequence position, the first is returned.
+-- | Given a sequence position and a sequence location relative to the
+-- same sequence, compute a new position representing the original
+-- position relative to the subsequence defined by the location.  If
+-- the sequence position lies outside of the sequence location,
+-- @Nothing@ is returned; thus, the offset of the new position will
+-- always be in the range @[0, length cloc - 1]@.
+--
+-- When the sequence positions in the location are not monotonic,
+-- there may be multiple possible posInto solutions.  That is, if the
+-- same outer sequence position is covered by two different contiguous
+-- blocks of the location, then it would have two possible sequence
+-- positions relative to the location. In this case, the position
+-- 5\'-most in the location orientation is returned.
 posInto :: Pos.Pos -> Loc -> Maybe Pos.Pos
 posInto seqpos (Loc contigs) = posIntoContigs seqpos contigs
 
@@ -76,9 +106,15 @@
         just@(Just _) -> just
         Nothing -> liftM (flip Pos.slide len) $ posIntoContigs seqpos rest
 
--- | For a 'Loc' region on a sequence and a 'Pos' relative to the
---   region, find the corresponding 'Pos' on the sequence, provided
---   that the position is within the bounds of the region.
+-- | Given a sequence location and a sequence position within that
+-- location, compute a new position representing the original position
+-- relative to the outer sequence.  If the sequence position lies
+-- outside the location, @Nothing@ is returned.  
+-- 
+-- This function inverts 'posInto' when the sequence position lies
+-- within the position is actually within the location.  Due to the
+-- possibility of redundant location-relative positions for a given
+-- absolute position, 'posInto' does not necessary invert 'posOutof'
 posOutof :: Pos.Pos -> Loc -> Maybe Pos.Pos
 posOutof pos (Loc contigs) = posOutofContigs pos contigs
 
@@ -89,8 +125,15 @@
         just@(Just _) -> just
         Nothing -> posOutofContigs (Pos.slide seqpos $ negate len) rest
 
--- | Extend a 'Loc' region by incorporating contigous nucleotide
---   regions of the specified lengths on the 5' and 3' ends
+-- | Returns a sequence location produced by extending the original
+-- location on each end, based on a pair of (/5\' extension/, /3\'
+-- extension/).  These add contiguous positions to the 5\' and 3\'
+-- ends of the original location.  The 5\' extension is applied to the
+-- 5\' end of the location on the location strand; if the location is
+-- on the 'RevCompl' strand, the 5\' end will have a higher offset
+-- than the 3\' end and this offset will increase by the amount of the
+-- 5\' extension.  Similarly, the 3\' extension is applied to the 3\'
+-- end of the location.
 extend :: (Offset, Offset) -- ^ (5' extension, 3' extension)
        -> Loc -> Loc
 extend _ (Loc []) = error "extendLoc on zero-contig Loc"
@@ -101,8 +144,8 @@
           extendContigs3 [clast] = [CLoc.extend (0, ext3) clast]
           extendContigs3 (contig:crest) = contig : extendContigs3 crest
 
--- | For a 'Pos' and a 'Loc' on the same sequence, does the position
---   fall within the 'Loc' region?
+-- | Returns @True@ when  a sequence position lies within a sequence
+-- location on the same sequence, and occupies the same strand.
 isWithin :: Pos.Pos -> Loc -> Bool
 isWithin seqpos (Loc contigs) = or $ map (CLoc.isWithin seqpos) contigs
 
@@ -110,10 +153,11 @@
 overlappingContigs (Loc contigs1) (Loc contigs2) 
     = filter (uncurry CLoc.overlaps) [(c1, c2) | c1 <- contigs1, c2 <- contigs2 ]
 
--- | For a pair of 'Loc' regions on the same sequence, do they overlap
---   at all?
+-- | Returns @True@ when two sequence locations overlap at any
+-- position.
 overlaps :: Loc -> Loc -> Bool
 overlaps l1 l2 = not $ null $ overlappingContigs l1 l2
 
+-- | Display a human-friendly, zero-based representation of a sequence location.
 display :: Loc -> String
 display (Loc contigs) = intercalate ";" $ map CLoc.display contigs
diff --git a/Bio/Location/OnSeq.hs b/Bio/Location/OnSeq.hs
--- a/Bio/Location/OnSeq.hs
+++ b/Bio/Location/OnSeq.hs
@@ -1,9 +1,21 @@
-module Bio.Location.OnSeq ( SeqName 
-                          , OnSeq(..), withSeqData, andSameSeq, onSameSeq
-                          , OnSeqs, perSeq, perSeqUpdate, withNameAndSeq
-                          )
-    where 
+{-| Data types for functorially lifting sequence positions and
+locations onto named sequences.  These are useful for taking functions
+that work with sequence positions and locations and associating them
+specific, named sequences.
 
+-}
+
+module Bio.Location.OnSeq ( 
+  -- * Data types
+  SeqName , OnSeq(..)
+
+  -- * Utility functions
+  , withSeqData, andSameSeq, onSameSeq
+
+  -- * Sequence collections indexed by name
+  , OnSeqs, perSeq, perSeqUpdate, withNameAndSeq
+  ) where 
+
 import Control.Monad.Error
 import qualified Data.ByteString.Lazy.Char8 as LBS
 import Data.List
@@ -12,8 +24,10 @@
 
 import Bio.Sequence.SeqData
 
+-- | Sequence name, as in a 'Sequence'
 type SeqName = SeqData
 
+-- | Data type for an object associated with a specific, named sequence
 data OnSeq a = OnSeq { onSeqName :: !SeqName
                      , onSeqObj :: !a
                      } deriving (Eq, Ord, Show)
@@ -21,27 +35,42 @@
 instance Functor OnSeq where
     fmap f (OnSeq seqname x) = OnSeq seqname (f x)
 
-withSeqData :: (Error e, MonadError e m) => (SeqData -> a -> m b) -> (SeqName -> m SeqData) -> OnSeq a -> m b
+-- | Looks up a sequence by name and applies a function to it
+withSeqData :: (Monad m) =>
+               (SeqData -> a -> m b) -- ^ Function using sequence data
+               -> (SeqName -> m SeqData) -- ^ Lookup sequence by name 
+               -> OnSeq a -- ^ Object with named sequence 
+               -> m b
 withSeqData f lookupSeq (OnSeq seqname x) = lookupSeq seqname >>= flip f x
 
+-- | Tests a predicate when two objects are on the same sequence,
+-- returning @False@ if they are on different sequences.
 andSameSeq :: (a -> b -> Bool) -> OnSeq a -> OnSeq b -> Bool
 andSameSeq f (OnSeq xname x) (OnSeq yname y) | xname == yname = f x y
                                              | otherwise = False
 
-onSameSeq :: (Monad m) => (a -> b -> m c) -> OnSeq a -> OnSeq b -> m c
+-- | Performs an action when two objects are on the same sequence and
+-- produces an error otherwise.
+onSameSeq :: (Error e, MonadError e m) => (a -> b -> m c) -> OnSeq a -> OnSeq b -> m c
 onSameSeq f (OnSeq xname x) (OnSeq yname y) 
     | xname == yname = f x y
-    | otherwise      = fail $ "onSameSeq: " ++ show (LBS.unpack xname) ++ " /= " ++ show (LBS.unpack yname)
-
+    | otherwise      = throwError $ strMsg $ "onSameSeq: " ++ show (LBS.unpack xname) ++ " /= " ++ show (LBS.unpack yname)
 
+-- | Data type for a collection of objects indexed by sequence name
 type OnSeqs a = M.Map SeqName a
 
+-- | Lifts a function on an underlying object to look up the sequence
+-- name in a name-indexed collection.
 perSeq :: (Monoid b) => (a -> b -> c) -> OnSeq a -> OnSeqs b -> c
 perSeq f (OnSeq seqname x) = f x . M.findWithDefault mempty seqname
 
+-- | Lifts a function that updates an underlying object to look up the
+-- named sequence and update a named-index collection.
 perSeqUpdate :: (Monoid b) => (a -> b -> b) -> OnSeq a -> OnSeqs b -> OnSeqs b
 perSeqUpdate upd onseq@(OnSeq seqname _) seqmap0 = M.insert seqname (perSeq upd onseq seqmap0) seqmap0
 
+-- | Lifts a function on underlying objects to look up a sequence in a
+-- name-indexed collection
 withNameAndSeq :: (Monad m) => (SeqName -> a -> b -> m c) -> OnSeq a -> OnSeqs b -> m c
 withNameAndSeq f (OnSeq seqname x) = mylookup seqname >=> f seqname x
     where mylookup k = maybe nameNotFound return . M.lookup k
diff --git a/Bio/Location/Position.hs b/Bio/Location/Position.hs
--- a/Bio/Location/Position.hs
+++ b/Bio/Location/Position.hs
@@ -1,8 +1,21 @@
--- | Positions on a sequence.  Zero-based 'Int64' indices are used
--- throughout, to facilitate direct use of indexing functions on
--- 'SeqData'.
+{-| Data type for a sequence position.
 
-module Bio.Location.Position ( Pos(..), slide, seqNt, seqNtPadded, display )
+Zero-based 'Offset' / 'Int64' indices are used throughout, to
+facilitate direct use of indexing functions on 'SeqData'.
+-}
+
+module Bio.Location.Position ( 
+  -- * Sequence positions
+  Pos(..)
+
+  -- * Manipulating positions
+  , slide
+
+  -- * Extracting sequences
+  , seqNt, seqNtPadded
+
+  -- * Displaying positions
+  , display )
     where 
 
 import Control.Monad.Error
@@ -15,25 +28,43 @@
 
 -- | Position in a sequence
 data Pos = Pos { offset :: !Offset -- ^ 0-based index of the position
-               , strand :: !Strand -- ^ Optional strand of the position
+               , strand :: !Strand -- ^ Strand of the position
                }
               deriving (Eq, Ord, Show, Read, Ix)
 
 instance Stranded Pos where
     revCompl (Pos off str) = Pos off $ revCompl str
 
--- | Slide a position by an offset
+-- | Returns a position resulting from sliding the original position
+-- along the sequence by a specified offset.  A positive offset will
+-- move the position away from the 5\' end of the forward stand of the
+-- sequence regardless of the strand of the position itself.  Thus,
+-- 
+-- > slide (revCompl pos) off == revCompl (slide pos off)
+
 slide :: Pos -> Offset -> Pos
 slide (Pos off str) doff = Pos (off + doff) str
 
+
+
+-- | Extract the nucleotide at a specific sequence position.  If the
+-- position lies outside the bounds of the sequence, an error results.
 seqNt :: (Error e, MonadError e m) => SeqData -> Pos -> m Char
 seqNt sequ (Pos off str) | off >= 0 && off < LBS.length sequ = return $ stranded str $ sequ `LBS.index` off
                          | otherwise = throwError $ strMsg $ "position " ++ show off ++ " out of SeqData bounds"
 
+-- | As 'seqNt', extract the nucleotide at a specific sequence
+-- position, but return @N@ when the position lies outside the
+-- bounds of the sequence.
+-- 
+-- > seqNtPadded sequ pos == (either 'N' id . seqNt sequ) pos
 seqNtPadded :: SeqData -> Pos -> Char
 seqNtPadded sequ (Pos off str) | off >= 0 && off < LBS.length sequ = stranded str $ sequ `LBS.index` off
                                | otherwise = 'N'
 
+
+
+-- | Display a human-friendly, zero-based representation of a sequence position.
 display :: Pos -> String
 display (Pos off str) = show off ++ displayStrand
     where displayStrand = case str of
diff --git a/Bio/Location/SeqLocMap.hs b/Bio/Location/SeqLocMap.hs
--- a/Bio/Location/SeqLocMap.hs
+++ b/Bio/Location/SeqLocMap.hs
@@ -1,6 +1,19 @@
-module Bio.Location.SeqLocMap ( SeqLocMap
-                              , empty
-                              , fromList, insert
+{-| 
+
+Efficient lookup of query positions in a collection of target sequence
+locations where positions and locations are associated with specific
+sequence names.  This is an extension of 'LocMap' to use locations and
+positions on named sequences as in 'SeqLocation'.
+
+-}
+
+module Bio.Location.SeqLocMap ( -- * Location lookup maps for named sequences
+                                SeqLocMap, empty, fromList
+
+                              -- * Modifying location lookup maps
+                              , insert
+                              
+                              -- * Searching for target locations
                               , lookupWithin, lookupOverlaps
                               )
     where 
@@ -13,21 +26,36 @@
 import qualified Bio.Location.LocMap as LM
 import qualified Bio.Location.SeqLocation as SeqLoc
 
+-- | A data structure for efficiently finding target sequence
+-- locations ('SeqLoc.Loc') that overlap query positions or locations.
+-- Each target location can be associated with an arbitrary additional
+-- value in the lookup map.
 type SeqLocMap a = OnSeqs (LM.LocMap a)
 
+-- | Empty lookup map.
 empty :: SeqLocMap a
 empty = M.empty
 
+-- | Creates a 'SeqLocMap' from a list of target locations and their
+-- associated objects
 fromList :: [(SeqLoc.SeqLoc, a)] -> SeqLocMap a
 fromList = foldl' (\lm (sl,x) -> insert sl x lm) M.empty
 
+-- | Inserts a new target location and associated object into the
+-- location lookup map.
 insert :: SeqLoc.SeqLoc -> a -> SeqLocMap a -> SeqLocMap a
 insert sloc x = perSeqUpdate (\loc locmap -> LM.insert loc x locmap) sloc
 
+-- | Find the (possibly empty) list of target locations and associated
+-- objects that contain a sequence position, in the sense of
+-- 'Loc.isWithin'.
 lookupWithin :: SeqLoc.SeqPos -> SeqLocMap a -> [(SeqLoc.SeqLoc, a)]
 lookupWithin = withNameAndSeq namedLookupWithin
     where namedLookupWithin seqname pos = map (first $ OnSeq seqname) . LM.lookupWithin pos
 
+-- | Find the (possibly empty) list of target locations and associated
+-- objects that overlap a sequence location, in the sense of
+-- 'Loc.overlaps'.
 lookupOverlaps :: SeqLoc.SeqLoc -> SeqLocMap a -> [(SeqLoc.SeqLoc, a)]
 lookupOverlaps = withNameAndSeq namedLookupOverlaps
     where namedLookupOverlaps seqname loc = map (first $ OnSeq seqname) . LM.lookupOverlaps loc
diff --git a/Bio/Location/SeqLocation.hs b/Bio/Location/SeqLocation.hs
--- a/Bio/Location/SeqLocation.hs
+++ b/Bio/Location/SeqLocation.hs
@@ -1,14 +1,32 @@
-module Bio.Location.SeqLocation ( SeqPos, displaySeqPos
-                                , ContigSeqLoc, withinContigSeqLoc, displayContigSeqLoc 
+{-| 
+
+Data types for sequence locations and sequence positions associated
+with specific, named sequences.
+
+-}
+
+module Bio.Location.SeqLocation ( -- * Positions on named sequences
+                                  SeqPos
+
+                                -- * Contiguous location spans on named sequences
+                                , ContigSeqLoc, withinContigSeqLoc
+
+                                -- * Arbitrary location spans on named sequences
                                 , SeqLoc
-                                , isWithin, overlaps, seqData
-                                , display
+
+                                -- * Testing for location intersection on named sequences
+                                , isWithin, overlaps
+
+                                -- * Extracting subsequences
+                                , seqData
+
+                                -- * Displaying locations on named sequences
+                                , displaySeqPos, displayContigSeqLoc, display
                                 )
     where 
 
 import Control.Monad.Error
 import qualified Data.ByteString.Lazy.Char8 as LBS
---import Data.List
 
 import qualified Bio.Location.ContigLocation as CLoc
 import qualified Bio.Location.Location as Loc
@@ -16,30 +34,50 @@
 import qualified Bio.Location.Position as Pos
 import Bio.Sequence.SeqData
 
+-- | A position on a named sequence
 type SeqPos = OnSeq Pos.Pos
 
+-- | Display a human-friendly representation of a 'SeqPos'
 displaySeqPos :: SeqPos -> String
 displaySeqPos (OnSeq refname pos) = LBS.unpack refname ++ "@" ++ Pos.display pos
 
+-- | A location consisting of a contiguous span of positions on a
+-- named sequence.
 type ContigSeqLoc = OnSeq CLoc.ContigLoc
 
+-- | Test whether a sequence position lies within a sequence location.
+-- This requires that the position lie within the location as per
+-- 'CLoc.isWithin' and have the same sequence name.
 withinContigSeqLoc :: SeqPos -> ContigSeqLoc -> Bool
 withinContigSeqLoc = andSameSeq CLoc.isWithin
 
+-- | Display a human-friendly representation of a 'ContigSeqLoc'
 displayContigSeqLoc :: ContigSeqLoc -> String
 displayContigSeqLoc (OnSeq refname cloc) = LBS.unpack refname ++ "@" ++ CLoc.display cloc
 
+-- | A general location, consisting of spans of sequence positions on
+-- a specific, named sequence.
 type SeqLoc = OnSeq Loc.Loc
 
+-- | Test whether a sequence position lies within a sequence location.
+-- This requires that the position lie within the location as per
+-- 'Loc.isWithin' and have the same sequence name.
 isWithin :: SeqPos -> SeqLoc -> Bool
 isWithin = andSameSeq Loc.isWithin
 
+-- | Test whether two sequence locations overlap in any position.
+-- This requires that the locations overlap as per 'Loc.overlaps' and
+-- have the same sequence name.
 overlaps :: SeqLoc -> SeqLoc -> Bool
 overlaps = andSameSeq Loc.overlaps
 
+-- | Extract the subsequence specified by a sequence location from a
+-- sequence database.  The sequence name is used to retrieve the full
+-- sequence and the subsequence is extracted as by 'Loc.seqData'
 seqData :: (Error e, MonadError e m) => (SeqName -> m SeqData) -> SeqLoc -> m SeqData
 seqData = withSeqData Loc.seqData
 
+-- | Display a human-friendly representation of a 'SeqLoc'
 display :: SeqLoc -> String
 display (OnSeq refname loc) = LBS.unpack refname ++ "@" ++ Loc.display loc
 
diff --git a/Bio/Location/Strand.hs b/Bio/Location/Strand.hs
--- a/Bio/Location/Strand.hs
+++ b/Bio/Location/Strand.hs
@@ -1,3 +1,9 @@
+{-| Utilities for manipulating nucleotide sequences and locations on
+nucleotide sequences that occur on a forward or a reverse-complement
+strand.
+
+-}
+
 module Bio.Location.Strand ( Strand(..)
                            , Stranded(..), stranded
                            )
@@ -11,11 +17,13 @@
 -- | Sequence strand
 data Strand = Fwd | RevCompl deriving (Eq, Ord, Show, Read, Bounded, Enum, Ix)
 
--- | Anything, such as a location or a sequence, which lies on a
--- strand and can thus be reverse complemented.
+-- | A nucleotide sequence or location on a nucleotide sequence that
+--   lies on a specific strand and has an orientation.
 class Stranded s where
     revCompl :: s -> s
 
+-- | Convert the orientation of a 'Stranded' thing based on a
+--   specified 'Strand'
 stranded :: (Stranded s) => Strand -> s -> s
 stranded Fwd      = id
 stranded RevCompl = revCompl
diff --git a/Bio/Sequence.hs b/Bio/Sequence.hs
--- a/Bio/Sequence.hs
+++ b/Bio/Sequence.hs
@@ -7,7 +7,7 @@
 module Bio.Sequence 
     (
     -- * Data structures etc ("Bio.Sequence.SeqData")
-      Sequence(..), Offset, SeqData, Qual, QualData
+      Sequence(..), Unknown, Offset, SeqData, Qual, QualData
     -- ** Accessor functions
     , seqlength, seqlabel, seqheader, seqdata, seqqual, (!)
     , appendHeader, setHeader
@@ -15,9 +15,9 @@
     -- ** Converting to and from String.
     , fromStr, toStr
     -- ** Nucleotide functionality.
-    , compl, revcompl
+    , compl, revcompl, revcompl', Nuc, castToNuc
     -- ** Protein sequence functionality
-    , Amino(..), translate, fromIUPAC, toIUPAC
+    , Amino(..), translate, fromIUPAC, toIUPAC, castToAmino
 
     -- * File formats
     -- ** The Fasta file format ("Bio.Sequence.Fasta")
diff --git a/Bio/Sequence/FastQ.hs b/Bio/Sequence/FastQ.hs
--- a/Bio/Sequence/FastQ.hs
+++ b/Bio/Sequence/FastQ.hs
@@ -15,6 +15,8 @@
 
    Currently, we only support the non-Solexa FastQ, adding\/subtracting 33 for 
    the quality values.
+
+   As far as I know, FastQ is only used for nucleotide sequences, never amino acid.
 -}
 
 module Bio.Sequence.FastQ 
@@ -32,18 +34,18 @@
 
 import Bio.Sequence.SeqData
 
-readFastQ :: FilePath -> IO [Sequence]
+readFastQ :: FilePath -> IO [Sequence Nuc]
 readFastQ f = (go . B.lines) `fmap` B.readFile f 
 
-hReadFastQ :: Handle -> IO [Sequence]
+hReadFastQ :: Handle -> IO [Sequence Nuc]
 hReadFastQ h = (go . B.lines) `fmap` B.hGetContents h
 
-go :: [B.ByteString] -> [Sequence]
+go :: [B.ByteString] -> [Sequence Nuc]
 go = map (either error id) . unfoldr parse
 
 -- | Parse one FastQ entry, suitable for using in 'unfoldr' over
 --   'B.lines' from a file
-parse :: [B.ByteString] -> Maybe (Either String Sequence, [B.ByteString])
+parse :: [B.ByteString] -> Maybe (Either String (Sequence Nuc), [B.ByteString])
 parse (h1:sd:h2:sq:rest) = 
     case (B.uncons h1,B.uncons h2) of
       (Just ('@',h1name), Just ('+',h2name))
@@ -56,14 +58,14 @@
                err = Left $ "Bio.Sequence.FastQ: illegal number of lines in FastQ format: " ++ showStanza
            in Just (err, [])
 
-writeFastQ :: FilePath -> [Sequence] -> IO ()
+writeFastQ :: FilePath -> [Sequence a] -> IO ()
 writeFastQ f = B.writeFile f . B.concat . map unparse
 
-hWriteFastQ :: Handle -> [Sequence] -> IO ()
+hWriteFastQ :: Handle -> [Sequence a] -> IO ()
 hWriteFastQ h = B.hPut h . B.concat . map unparse
 
 -- helper function for writing
-unparse :: Sequence -> B.ByteString
+unparse :: Sequence a -> B.ByteString
 unparse (Seq h sd (Just sq)) = 
     B.unlines [B.cons '@' h, sd, B.cons '+' h, BB.map (+33) sq]
 unparse (Seq h _ Nothing) = error ("Bio.Sequence.FastQ: sequence " ++ show (B.unpack h) 
diff --git a/Bio/Sequence/Fasta.hs b/Bio/Sequence/Fasta.hs
--- a/Bio/Sequence/Fasta.hs
+++ b/Bio/Sequence/Fasta.hs
@@ -5,6 +5,10 @@
    sequence data in the Fasta format.
    Each sequence consists of a header (with a '>' prefix)
    and a set of lines containing the sequence data.
+
+   As Fasta is used for both amino acids and nucleotides, the
+   resulting 'Sequence's are type-tagged with 'Unknown'.  If you know the 
+   type of sequence you are reading, use 'castToAmino' or 'castToNuc'.
 -}
 
 module Bio.Sequence.Fasta
@@ -57,23 +61,23 @@
 -}
 
 -- | Lazily read sequences from a FASTA-formatted file
-readFasta :: FilePath -> IO [Sequence]
+readFasta :: FilePath -> IO [Sequence Unknown]
 readFasta f = (mkSeqs . B.lines) `fmap` B.readFile f
 
 -- | Write sequences to a FASTA-formatted file.
 --   Line length is 60.
-writeFasta :: FilePath -> [Sequence] -> IO ()
+writeFasta :: FilePath -> [Sequence a] -> IO ()
 writeFasta f ss = do
   h <- openFile f WriteMode
   hWriteFasta h ss
   hClose h
 
 -- | Read quality data for sequences to a file.
-readQual :: FilePath -> IO [Sequence]
+readQual :: FilePath -> IO [Sequence Unknown]
 readQual f = (mkQual . B.lines) `fmap` B.readFile f
 
 -- | Write quality data for sequences to a file.
-writeQual :: FilePath -> [Sequence] -> IO ()
+writeQual :: FilePath -> [Sequence a] -> IO ()
 writeQual f ss = do
   h <- openFile f WriteMode
   hWriteQual h ss
@@ -81,7 +85,7 @@
 
 -- | Read sequence and associated quality.  Will error if
 --   the sequences and qualites do not match one-to-one in sequence.
-readFastaQual :: FilePath -> FilePath -> IO [Sequence]
+readFastaQual :: FilePath -> FilePath -> IO [Sequence Unknown]
 readFastaQual  s q = do
   ss <- readFasta s
   qs <- readQual q
@@ -100,7 +104,7 @@
 
 -- | Write sequence and quality data simulatnously
 --   This may be more laziness-friendly.
-writeFastaQual :: FilePath -> FilePath -> [Sequence] -> IO ()
+writeFastaQual :: FilePath -> FilePath -> [Sequence a] -> IO ()
 writeFastaQual f q ss = do
   hf <- openFile f WriteMode
   hq <- openFile q WriteMode
@@ -108,16 +112,16 @@
   hClose hq
   hClose hf
 
-hWriteFastaQual :: Handle -> Handle -> [Sequence] -> IO ()
+hWriteFastaQual :: Handle -> Handle -> [Sequence a] -> IO ()
 hWriteFastaQual hf hq = mapM_ wFQ
     where wFQ s = wFasta hf s >> wQual hq s
 
 -- | Lazily read sequence from handle
-hReadFasta :: Handle -> IO [Sequence]
+hReadFasta :: Handle -> IO [Sequence Unknown]
 hReadFasta h = (mkSeqs . B.lines) `fmap` B.hGetContents h
 
 -- | Write sequences in FASTA format to a handle.
-hWriteFasta :: Handle -> [Sequence] -> IO ()
+hWriteFasta :: Handle -> [Sequence a] -> IO ()
 hWriteFasta h = mapM_ (wFasta h)
 
 wHead :: Handle -> SeqData -> IO ()
@@ -126,17 +130,17 @@
    B.hPut h l
    B.hPut h $ B.pack "\n"
 
-wFasta :: Handle -> Sequence -> IO ()
+wFasta :: Handle -> Sequence a -> IO ()
 wFasta h (Seq l d _) = do
   wHead h l
   let ls = splitsAt 60 d
   mapM_ (B.hPut h) $ intersperse (B.pack "\n") ls
   B.hPut h $ B.pack "\n"
 
-hWriteQual :: Handle -> [Sequence] -> IO ()
+hWriteQual :: Handle -> [Sequence a] -> IO ()
 hWriteQual h = mapM_ (wQual h)
 
-wQual :: Handle -> Sequence -> IO ()
+wQual :: Handle -> Sequence a -> IO ()
 wQual h (Seq l _ (Just q)) = do
   wHead h l
   let qls = splitsAt 20 q
@@ -148,15 +152,15 @@
 --   Blank lines are ignored.
 --   Comment lines start with "#" are allowed between sequences (and ignored).
 --   Lines starting with ">" initiate a new sequence.
-mkSeqs :: [ByteString] -> [Sequence]
+mkSeqs :: [ByteString] -> [Sequence Unknown]
 mkSeqs = map mkSeq . blocks
 
-mkSeq :: [ByteString] -> Sequence
+mkSeq :: [ByteString] -> Sequence Unknown
 mkSeq (l:ls) = Seq (B.drop 1 l) (B.concat $ takeWhile isSeq ls) Nothing
     where isSeq s = (not . B.null) s && ((flip elem) (['A'..'Z']++['a'..'z']) . B.head) s
 mkSeq [] = error "empty input to mkSeq"
 
-mkQual :: [ByteString] -> [Sequence]
+mkQual :: [ByteString] -> [Sequence Unknown]
 mkQual = map f . blocks
     where f (l:ls) = Seq (B.drop 1 l) B.empty
                      (Just $ BB.pack (map int (B.words $ B.unlines ls)))
diff --git a/Bio/Sequence/Phd.hs b/Bio/Sequence/Phd.hs
--- a/Bio/Sequence/Phd.hs
+++ b/Bio/Sequence/Phd.hs
@@ -15,32 +15,33 @@
 import System.IO
 
 -- | Parse a .phd file, extracting the contents as a Sequence
-readPhd :: FilePath -> IO Sequence
+readPhd :: FilePath -> IO (Sequence Nuc)
 readPhd f = return . mkPhd =<< B.readFile f
 
 -- | Parse .phd contents from a handle
-hReadPhd :: Handle -> IO Sequence
+hReadPhd :: Handle -> IO (Sequence Nuc)
 hReadPhd h = return . mkPhd =<< B.hGetContents h
 
 -- | The actual phd parser.
 
 --  aesthetics is not a major design goal...
 --  but error checking really should have been.  Sigh.
-mkPhd :: B.ByteString -> Sequence
+mkPhd :: B.ByteString -> (Sequence Nuc)
 mkPhd inp = 
   let (hd:fs) = filter (not . B.null) . B.lines $ inp
       (comment,sd) = break (==B.pack "BEGIN_DNA") fs
-      label = B.drop 15 hd
+      (magic,label) = B.splitAt 15 hd
+      more_magic = magic == B.pack "BEGIN_SEQUENCE "
       fields = B.words . B.unlines
                . filter (not . isSubstr (B.pack "_COMMENT")) $ comment
       sdata = filter ((==3).length) . map B.words $ sd
       err = error "failed to parse quality value"
       qual = BB.fromChunks [BBB.pack . map (maybe err (fromIntegral . fst) . B.readInt . (!!1)) $ sdata]
-  in qual `seq` (Seq (compact $ B.unwords (label:fields)) 
-      (compact $ B.concat $ map head sdata)
-      (Just qual))
---   Todo: check that we start with "BEGIN_SEQUENCE", and that we
---   have a BEGIN_DNA/END_DNA region there.
+  in if more_magic then qual `seq` (Seq (compact $ B.unwords (label:fields)) 
+                                    (compact $ B.concat $ map head sdata)
+                                    (Just qual))
+     else error "Incorrectly formatted PHD file - missing BEGIN_SEQUENCE"
+     --   Todo: also check that we have a BEGIN_DNA/END_DNA region there.
 
 isSubstr :: B.ByteString -> B.ByteString -> Bool
 isSubstr s = any (B.isPrefixOf s) . B.tails
diff --git a/Bio/Sequence/SFF.hs b/Bio/Sequence/SFF.hs
--- a/Bio/Sequence/SFF.hs
+++ b/Bio/Sequence/SFF.hs
@@ -13,16 +13,19 @@
 
 module Bio.Sequence.SFF ( SFF(..), CommonHeader(..)
                         , ReadHeader(..), ReadBlock(..)
-                        , readSFF, writeSFF
+                        , readSFF, writeSFF, writeSFF'
                         , sffToSequence
                         , test, convert
-                        , Flow, Qual, Index
+                        , Flow, Qual, Index, SeqData, QualData
+                        , ReadName (..), decodeReadName, encodeReadName
                         ) where
 
 import Bio.Sequence.SeqData
+import Bio.Sequence.SFF_name
 
 import Data.Int
 import qualified Data.ByteString.Lazy as LB
+import qualified Data.ByteString.Lazy.Char8 as LBC
 import qualified Data.ByteString as B
 import qualified Data.ByteString.Char8 as BC
 import Data.ByteString (ByteString)
@@ -30,11 +33,12 @@
 
 import Data.List (intersperse)
 import Data.Binary
-import Data.Binary.Get (getByteString)
+import Data.Binary.Get (getByteString,getLazyByteString)
 import qualified Data.Binary.Get as G
-import Data.Binary.Put (putByteString)
+import Data.Binary.Put (putByteString,putLazyByteString)
 
 import Text.Printf (printf)
+import System.IO
 
 -- | The type of flowgram value
 type Flow = Int16
@@ -51,17 +55,47 @@
 readSFF :: FilePath -> IO SFF
 readSFF f = return . decode =<< LB.readFile f
 
-sffToSequence :: SFF -> [Sequence]
+sffToSequence :: SFF -> [Sequence Nuc]
 sffToSequence (SFF ch rs) = map r2s rs
     where r2s r = clip (read_name $ read_header r, bases r, quality r)
           lb x = LB.fromChunks [x]
-          clip (n,s,q) = let (k,s2) = B.splitAt (fromIntegral $ key_length ch) s
-                         in if k==key ch then Seq (lb n) (lb s2) (Just $ lb $ B.drop (fromIntegral $ key_length ch) q)
-                                else error ("Couldn't match key in sequence "++BC.unpack n++" ("++BC.unpack k++" vs. "++BC.unpack (key ch)++")!")
+          clip (n,s,q) = let (k,s2) = LB.splitAt (fromIntegral $ key_length ch) s
+                         in if k==lb (key ch) then Seq (lb n) s2 (Just $ LB.drop (fromIntegral $ key_length ch) q)
+                                else error ("Couldn't match key in sequence "++BC.unpack n++" ("++LBC.unpack k++" vs. "++BC.unpack (key ch)++")!")
 
+-- | Write an 'SFF' to the specified file name
 writeSFF :: FilePath -> SFF -> IO ()
 writeSFF = encodeFile
 
+
+-- | Write an 'SFF' to the specified file name, but go back and
+--   update the read count.  Useful if you want to output a lazy
+--   stream of 'ReadBlock's.
+writeSFF' :: FilePath -> SFF -> IO ()
+writeSFF' f (SFF hs rs) = do
+  h <- openFile f WriteMode
+  LBC.hPut h $ encode hs
+  c <- writeReads h (fromIntegral $ flow_length hs) rs
+  putStrLn ("Count: "++show c)
+  hSeek h AbsoluteSeek 20
+  LBC.hPut h $ encode c
+  hClose h
+
+writeReads :: Handle -> Int -> [ReadBlock] -> IO Int32
+writeReads _ _ [] = return 0
+writeReads h i (r:rs) = do
+  LBC.hPut h $ encode (RBI i r)
+  c <- writeReads h i rs
+  return (c+1)
+
+data RBI = RBI Int ReadBlock
+
+-- | Wrapper for ReadBlocks since they need additional information
+instance Binary RBI where 
+    put (RBI c r) = do
+      putRB c r
+    get = undefined
+      
 -- --------------------------------------------------
 -- | test serialization by output'ing the header and first two reads 
 --   in an SFF, and the same after a decode + encode cycle.
@@ -102,11 +136,12 @@
                                    (do 
                                       rh <- get :: Get ReadHeader
                                       let nb = fromIntegral $ num_bases rh
+                                          nb' = fromIntegral $ num_bases rh
                                           fl = fromIntegral $ flow_length chead
                                       fg <- getByteString (2*fl)
                                       fi <- getByteString nb
-                                      bs <- getByteString nb
-                                      qty <- getByteString nb
+                                      bs <- getLazyByteString nb'
+                                      qty <- getLazyByteString nb'
                                       let l = (fl*2+nb*3) `mod` 8
                                       when (l > 0) (skip (8-l))
                                       return (ReadBlock rh (decodeArray fg) fi bs qty)
@@ -115,7 +150,7 @@
 
     put (SFF hd rds) = do
       put hd
-      mapM_ (putRB (fromIntegral $ flow_length hd)) rds
+      mapM_ (put . RBI (fromIntegral $ flow_length hd)) rds
 
 -- | A ReadBlock can't be an instance of Binary directly, since it depends on
 --   information from the CommonHeader.
@@ -123,9 +158,11 @@
 putRB fl rb = do
   put (read_header rb)
   mapM_ put (flowgram rb)
+  -- ensure that flowgram has correct lenght
+  replicateM (2*(fl-length (flowgram rb))) (put (0::Word8))
   putByteString (flow_index rb)
-  putByteString (bases rb)
-  putByteString (quality rb)
+  putLazyByteString (bases rb)
+  putLazyByteString (quality rb)
   let nb = fromIntegral $ num_bases $ read_header rb
       l = (fl*2+nb*3) `mod` 8
   when (l > 0) (pad (8-l))
@@ -233,7 +270,9 @@
       read_header                :: ReadHeader
     -- The data block
     , flowgram                   :: [Flow]
-    , flow_index, bases, quality :: ByteString
+    , flow_index                 :: ByteString
+    , bases                      :: SeqData
+    , quality                    :: QualData
     }
 
 instance Show ReadBlock where
@@ -241,7 +280,7 @@
         show h ++ unlines (map ("     "++) 
             ["flowgram:\t"++show f
             , "index:\t"++(concat . intersperse " " . map show . B.unpack) i
-            , "bases:\t"++BC.unpack b
-            , "quality:\t"++(concat . intersperse " " . map show . B.unpack) q
+            , "bases:\t"++LBC.unpack b
+            , "quality:\t"++(concat . intersperse " " . map show . LB.unpack) q
             , ""
             ])
diff --git a/Bio/Sequence/SFF_name.hs b/Bio/Sequence/SFF_name.hs
new file mode 100644
--- /dev/null
+++ b/Bio/Sequence/SFF_name.hs
@@ -0,0 +1,86 @@
+module Bio.Sequence.SFF_name where
+
+import qualified Data.ByteString.Char8 as B
+import Data.ByteString.Char8 (ByteString, pack)
+import Data.Array.Unboxed
+import Data.Char (ord)
+
+-- | Read names encode various information, as per this struct.
+data ReadName = ReadName { date :: (Int,Int,Int)
+                         , time :: (Int,Int,Int)
+                         , region :: Int
+                         , x_loc, y_loc :: Int } deriving Show
+
+-- ----------------------------------------------------------
+-- Decoding
+
+decodeReadName :: ByteString -> Maybe ReadName
+decodeReadName b = do t <- decodeDate $ B.take 6 b
+                      r <- fst `fmap` (B.readInt $ B.take 2 $ B.drop 7 b)
+                      l <- decodeLocation $ B.drop 9 b
+                      return $ ReadName { date = (\[y,m,d] -> (y,m,d)) (take 3 t)
+                               , time = (\[hh,mm,ss] -> (hh,mm,ss)) (drop 3 t)
+                               , region = r
+                               , x_loc = fst l, y_loc = snd l }
+
+decodeLocation :: ByteString -> Maybe (Int,Int)
+decodeLocation l = (`divMod` 4096) `fmap` decode36 l
+
+decodeDate :: ByteString -> Maybe [Int]
+decodeDate d    = (fixyear . reverse . (`divMods` [60,60,24,32,13])) =<< decode36 d
+    where fixyear (i:is) = Just (2000+i:is)
+          fixyear []     = Nothing
+
+-- ----------------------------------------------------------
+-- Encoding
+
+encodeReadName :: ReadName -> ByteString
+encodeReadName r =  B.concat [ encodeDate (date r) (time r) 
+                             , encodeRegion (region r)
+                             , encodeLocation (x_loc r) (y_loc r)]
+
+encodeLocation :: Int -> Int -> ByteString
+encodeLocation = undefined
+
+encodeRegion :: Int -> ByteString
+encodeRegion = undefined
+
+encodeDate :: (Int,Int,Int) -> (Int,Int,Int) -> ByteString
+encodeDate = undefined
+
+-- ----------------------------------------------------------
+
+divMods :: Int -> [Int] -> [Int]
+divMods x (i:is) = let (a,b) = x `divMod` i
+                   in b : divMods a is
+divMods x [] = [x]
+
+-- ----------------------------------------------------------
+-- Decoding base36 strings
+
+decode36 :: ByteString -> Maybe Int
+decode36 s = (foldr1 (\a b -> b*36+a) . reverse) `fmap` (mapM decCh . B.unpack $ s)
+
+{-
+decode36' = dec 0
+    where dec i b = case uncons b of Just (c,rest) -> dec (i*36+fromJust (decCh c)) rest
+                                     Nothing       -> i
+          fromJust (Just z) = z
+-}
+
+decCh :: Char -> Maybe Int
+decCh x | x >= 'A' && x <= 'Z' = Just (ord x - ord 'A')
+        | x >= '0' && x <= '9' = Just (26 + ord x - ord '0')
+        | otherwise            = Nothing -- error ("decode36: can't decode "++show x)
+
+encode36 :: Int -> ByteString
+encode36 = pack . map (b36!) . reverse . enc
+    where
+      enc 0 = []
+      enc i = let (a,b) = i `divMod` 36
+              in b : enc a
+
+b36 :: UArray Int Char
+b36 = listArray (0,35) (['A'..'Z']++['0'..'9'])
+
+
diff --git a/Bio/Sequence/SeqData.hs b/Bio/Sequence/SeqData.hs
--- a/Bio/Sequence/SeqData.hs
+++ b/Bio/Sequence/SeqData.hs
@@ -5,10 +5,14 @@
    like @revcompl@, only makes sense for nucleotides.
 -}
 
+{-# LANGUAGE EmptyDataDecls #-}
+
 module Bio.Sequence.SeqData 
     (
       -- * Data structure
-      -- | A sequence is a header, sequence data itself, and optional quality data.
+      -- | A sequence is a header, sequence data itself, and optional quality data. 
+      --   Sequences are type-tagged to identify them as nucleotide, amino acids,
+      --   or unknown type. 
       --   All items are lazy bytestrings.  The Offset type can be used for indexing.
       Sequence(..), Offset, SeqData,
 
@@ -25,25 +29,37 @@
       -- * Converting to and from [Char]
     , fromStr, toStr
 
+      -- * Sequence utilities
+    , defragSeq
+
       -- * Nucleotide functionality
       -- | Nucleotide sequences contain the alphabet [A,C,G,T].
       -- IUPAC specifies an extended nucleotide alphabet with wildcards, but
       -- it is not supported at this point.
-    , compl, revcompl
+    , compl, revcompl, revcompl', Nuc, castToNuc
 
       -- * Protein functionality
       -- | Proteins are chains of amino acids, represented by the IUPAC alphabet.
     , Amino(..), translate, fromIUPAC, toIUPAC -- amino acids
+    , castToAmino
+    
+      -- * Display a nicely formated sequence.
+    , putSeqLn, seqToStr
+
+      -- * Default type for sequences
+    , Unknown
     ) where
 
-import Data.List (unfoldr)
-import Data.Char (toUpper)
+import Data.List (unfoldr, intercalate, isPrefixOf)
+import Data.Char (toUpper, isNumber)
 import Data.Int
 import Data.Word
 import Data.Maybe (fromJust, isJust)
 import qualified Data.ByteString.Lazy.Char8 as B
 import qualified Data.ByteString.Lazy as BB
-
+import qualified Data.ByteString.Char8 as BS
+import qualified Data.ByteString as BBS
+import Text.Printf (printf)
 
 -- | An offset, index, or length of a 'SeqData'
 type Offset  = Int64
@@ -59,9 +75,186 @@
 type QualData = BB.ByteString -- mild abuse
 
 -- | A sequence consists of a header, the sequence data itself, and optional quality data.
-data Sequence = Seq !SeqData !SeqData !(Maybe QualData) -- ^ header and actual sequence
-                deriving (Show,Eq)
+--   The type parameter is a phantom type to separate nucleotide and amino acid sequences
+data Sequence t = Seq !SeqData !SeqData !(Maybe QualData) -- ^ header and actual sequence
+                  deriving Eq
 
+-- | For type tagging sequences (protein sequences use 'Amino' below)
+data Nuc
+data Unknown
+
+-- | Phantom type functionality
+castToNuc :: Sequence a -> Sequence Nuc
+castToNuc (Seq h d mq) = Seq h d mq
+
+castToAmino :: Sequence a -> Sequence Amino
+castToAmino (Seq h d mq) = Seq h d mq
+
+
+-- | A more arranged show instance for Sequences reassembling the display of
+--   the fasta-format
+instance Show (Sequence a) where
+    show s = seqToStr s 6 10 []
+
+-- | Returns a properly formatted and probably highlighted string
+-- | representation of a sequence. Highlighting is done using ANSI-Escape
+-- | sequences.
+seqToStr :: Sequence a -> Int -> Int -> [(Int, Int)] -> [Char]
+seqToStr (Seq header sd _) chunks len parts =
+    let strSeq     = toStr sd
+        strHeader  = toStr header
+        identifier = head (words strHeader)
+        comment    = intercalate " " $ tail (words strHeader)
+
+        numLength  = (length . show . length) strSeq
+        width      = chunks * len + numLength + chunks
+        format     = intercalate "\n" . splits (width)
+    in
+    (line "ID" width)      ++ (format identifier) ++ "\n\n" ++
+    (line "COMMENT" width) ++ (format comment)    ++ "\n\n" ++
+    (line "DATA" width)    ++ (showDNA strSeq len chunks parts)
+ 
+-- | A simple function to display a sequence: we generate the sequence string and
+-- | call putStrLn
+putSeqLn :: Sequence a -> Int -> Int -> [(Int, Int)] -> IO ()
+putSeqLn s chunks len parts = putStrLn $ seqToStr s chunks len parts
+ 
+-- Creates a header line of length n
+line :: [Char] -> Int -> [Char]
+line header width =
+  header ++ " " ++ (take (width - (length header) - 1) (repeat '-')) ++ "\n"
+ 
+-- | Splits a string into parts of size width. The last element can be shorter.
+splits :: Int -> [t] -> [[t]]
+splits _     [] = []
+splits width xs =
+    let (a, b) = splitAt width xs in
+    [a] ++ (splits width b)
+
+-- Highlighting and Dehighlighting ANSI-Escape sequences
+beginHighlight, endHighlight :: String
+beginHighlight = "\ESC[7m"
+endHighlight   = "\ESC[0m"
+
+-- Creates a string representation of the dna-part of a sequence. We can either
+-- highlight specific areas, or, if this is empty, numberize the rows. At the
+-- moment, both at the same time is not possible (see comment in the numberize-
+-- function). If you want numberized output, take a look at the showDNANumber-
+-- function.
+showDNA :: [Char] -> Int -> Int -> [(Int, Int)] -> [Char]
+showDNA xs len chunks highs =
+    let hs  = highlight xs highs
+        ks  = splitsSkip hs len [beginHighlight, endHighlight]
+        ch  = chunkify ks chunks       
+        ns  = numberize ch chunks len
+        ns' = clean ns False
+    in ns'
+    
+-- Add orientation numbers to each line.
+numberize :: [[Char]] -> Int -> Int -> [Char]
+numberize ls chunks len = 
+    let -- The initial numbering of the single elements in lines
+        numbered = zip [0, (chunks * len)..] ls 
+        
+        -- Create the format string with the right alignment, depending on the 
+        -- number of elements in lines
+        numLength   = (show . length . show ) (chunks * len * (length ls))
+        formatStr   = "%" ++ numLength ++ "d "
+        toNum (a,b) = (printf formatStr a) ++ b        
+    in intercalate "\n" (map toNum numbered)
+  
+-- Concatenates each n elements of a list.    
+chunkify :: [[Char]] -> Int -> [[Char]]
+chunkify [] _n = []
+chunkify xs  n =  
+    [intercalate " " (take n xs)] ++ (chunkify (drop n xs) n)
+      
+-- Highlighting should not occur on spaces or over newline borders. We walk
+-- through a generated string and add appropiate stop-highlight and start-
+-- highlight markers.
+--
+-- This is one of the functions which can be definitely optimized by a seasoned
+-- Haskell programmer.
+clean :: [Char] -> Bool -> [Char]
+clean []         _      = []
+
+-- Handling spaces
+clean (' ':xs)   True   = endHighlight ++ " " ++ beginHighlight ++ clean xs True
+clean (' ':xs)   False  = " " ++ clean xs False
+
+-- Newlines and numbering
+clean ('\n':xs)   True   = 
+    let (num, dna) = span (\c -> isNumber c) xs
+    in endHighlight ++ "\n" ++ num ++ beginHighlight ++ clean dna True
+clean ('\n':xs)   False  = "\n" ++ clean xs False
+
+-- Checking for highlighting at all
+clean str@(x:xs) isHigh =
+    -- Special cases: check for the beginning or ending of highlighting 
+    let a = isPrefixOf beginHighlight str
+        b = isPrefixOf endHighlight str
+    in if a then x : clean xs True
+            else if b then x : clean xs False
+                      else x : clean xs isHigh
+
+-- Highlights a range of strings. Each start and end to be highlighted is 
+-- described by a position-tuple (start, end).
+highlight :: [Char] -> [(Int, Int)] -> [Char]
+highlight ss [] = ss 
+highlight ss ((start, end):rest) = 
+    let sl = [beginHighlight, endHighlight]
+        a  = insert ss start   beginHighlight sl
+        b  = insert a  (end+1) endHighlight   sl
+    in highlight b rest
+
+-- Generates lists of length len from xs, skipping sublist in skips while
+-- counting list elements.
+splitsSkip :: (Eq t) => [t] -> Int -> [[t]] -> [[t]]
+splitsSkip [] _   _     = []
+splitsSkip xs len skips =
+    (takeSkip xs len skips) : splitsSkip (dropSkip xs len skips) len skips
+    
+-- Emulates the normal take behaviour, skipping over sublists in skips.
+takeSkip :: (Eq a) => [a] -> Int -> [[a]] -> [a]
+takeSkip []         _ _     = []
+takeSkip _          0 _     = []
+takeSkip str@(x:xs) n skips =
+    let n' = case (testPrefixes str skips) of
+               Nothing  -> n-1
+               Just len -> n + len - 1
+    in x : takeSkip xs n' skips
+
+-- Emulates the normal drop behaviour, skipping over sublists in skips.
+dropSkip :: (Eq a) => [a] -> Int -> [[a]] -> [a]
+dropSkip []         _ _     = []
+dropSkip xs         0 _     = xs
+dropSkip str@(x:xs) n skips =
+    let n' = case (testPrefixes str skips) of
+               Nothing  -> n-1
+               Just len -> n + len - 1
+    in dropSkip xs n' skips
+
+-- Another definition of insert, which provides support for handling skipping
+-- of strings. Inserts elem *before* pos into xs, but skips all strings in
+-- skips while calculating the specified position
+insert :: (Eq a) => [a] -> Int -> [a] -> [[a]] -> [a]
+insert []         _   _    _     = []
+insert xs         0   elem _     = elem ++ xs
+insert str@(x:xs) pos elem skips = 
+    let pos' = case (testPrefixes str skips) of
+               Nothing  -> pos-1
+               Just len -> pos + len - 1
+    in x : insert xs pos' elem skips
+        
+-- Test, if one of the prefixes is the prefix of xs. If yes, return it's length
+-- else Nothing  
+testPrefixes :: (Eq a) => [a] -> [[a]] -> Maybe Int
+testPrefixes _  []     = Nothing
+testPrefixes xs (p:ps) = 
+    if (isPrefixOf p xs)
+        then return (length p)
+        else testPrefixes xs ps
+        
 -- | Convert a String to 'SeqData'
 fromStr :: String -> SeqData
 fromStr = B.pack
@@ -72,47 +265,63 @@
 
 -- | Read the character at the specified position in the sequence.
 {-# INLINE (!) #-}
-(!) :: Sequence -> Offset -> Char
+(!) :: Sequence a -> Offset -> Char
 (!) (Seq _ bs _) = B.index bs
 
 {-# INLINE (?) #-}
-(?) :: Sequence -> Offset -> Qual
+(?) :: Sequence a -> Offset -> Qual
 (?) (Seq _ _ (Just qs)) = BB.index qs
 (?) (Seq _ _ Nothing) = 
     error "Attempting to use 'seqqual' on sequence without associated quality data"
 
 -- | Return sequence length.
-seqlength :: Sequence -> Offset
+seqlength :: Sequence a -> Offset
 seqlength (Seq _ bs _) = B.length bs
 
 -- | Return sequence label (first word of header)
-seqlabel :: Sequence -> SeqData
+seqlabel :: Sequence a -> SeqData
 seqlabel (Seq l _ _) = head (B.words l)
 
 -- | Return full header.
-seqheader :: Sequence -> SeqData
+seqheader :: Sequence a -> SeqData
 seqheader (Seq l _ _) = l
 
 -- | Return the sequence data.
-seqdata :: Sequence -> SeqData
+seqdata :: Sequence a -> SeqData
 seqdata (Seq _ d _) = d
 
 -- | Return the quality data, or error if none exist.  Use hasqual if in doubt.
-seqqual :: Sequence -> QualData
+seqqual :: Sequence a -> QualData
 seqqual (Seq _ _ (Just q)) = q
 seqqual (Seq _ _ Nothing) = 
     error "Attempting to use 'seqqual' on sequence without associated quality data"
 
 -- | Check whether the sequence has associated quality data.
-hasqual :: Sequence -> Bool
+hasqual :: Sequence a -> Bool
 hasqual (Seq _ _ q) = isJust q
 
 -- | Modify the header by appending text, or by replacing
 --   all but the sequence label (i.e. first word).
-appendHeader, setHeader :: Sequence -> String -> Sequence
+appendHeader, setHeader :: Sequence a -> String -> Sequence a
 appendHeader (Seq h d q) t = (Seq (B.append h (B.pack (" "++t))) d q)
 setHeader (Seq h d q) t = (Seq (B.append (head $ B.words h) (B.pack (" "++t))) d q)
 
+-- | Returns a sequence with all internal storage freshly copied and
+-- with sequence and quality data present as a single chunk.  
+--  
+-- By freshly copying internal storage, 'defragSeq' allows garbage
+-- collection of the original data source whence the sequence was
+-- read; otherwise, use of just a short sequence name can cause an
+-- entire sequence file buffer to be retained.
+-- 
+-- By compacting sequence data into a single chunk, 'defragSeq' avoids
+-- linear-time traversal of sequence chunks during random access into
+-- sequence data.
+defragSeq :: Sequence t -> Sequence t
+defragSeq (Seq name sequ qual) = Seq (defragB name) (defragB sequ) (fmap defragBB qual)
+    where defragB = B.fromChunks . (: []) . BS.concat . B.toChunks
+          defragBB = BB.fromChunks . (: []) . BBS.concat . BB.toChunks
+
 ------------------------------------------------------------
 -- Nucleotide (DNA, RNA) specific stuff
 ------------------------------------------------------------
@@ -120,10 +329,14 @@
 -- | Calculate the reverse complement.
 --   This is only relevant for the nucleotide alphabet, 
 --   and it leaves other characters unmodified.  
-revcompl :: Sequence -> Sequence
-revcompl (Seq l bs qs) = Seq l (B.reverse $ B.map compl bs)
+revcompl :: Sequence Nuc -> Sequence Nuc
+revcompl (Seq l bs qs) = Seq l (revcompl' bs)
                         (maybe Nothing (Just . BB.reverse) qs)
 
+-- | Calculate the reverse complent for SeqData only.
+revcompl' :: SeqData -> SeqData
+revcompl' = B.reverse . B.map compl
+
 -- | Complement a single character.  I.e. identify the nucleotide it 
 --   can hybridize with.  Note that for multiple nucleotides, you usually
 --   want the reverse complement (see 'revcompl' for that).
@@ -145,7 +358,7 @@
 -- | Translate a nucleotide sequence into the corresponding protein
 --   sequence.  This works rather blindly, with no attempt to identify ORFs
 --   or otherwise QA the result.
-translate :: Sequence -> Offset -> [Amino]
+translate :: Sequence Nuc -> Offset -> [Amino]
 translate s' o' = unfoldr triples (s',o')
    where triples (s,o) = 
              if o > seqlength s - 3 then Nothing
diff --git a/Bio/Sequence/TwoBit.hs b/Bio/Sequence/TwoBit.hs
--- a/Bio/Sequence/TwoBit.hs
+++ b/Bio/Sequence/TwoBit.hs
@@ -1,4 +1,5 @@
 {- |
+
    This module implements the 2bit format for sequences.
 
    Based on: <http://genome.ucsc.edu/FAQ/FAQformat#format7>
@@ -11,9 +12,15 @@
 
 
 module Bio.Sequence.TwoBit
-   ( decode2Bit, read2Bit, hRead2Bit
+   ( decode2Bit, 
+     read2Bit, 
+     hRead2Bit,
+     encode2Bit,
+     write2Bit,
+     hWrite2Bit
    ) where
 
+
 import Bio.Sequence.SeqData
 import qualified Data.ByteString.Lazy as BB
 import qualified Data.ByteString.Lazy.Char8 as B
@@ -30,13 +37,13 @@
 import Test.QuickCheck hiding (check)    -- QC 1.0
 -- import Test.QuickCheck hiding ((.&.)) -- QC 2.0
 
-type ByteString = B.ByteString
 
 -- constants
 default_magic, default_version :: Word32
 default_magic   = 0x1A412743
 default_version = 0
 
+
 -- binary extras
 check :: Monad m => (a -> Bool) ->  a -> m a
 check p x = if (p x) then return x else fail "check failed"
@@ -51,6 +58,8 @@
 unbytes :: Integral a => [Word8] -> a
 unbytes = Data.List.foldr (\x y -> y*256+x) 0 . map fromIntegral
 
+
+
 -- Conflicts with Bio.Util.TestBase
 -- instance Arbitrary Word8 where
 --    arbitrary = choose (0,255::Int) >>= return . fromIntegral
@@ -58,166 +67,391 @@
 -- prop_bswap :: Word8 -> Word8 -> Word8 -> Word8 -> Bool
 -- prop_bswap a1 a2 a3 a4 = (bswap 4 . decode . BB.pack) [a1,a2,a3,a4] == ((decode . BB.pack) [a4,a3,a2,a1] :: Word32)
 
--- basic data types
+
+
+
+
+-- "in-core" representation of 2Bit data types
+
 data Header = Header { swap :: Bool,  version, count, reserved :: Word32 }
 instance Show Header where show (Header _ v c r) = "H "++show (v,c,r)
 
 instance Binary Header where
     get = do
-      m <- get
-      v <- get >>= check (==default_version)
-      c <- get
-      r <- get
-      let s = if m==default_magic then id else if m==bswap 4 default_magic then bswap 4
-                                               else error "2bit decode: incorrect magic number"
-      return (Header (m/=default_magic) (s v) (s c) (s r))
+       m <- get
+       v <- get >>= check (==default_version)
+       c <- get
+       r <- get
+       let s = if m == default_magic 
+                   then id 
+                   else if m == bswap 4 default_magic 
+                            then bswap 4
+                            else error "2bit decode: incorrect magic number"
+       return (Header (m /= default_magic) (s v) (s c) (s r))
     put (Header m v c r) = do
-      put default_magic
-      put default_version
-      put c
-      put (0 :: Word32)
+       put default_magic
+       put default_version
+       put c
+       put (0 :: Word32)
 
-data Entry  = Entry  { name :: ByteString, offset :: Word32 }
+
+
+data Entry  = Entry  { name :: B.ByteString, offset :: Word32 }
               deriving Show
 
-swapEntry e = e { offset = bswap 4 (offset e) }
+-- Byte swap an Entry's offset!
+swapEntry :: Entry -> Entry
+swapEntry entry = entry { offset = bswap 4 (offset entry) }
 
 instance Binary Entry where
     get = do
-      len <- getWord8
-      name <- replicateM (fromIntegral len) getWord8
-      offset <- get
-      return (Entry (BB.pack name) offset)
-    put (Entry bs offset) = do
-      let l = fromIntegral $ B.length bs :: Word8
-      put l
-      mapM_ put $ BB.unpack bs
-      put offset
+       len <- getWord8
+       name <- replicateM (fromIntegral len) getWord8
+       offset <- get
+       return (Entry (BB.pack name) offset)
+    put (Entry byteString offset) = do
+       let len = fromIntegral $ B.length byteString :: Word8
+       put len  
+       mapM_ put $ BB.unpack byteString
+       put offset
 
+
+
 data Entries = Entries Header [Entry]
-instance Show Entries where show (Entries h es) = unlines (show h : map show es)
 
+instance Show Entries 
+   where show (Entries h es) = unlines (show h : map show es)
+
 instance Binary Entries where
     get = do
-      h <- get
-      es <- replicateM (fromIntegral $ count h) get
-      return (Entries h $ if swap h then map swapEntry es else es)
+       h <- get
+       es <- replicateM (fromIntegral $ count h) get
+       return (Entries h $ if swap h then map swapEntry es else es)
+    put (Entries h es) = do
+       put h
+       put es   
 
+
+
 {-
- dnaSize - number of bases of DNA in the sequence
- nBlockCount - the number of blocks of Ns in the file (representing unknown sequence)
- nBlockStarts - the starting position for each block of Ns
- nBlockSizes - the size of each block of Ns
- maskBlockCount - the number of masked (lower-case) blocks
+
+   Sequence Record definition
+
+ dnaSize         - number of bases of DNA in the sequence
+ nBlockCount     - the number of blocks of Ns in the file (representing unknown sequence)
+ nBlockStarts    - the starting position for each block of Ns
+ nBlockSizes     - the size of each block of Ns
+ maskBlockCount  - the number of masked (lower-case) blocks
  maskBlockStarts - the starting position for each masked block
- maskBlockSizes - the size of each masked block
- packedDna - the DNA packed to two bits per base
+ maskBlockSizes  - the size of each masked block
+ packedDna       - the DNA packed to two bits per base
+
 -}
 
-data SR = SR { dnaSize, nBlockCount :: Word32
-                         , nBlockStarts, nBlockSizes :: [Word32]
-                         , maskBlockCount :: Word32
-                         , maskBlockStarts, maskBlockSizes :: [Word32]
-                         , packedDna :: [Word8]
-                         , reserved2 :: Word32 }
+
+data SR = SR { dnaSize         :: Word32,
+               nBlockCount     :: Word32,
+               nBlockStarts    :: [Word32],
+               nBlockSizes     :: [Word32],
+               maskBlockCount  :: Word32,
+               maskBlockStarts :: [Word32],
+               maskBlockSizes  :: [Word32],
+               packedDna       :: [Word8],
+               reserved2       :: Word32 }
         deriving Show
 
+
+
 -- big- and little-endian variants (what a mess)
 newtype SRBE = SRBE SR deriving Show
 newtype SRLE = SRLE SR deriving Show
 
+
 instance Binary SRBE where
-    get = undefined
+    get = do
+       dz <- get :: Get Word32
+       nc <- get :: Get Word32
+       let n = fromIntegral nc
+       nbs <- replicateM n get
+       nbsz <- replicateM n get
+       mc <- get :: Get Word32
+       let m = fromIntegral mc
+       mbs <- replicateM m get
+       mbsz <- replicateM m get
+       _reserved <- get :: Get Word32 -- !!!! oops?
+       let d = fromIntegral dz
+       pdna <- replicateM ((d+3) `div` 4)  get
+       return (SRBE $ SR dz nc nbs nbsz
+              mc mbs mbsz pdna _reserved)
+    -- should this happen?  Why not just write default format?
+    put (SRBE sr) = do
+       put $ dnaSize sr
+       put $ nBlockCount sr
+       mapM_ put (nBlockStarts sr)
+       mapM_ put (nBlockSizes sr)
+       put $ maskBlockCount sr
+       mapM_ put (maskBlockStarts sr)
+       mapM_ put (maskBlockSizes sr)
+       put (0::Word32)
+       mapM_ put (packedDna sr)
 
+
 instance Binary SRLE where
     get = do
-      dz <- get :: Get Word32
-      nc <- get :: Get Word32
-      let n = fromIntegral $ bswap 4 nc
-      nbs <- replicateM n get
-      nbsz <- replicateM n get
-      mc <- get :: Get Word32
-      let m = fromIntegral $ bswap 4 mc
-      mbs <- replicateM m get
-      mbsz <- replicateM m get
-      _reserved <- get :: Get Word32 -- !!!! oops?
-      let d = fromIntegral $ bswap 4 dz
-      pdna <- replicateM ((d+3) `div` 4)  get
-      return (SRLE $ SR (bswap 4 dz) (bswap 4 nc)
+       dz <- get :: Get Word32
+       nc <- get :: Get Word32
+       let n = fromIntegral $ bswap 4 nc
+       nbs <- replicateM n get
+       nbsz <- replicateM n get
+       mc <- get :: Get Word32
+       let m = fromIntegral $ bswap 4 mc
+       mbs <- replicateM m get
+       mbsz <- replicateM m get
+       _reserved <- get :: Get Word32 -- !!!! oops?
+       let d = fromIntegral $ bswap 4 dz
+       pdna <- replicateM ((d+3) `div` 4)  get
+       return (SRLE $ SR (bswap 4 dz) (bswap 4 nc)
               (map (bswap 4) nbs) (map (bswap 4) nbsz)
               (bswap 4 mc) (map (bswap 4) mbs) (map (bswap 4) mbsz) pdna (bswap 4 _reserved))
     -- should this happen?  Why not just write default format?
     put (SRLE sr) = do
-      put (bswap 4 $ dnaSize sr)
-      put (bswap 4 $ nBlockCount sr)
-      mapM_ (put . bswap 4) (nBlockStarts sr)
-      mapM_ (put . bswap 4) (nBlockSizes sr)
-      put (bswap 4 $ maskBlockCount sr)
-      mapM_ (put . bswap 4) (maskBlockStarts sr)
-      mapM_ (put . bswap 4) (maskBlockSizes sr)
-      put (0::Word32)
-      mapM_ put (packedDna sr)
+       put (bswap 4 $ dnaSize sr)
+       put (bswap 4 $ nBlockCount sr)
+       mapM_ (put . bswap 4) (nBlockStarts sr)
+       mapM_ (put . bswap 4) (nBlockSizes sr)
+       put (bswap 4 $ maskBlockCount sr)
+       mapM_ (put . bswap 4) (maskBlockStarts sr)
+       mapM_ (put . bswap 4) (maskBlockSizes sr)
+       put (0::Word32)
+       mapM_ put (packedDna sr)
 
--- used to convert to/from sequence data in the Sequence data structure
-fromSR :: SR -> ByteString
+
+
+-- Used to convert from sequence data in the Sequence data structure to ByteString
+fromSR :: SR -> B.ByteString
 fromSR sr = B.unfoldr go (0,low,ns,take (fromIntegral $ dnaSize sr) dna)
     where
-    low = combine (maskBlockStarts sr) (maskBlockSizes sr)
-    ns  = combine (nBlockStarts sr) (nBlockSizes sr)
-    combine starts lengths = concatMap (\(p,l) -> [p..p+l-1]) $ zip starts lengths
+       low = combine (maskBlockStarts sr) (maskBlockSizes sr)
+       ns  = combine (nBlockStarts sr) (nBlockSizes sr)
 
-    dna = decodeDNA $ packedDna sr
-    decodeDNA = concatMap (\x -> [shiftR (x .&. 0xC0) 6, shiftR (x .&. 0x30) 4, shiftR (x .&. 0x0C) 2, x .&. 0x03])
-    dec1 x = case x of 0 -> 'T'; 1 -> 'C'; 2 -> 'A'; 3 -> 'G'; _ -> error ("can't decode value "++show x)
+       combine :: (Num t, Enum t) => [t] -> [t] -> [t]
+       combine starts lengths = concatMap (\(p,l) -> [p..p+l-1]) $ zip starts lengths
 
-    go (_,_,_,[]) = Nothing
-    go (pos,(l:ls),(n:ns),(d:ds))
-        | pos == l && pos == n = Just ('n',(pos+1,ls,ns,ds))
-        | pos == l             = Just (toLower (dec1 d),(pos+1,ls,n:ns,ds))
-        |             pos == n = Just ('N',(pos+1,l:ls,ns,ds))
-        | otherwise            = Just (dec1 d, (pos+1,l:ls,n:ns,ds))
-    go (pos,[],n:ns,d:ds)
-        |             pos == n = Just ('N',(pos+1,[],ns,ds))
-        | otherwise            = Just (dec1 d, (pos+1,[],n:ns,ds))
-    go (pos,l:ls,[],d:ds)
-        | pos == l             = Just (toLower (dec1 d),(pos+1,ls,[],ds))
-        | otherwise            = Just (dec1 d, (pos+1,l:ls,[],ds))
-    go (pos,[],[],d:ds)        = Just (dec1 d, (pos+1,[],[],ds))
+       -- Unpack a 2Bit packed DNA sequence into a list of 2 bit values 0-3 
+       dna = decodeDNA $ packedDna sr
+       decodeDNA = concatMap (\x -> [shiftR (x .&. 0xC0) 6, shiftR (x .&. 0x30) 4, shiftR (x .&. 0x0C) 2, x .&. 0x03])
+
+
+       -- Map 2Bit nucleotide encodings to their character equivalents
+       dec1 :: (Num t) => t -> Char
+       dec1 x = case x of 
+                   0 -> 'T'; 
+                   1 -> 'C'; 
+                   2 -> 'A'; 
+                   3 -> 'G'; 
+                   _ -> error ("can't decode value "++show x)
+
+
+       go :: (Num a, Num t) => (a, [a], [a], [t]) -> Maybe (Char, (a, [a], [a], [t]))
+       go (_,_,_,[])              = Nothing
+       go (pos,(l:ls),(n:ns),(d:ds))
+           | pos == l && pos == n = Just ('n',(pos+1,ls,ns,ds))
+           | pos == l             = Just (toLower (dec1 d),(pos+1,ls,n:ns,ds))
+           |             pos == n = Just ('N',(pos+1,l:ls,ns,ds))
+           | otherwise            = Just (dec1 d, (pos+1,l:ls,n:ns,ds))
+       go (pos,[],n:ns,d:ds)
+           |             pos == n = Just ('N',(pos+1,[],ns,ds))
+           | otherwise            = Just (dec1 d, (pos+1,[],n:ns,ds))
+       go (pos,l:ls,[],d:ds)
+           | pos == l             = Just (toLower (dec1 d),(pos+1,ls,[],ds))
+           | otherwise            = Just (dec1 d, (pos+1,l:ls,[],ds))
+       go (pos,[],[],d:ds)        = Just (dec1 d, (pos+1,[],[],ds))
 --    go x = error (show x)
 
+
+
+toSR :: B.ByteString -> SR
+toSR bs = undefined
+
+
+
+splits :: [Int64] -> B.ByteString -> [B.ByteString] 
+splits [] cs = [cs]
+splits (e:es) cs = let (this,rest) = B.splitAt e cs
+                   in this : splits es rest
+
+
+
+
 -- | Parse a (lazy) ByteString as sequences in the 2bit format.
-decode2Bit :: B.ByteString -> [Sequence]
-decode2Bit cs = let (Entries h es) = decode cs :: Entries
+decode2Bit :: B.ByteString -> [Sequence Unknown]
+decode2Bit cs = let 
+                    -- decode to (Header, [Entry]) from ByteString
+                    (Entries h es) = decode cs :: Entries
+
+                    -- map to [64-bit offset] from [32-bit offset]!!! 
                     ms = map (fromIntegral . offset) es
+
+                    -- break the ByteString up into chunks of raw 2Bit "sequences"
                     (c:chunks) = zipWith (-) ms (0:ms)
+
+                    --  build raw undecoded [Sequence]!!!
+                    --  chunks :: [Int64]
+                    --          "drop c cs" .. chop off all cruft in the beginning!!!
                     ss = splits chunks $ B.drop c cs
-               in map ($ Nothing) $ zipWith Seq (map name es)
-                      $ map fromSR $ case swap h of
-                                       True -> map (unSRLE.decode) ss
-                                       False -> map (unSRBE.decode) ss
 
-splits [] cs = [cs]
-splits (e:es) cs = let (this,rest) = B.splitAt e cs
-                   in this : splits es rest
 
+                --                              ($ Nothing) :: (Maybe a -> b) -> b
+                --                              zipWith :: (a -> b -> c) -> [a] -> [b] -> [c]
+                --                              Seq :: SeqData -> SeqData -> Maybe QualData -> Sequence
+                --                              map fromSR :: [SR] -> [B.ByteString] 
+                --                              map (unSRLE.decode) :: [B.ByteString] -> [SR]
+                --                              map (unSRBE.decode) :: [B.ByteString] -> [SR]
+                in map ($ Nothing) $ zipWith Seq (map name es) $ map fromSR $ case swap h of
+                                                                                 True -> map (unSRLE.decode) ss
+                                                                                 False -> map (unSRBE.decode) ss
+                   
+
+
+
+type SequenceLabel = SeqData
+type SequenceSize = Offset
+type SequenceData = SeqData 
+type TwoBitData = (SequenceLabel, SequenceSize, SequenceData)
+
+-- | Marshall from neutral representation to the 2Bit ByteString rep
+encode2Bit :: [Sequence a] -> B.ByteString
+encode2Bit ss = let
+
+                   buildHeader :: [Sequence a] -> Header
+                   buildHeader ss = Header {swap = True, 
+                                            version = default_version, 
+                                            count = sequenceListLength ss, 
+                                            reserved = 0}
+
+
+                   -- Build the list of 2Bit Entries 
+                   buildEntries :: [TwoBitData] -> Offset -> [Entry]
+                   buildEntries [] _ = []
+                   buildEntries ((label, length, _):xs) currentOffset = Entry {name=label, 
+                                                                               offset=fromIntegral $ currentOffset} : buildEntries xs (currentOffset+length)   
+
+                   -- Build a 2Bit Sequence Record!!
+                   buildSR :: TwoBitData -> SR
+                   buildSR (_, size, dnaData) = SR {dnaSize = (fromIntegral size), 
+                                                    nBlockCount = 0,
+                                                    nBlockStarts = [],
+                                                    nBlockSizes = [],
+                                                    maskBlockCount = 0,
+                                                    maskBlockStarts = [], 
+                                                    maskBlockSizes = [],
+                                                    packedDna = encodeDNA $ splitWord8 $ explode dnaData, 
+                                                    reserved2 = 0}
+
+
+                   -- Total # of sequences present
+                   sequenceListLength :: [Sequence a] -> Word32
+                   sequenceListLength [] = 0
+                   sequenceListLength (s:ss) = 1 + sequenceListLength ss
+  
+
+                   -- Build a list of vital data to do 2Bit marshalling
+                   sequenceListExtract :: [Sequence a] -> [TwoBitData]
+                   sequenceListExtract ss = map (\seq -> (seqlabel seq, seqlength seq, seqdata seq)) ss
+
+  
+                   -- Map a nucleotide to its respective 2Bit encoding
+                   enc1 :: (Num t) => Char -> t 
+                   enc1 c = case c of 
+                               'T' -> 0; 
+                               'C' -> 1; 
+                               'A' -> 2; 
+                               'G' -> 3;
+
+
+                   -- Take a ByteString of nucleotides in character encoding and map to a list of their corresponding TwoBit encodings
+                   explode :: SeqData -> [Word8]
+                   explode seq = map enc1 $ B.unpack seq 
+
+
+                   -- Split a list of Word8's into 4 character sublists (i.e. quads) 
+                   splitWord8 ::  [Word8] -> [[Word8]]
+                   splitWord8 [] = []
+                   splitWord8 cs = let (this,rest) = splitAt 4 cs
+                                      in this : splitWord8 rest
+
+
+                   -- Build a 2Bit packed DNA sequence
+                   encodeDNA :: [[Word8]] -> [Word8]
+                   encodeDNA [] = []
+                   encodeDNA (w:ws) = (pack2Bit $ getQuad) : encodeDNA ws
+
+                                     where
+
+                                        -- Length of a sequence
+                                        len = length w
+
+
+                                        -- Build a 2Bit encoded Word8 from "exploded" encoding 
+                                        pack2Bit :: (Word8, Word8, Word8, Word8) -> Word8
+                                        pack2Bit (d1, d2, d3, d4) = shiftL (d1) 6 .|. shiftL (d2) 4 .|. shiftL (d3) 2 .|. d4
+
+                                        -- Build a quad, i.e. a 4 tuple of Word8
+                                        getQuad :: (Word8, Word8, Word8, Word8)
+                                        getQuad = case len of
+                                                        1 -> (w !! 0, 0, 0, 0);
+                                                        2 -> (w !! 0, w !! 1, 0, 0);
+                                                        3 -> (w !! 0, w !! 1, w !! 2, 0);
+                                                        4 -> (w !! 0, w !! 1, w !! 2, w !! 3);
+
+                          
+
+               in 
+
+                   -- Serialize/marshall into ByteString representation
+                   B.append (encode (buildHeader ss))
+                      (B.append
+                         (B.concat (map encode (buildEntries (sequenceListExtract ss) 56)))  -- TEMP WNH
+                         (B.concat (map (encode . SRBE) (map buildSR (sequenceListExtract ss)))))      -- TEMP WNH
+                         -- WNH (B.concat (map (encode . SRBE) (map buildSR (sequenceListExtract ss)))))
+
+
+
+unSRBE :: SRBE -> SR
 unSRBE (SRBE x) = x
+
+unSRLE :: SRLE -> SR
 unSRLE (SRLE x) = x
 
-toSR :: ByteString -> SR
-toSR bs = undefined
 
--- | Extract sequences from a file in 2bit format.
-read2Bit  :: FilePath -> IO [Sequence]
+
+
+-- | Read sequences from a file in 2bit format and 
+-- | unmarshall/deserialize into Sequence format.
+read2Bit  :: FilePath -> IO [Sequence Unknown]
 read2Bit f = B.readFile f >>= return . decode2Bit
 
--- | Extract sequences in the 2bit format from a handle.
-hRead2Bit :: Handle   -> IO [Sequence]
+-- | Read sequences from a file handle in the 2bit format and
+-- | unmarshall/deserialze into Sequence format.
+hRead2Bit :: Handle   -> IO [Sequence Unknown]
 hRead2Bit h = B.hGetContents h >>= return . decode2Bit
 
--- | Write sequences to file in the 2bit format.
-write2Bit  :: FilePath -> [Sequence] -> IO ()
-write2Bit = undefined
 
--- | Write sequences to a handle in the 2bit format.
-hWrite2Bit :: Handle   -> [Sequence] -> IO ()
-hWrite2Bit = undefined
+
+
+-- | Marshall/serialize [Sequence] into 2Bit format and write to a file. 
+write2Bit  :: FilePath -> [Sequence a] -> IO ()
+write2Bit f seq = do
+                   let byteString = encode2Bit seq    
+                   B.writeFile f byteString
+                   return ()
+
+-- | Marshall/serialize [Sequence] into 2Bit format and write to a file using handle. 
+hWrite2Bit :: Handle   -> [Sequence a] -> IO ()
+hWrite2Bit h seq = do
+                   let byteString = encode2Bit seq
+                   B.hPut h byteString
+                   return ()
+
+
+
diff --git a/Bio/Util/TestBase.hs b/Bio/Util/TestBase.hs
--- a/Bio/Util/TestBase.hs
+++ b/Bio/Util/TestBase.hs
@@ -25,14 +25,14 @@
 fromQ (Q c) = c
 
 -- | For testing, variable lengths
-newtype EST = E Sequence deriving Show
-newtype ESTq = Eq Sequence deriving Show
-newtype Protein = P Sequence deriving Show
+newtype EST = E (Sequence Nuc) deriving Show
+newtype ESTq = Eq (Sequence Nuc) deriving Show
+newtype Protein = P (Sequence Amino) deriving Show
 
 -- | For benchmarking, fixed lengths
-newtype EST_short = ES Sequence deriving Show
-newtype EST_long  = EL Sequence deriving Show
-newtype EST_set  = ESet [Sequence] deriving Show
+newtype EST_short = ES (Sequence Nuc) deriving Show
+newtype EST_long  = EL (Sequence Nuc) deriving Show
+newtype EST_set  = ESet [Sequence Nuc] deriving Show
 
 -- | Take time (CPU and wall clock) and report it
 time :: String -> IO () -> IO ()
diff --git a/bio.cabal b/bio.cabal
--- a/bio.cabal
+++ b/bio.cabal
@@ -1,5 +1,5 @@
 Name:                bio
-Version:             0.3.5
+Version:             0.4
 License:             LGPL
 License-file:        LICENSE
 Author:              Ketil Malde
@@ -8,15 +8,16 @@
 Category:            Bioinformatics
 Synopsis:            A bioinformatics library
 Description:         This is a collection of data structures and algorithms
-                     I've found useful when building various bioinformatics-related tools
+                     useful for building bioinformatics-related tools
                      and utilities.
                      .
                      Current list of features includes: a Sequence data type supporting
-                     protein and nucleotide sequences and conversion between them, quality
+                     protein and nucleotide sequences and conversion between them.  As of version
+		     0.4, different kinds of sequence have different types.  Support for quality
                      data, reading and writing Fasta formatted files, reading TwoBit and
-                     phd formats.  Rudimentary support for doing alignments - including
-                     dynamic adjustment of scores based on sequence quality - and Blast
-                     output parsing.  Partly implemented single linkage clustering, and
+                     phd formats, and Roche/454 SFF files.  Rudimentary (i.e. unoptimized) support
+		     for doing alignments - including dynamic adjustment of scores based on sequence quality. 
+		     Also Blast output parsing.  Partly implemented single linkage clustering, and
 		     multiple alignment.  Reading Gene Ontology (GO) annotations (GOA) and
 		     definitions\/hierarchy.
 		     .
@@ -25,7 +26,7 @@
 
 Tested-With:         GHC==6.8.2
 Build-Type:          Simple
-Build-Depends:       base>3, QuickCheck<2, binary, tagsoup>=0.4, bytestring >= 0.9.1,
+Build-Depends:       base>=3 && <4, QuickCheck<2, binary, tagsoup>=0.4, bytestring >= 0.9.1,
                      containers, array, parallel, parsec, random, old-time, mtl
 -- add fps for ghc 6.4.2; change imports in Bio/Sequence/TwoBit.hs if you want QC 2
 
@@ -40,12 +41,13 @@
                      Bio.Sequence.GOA,
                      Bio.Sequence.GeneOntology,
 		     Bio.Sequence.KEGG,
-		     Bio.Sequence.SFF,
+		     Bio.Sequence.SFF, Bio.Sequence.SFF_name
                      Bio.Alignment.BlastData, Bio.Alignment.BlastFlat,
                      Bio.Alignment.Blast, Bio.Alignment.BlastXML,
                      Bio.Alignment.AlignData, Bio.Alignment.Matrices,
                      Bio.Alignment.SAlign, Bio.Alignment.AAlign, Bio.Alignment.QAlign
                      Bio.Alignment.Multiple, Bio.Alignment.ACE,
+                     Bio.Alignment.Bowtie,
                      Bio.Alignment.Soap,
                      Bio.Clustering,
                      Bio.Util, Bio.Util.Parsex, Bio.Util.TestBase