| This module contains a collection of functions for calculate hash values from nucleotide words. (Actually, it mostly packs nucleotide words from into integers, so it's a one-to-one relationship.) It is intended to be used for indexing large sequence collections. \begin{code}

    module Bio.Sequence.HashWord where

    import Bio.Sequence.SeqData
    import Data.List (sort)
    import Data.Char (toUpper)
    import Data.Int

    import qualified Data.ByteString.Lazy.Char8 as B

    -- | This is a struct for containing a set of hashing functions
    data HashF k = HF { hash   :: SeqData -> Offset -> Maybe k -- ^ calculates the hash at a given offset in the sequence
                     , hashes  :: SeqData -> [(k,Offset)]      -- ^ calculate all hashes from a sequence, and their indices
                     , ksort :: [k] -> [k]                     -- ^ for sorting hashes
                     }

    -- | Adds a default "hashes" function to a @HashF@, when "hash" is defined.
    genkeys :: HashF k -> HashF k
    genkeys kf = kf { hashes = gkeys }
        where gkeys s = mkkeys (hash kf $ s)  [0..B.length s-1]
              mkkeys f (p:ps) = case f p of Nothing -> mkkeys f ps
                                            Just k -> (k,p) : mkkeys f ps
              mkkeys _ [] = []

\end{code} \subsection{Generating hashes} Contigous hashes. \begin{code}

    -- | Contigous constructs an int\/eger from a contigous k-word.
    contigous :: Integral k => Int -> HashF k
    contigous k' = let
       k = fromIntegral k'

       c_key s i | k+i > B.length s = Nothing
                 | otherwise = let s' = B.take k $ B.drop i s
                               in if B.find isN s' == Nothing then Just (n2k k' s')
                                  else Nothing
       c_keys s = c_keys' 0 s
       c_keys' c s | k > B.length s = []
                   | otherwise = case B.findIndex isN w0 of
                                   Nothing -> let cur = n2k (fromIntegral k) s
                                              in (cur,c) : c_keys'' cur c (B.drop k s)
                                   Just i  -> c_keys' (c+i+1) (B.drop (i+1) s)
                where w0 = B.take k s

       c_keys'' cur p s | B.null s = []
                        | isN (B.head s) = c_keys' (p+k+1) (B.tail s)
                        | otherwise = let new = (cur `mod` 4^(k-1))*4+val (B.head s)
                                      in (new,p+1) : c_keys'' new (p+1) (B.tail s)

                  in HF { hash = c_key
                        , hashes  = c_keys
                        , ksort = Data.List.sort
                        }

\end{code} \begin{code}

    -- | Like 'contigous', but returns the same hash for a word and its reverse complement.
    {-# SPECIALIZE rcontig :: Int -> HashF Int #-}
    rcontig :: Integral k => Int -> HashF k
    rcontig k' = let
       k = fromIntegral k'

       c_key s i | k+i > B.length s = Nothing
                 | B.find isN s' /= Nothing = Nothing
                 | otherwise = Just $ min (n2k k' s') (n2k k' rs')
                 where s' = B.take k $ B.drop i s
                       rs' = B.reverse $ B.map complement s'

       c_keys s = c_keys' 0 s
       c_keys' c s | k > B.length s = []
                   | otherwise = case B.findIndex isN w0 of
                                   Nothing -> let c1 = n2k (fromIntegral k) w0
                                                  c2 = n2k (fromIntegral k) w0r
                                              in (min c1 c2,c) : c_keys'' (c1,c2) c (B.drop k s)
                                   Just i  -> c_keys' (c+i+1) (B.drop (i+1) s)
                where w0  = B.take k s
                      w0r = B.reverse $ B.map complement w0

       c_keys'' (c1,c2) p s | B.null s = []
                            | isN (B.head s) = c_keys' (p+k+1) (B.tail s)
                            | otherwise = let n1 = (c1 `mod` 4^(k-1))*4+val (B.head s)
                                              n2 = c2 `div` 4+4^(k-1)*(val . complement) (B.head s)
                                          in (min n1 n2,p+1) : c_keys'' (n1,n2) (p+1) (B.tail s)

                  in HF { hash = c_key
                        , hashes  = c_keys
                        , ksort = Data.List.sort
                        }

\end{code} \section{454-hashes} 454 sequencing results in relatively few substitution errors, but has problems with mulitplicity of mononucleotides. One option is to ignore multiplicity, and simply hash nucleotide changes. Note that genkeys won't work on this one. \begin{code}

    -- what about Ns?
    compact :: SeqData -> [SeqData]
    compact bs | B.null bs = []
               | otherwise = let (h,t) = B.break (/=B.head bs) bs in h : compact t

    -- | Like @rcontig@, but ignoring monomers (i.e. arbitrarily long runs of a single nucelotide
    --   are treated the same a single nucleotide.
    rcpacked :: Integral k => Int -> HashF k
    rcpacked k' = let
        k = fromIntegral k'

        c_key s i | k > B.length s' = Nothing
                  | B.any isN s'    = Nothing
                  | otherwise = Just $ min (n2k k' s') (n2k k' rs')
                  where s1 = take k' $ compact $ B.drop i s
                        s' = fromStr $ map B.head s1
                        rs' = B.reverse $ B.map complement s'

        c_keys = c_keys' 0 . compact -- :: Integral k => SeqData -> [(k,Offset)]
    --    c_keys' :: Offset -> [SeqData] -> [(k,Offset)]
        c_keys' i ss | k' > length ss = []
                     | any isN $ map B.head $ take k' ss = c_keys' (i+(B.length . head) ss) (tail ss)
                     | otherwise = let s' = fromStr $ map B.head $ take k' ss
                                       key1 = n2k k' s'
                                       key2 = n2k k' (B.reverse $ B.map complement s')
                                   in  (min key1 key2,i) : c_keys'' (key1,key2) (i+(B.length . head) ss) (tail ss)

        c_keys'' (c1,c2) i [] = []
        c_keys'' (c1,c2) i (s:ss)
            | isN $ B.head s = c_keys' (i + (sum . map B.length) (s:take k' ss)) ss
            | otherwise = let s1 = B.head s
                              n1 = (c1 `mod` 4^(k-1))*4+val s1
                              n2 = c2 `div` 4+4^(k-1)*(val . complement) s1
                          in (min n1 n2,i) : c_keys'' (n1,n2) (i+B.length s) ss

                 in HF { hash = c_key
                       , hashes = c_keys
                       , ksort = Data.List.sort
                       }

\end{code} \section{Gapped words} Shape uses a gapped shape specified via a string Note that this will be challenging to get symmetric, since an 'N' character may appear in the mask in only one direction. \begin{code}

    type Shape = String
    gapped :: Integral k => Shape -> HashF k
    gapped = error "gapped"

\end{code} \subsection{Key arithmetic} A hash is an integer representing packed k-words. This is currently limited to the nucleotide alphabet, but could easily(?) be extended to arbitrary alphabets. Ideally, it should also handle wildcards (i.e. inserting a position with multiple hashes) \begin{code}

    isN :: Char -> Bool
    isN = not . flip elem "ACGT" . toUpper

    -- conversion between integers (hashes) and strings
    -- n2k ignores contents beyond k chars
    n2k :: Integral k => Int -> SeqData -> k
    n2k k = n2i' 0 . B.take (fromIntegral k)

    n2i' :: (Num a) => a -> SeqData -> a
    n2i' acc gs = if B.null gs then acc else n2i' (4*acc + val (B.head gs)) (B.tail gs)

    -- inefficient, but not used much in critical code anyway
    k2n :: Integral k => Int -> k -> SeqData
    k2n k = fromStr . k2n' k

    k2n' k i = if k==1 then [unval i]
               else let (q,r) = i `divMod` (4^(k-1)) in unval q : k2n' (k-1) r

    {-
    shiftRight, shiftLeft :: Integral k => Int -> Int -> k -> String -> k
    shiftRight k l pk str = (pk `mod` 4^l)*4^(k-l) + n2k 0 str
    shiftLeft = error "Keys.shiftLeft not implemented"
    -}

    val :: (Num t) => Char -> t
    val g = case toUpper g of {'A' -> 0; 'C' -> 1; 'G' -> 2; 'T' -> 3;
                               _ -> error ("val: illegal character" ++ show g)}

    unval :: (Num a) => a -> Char
    unval i = case i of {0 -> 'A'; 1 -> 'C'; 2 -> 'G'; 3 -> 'T';
                              _ -> error ("unval: illegal value "++show i)}

    complement :: Char -> Char
    complement g = case g of {'A' -> 'T'; 'T' -> 'A'; 'a' -> 't'; 't' -> 'A';
                              'C' -> 'G'; 'G' -> 'C'; 'c' -> 'g'; 'g' -> 'c';
                              _ -> error ("Can't complement "++show g)}

\end{code}