diff --git a/Bio/Alignment/AAlign.hs b/Bio/Alignment/AAlign.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/AAlign.hs @@ -0,0 +1,74 @@ +{- | + Implement alignments\/edit distance with affine gap penalties + + I've seen g = (-10,-1) as the suggested price to pay for a gaps + using BLOSUM62. Good choice as any, I guess. +-} + +module Bio.Alignment.AAlign ( + -- * Smith-Waterman, or locally optimal alignment with affine gaps + local_score, local_align + + -- * Needleman-Wunsch, or globally optimal alignment with affine gaps + , global_score, global_align + + ) where + +import Data.List (partition,maximumBy) +import Bio.Sequence.SeqData +import Bio.Alignment.AlignData + +-- | Minus infinity (or an approximation thereof) +minf :: Num a => a +minf = -100000000 + +-- ------------------------------------------------------------ +-- 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 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 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 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] + g = maximum $ minf:[max (sub+go) (gap+ge) | ((sub,gap),_) <- ids] + in (s,g) + +-- ------------------------------------------------------------ +-- Alignments + +-- maximum and addition for compound values +max' (x,ax) (y,yx) = if x>=y then (x,ax) else (y,yx) +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,Alignment) +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,Alignment) +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 + +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,Alignment),(a,Alignment)) +align_select minf mx (go,ge) cds = + let (reps,ids) = partition (isRepl.snd) cds + s = maximumBy (compare `on` fst) + $ (minf,[]):[max' sub gap `fp` (mx (x,y),e) | ((sub,gap),e@(Repl x y)) <- reps] + g = maximumBy (compare `on` fst) + $ (minf,[]):[max' (sub `fp` (go,e)) (gap `fp` (ge,e)) | ((sub,gap),e) <- ids] + in (s,g) diff --git a/Bio/Alignment/ACE.hs b/Bio/Alignment/ACE.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/ACE.hs @@ -0,0 +1,259 @@ +{- | + Read ACE format assembly files + + These are typically output by sequence assembly tools, + like CAP3 or Phrap. + + Documented in the section labelled \"ACE FILE FORMAT\" at + + + Briefly: each field is a line starting with a two letter code, + in some cases followed by data lines termintated by a blank line. + -- AS contigs reads + -- CO contig_name bases reads segments compl (CAP3: segments=0) + -- sequence + -- BQ + -- base_qualities + -- AF read1 compl padded_start_consensus (negatives meaning?) + -- AF read2 .. + -- BS segments + -- RD read1 bases info_items info_tags (latter two set to 0 by CAP3) + -- sequence + -- QA read1 qual_start qual_end align_start align_end + -- DS (phred header? left empty by CAP3) + -- RD read2 ... +-} + +{-# LANGUAGE CPP #-} + +module Bio.Alignment.ACE (readACE, writeACE, Assembly(..) + ,Gaps,extractGaps, insertGaps, ptest + ) +where + +#include +import Interlude hiding (lines,words,readFile,unwords -- ByteString clashes + ,reads -- Assembly clash + ) + +import Bio.Sequence.SeqData + +import qualified Data.ByteString.Lazy.Char8 as B +import Data.ByteString.Lazy.Char8 (ByteString,words,pack,unpack,readFile,unwords) + +import Text.ParserCombinators.Parsec +import Text.ParserCombinators.Parsec.Pos (newPos) + +import Control.Monad (liftM) -- ,when? +import Data.Char (chr) + +data Dir = Fwd | Rev deriving (Eq,Show) +data Assembly = Asm { contig :: (Sequence,Gaps), reads :: [(Offset,Dir,Sequence,Gaps)] } + deriving Show + +type Gaps = [Offset] +type Str = ByteString + +-- | ACE header lines with parameters +-- The tokenizer (scanner) should convert input into a list of these, +-- which in turn can be parsed by Parsec +data ACE = AS Str Str | CO Str Str Str Str Str + | BQ | AF Str Str Str | BS [Str] + | RD Str Str Str Str + | QA Str Str Str Str Str + | DS [Str] | Other [Str] | Empty + deriving (Eq) + +instance Show ACE where + show (AS x y) = "AS "++uw [x,y] + show (CO a b c d e) = "CO "++uw [a,b,c,d,e] + show BQ = "BQ" + show (AF a b c) = "AF "++uw [a,b,c] + show (RD a b c d) = "RD "++uw [a,b,c,d] + show (QA a b c d e) = "QA "++uw [a,b,c,d,e] + show (DS ss) = "DS "++uw ss + show (Other ss) = uw ss + show Empty = "(blank)" + show _ = "unknown ACE string" + +uw = unpack . unwords + + +-- | The Parsec parser type +type AceParser a = GenParser (SourcePos,ACE) () a + +-- | Parse a single token, primitive parser +parse1 :: (ACE -> Maybe a) -> AceParser a +parse1 p = token sho pos pred + where sho (_,t) = show t + pos (n,_) = n + pred (_,t) = p t + +-- | Test parser p on a list of ACE elements +ptest :: Show a => String -> AceParser a -> [ACE] -> IO () +ptest m p = parseTest p . source m + +-- | Add SourcePoses to a stream of ACEs. +source :: String -> [ACE] -> [(SourcePos,ACE)] +source m = zip (iterate (\sp -> incSourceLine sp 1) (newPos m 1 0)) + +-- | Parse a complete ACE file as a set of assemblies. +ace :: AceParser [[Assembly]] +ace = many1 ace1 + +ace1 :: AceParser [Assembly] +ace1 = do + as + many blank + many (do ctg >>= asm) -- apparently, CAP3 outputs empty assemblies (AS 0 0) + +-- | parse the initial header +as :: AceParser (Int,Int) +as = parse1 (\t -> case t of AS cs rs -> do c <- B.readInt cs + r <- B.readInt rs + return (fst c,fst r) + _ -> Nothing) "AS " + +blank = parse1 (\t -> case t of Empty -> Just () + _ -> Nothing) "empty line" + +-- | Gaps are coded as '*'s, this function removes them, and returns +-- the sequence along with the list of gap positions. +extractGaps :: Str -> (Str,Gaps) +extractGaps str = (B.filter (/='*') str,B.elemIndices '*' str) + +-- todo: faster to lift concat out of the inner loop? +insertGaps :: Char -> (Str,Gaps) -> Str +insertGaps c (str,gaps) = go str B.empty 0 gaps + where go str acc p (next:rest) = let (a,b) = B.splitAt (next-p) str + in go b (B.concat [acc,a,pack [c]]) (next+1) rest + go str acc _ [] = B.append acc str + +-- | parse the contig and quality information (CO, BQ) +ctg :: AceParser (Sequence,Gaps) +ctg = do + name <- co + sd <- sdata + let (sd',gaps) = extractGaps sd + many blank + bq + sq <- qdata + many blank + -- todo: gaps? + return (Seq name sd' (Just sq),gaps) + +co, sdata, qdata :: AceParser Str +co = parse1 (\t -> case t of CO name a b c _compl -> do + _bs <- B.readInt a + _rds <- B.readInt b + _seg <- B.readInt c + return name + _ -> Nothing) "CO name bool" + +sdata = do return . B.concat =<< many1 sdata1 + where sdata1 = parse1 (\t -> case t of Other sd -> Just (unwords sd) + _ -> Nothing) "sequence data" + +qdata = do return . B.concat =<< many1 qdata1 + where qdata1 = parse1 (\t -> case t of Other sd -> liftM (pack . map chr) (readInts sd) + _ -> Nothing) "sequence data" + +-- | Read a list of Ints in the Maybe monad +readInts :: [ByteString] -> Maybe [Int] +readInts [] = Just [] +readInts (x:xs) = do (a,_) <- B.readInt x + as <- readInts xs + return $ (a:as) + +bq :: AceParser () +bq = parse1 (\t -> case t of BQ -> Just (); _ -> Nothing) "BQ" + +-- | Given the CO info, get the AFS'es +asm :: (Sequence,Gaps) -> AceParser Assembly +asm cg = do + many blank + afs <- many1 af + _bss <- many bs + many blank + rds cg afs + +-- | Parse a list of AFS, followed by actual read, and merge them +-- afs :: Sequence -> AceParser [Sequence] -- plus some auxiliary info? +-- todo: better error handling! +af :: AceParser (Str,Dir,Offset) +af = parse1 (\t -> case t of AF a b c -> Just (a,f b,readInt' c) + _ -> Nothing) "AF name compl pad_start" + where f b = case unpack b of "U" -> Fwd; "C" -> Rev + readInt' x = case B.readInt x of Just (a,_) -> fromIntegral a + +bs :: AceParser (Int,Int,Str) +bs = parse1 (\t -> case t of BS [x,y,n] -> Just (readInt' x, readInt' y,n) + _ -> Nothing) "BS x y name" + where readInt' i = case B.readInt i of Just (a,_) -> fromIntegral a + +rds :: (Sequence,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, reads = zipWith f xs r } + +rseq :: AceParser (Sequence,Gaps) +rseq = do + (rn,_len,_,_) <- rd + (s,gaps) <- return . extractGaps =<< sdata + -- when (B.length s == fromIntegral len) $ (fail "Incorrect sequence length!") + -- todo: fix gaps! + many1 blank + qa + ds + many blank + return (Seq rn s Nothing,gaps) -- huh? + +-- | parse each read (RD, QA, DS) +rd :: AceParser (Str,Int,Int,Int) +rd = parse1 (\t -> case t of RD a b c d -> do [x,y,z] <- readInts [b,c,d] + return (a,x,y,z) + _ -> Nothing) "RD " + +qa :: AceParser () +qa = parse1 (\_ -> Just ()) + +ds :: AceParser () +ds = parse1 (\t -> case t of DS _ -> Just (); _ -> Nothing) "DS" + +-- ---------------------------------------------------------- +-- Convert lines into tokens +tokenize :: ByteString -> [ACE] +tokenize = map tokenize1 . B.lines + +-- Tokenise a single line +-- todo: error on incorrect (partial) format, error reports with line number +tokenize1 :: ByteString -> ACE +tokenize1 l | B.null l = Empty + | otherwise = let (h:ws) = words l + in case (unpack h,ws) of + ("AS",[cs,rs]) -> AS cs rs + ("CO",[nm,bs,rs,segs,compl]) -> CO nm bs rs segs compl + ("BQ",[]) -> BQ + ("AF",[a,b,c]) -> AF a b c + ("BS",_) -> BS ws + ("RD",[a,b,c,d]) -> RD a b c d + ("QA",[a,b,c,d,e]) -> QA a b c d e + ("DS",_) -> DS ws + _ -> Other (h:ws) + +-- | Reading an ACE file. +readACE :: FilePath -> IO [[Assembly]] +readACE f = parseit =<< B.readFile f + where parseit = \s -> case (parse ace f . source f . tokenize) s of + Left e -> fail (show e) + Right a -> return a + +-- formatError msg = error ("readACE: incorrect format in "++msg) + +writeACE :: FilePath -> [Assembly] -> IO () +writeACE = undefined + +-- todo: hWrite etc + diff --git a/Bio/Alignment/AlignData.hs b/Bio/Alignment/AlignData.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/AlignData.hs @@ -0,0 +1,80 @@ +-- | Data structures and helper functions for calculating alignments + +module Bio.Alignment.AlignData ( + -- * Data types for alignments + Edit(..), Alignment, SubstMx, Selector, Chr + -- * Helper functions + , columns, eval, isRepl, on + , showalign, toStrings + ) where + +import qualified Data.ByteString.Lazy as B +import Bio.Sequence.SeqData +import Data.List (unfoldr) +import Data.Word +import Data.Char (chr) + +-- Q&D helper function +showalign a = let (s1,s2) = toStrings a in s1++"\n"++s2 + +-- | turn an alignment into sequences with '-' representing gaps +-- (for checking, filtering out the '-' characters should return +-- the original sequences, provided '-' isn't part of the sequence +-- alphabet) +toStrings :: Alignment -> (String,String) +toStrings [] = ("","") +toStrings (x:xs) = let (a1',a2') = toStrings xs + chr' = chr . fromIntegral + in case x of Ins c -> ('-':a1', chr' c:a2') + Del c -> (chr' c:a1', '-':a2') + Repl c1 c2 -> (chr' c1:a1', chr' c2:a2') + +-- | The sequence element type, used in alignments. +type Chr = Word8 + +-- | An Edit is either the insertion, the deletion, +-- or the replacement of a character. +data Edit = Ins Chr | Del Chr | Repl Chr Chr deriving (Show,Eq) + +-- | An alignment is a sequence of edits. +type Alignment = [Edit] + +-- | True if the Edit is a Repl. +isRepl :: Edit -> Bool +isRepl (Repl _ _) = True +isRepl _ = False + +-- | A substitution matrix gives scores for replacing a character with another. +-- Typically, it will be symmetric. +type SubstMx a = (Chr,Chr) -> a + +-- | Evaluate an Edit based on SubstMx and gap penalty +eval :: SubstMx 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 +-- that chooses a possible Edit operation, and generates an updated result. +type Selector a = [(a,Edit)] -> a + +-- ------------------------------------------------------------ +-- | 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 f z (Seq _ s1 _) (Seq _ s2 _) = columns' f z s1 s2 + +columns' :: Selector a -> a -> SeqData -> SeqData -> [[a]] +columns' f zero s1 s2 = let + -- the first column consists of increasing numbers of insertions + c0 = zero : map (f.return) (zip c0 (map Ins (B.unpack s2))) + -- given the previous column, and the remains of s2, calculate the next column + mkcol (p0:prev,x) = if B.null x then Nothing + else let xi = B.head x + ys = B.unpack s2 + c = f [(p0,Del xi)] : [f [del,ins,rep] | del <- zip prev $ repeat (Del xi) + | ins <- zip c $ map Ins ys + | rep <- zip (p0:prev) $ map (Repl xi) ys] + in Just (c,(c,B.tail x)) + in c0 : unfoldr mkcol (c0,s1) + +on c f x y = c (f x) (f y) diff --git a/Bio/Alignment/Blast.hs b/Bio/Alignment/Blast.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/Blast.hs @@ -0,0 +1,105 @@ +{- | + This module implements a parser for BLAST results. + + This module is DEPRECATED. + It is *very* recommended that you run blast with XML output + instaed, and use the BlastXML module to parse it. + Don't say I didn't warn you! + + BLAST is a tool for searching in (biological) sequences for + similarity. This library is tested against NCBI-blast version + 2.2.14. There exist several independent versions, so expect some + incompatbilities if you're using a different BLAST version. + + The format is straightforward (and non-recursive), and this implementation + uses a simple line-based, hierarchical parser. + + For more information on BLAST, check + +-} + +{-# LANGUAGE CPP #-} + +module Bio.Alignment.Blast + {-# DEPRECATED "Use XML output and the BlastXML module" #-} +(parse) +where + +import Prelude hiding (lines,words) +import qualified Data.ByteString.Lazy.Char8 as B +import Data.ByteString.Lazy.Char8 (ByteString,isPrefixOf,lines,words) + +import Data.List (groupBy) +import Data.Char (isDigit) + +import Bio.Alignment.BlastData +import Bio.Util (splitWhen) + +-- String constant used in parsing +str_query = B.pack "Query= " +str_gt = B.pack ">" +str_score = B.pack " Score = " +str_refer = B.pack "Reference: " +str_datab = B.pack "Database: " +str_search= B.pack "Searching" + +#define BUG(C_,M_) (error ("Program error - '"++C_++"' failed: "++M_++". Location: "++__FILE__++" line: "++ show __LINE__)) +#define readX (\s -> case [ x | (x,t) <- reads (B.unpack s), ("","") <- lex t] of { [x] -> x ; [] -> BUG("read",("no parse for "++B.unpack s)); _ -> BUG("read","ambigous parse")}) +-------------------------------------------------- +-- Splitting up the input in its relevant parts -- +-------------------------------------------------- + +queries :: [ByteString] -> [[ByteString]] +queries = splitWhen (isPrefixOf str_query) + +qhits :: [ByteString] -> [[ByteString]] +qhits = splitWhen (isPrefixOf str_gt) + +hmatches :: [ByteString] -> [[ByteString]] +hmatches = splitWhen (isPrefixOf str_score) + +-------------------------------------------------- +-- parsing each part -- +-------------------------------------------------- + +-- top level parsing function +parse :: ByteString -> BlastResult +parse s = let (p:qs) = queries $ lines s + br = parse_preamble p + rs = map parse_query qs + in br { results = rs } + + +-- parse metadata from the preamble +parse_preamble (l:ls) = BlastResult { blastprogram = w1, blastversion = w2, blastdate = w3 + , blastreferences = B.concat $ map (B.drop (B.length str_refer)) rs + , database = B.drop (B.length str_datab) d + , dbsequences = -1, dbchars = -1, results = [] } + where [w1,w2,w3] = words l + records = map B.unlines $ splitWhen (\p -> or $ map (\x -> isPrefixOf x p) [str_refer,str_datab,str_search]) ls + rs = filter (isPrefixOf str_refer) records + d = case filter (isPrefixOf str_datab) records of [d1] -> d1; _ -> B.pack "" + +-- parse the blast record for one query sequence +parse_query ls = let ((s1:s2:_):hs) = qhits ls + in BlastRecord { query = B.drop (B.length str_query) s1 + , qlength = readX $ B.takeWhile isDigit $ (!!1) $ B.split '(' s2 + , hits = map parse_hit hs } + +-- parse a hit against a sequence in the data base +parse_hit ls = let (h:ms) = hmatches ls + in BlastHit { subject = B.concat h, slength = -1, matches = map parse_match ms} + +-- parse a match between a query and a subject +-- the format is a bit variable, it seems +parse_match (l1:l2:l3:_) = let + hdr = words $ B.concat [l1,l2,l3] + [bs,ev,ident,str1,str2] = map (hdr!!) [2,7,10] ++ map (reverse hdr!!) [2,0] -- Strand is last, ignore optional Gap + identX = let [nom,_:denom] = groupBy (const isDigit) $ B.unpack ident in (read nom, read denom) + in BlastMatch { bits = readX bs, + e_val = readX ev, + identity = identX, + aux = Strands (readX str1) (readX str2), + q_from = ee, q_to = ee, h_from = ee, h_to = ee + } + where ee = error "hit coordinates is not supported - use XML format" diff --git a/Bio/Alignment/BlastData.hs b/Bio/Alignment/BlastData.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/BlastData.hs @@ -0,0 +1,61 @@ +{- | + This module implements a hierarchical data structure for BLAST results, + there is an alternative flat structure in the "Bio.Alignment.BlastFlat" module. + + BLAST is a tool for searching in (biological) sequences for + similarity. This library is tested against NCBI-blast version + 2.2.14. There exist several independent versions of BLAST, so expect some + incompatbilities if you're using a different BLAST version. + + For parsing BLAST results, the XML format (blastall -m 7) is by far the most + robust choice, and is implemented in the "Bio.Alignment.BlastXML" module. + + The format is straightforward (and non-recursive). + For more information on BLAST, check + +-} + +module Bio.Alignment.BlastData where + +import Data.ByteString.Lazy.Char8 (ByteString) + +-- ------------------------------------------------------------ + +-- | The sequence id, i.e. the first word of the header field. +type SeqId = ByteString + +-- | The 'Strand' indicates the direction of the match, i.e. the plain sequence or +-- its reverse complement. +data Strand = Plus | Minus deriving (Read,Show,Eq) + +-- | The Aux field in the BLAST output includes match information that depends +-- on the BLAST flavor (blastn, blastx, or blastp). This data structure captures +-- those variations. +data Aux = Strands !Strand !Strand -- ^ blastn + | Frame !Strand !Int -- ^ blastx + deriving (Show,Eq) + +-- | A 'BlastResult' is the root of the hierarchy. +data BlastResult = BlastResult + { blastprogram, blastversion, blastdate :: !ByteString + , blastreferences :: !ByteString + , database :: !ByteString + , dbsequences, dbchars :: !Integer + , results :: [BlastRecord] } + deriving Show + +-- | Each query sequence generates a 'BlastRecord' +data BlastRecord = BlastRecord { query :: !SeqId, qlength :: !Int + , hits :: [BlastHit] } deriving Show + +-- | Each match between a query and a target sequence (or subject) +-- is a 'BlastHit'. +data BlastHit = BlastHit { subject :: !SeqId, slength :: !Int + , matches :: [BlastMatch] } deriving Show +-- | A 'BlastHit' may contain multiple separate matches (typcially when +-- an indel causes a frameshift that blastx is unable to bridge). +data BlastMatch = BlastMatch { bits :: !Double, e_val :: !Double + , identity :: (Int,Int) + , q_from, q_to, h_from, h_to :: !Int + , aux :: !Aux } deriving Show + diff --git a/Bio/Alignment/BlastFlat.hs b/Bio/Alignment/BlastFlat.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/BlastFlat.hs @@ -0,0 +1,55 @@ +{- | + This module implements a \"flattened\" data structure for Blast hits, + as opposed to the hierarchical structure in "Bio.Alignment.BlastData". + + The flat data type is useful in many cases where it is more natural + to see the result as a set of rows (e.g. for insertaion in a database). + + It would probably be more (memory-) efficient to go the other way + (i.e. from flat to hierarchical), as passing the current, partially + built "BlastFlat" object down the stream of results and stamping + out a stream of completed ones. (See "Bio.Alignment.BlastXML.breaks" + for this week's most cumbersome use of parallelism to avoid the + memory issue.) +-} + +module Bio.Alignment.BlastFlat + ( + -- * The BlastFlat data type + BlastFlat(..) + -- * Convert from hierarchical to flat structure + , flatten + -- * Re-exports from the hierarchical module ("Bio.Alignment.BlastData") + , B.BlastRecord + , B.blastprogram, B.blastversion, B.blastdate, B.blastreferences + , B.database, B.dbsequences, B.dbchars, B.results + , B.Aux(..), B.Strand(..) + )where + +import qualified Bio.Alignment.BlastData as B +import Data.ByteString.Lazy.Char8 (empty) + +-- | The BlastFlat data structure contains information about a single match +data BlastFlat = BlastFlat { + query :: !B.SeqId, qlength :: !Int -- BlastRecord + , subject :: !B.SeqId, slength :: !Int -- BlastHit + , bits :: !Double, e_val :: !Double -- BlastMatch + , identity :: (Int,Int) + , q_from, q_to, h_from, h_to :: !Int + , aux :: !B.Aux + } + +-- | Convert BlastRecords into BlastFlats (representing a depth-first traversal of the +-- BlastRecord structure.) +flatten :: [B.BlastRecord] -> [BlastFlat] +flatten = concatMap frecord + where frecord r = + concatMap (fhit (bf0 { query = B.query r, qlength = B.qlength r })) $ B.hits r + fhit f h = + map (fmatch f { subject = B.subject h, slength = B.slength h }) $ B.matches h + fmatch f m = + f { bits = B.bits m, e_val = B.e_val m, identity = B.identity m + , q_from = B.q_from m, q_to = B.q_to m + , h_from = B.h_from m, h_to = B.h_to m, aux = B.aux m} + bf0 = BlastFlat e 0 e 0 0 0 (0,0) 0 0 0 0 (B.Frame B.Plus 0) + e = empty diff --git a/Bio/Alignment/BlastXML.hs b/Bio/Alignment/BlastXML.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/BlastXML.hs @@ -0,0 +1,90 @@ +-- Parse blast XML output using tagsoup + +module Bio.Alignment.BlastXML where + +import Bio.Alignment.BlastData + +import qualified Data.ByteString.Lazy.Char8 as B +import Text.HTML.TagSoup +import Control.Monad + +import Control.Parallel + +readXML :: FilePath -> IO [BlastResult] +readXML fp = do + ts <- return . parseTags =<< readFile fp + let (h:iters) = breaks (\t -> isTagOpenName "Iteration" t || isTagOpenName "Hit" t) ts + return [xml2br h iters] + +-- | breaks p = groupBy (const (not.p)) +breaks :: (a -> Bool) -> [a] -> [[a]] +breaks p (x:xs) = let first = x : takeWhile (not.p) xs + rest = dropWhile (not.p) xs + in rest `par` first : if null rest then [] else breaks p rest +breaks _ [] = [] + +getFrom list tag = let xs = sections (isTagOpenName tag) list + in if null xs || null (head xs) || (null . drop 1 . head) xs + then error ("Couldn't find tag '"++show tag++"' in\n"++showSome list) + else case xs !! 0 !! 1 of + TagText s -> s + x -> error ("Unexpeced tag: "++ show x) + +-- Use pattern match since 'length' is strict, defeating the purpose. +showSome a@(_:_:_:_:_:_:_) = (init . show . take 5 $ a)++" ... ]" +showSome a = show a + +xml2br :: [Tag] -> [[Tag]] -> BlastResult +xml2br h is = BlastResult { blastprogram = get "BlastOutput_program" + , blastversion = bv + , blastdate = bd + , blastreferences = get "BlastOutput_reference" + , database = get "BlastOutput_db" + , dbsequences = 0 + , dbchars = 0 + , results = map iter2rec $ breaks (isTagOpenName "Iteration" . head) is + } + where (bv,bd) = B.break (=='[') $ get "BlastOutput_version" + get = B.pack . getFrom h + +iter2rec :: [[Tag]] -> BlastRecord +iter2rec (i:hs) = BlastRecord + { query = B.pack $ get "Iteration_query-def" + , qlength = read $ get "Iteration_query-len" + , hits = map hit2hit hs + } + where get = getFrom i +iter2rec [] = error "iter2rec: got empty list of sections!" + +hit2hit :: [Tag] -> BlastHit +hit2hit hs = BlastHit + { subject = B.pack $ get "Hit_def" + , slength = read $ get "Hit_len" + , matches = map hsp2match $ partitions (isTagOpenName "Hsp") hs + } + where get = getFrom hs + +hsp2match :: [Tag] -> BlastMatch +hsp2match ms = BlastMatch + { bits = read $ get "Hsp_bit-score" + , e_val = read $ get "Hsp_evalue" + , q_from = read $ get "Hsp_query-from" + , q_to = read $ get "Hsp_query-to" + , h_from = read $ get "Hsp_hit-from" + , h_to = read $ get "Hsp_hit-to" + , identity = (read $ get "Hsp_identity", read $ get "Hsp_align-len") + -- blastx has query-frame, tblastn has hit-frame, blastn has both hit and query + , aux = case sections (isTagOpenName "Hsp_hit-frame") ms of + [] -> mkFrame $ get "Hsp_query-frame" + [(_o:TagText hf:_c)] -> case sections (isTagOpenName "Hsp_query-frame") ms of + [] -> mkFrame hf + [(_o:TagText qf:_c)] -> mkStrands hf qf + e -> error ("hsp2match: failed to determine frame:\n"++show e) + } + where get = getFrom ms + mkFrame f = Frame (strand' $ signum $ read f) (abs $ read f) + mkStrands h q = Strands (strand' $ read h) (strand' $ read q) + strand' s = case s of 1 -> Plus; -1 -> Minus + _ -> error ("Strand must be +1 or -1, but was"++show s) + + diff --git a/Bio/Alignment/Matrices.hs b/Bio/Alignment/Matrices.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/Matrices.hs @@ -0,0 +1,724 @@ +{- | + Common substitution matrices for alignments. + + When in doubt, use BLOSUM62. + Consult + for some hints on good parameters for nucleotide alignments. +-} + +module Bio.Alignment.Matrices + ( + -- * BLOSUM matrices (Henikoff and Henikoff) + blosum45, blosum62, blosum80 + + -- * PAM matrices (Dayhoff et al) + , pam30, pam70 + + -- * BLASTn defaults, for nucleotide sequences + , blastn_default + + -- * Generic and simple "matrix" generator + , simpleMx + ) where + +import Bio.Alignment.AlignData (Chr) +import qualified Data.Map as M + +-- | The standard BLOSUM45 matrix. +blosum45, blosum62, blosum80, pam30, pam70 :: (Char,Char) -> Int +blosum45 m = M.findWithDefault (-5) m $ M.fromList [(('A','A'),5), + (('A','B'),-1), (('A','C'),-1), (('A','D'),-2), (('A','E'),-1), + (('A','F'),-2), (('A','G'),0), (('A','H'),-2), (('A','I'),-1), + (('A','K'),-1), (('A','L'),-1), (('A','M'),-1), (('A','N'),-1), + (('A','P'),-1), (('A','Q'),-1), (('A','R'),-2), (('A','S'),1), + (('A','T'),0), (('A','V'),0), (('A','W'),-2), (('A','X'),-1), + (('A','Y'),-2), (('A','Z'),-1), (('B','A'),-1), (('B','B'),4), + (('B','C'),-2), (('B','D'),5), (('B','E'),1), (('B','F'),-3), + (('B','G'),-1), (('B','H'),0), (('B','I'),-3), (('B','K'),0), + (('B','L'),-3), (('B','M'),-2), (('B','N'),4), (('B','P'),-2), + (('B','Q'),0), (('B','R'),-1), (('B','S'),0), (('B','T'),0), + (('B','V'),-3), (('B','W'),-4), (('B','X'),-1), (('B','Y'),-2), + (('B','Z'),2), (('C','A'),-1), (('C','B'),-2), (('C','C'),12), + (('C','D'),-3), (('C','E'),-3), (('C','F'),-2), (('C','G'),-3), + (('C','H'),-3), (('C','I'),-3), (('C','K'),-3), (('C','L'),-2), + (('C','M'),-2), (('C','N'),-2), (('C','P'),-4), (('C','Q'),-3), + (('C','R'),-3), (('C','S'),-1), (('C','T'),-1), (('C','V'),-1), + (('C','W'),-5), (('C','X'),-1), (('C','Y'),-3), (('C','Z'),-3), + (('D','A'),-2), (('D','B'),5), (('D','C'),-3), (('D','D'),7), + (('D','E'),2), (('D','F'),-4), (('D','G'),-1), (('D','H'),0), + (('D','I'),-4), (('D','K'),0), (('D','L'),-3), (('D','M'),-3), + (('D','N'),2), (('D','P'),-1), (('D','Q'),0), (('D','R'),-1), + (('D','S'),0), (('D','T'),-1), (('D','V'),-3), (('D','W'),-4), + (('D','X'),-1), (('D','Y'),-2), (('D','Z'),1), (('E','A'),-1), + (('E','B'),1), (('E','C'),-3), (('E','D'),2), (('E','E'),6), + (('E','F'),-3), (('E','G'),-2), (('E','H'),0), (('E','I'),-3), + (('E','K'),1), (('E','L'),-2), (('E','M'),-2), (('E','N'),0), + (('E','P'),0), (('E','Q'),2), (('E','R'),0), (('E','S'),0), + (('E','T'),-1), (('E','V'),-3), (('E','W'),-3), (('E','X'),-1), + (('E','Y'),-2), (('E','Z'),4), (('F','A'),-2), (('F','B'),-3), + (('F','C'),-2), (('F','D'),-4), (('F','E'),-3), (('F','F'),8), + (('F','G'),-3), (('F','H'),-2), (('F','I'),0), (('F','K'),-3), + (('F','L'),1), (('F','M'),0), (('F','N'),-2), (('F','P'),-3), + (('F','Q'),-4), (('F','R'),-2), (('F','S'),-2), (('F','T'),-1), + (('F','V'),0), (('F','W'),1), (('F','X'),-1), (('F','Y'),3), + (('F','Z'),-3), (('G','A'),0), (('G','B'),-1), (('G','C'),-3), + (('G','D'),-1), (('G','E'),-2), (('G','F'),-3), (('G','G'),7), + (('G','H'),-2), (('G','I'),-4), (('G','K'),-2), (('G','L'),-3), + (('G','M'),-2), (('G','N'),0), (('G','P'),-2), (('G','Q'),-2), + (('G','R'),-2), (('G','S'),0), (('G','T'),-2), (('G','V'),-3), + (('G','W'),-2), (('G','X'),-1), (('G','Y'),-3), (('G','Z'),-2), + (('H','A'),-2), (('H','B'),0), (('H','C'),-3), (('H','D'),0), + (('H','E'),0), (('H','F'),-2), (('H','G'),-2), (('H','H'),10), + (('H','I'),-3), (('H','K'),-1), (('H','L'),-2), (('H','M'),0), + (('H','N'),1), (('H','P'),-2), (('H','Q'),1), (('H','R'),0), + (('H','S'),-1), (('H','T'),-2), (('H','V'),-3), (('H','W'),-3), + (('H','X'),-1), (('H','Y'),2), (('H','Z'),0), (('I','A'),-1), + (('I','B'),-3), (('I','C'),-3), (('I','D'),-4), (('I','E'),-3), + (('I','F'),0), (('I','G'),-4), (('I','H'),-3), (('I','I'),5), + (('I','K'),-3), (('I','L'),2), (('I','M'),2), (('I','N'),-2), + (('I','P'),-2), (('I','Q'),-2), (('I','R'),-3), (('I','S'),-2), + (('I','T'),-1), (('I','V'),3), (('I','W'),-2), (('I','X'),-1), + (('I','Y'),0), (('I','Z'),-3), (('K','A'),-1), (('K','B'),0), + (('K','C'),-3), (('K','D'),0), (('K','E'),1), (('K','F'),-3), + (('K','G'),-2), (('K','H'),-1), (('K','I'),-3), (('K','K'),5), + (('K','L'),-3), (('K','M'),-1), (('K','N'),0), (('K','P'),-1), + (('K','Q'),1), (('K','R'),3), (('K','S'),-1), (('K','T'),-1), + (('K','V'),-2), (('K','W'),-2), (('K','X'),-1), (('K','Y'),-1), + (('K','Z'),1), (('L','A'),-1), (('L','B'),-3), (('L','C'),-2), + (('L','D'),-3), (('L','E'),-2), (('L','F'),1), (('L','G'),-3), + (('L','H'),-2), (('L','I'),2), (('L','K'),-3), (('L','L'),5), + (('L','M'),2), (('L','N'),-3), (('L','P'),-3), (('L','Q'),-2), + (('L','R'),-2), (('L','S'),-3), (('L','T'),-1), (('L','V'),1), + (('L','W'),-2), (('L','X'),-1), (('L','Y'),0), (('L','Z'),-2), + (('M','A'),-1), (('M','B'),-2), (('M','C'),-2), (('M','D'),-3), + (('M','E'),-2), (('M','F'),0), (('M','G'),-2), (('M','H'),0), + (('M','I'),2), (('M','K'),-1), (('M','L'),2), (('M','M'),6), + (('M','N'),-2), (('M','P'),-2), (('M','Q'),0), (('M','R'),-1), + (('M','S'),-2), (('M','T'),-1), (('M','V'),1), (('M','W'),-2), + (('M','X'),-1), (('M','Y'),0), (('M','Z'),-1), (('N','A'),-1), + (('N','B'),4), (('N','C'),-2), (('N','D'),2), (('N','E'),0), + (('N','F'),-2), (('N','G'),0), (('N','H'),1), (('N','I'),-2), + (('N','K'),0), (('N','L'),-3), (('N','M'),-2), (('N','N'),6), + (('N','P'),-2), (('N','Q'),0), (('N','R'),0), (('N','S'),1), + (('N','T'),0), (('N','V'),-3), (('N','W'),-4), (('N','X'),-1), + (('N','Y'),-2), (('N','Z'),0), (('P','A'),-1), (('P','B'),-2), + (('P','C'),-4), (('P','D'),-1), (('P','E'),0), (('P','F'),-3), + (('P','G'),-2), (('P','H'),-2), (('P','I'),-2), (('P','K'),-1), + (('P','L'),-3), (('P','M'),-2), (('P','N'),-2), (('P','P'),9), + (('P','Q'),-1), (('P','R'),-2), (('P','S'),-1), (('P','T'),-1), + (('P','V'),-3), (('P','W'),-3), (('P','X'),-1), (('P','Y'),-3), + (('P','Z'),-1), (('Q','A'),-1), (('Q','B'),0), (('Q','C'),-3), + (('Q','D'),0), (('Q','E'),2), (('Q','F'),-4), (('Q','G'),-2), + (('Q','H'),1), (('Q','I'),-2), (('Q','K'),1), (('Q','L'),-2), + (('Q','M'),0), (('Q','N'),0), (('Q','P'),-1), (('Q','Q'),6), + (('Q','R'),1), (('Q','S'),0), (('Q','T'),-1), (('Q','V'),-3), + (('Q','W'),-2), (('Q','X'),-1), (('Q','Y'),-1), (('Q','Z'),4), + (('R','A'),-2), (('R','B'),-1), (('R','C'),-3), (('R','D'),-1), + (('R','E'),0), (('R','F'),-2), (('R','G'),-2), (('R','H'),0), + (('R','I'),-3), (('R','K'),3), (('R','L'),-2), (('R','M'),-1), + (('R','N'),0), (('R','P'),-2), (('R','Q'),1), (('R','R'),7), + (('R','S'),-1), (('R','T'),-1), (('R','V'),-2), (('R','W'),-2), + (('R','X'),-1), (('R','Y'),-1), (('R','Z'),0), (('S','A'),1), + (('S','B'),0), (('S','C'),-1), (('S','D'),0), (('S','E'),0), + (('S','F'),-2), (('S','G'),0), (('S','H'),-1), (('S','I'),-2), + (('S','K'),-1), (('S','L'),-3), (('S','M'),-2), (('S','N'),1), + (('S','P'),-1), (('S','Q'),0), (('S','R'),-1), (('S','S'),4), + (('S','T'),2), (('S','V'),-1), (('S','W'),-4), (('S','X'),-1), + (('S','Y'),-2), (('S','Z'),0), (('T','A'),0), (('T','B'),0), + (('T','C'),-1), (('T','D'),-1), (('T','E'),-1), (('T','F'),-1), + (('T','G'),-2), (('T','H'),-2), (('T','I'),-1), (('T','K'),-1), + (('T','L'),-1), (('T','M'),-1), (('T','N'),0), (('T','P'),-1), + (('T','Q'),-1), (('T','R'),-1), (('T','S'),2), (('T','T'),5), + (('T','V'),0), (('T','W'),-3), (('T','X'),-1), (('T','Y'),-1), + (('T','Z'),-1), (('V','A'),0), (('V','B'),-3), (('V','C'),-1), + (('V','D'),-3), (('V','E'),-3), (('V','F'),0), (('V','G'),-3), + (('V','H'),-3), (('V','I'),3), (('V','K'),-2), (('V','L'),1), + (('V','M'),1), (('V','N'),-3), (('V','P'),-3), (('V','Q'),-3), + (('V','R'),-2), (('V','S'),-1), (('V','T'),0), (('V','V'),5), + (('V','W'),-3), (('V','X'),-1), (('V','Y'),-1), (('V','Z'),-3), + (('W','A'),-2), (('W','B'),-4), (('W','C'),-5), (('W','D'),-4), + (('W','E'),-3), (('W','F'),1), (('W','G'),-2), (('W','H'),-3), + (('W','I'),-2), (('W','K'),-2), (('W','L'),-2), (('W','M'),-2), + (('W','N'),-4), (('W','P'),-3), (('W','Q'),-2), (('W','R'),-2), + (('W','S'),-4), (('W','T'),-3), (('W','V'),-3), (('W','W'),15), + (('W','X'),-1), (('W','Y'),3), (('W','Z'),-2), (('X','A'),-1), + (('X','B'),-1), (('X','C'),-1), (('X','D'),-1), (('X','E'),-1), + (('X','F'),-1), (('X','G'),-1), (('X','H'),-1), (('X','I'),-1), + (('X','K'),-1), (('X','L'),-1), (('X','M'),-1), (('X','N'),-1), + (('X','P'),-1), (('X','Q'),-1), (('X','R'),-1), (('X','S'),-1), + (('X','T'),-1), (('X','V'),-1), (('X','W'),-1), (('X','X'),-1), + (('X','Y'),-1), (('X','Z'),-1), (('Y','A'),-2), (('Y','B'),-2), + (('Y','C'),-3), (('Y','D'),-2), (('Y','E'),-2), (('Y','F'),3), + (('Y','G'),-3), (('Y','H'),2), (('Y','I'),0), (('Y','K'),-1), + (('Y','L'),0), (('Y','M'),0), (('Y','N'),-2), (('Y','P'),-3), + (('Y','Q'),-1), (('Y','R'),-1), (('Y','S'),-2), (('Y','T'),-1), + (('Y','V'),-1), (('Y','W'),3), (('Y','X'),-1), (('Y','Y'),8), + (('Y','Z'),-2), (('Z','A'),-1), (('Z','B'),2), (('Z','C'),-3), + (('Z','D'),1), (('Z','E'),4), (('Z','F'),-3), (('Z','G'),-2), + (('Z','H'),0), (('Z','I'),-3), (('Z','K'),1), (('Z','L'),-2), + (('Z','M'),-1), (('Z','N'),0), (('Z','P'),-1), (('Z','Q'),4), + (('Z','R'),0), (('Z','S'),0), (('Z','T'),-1), (('Z','V'),-3), + (('Z','W'),-2), (('Z','X'),-1), (('Z','Y'),-2), (('Z','Z'),4)] + +-- | The standard BLOSUM62 matrix. +blosum62 m = M.findWithDefault (-4) m $ M.fromList + [(('A','A'),4), (('A','B'),-2), (('A','C'),0), (('A','D'),-2), + (('A','E'),-1), (('A','F'),-2), (('A','G'),0), (('A','H'),-2), + (('A','I'),-1), (('A','K'),-1), (('A','L'),-1), (('A','M'),-1), + (('A','N'),-2), (('A','P'),-1), (('A','Q'),-1), (('A','R'),-1), + (('A','S'),1), (('A','T'),0), (('A','V'),0), (('A','W'),-3), + (('A','X'),-1), (('A','Y'),-2), (('A','Z'),-1), (('B','A'),-2), + (('B','B'),4), (('B','C'),-3), (('B','D'),4), (('B','E'),1), + (('B','F'),-3), (('B','G'),-1), (('B','H'),0), (('B','I'),-3), + (('B','K'),0), (('B','L'),-4), (('B','M'),-3), (('B','N'),3), + (('B','P'),-2), (('B','Q'),0), (('B','R'),-1), (('B','S'),0), + (('B','T'),-1), (('B','V'),-3), (('B','W'),-4), (('B','X'),-1), + (('B','Y'),-3), (('B','Z'),1), (('C','A'),0), (('C','B'),-3), + (('C','C'),9), (('C','D'),-3), (('C','E'),-4), (('C','F'),-2), + (('C','G'),-3), (('C','H'),-3), (('C','I'),-1), (('C','K'),-3), + (('C','L'),-1), (('C','M'),-1), (('C','N'),-3), (('C','P'),-3), + (('C','Q'),-3), (('C','R'),-3), (('C','S'),-1), (('C','T'),-1), + (('C','V'),-1), (('C','W'),-2), (('C','X'),-1), (('C','Y'),-2), + (('C','Z'),-3), (('D','A'),-2), (('D','B'),4), (('D','C'),-3), + (('D','D'),6), (('D','E'),2), (('D','F'),-3), (('D','G'),-1), + (('D','H'),-1), (('D','I'),-3), (('D','K'),-1), (('D','L'),-4), + (('D','M'),-3), (('D','N'),1), (('D','P'),-1), (('D','Q'),0), + (('D','R'),-2), (('D','S'),0), (('D','T'),-1), (('D','V'),-3), + (('D','W'),-4), (('D','X'),-1), (('D','Y'),-3), (('D','Z'),1), + (('E','A'),-1), (('E','B'),1), (('E','C'),-4), (('E','D'),2), + (('E','E'),5), (('E','F'),-3), (('E','G'),-2), (('E','H'),0), + (('E','I'),-3), (('E','K'),1), (('E','L'),-3), (('E','M'),-2), + (('E','N'),0), (('E','P'),-1), (('E','Q'),2), (('E','R'),0), + (('E','S'),0), (('E','T'),-1), (('E','V'),-2), (('E','W'),-3), + (('E','X'),-1), (('E','Y'),-2), (('E','Z'),4), (('F','A'),-2), + (('F','B'),-3), (('F','C'),-2), (('F','D'),-3), (('F','E'),-3), + (('F','F'),6), (('F','G'),-3), (('F','H'),-1), (('F','I'),0), + (('F','K'),-3), (('F','L'),0), (('F','M'),0), (('F','N'),-3), + (('F','P'),-4), (('F','Q'),-3), (('F','R'),-3), (('F','S'),-2), + (('F','T'),-2), (('F','V'),-1), (('F','W'),1), (('F','X'),-1), + (('F','Y'),3), (('F','Z'),-3), (('G','A'),0), (('G','B'),-1), + (('G','C'),-3), (('G','D'),-1), (('G','E'),-2), (('G','F'),-3), + (('G','G'),6), (('G','H'),-2), (('G','I'),-4), (('G','K'),-2), + (('G','L'),-4), (('G','M'),-3), (('G','N'),0), (('G','P'),-2), + (('G','Q'),-2), (('G','R'),-2), (('G','S'),0), (('G','T'),-2), + (('G','V'),-3), (('G','W'),-2), (('G','X'),-1), (('G','Y'),-3), + (('G','Z'),-2), (('H','A'),-2), (('H','B'),0), (('H','C'),-3), + (('H','D'),-1), (('H','E'),0), (('H','F'),-1), (('H','G'),-2), + (('H','H'),8), (('H','I'),-3), (('H','K'),-1), (('H','L'),-3), + (('H','M'),-2), (('H','N'),1), (('H','P'),-2), (('H','Q'),0), + (('H','R'),0), (('H','S'),-1), (('H','T'),-2), (('H','V'),-3), + (('H','W'),-2), (('H','X'),-1), (('H','Y'),2), (('H','Z'),0), + (('I','A'),-1), (('I','B'),-3), (('I','C'),-1), (('I','D'),-3), + (('I','E'),-3), (('I','F'),0), (('I','G'),-4), (('I','H'),-3), + (('I','I'),4), (('I','K'),-3), (('I','L'),2), (('I','M'),1), + (('I','N'),-3), (('I','P'),-3), (('I','Q'),-3), (('I','R'),-3), + (('I','S'),-2), (('I','T'),-1), (('I','V'),3), (('I','W'),-3), + (('I','X'),-1), (('I','Y'),-1), (('I','Z'),-3), (('K','A'),-1), + (('K','B'),0), (('K','C'),-3), (('K','D'),-1), (('K','E'),1), + (('K','F'),-3), (('K','G'),-2), (('K','H'),-1), (('K','I'),-3), + (('K','K'),5), (('K','L'),-2), (('K','M'),-1), (('K','N'),0), + (('K','P'),-1), (('K','Q'),1), (('K','R'),2), (('K','S'),0), + (('K','T'),-1), (('K','V'),-2), (('K','W'),-3), (('K','X'),-1), + (('K','Y'),-2), (('K','Z'),1), (('L','A'),-1), (('L','B'),-4), + (('L','C'),-1), (('L','D'),-4), (('L','E'),-3), (('L','F'),0), + (('L','G'),-4), (('L','H'),-3), (('L','I'),2), (('L','K'),-2), + (('L','L'),4), (('L','M'),2), (('L','N'),-3), (('L','P'),-3), + (('L','Q'),-2), (('L','R'),-2), (('L','S'),-2), (('L','T'),-1), + (('L','V'),1), (('L','W'),-2), (('L','X'),-1), (('L','Y'),-1), + (('L','Z'),-3), (('M','A'),-1), (('M','B'),-3), (('M','C'),-1), + (('M','D'),-3), (('M','E'),-2), (('M','F'),0), (('M','G'),-3), + (('M','H'),-2), (('M','I'),1), (('M','K'),-1), (('M','L'),2), + (('M','M'),5), (('M','N'),-2), (('M','P'),-2), (('M','Q'),0), + (('M','R'),-1), (('M','S'),-1), (('M','T'),-1), (('M','V'),1), + (('M','W'),-1), (('M','X'),-1), (('M','Y'),-1), (('M','Z'),-1), + (('N','A'),-2), (('N','B'),3), (('N','C'),-3), (('N','D'),1), + (('N','E'),0), (('N','F'),-3), (('N','G'),0), (('N','H'),1), + (('N','I'),-3), (('N','K'),0), (('N','L'),-3), (('N','M'),-2), + (('N','N'),6), (('N','P'),-2), (('N','Q'),0), (('N','R'),0), + (('N','S'),1), (('N','T'),0), (('N','V'),-3), (('N','W'),-4), + (('N','X'),-1), (('N','Y'),-2), (('N','Z'),0), (('P','A'),-1), + (('P','B'),-2), (('P','C'),-3), (('P','D'),-1), (('P','E'),-1), + (('P','F'),-4), (('P','G'),-2), (('P','H'),-2), (('P','I'),-3), + (('P','K'),-1), (('P','L'),-3), (('P','M'),-2), (('P','N'),-2), + (('P','P'),7), (('P','Q'),-1), (('P','R'),-2), (('P','S'),-1), + (('P','T'),-1), (('P','V'),-2), (('P','W'),-4), (('P','X'),-1), + (('P','Y'),-3), (('P','Z'),-1), (('Q','A'),-1), (('Q','B'),0), + (('Q','C'),-3), (('Q','D'),0), (('Q','E'),2), (('Q','F'),-3), + (('Q','G'),-2), (('Q','H'),0), (('Q','I'),-3), (('Q','K'),1), + (('Q','L'),-2), (('Q','M'),0), (('Q','N'),0), (('Q','P'),-1), + (('Q','Q'),5), (('Q','R'),1), (('Q','S'),0), (('Q','T'),-1), + (('Q','V'),-2), (('Q','W'),-2), (('Q','X'),-1), (('Q','Y'),-1), + (('Q','Z'),3), (('R','A'),-1), (('R','B'),-1), (('R','C'),-3), + (('R','D'),-2), (('R','E'),0), (('R','F'),-3), (('R','G'),-2), + (('R','H'),0), (('R','I'),-3), (('R','K'),2), (('R','L'),-2), + (('R','M'),-1), (('R','N'),0), (('R','P'),-2), (('R','Q'),1), + (('R','R'),5), (('R','S'),-1), (('R','T'),-1), (('R','V'),-3), + (('R','W'),-3), (('R','X'),-1), (('R','Y'),-2), (('R','Z'),0), + (('S','A'),1), (('S','B'),0), (('S','C'),-1), (('S','D'),0), + (('S','E'),0), (('S','F'),-2), (('S','G'),0), (('S','H'),-1), + (('S','I'),-2), (('S','K'),0), (('S','L'),-2), (('S','M'),-1), + (('S','N'),1), (('S','P'),-1), (('S','Q'),0), (('S','R'),-1), + (('S','S'),4), (('S','T'),1), (('S','V'),-2), (('S','W'),-3), + (('S','X'),-1), (('S','Y'),-2), (('S','Z'),0), (('T','A'),0), + (('T','B'),-1), (('T','C'),-1), (('T','D'),-1), (('T','E'),-1), + (('T','F'),-2), (('T','G'),-2), (('T','H'),-2), (('T','I'),-1), + (('T','K'),-1), (('T','L'),-1), (('T','M'),-1), (('T','N'),0), + (('T','P'),-1), (('T','Q'),-1), (('T','R'),-1), (('T','S'),1), + (('T','T'),5), (('T','V'),0), (('T','W'),-2), (('T','X'),-1), + (('T','Y'),-2), (('T','Z'),-1), (('V','A'),0), (('V','B'),-3), + (('V','C'),-1), (('V','D'),-3), (('V','E'),-2), (('V','F'),-1), + (('V','G'),-3), (('V','H'),-3), (('V','I'),3), (('V','K'),-2), + (('V','L'),1), (('V','M'),1), (('V','N'),-3), (('V','P'),-2), + (('V','Q'),-2), (('V','R'),-3), (('V','S'),-2), (('V','T'),0), + (('V','V'),4), (('V','W'),-3), (('V','X'),-1), (('V','Y'),-1), + (('V','Z'),-2), (('W','A'),-3), (('W','B'),-4), (('W','C'),-2), + (('W','D'),-4), (('W','E'),-3), (('W','F'),1), (('W','G'),-2), + (('W','H'),-2), (('W','I'),-3), (('W','K'),-3), (('W','L'),-2), + (('W','M'),-1), (('W','N'),-4), (('W','P'),-4), (('W','Q'),-2), + (('W','R'),-3), (('W','S'),-3), (('W','T'),-2), (('W','V'),-3), + (('W','W'),11), (('W','X'),-1), (('W','Y'),2), (('W','Z'),-3), + (('X','A'),-1), (('X','B'),-1), (('X','C'),-1), (('X','D'),-1), + (('X','E'),-1), (('X','F'),-1), (('X','G'),-1), (('X','H'),-1), + (('X','I'),-1), (('X','K'),-1), (('X','L'),-1), (('X','M'),-1), + (('X','N'),-1), (('X','P'),-1), (('X','Q'),-1), (('X','R'),-1), + (('X','S'),-1), (('X','T'),-1), (('X','V'),-1), (('X','W'),-1), + (('X','X'),-1), (('X','Y'),-1), (('X','Z'),-1), (('Y','A'),-2), + (('Y','B'),-3), (('Y','C'),-2), (('Y','D'),-3), (('Y','E'),-2), + (('Y','F'),3), (('Y','G'),-3), (('Y','H'),2), (('Y','I'),-1), + (('Y','K'),-2), (('Y','L'),-1), (('Y','M'),-1), (('Y','N'),-2), + (('Y','P'),-3), (('Y','Q'),-1), (('Y','R'),-2), (('Y','S'),-2), + (('Y','T'),-2), (('Y','V'),-1), (('Y','W'),2), (('Y','X'),-1), + (('Y','Y'),7), (('Y','Z'),-2), (('Z','A'),-1), (('Z','B'),1), + (('Z','C'),-3), (('Z','D'),1), (('Z','E'),4), (('Z','F'),-3), + (('Z','G'),-2), (('Z','H'),0), (('Z','I'),-3), (('Z','K'),1), + (('Z','L'),-3), (('Z','M'),-1), (('Z','N'),0), (('Z','P'),-1), + (('Z','Q'),3), (('Z','R'),0), (('Z','S'),0), (('Z','T'),-1), + (('Z','V'),-2), (('Z','W'),-3), (('Z','X'),-1), (('Z','Y'),-2), + (('Z','Z'),4)] + +-- | The standard BLOSUM80 matrix. +blosum80 m = M.findWithDefault (-6) m $ M.fromList [(('A','A'),5), + (('A','B'),-2), (('A','C'),-1), (('A','D'),-2), (('A','E'),-1), + (('A','F'),-3), (('A','G'),0), (('A','H'),-2), (('A','I'),-2), + (('A','K'),-1), (('A','L'),-2), (('A','M'),-1), (('A','N'),-2), + (('A','P'),-1), (('A','Q'),-1), (('A','R'),-2), (('A','S'),1), + (('A','T'),0), (('A','V'),0), (('A','W'),-3), (('A','X'),-1), + (('A','Y'),-2), (('A','Z'),-1), (('B','A'),-2), (('B','B'),4), + (('B','C'),-4), (('B','D'),4), (('B','E'),1), (('B','F'),-4), + (('B','G'),-1), (('B','H'),-1), (('B','I'),-4), (('B','K'),-1), + (('B','L'),-4), (('B','M'),-3), (('B','N'),4), (('B','P'),-2), + (('B','Q'),0), (('B','R'),-2), (('B','S'),0), (('B','T'),-1), + (('B','V'),-4), (('B','W'),-5), (('B','X'),-1), (('B','Y'),-3), + (('B','Z'),0), (('C','A'),-1), (('C','B'),-4), (('C','C'),9), + (('C','D'),-4), (('C','E'),-5), (('C','F'),-3), (('C','G'),-4), + (('C','H'),-4), (('C','I'),-2), (('C','K'),-4), (('C','L'),-2), + (('C','M'),-2), (('C','N'),-3), (('C','P'),-4), (('C','Q'),-4), + (('C','R'),-4), (('C','S'),-2), (('C','T'),-1), (('C','V'),-1), + (('C','W'),-3), (('C','X'),-1), (('C','Y'),-3), (('C','Z'),-4), + (('D','A'),-2), (('D','B'),4), (('D','C'),-4), (('D','D'),6), + (('D','E'),1), (('D','F'),-4), (('D','G'),-2), (('D','H'),-2), + (('D','I'),-4), (('D','K'),-1), (('D','L'),-5), (('D','M'),-4), + (('D','N'),1), (('D','P'),-2), (('D','Q'),-1), (('D','R'),-2), + (('D','S'),-1), (('D','T'),-1), (('D','V'),-4), (('D','W'),-6), + (('D','X'),-1), (('D','Y'),-4), (('D','Z'),1), (('E','A'),-1), + (('E','B'),1), (('E','C'),-5), (('E','D'),1), (('E','E'),6), + (('E','F'),-4), (('E','G'),-3), (('E','H'),0), (('E','I'),-4), + (('E','K'),1), (('E','L'),-4), (('E','M'),-2), (('E','N'),-1), + (('E','P'),-2), (('E','Q'),2), (('E','R'),-1), (('E','S'),0), + (('E','T'),-1), (('E','V'),-3), (('E','W'),-4), (('E','X'),-1), + (('E','Y'),-3), (('E','Z'),4), (('F','A'),-3), (('F','B'),-4), + (('F','C'),-3), (('F','D'),-4), (('F','E'),-4), (('F','F'),6), + (('F','G'),-4), (('F','H'),-2), (('F','I'),-1), (('F','K'),-4), + (('F','L'),0), (('F','M'),0), (('F','N'),-4), (('F','P'),-4), + (('F','Q'),-4), (('F','R'),-4), (('F','S'),-3), (('F','T'),-2), + (('F','V'),-1), (('F','W'),0), (('F','X'),-1), (('F','Y'),3), + (('F','Z'),-4), (('G','A'),0), (('G','B'),-1), (('G','C'),-4), + (('G','D'),-2), (('G','E'),-3), (('G','F'),-4), (('G','G'),6), + (('G','H'),-3), (('G','I'),-5), (('G','K'),-2), (('G','L'),-4), + (('G','M'),-4), (('G','N'),-1), (('G','P'),-3), (('G','Q'),-2), + (('G','R'),-3), (('G','S'),-1), (('G','T'),-2), (('G','V'),-4), + (('G','W'),-4), (('G','X'),-1), (('G','Y'),-4), (('G','Z'),-3), + (('H','A'),-2), (('H','B'),-1), (('H','C'),-4), (('H','D'),-2), + (('H','E'),0), (('H','F'),-2), (('H','G'),-3), (('H','H'),8), + (('H','I'),-4), (('H','K'),-1), (('H','L'),-3), (('H','M'),-2), + (('H','N'),0), (('H','P'),-3), (('H','Q'),1), (('H','R'),0), + (('H','S'),-1), (('H','T'),-2), (('H','V'),-4), (('H','W'),-3), + (('H','X'),-1), (('H','Y'),2), (('H','Z'),0), (('I','A'),-2), + (('I','B'),-4), (('I','C'),-2), (('I','D'),-4), (('I','E'),-4), + (('I','F'),-1), (('I','G'),-5), (('I','H'),-4), (('I','I'),5), + (('I','K'),-3), (('I','L'),1), (('I','M'),1), (('I','N'),-4), + (('I','P'),-4), (('I','Q'),-3), (('I','R'),-3), (('I','S'),-3), + (('I','T'),-1), (('I','V'),3), (('I','W'),-3), (('I','X'),-1), + (('I','Y'),-2), (('I','Z'),-4), (('K','A'),-1), (('K','B'),-1), + (('K','C'),-4), (('K','D'),-1), (('K','E'),1), (('K','F'),-4), + (('K','G'),-2), (('K','H'),-1), (('K','I'),-3), (('K','K'),5), + (('K','L'),-3), (('K','M'),-2), (('K','N'),0), (('K','P'),-1), + (('K','Q'),1), (('K','R'),2), (('K','S'),-1), (('K','T'),-1), + (('K','V'),-3), (('K','W'),-4), (('K','X'),-1), (('K','Y'),-3), + (('K','Z'),1), (('L','A'),-2), (('L','B'),-4), (('L','C'),-2), + (('L','D'),-5), (('L','E'),-4), (('L','F'),0), (('L','G'),-4), + (('L','H'),-3), (('L','I'),1), (('L','K'),-3), (('L','L'),4), + (('L','M'),2), (('L','N'),-4), (('L','P'),-3), (('L','Q'),-3), + (('L','R'),-3), (('L','S'),-3), (('L','T'),-2), (('L','V'),1), + (('L','W'),-2), (('L','X'),-1), (('L','Y'),-2), (('L','Z'),-3), + (('M','A'),-1), (('M','B'),-3), (('M','C'),-2), (('M','D'),-4), + (('M','E'),-2), (('M','F'),0), (('M','G'),-4), (('M','H'),-2), + (('M','I'),1), (('M','K'),-2), (('M','L'),2), (('M','M'),6), + (('M','N'),-3), (('M','P'),-3), (('M','Q'),0), (('M','R'),-2), + (('M','S'),-2), (('M','T'),-1), (('M','V'),1), (('M','W'),-2), + (('M','X'),-1), (('M','Y'),-2), (('M','Z'),-2), (('N','A'),-2), + (('N','B'),4), (('N','C'),-3), (('N','D'),1), (('N','E'),-1), + (('N','F'),-4), (('N','G'),-1), (('N','H'),0), (('N','I'),-4), + (('N','K'),0), (('N','L'),-4), (('N','M'),-3), (('N','N'),6), + (('N','P'),-3), (('N','Q'),0), (('N','R'),-1), (('N','S'),0), + (('N','T'),0), (('N','V'),-4), (('N','W'),-4), (('N','X'),-1), + (('N','Y'),-3), (('N','Z'),0), (('P','A'),-1), (('P','B'),-2), + (('P','C'),-4), (('P','D'),-2), (('P','E'),-2), (('P','F'),-4), + (('P','G'),-3), (('P','H'),-3), (('P','I'),-4), (('P','K'),-1), + (('P','L'),-3), (('P','M'),-3), (('P','N'),-3), (('P','P'),8), + (('P','Q'),-2), (('P','R'),-2), (('P','S'),-1), (('P','T'),-2), + (('P','V'),-3), (('P','W'),-5), (('P','X'),-1), (('P','Y'),-4), + (('P','Z'),-2), (('Q','A'),-1), (('Q','B'),0), (('Q','C'),-4), + (('Q','D'),-1), (('Q','E'),2), (('Q','F'),-4), (('Q','G'),-2), + (('Q','H'),1), (('Q','I'),-3), (('Q','K'),1), (('Q','L'),-3), + (('Q','M'),0), (('Q','N'),0), (('Q','P'),-2), (('Q','Q'),6), + (('Q','R'),1), (('Q','S'),0), (('Q','T'),-1), (('Q','V'),-3), + (('Q','W'),-3), (('Q','X'),-1), (('Q','Y'),-2), (('Q','Z'),3), + (('R','A'),-2), (('R','B'),-2), (('R','C'),-4), (('R','D'),-2), + (('R','E'),-1), (('R','F'),-4), (('R','G'),-3), (('R','H'),0), + (('R','I'),-3), (('R','K'),2), (('R','L'),-3), (('R','M'),-2), + (('R','N'),-1), (('R','P'),-2), (('R','Q'),1), (('R','R'),6), + (('R','S'),-1), (('R','T'),-1), (('R','V'),-3), (('R','W'),-4), + (('R','X'),-1), (('R','Y'),-3), (('R','Z'),0), (('S','A'),1), + (('S','B'),0), (('S','C'),-2), (('S','D'),-1), (('S','E'),0), + (('S','F'),-3), (('S','G'),-1), (('S','H'),-1), (('S','I'),-3), + (('S','K'),-1), (('S','L'),-3), (('S','M'),-2), (('S','N'),0), + (('S','P'),-1), (('S','Q'),0), (('S','R'),-1), (('S','S'),5), + (('S','T'),1), (('S','V'),-2), (('S','W'),-4), (('S','X'),-1), + (('S','Y'),-2), (('S','Z'),0), (('T','A'),0), (('T','B'),-1), + (('T','C'),-1), (('T','D'),-1), (('T','E'),-1), (('T','F'),-2), + (('T','G'),-2), (('T','H'),-2), (('T','I'),-1), (('T','K'),-1), + (('T','L'),-2), (('T','M'),-1), (('T','N'),0), (('T','P'),-2), + (('T','Q'),-1), (('T','R'),-1), (('T','S'),1), (('T','T'),5), + (('T','V'),0), (('T','W'),-4), (('T','X'),-1), (('T','Y'),-2), + (('T','Z'),-1), (('V','A'),0), (('V','B'),-4), (('V','C'),-1), + (('V','D'),-4), (('V','E'),-3), (('V','F'),-1), (('V','G'),-4), + (('V','H'),-4), (('V','I'),3), (('V','K'),-3), (('V','L'),1), + (('V','M'),1), (('V','N'),-4), (('V','P'),-3), (('V','Q'),-3), + (('V','R'),-3), (('V','S'),-2), (('V','T'),0), (('V','V'),4), + (('V','W'),-3), (('V','X'),-1), (('V','Y'),-2), (('V','Z'),-3), + (('W','A'),-3), (('W','B'),-5), (('W','C'),-3), (('W','D'),-6), + (('W','E'),-4), (('W','F'),0), (('W','G'),-4), (('W','H'),-3), + (('W','I'),-3), (('W','K'),-4), (('W','L'),-2), (('W','M'),-2), + (('W','N'),-4), (('W','P'),-5), (('W','Q'),-3), (('W','R'),-4), + (('W','S'),-4), (('W','T'),-4), (('W','V'),-3), (('W','W'),11), + (('W','X'),-1), (('W','Y'),2), (('W','Z'),-4), (('X','A'),-1), + (('X','B'),-1), (('X','C'),-1), (('X','D'),-1), (('X','E'),-1), + (('X','F'),-1), (('X','G'),-1), (('X','H'),-1), (('X','I'),-1), + (('X','K'),-1), (('X','L'),-1), (('X','M'),-1), (('X','N'),-1), + (('X','P'),-1), (('X','Q'),-1), (('X','R'),-1), (('X','S'),-1), + (('X','T'),-1), (('X','V'),-1), (('X','W'),-1), (('X','X'),-1), + (('X','Y'),-1), (('X','Z'),-1), (('Y','A'),-2), (('Y','B'),-3), + (('Y','C'),-3), (('Y','D'),-4), (('Y','E'),-3), (('Y','F'),3), + (('Y','G'),-4), (('Y','H'),2), (('Y','I'),-2), (('Y','K'),-3), + (('Y','L'),-2), (('Y','M'),-2), (('Y','N'),-3), (('Y','P'),-4), + (('Y','Q'),-2), (('Y','R'),-3), (('Y','S'),-2), (('Y','T'),-2), + (('Y','V'),-2), (('Y','W'),2), (('Y','X'),-1), (('Y','Y'),7), + (('Y','Z'),-3), (('Z','A'),-1), (('Z','B'),0), (('Z','C'),-4), + (('Z','D'),1), (('Z','E'),4), (('Z','F'),-4), (('Z','G'),-3), + (('Z','H'),0), (('Z','I'),-4), (('Z','K'),1), (('Z','L'),-3), + (('Z','M'),-2), (('Z','N'),0), (('Z','P'),-2), (('Z','Q'),3), + (('Z','R'),0), (('Z','S'),0), (('Z','T'),-1), (('Z','V'),-3), + (('Z','W'),-4), (('Z','X'),-1), (('Z','Y'),-3), (('Z','Z'),4)] + +-- | The standard PAM30 matrix +pam30 m = M.findWithDefault (-17) m $ M.fromList [(('A','A'),6), + (('A','B'),-3), (('A','C'),-6), (('A','D'),-3), (('A','E'),-2), + (('A','F'),-8), (('A','G'),-2), (('A','H'),-7), (('A','I'),-5), + (('A','K'),-7), (('A','L'),-6), (('A','M'),-5), (('A','N'),-4), + (('A','P'),-2), (('A','Q'),-4), (('A','R'),-7), (('A','S'),0), + (('A','T'),-1), (('A','V'),-2), (('A','W'),-13), (('A','X'),-3), + (('A','Y'),-8), (('A','Z'),-3), (('B','A'),-3), (('B','B'),6), + (('B','C'),-12), (('B','D'),6), (('B','E'),1), (('B','F'),-10), + (('B','G'),-3), (('B','H'),-1), (('B','I'),-6), (('B','K'),-2), + (('B','L'),-9), (('B','M'),-10), (('B','N'),6), (('B','P'),-7), + (('B','Q'),-3), (('B','R'),-7), (('B','S'),-1), (('B','T'),-3), + (('B','V'),-8), (('B','W'),-10), (('B','X'),-5), (('B','Y'),-6), + (('B','Z'),0), (('C','A'),-6), (('C','B'),-12), (('C','C'),10), + (('C','D'),-14), (('C','E'),-14), (('C','F'),-13), (('C','G'),-9), + (('C','H'),-7), (('C','I'),-6), (('C','K'),-14), (('C','L'),-15), + (('C','M'),-13), (('C','N'),-11), (('C','P'),-8), (('C','Q'),-14), + (('C','R'),-8), (('C','S'),-3), (('C','T'),-8), (('C','V'),-6), + (('C','W'),-15), (('C','X'),-9), (('C','Y'),-4), (('C','Z'),-14), + (('D','A'),-3), (('D','B'),6), (('D','C'),-14), (('D','D'),8), + (('D','E'),2), (('D','F'),-15), (('D','G'),-3), (('D','H'),-4), + (('D','I'),-7), (('D','K'),-4), (('D','L'),-12), (('D','M'),-11), + (('D','N'),2), (('D','P'),-8), (('D','Q'),-2), (('D','R'),-10), + (('D','S'),-4), (('D','T'),-5), (('D','V'),-8), (('D','W'),-15), + (('D','X'),-5), (('D','Y'),-11), (('D','Z'),1), (('E','A'),-2), + (('E','B'),1), (('E','C'),-14), (('E','D'),2), (('E','E'),8), + (('E','F'),-14), (('E','G'),-4), (('E','H'),-5), (('E','I'),-5), + (('E','K'),-4), (('E','L'),-9), (('E','M'),-7), (('E','N'),-2), + (('E','P'),-5), (('E','Q'),1), (('E','R'),-9), (('E','S'),-4), + (('E','T'),-6), (('E','V'),-6), (('E','W'),-17), (('E','X'),-5), + (('E','Y'),-8), (('E','Z'),6), (('F','A'),-8), (('F','B'),-10), + (('F','C'),-13), (('F','D'),-15), (('F','E'),-14), (('F','F'),9), + (('F','G'),-9), (('F','H'),-6), (('F','I'),-2), (('F','K'),-14), + (('F','L'),-3), (('F','M'),-4), (('F','N'),-9), (('F','P'),-10), + (('F','Q'),-13), (('F','R'),-9), (('F','S'),-6), (('F','T'),-9), + (('F','V'),-8), (('F','W'),-4), (('F','X'),-8), (('F','Y'),2), + (('F','Z'),-13), (('G','A'),-2), (('G','B'),-3), (('G','C'),-9), + (('G','D'),-3), (('G','E'),-4), (('G','F'),-9), (('G','G'),6), + (('G','H'),-9), (('G','I'),-11), (('G','K'),-7), (('G','L'),-10), + (('G','M'),-8), (('G','N'),-3), (('G','P'),-6), (('G','Q'),-7), + (('G','R'),-9), (('G','S'),-2), (('G','T'),-6), (('G','V'),-5), + (('G','W'),-15), (('G','X'),-5), (('G','Y'),-14), (('G','Z'),-5), + (('H','A'),-7), (('H','B'),-1), (('H','C'),-7), (('H','D'),-4), + (('H','E'),-5), (('H','F'),-6), (('H','G'),-9), (('H','H'),9), + (('H','I'),-9), (('H','K'),-6), (('H','L'),-6), (('H','M'),-10), + (('H','N'),0), (('H','P'),-4), (('H','Q'),1), (('H','R'),-2), + (('H','S'),-6), (('H','T'),-7), (('H','V'),-6), (('H','W'),-7), + (('H','X'),-5), (('H','Y'),-3), (('H','Z'),-1), (('I','A'),-5), + (('I','B'),-6), (('I','C'),-6), (('I','D'),-7), (('I','E'),-5), + (('I','F'),-2), (('I','G'),-11), (('I','H'),-9), (('I','I'),8), + (('I','K'),-6), (('I','L'),-1), (('I','M'),-1), (('I','N'),-5), + (('I','P'),-8), (('I','Q'),-8), (('I','R'),-5), (('I','S'),-7), + (('I','T'),-2), (('I','V'),2), (('I','W'),-14), (('I','X'),-5), + (('I','Y'),-6), (('I','Z'),-6), (('K','A'),-7), (('K','B'),-2), + (('K','C'),-14), (('K','D'),-4), (('K','E'),-4), (('K','F'),-14), + (('K','G'),-7), (('K','H'),-6), (('K','I'),-6), (('K','K'),7), + (('K','L'),-8), (('K','M'),-2), (('K','N'),-1), (('K','P'),-6), + (('K','Q'),-3), (('K','R'),0), (('K','S'),-4), (('K','T'),-3), + (('K','V'),-9), (('K','W'),-12), (('K','X'),-5), (('K','Y'),-9), + (('K','Z'),-4), (('L','A'),-6), (('L','B'),-9), (('L','C'),-15), + (('L','D'),-12), (('L','E'),-9), (('L','F'),-3), (('L','G'),-10), + (('L','H'),-6), (('L','I'),-1), (('L','K'),-8), (('L','L'),7), + (('L','M'),1), (('L','N'),-7), (('L','P'),-7), (('L','Q'),-5), + (('L','R'),-8), (('L','S'),-8), (('L','T'),-7), (('L','V'),-2), + (('L','W'),-6), (('L','X'),-6), (('L','Y'),-7), (('L','Z'),-7), + (('M','A'),-5), (('M','B'),-10), (('M','C'),-13), (('M','D'),-11), + (('M','E'),-7), (('M','F'),-4), (('M','G'),-8), (('M','H'),-10), + (('M','I'),-1), (('M','K'),-2), (('M','L'),1), (('M','M'),11), + (('M','N'),-9), (('M','P'),-8), (('M','Q'),-4), (('M','R'),-4), + (('M','S'),-5), (('M','T'),-4), (('M','V'),-1), (('M','W'),-13), + (('M','X'),-5), (('M','Y'),-11), (('M','Z'),-5), (('N','A'),-4), + (('N','B'),6), (('N','C'),-11), (('N','D'),2), (('N','E'),-2), + (('N','F'),-9), (('N','G'),-3), (('N','H'),0), (('N','I'),-5), + (('N','K'),-1), (('N','L'),-7), (('N','M'),-9), (('N','N'),8), + (('N','P'),-6), (('N','Q'),-3), (('N','R'),-6), (('N','S'),0), + (('N','T'),-2), (('N','V'),-8), (('N','W'),-8), (('N','X'),-3), + (('N','Y'),-4), (('N','Z'),-3), (('P','A'),-2), (('P','B'),-7), + (('P','C'),-8), (('P','D'),-8), (('P','E'),-5), (('P','F'),-10), + (('P','G'),-6), (('P','H'),-4), (('P','I'),-8), (('P','K'),-6), + (('P','L'),-7), (('P','M'),-8), (('P','N'),-6), (('P','P'),8), + (('P','Q'),-3), (('P','R'),-4), (('P','S'),-2), (('P','T'),-4), + (('P','V'),-6), (('P','W'),-14), (('P','X'),-5), (('P','Y'),-13), + (('P','Z'),-4), (('Q','A'),-4), (('Q','B'),-3), (('Q','C'),-14), + (('Q','D'),-2), (('Q','E'),1), (('Q','F'),-13), (('Q','G'),-7), + (('Q','H'),1), (('Q','I'),-8), (('Q','K'),-3), (('Q','L'),-5), + (('Q','M'),-4), (('Q','N'),-3), (('Q','P'),-3), (('Q','Q'),8), + (('Q','R'),-2), (('Q','S'),-5), (('Q','T'),-5), (('Q','V'),-7), + (('Q','W'),-13), (('Q','X'),-5), (('Q','Y'),-12), (('Q','Z'),6), + (('R','A'),-7), (('R','B'),-7), (('R','C'),-8), (('R','D'),-10), + (('R','E'),-9), (('R','F'),-9), (('R','G'),-9), (('R','H'),-2), + (('R','I'),-5), (('R','K'),0), (('R','L'),-8), (('R','M'),-4), + (('R','N'),-6), (('R','P'),-4), (('R','Q'),-2), (('R','R'),8), + (('R','S'),-3), (('R','T'),-6), (('R','V'),-8), (('R','W'),-2), + (('R','X'),-6), (('R','Y'),-10), (('R','Z'),-4), (('S','A'),0), + (('S','B'),-1), (('S','C'),-3), (('S','D'),-4), (('S','E'),-4), + (('S','F'),-6), (('S','G'),-2), (('S','H'),-6), (('S','I'),-7), + (('S','K'),-4), (('S','L'),-8), (('S','M'),-5), (('S','N'),0), + (('S','P'),-2), (('S','Q'),-5), (('S','R'),-3), (('S','S'),6), + (('S','T'),0), (('S','V'),-6), (('S','W'),-5), (('S','X'),-3), + (('S','Y'),-7), (('S','Z'),-5), (('T','A'),-1), (('T','B'),-3), + (('T','C'),-8), (('T','D'),-5), (('T','E'),-6), (('T','F'),-9), + (('T','G'),-6), (('T','H'),-7), (('T','I'),-2), (('T','K'),-3), + (('T','L'),-7), (('T','M'),-4), (('T','N'),-2), (('T','P'),-4), + (('T','Q'),-5), (('T','R'),-6), (('T','S'),0), (('T','T'),7), + (('T','V'),-3), (('T','W'),-13), (('T','X'),-4), (('T','Y'),-6), + (('T','Z'),-6), (('V','A'),-2), (('V','B'),-8), (('V','C'),-6), + (('V','D'),-8), (('V','E'),-6), (('V','F'),-8), (('V','G'),-5), + (('V','H'),-6), (('V','I'),2), (('V','K'),-9), (('V','L'),-2), + (('V','M'),-1), (('V','N'),-8), (('V','P'),-6), (('V','Q'),-7), + (('V','R'),-8), (('V','S'),-6), (('V','T'),-3), (('V','V'),7), + (('V','W'),-15), (('V','X'),-5), (('V','Y'),-7), (('V','Z'),-6), + (('W','A'),-13), (('W','B'),-10), (('W','C'),-15), + (('W','D'),-15), (('W','E'),-17), (('W','F'),-4), (('W','G'),-15), + (('W','H'),-7), (('W','I'),-14), (('W','K'),-12), (('W','L'),-6), + (('W','M'),-13), (('W','N'),-8), (('W','P'),-14), (('W','Q'),-13), + (('W','R'),-2), (('W','S'),-5), (('W','T'),-13), (('W','V'),-15), + (('W','W'),13), (('W','X'),-11), (('W','Y'),-5), (('W','Z'),-14), + (('X','A'),-3), (('X','B'),-5), (('X','C'),-9), (('X','D'),-5), + (('X','E'),-5), (('X','F'),-8), (('X','G'),-5), (('X','H'),-5), + (('X','I'),-5), (('X','K'),-5), (('X','L'),-6), (('X','M'),-5), + (('X','N'),-3), (('X','P'),-5), (('X','Q'),-5), (('X','R'),-6), + (('X','S'),-3), (('X','T'),-4), (('X','V'),-5), (('X','W'),-11), + (('X','X'),-5), (('X','Y'),-7), (('X','Z'),-5), (('Y','A'),-8), + (('Y','B'),-6), (('Y','C'),-4), (('Y','D'),-11), (('Y','E'),-8), + (('Y','F'),2), (('Y','G'),-14), (('Y','H'),-3), (('Y','I'),-6), + (('Y','K'),-9), (('Y','L'),-7), (('Y','M'),-11), (('Y','N'),-4), + (('Y','P'),-13), (('Y','Q'),-12), (('Y','R'),-10), (('Y','S'),-7), + (('Y','T'),-6), (('Y','V'),-7), (('Y','W'),-5), (('Y','X'),-7), + (('Y','Y'),10), (('Y','Z'),-9), (('Z','A'),-3), (('Z','B'),0), + (('Z','C'),-14), (('Z','D'),1), (('Z','E'),6), (('Z','F'),-13), + (('Z','G'),-5), (('Z','H'),-1), (('Z','I'),-6), (('Z','K'),-4), + (('Z','L'),-7), (('Z','M'),-5), (('Z','N'),-3), (('Z','P'),-4), + (('Z','Q'),6), (('Z','R'),-4), (('Z','S'),-5), (('Z','T'),-6), + (('Z','V'),-6), (('Z','W'),-14), (('Z','X'),-5), (('Z','Y'),-9), + (('Z','Z'),6)] + +-- | The standard PAM70 matrix. +pam70 m = M.findWithDefault (-11) m $ M.fromList [(('A','A'),5), + (('A','B'),-1), (('A','C'),-4), (('A','D'),-1), (('A','E'),-1), + (('A','F'),-6), (('A','G'),0), (('A','H'),-4), (('A','I'),-2), + (('A','K'),-4), (('A','L'),-4), (('A','M'),-3), (('A','N'),-2), + (('A','P'),0), (('A','Q'),-2), (('A','R'),-4), (('A','S'),1), + (('A','T'),1), (('A','V'),-1), (('A','W'),-9), (('A','X'),-2), + (('A','Y'),-5), (('A','Z'),-1), (('B','A'),-1), (('B','B'),5), + (('B','C'),-8), (('B','D'),5), (('B','E'),2), (('B','F'),-7), + (('B','G'),-1), (('B','H'),0), (('B','I'),-4), (('B','K'),-1), + (('B','L'),-6), (('B','M'),-6), (('B','N'),5), (('B','P'),-4), + (('B','Q'),-1), (('B','R'),-4), (('B','S'),0), (('B','T'),-1), + (('B','V'),-5), (('B','W'),-7), (('B','X'),-2), (('B','Y'),-4), + (('B','Z'),1), (('C','A'),-4), (('C','B'),-8), (('C','C'),9), + (('C','D'),-9), (('C','E'),-9), (('C','F'),-8), (('C','G'),-6), + (('C','H'),-5), (('C','I'),-4), (('C','K'),-9), (('C','L'),-10), + (('C','M'),-9), (('C','N'),-7), (('C','P'),-5), (('C','Q'),-9), + (('C','R'),-5), (('C','S'),-1), (('C','T'),-5), (('C','V'),-4), + (('C','W'),-11), (('C','X'),-6), (('C','Y'),-2), (('C','Z'),-9), + (('D','A'),-1), (('D','B'),5), (('D','C'),-9), (('D','D'),6), + (('D','E'),3), (('D','F'),-10), (('D','G'),-1), (('D','H'),-1), + (('D','I'),-5), (('D','K'),-2), (('D','L'),-8), (('D','M'),-7), + (('D','N'),3), (('D','P'),-4), (('D','Q'),0), (('D','R'),-6), + (('D','S'),-1), (('D','T'),-2), (('D','V'),-5), (('D','W'),-10), + (('D','X'),-3), (('D','Y'),-7), (('D','Z'),2), (('E','A'),-1), + (('E','B'),2), (('E','C'),-9), (('E','D'),3), (('E','E'),6), + (('E','F'),-9), (('E','G'),-2), (('E','H'),-2), (('E','I'),-4), + (('E','K'),-2), (('E','L'),-6), (('E','M'),-4), (('E','N'),0), + (('E','P'),-3), (('E','Q'),2), (('E','R'),-5), (('E','S'),-2), + (('E','T'),-3), (('E','V'),-4), (('E','W'),-11), (('E','X'),-3), + (('E','Y'),-6), (('E','Z'),5), (('F','A'),-6), (('F','B'),-7), + (('F','C'),-8), (('F','D'),-10), (('F','E'),-9), (('F','F'),8), + (('F','G'),-7), (('F','H'),-4), (('F','I'),0), (('F','K'),-9), + (('F','L'),-1), (('F','M'),-2), (('F','N'),-6), (('F','P'),-7), + (('F','Q'),-9), (('F','R'),-7), (('F','S'),-4), (('F','T'),-6), + (('F','V'),-5), (('F','W'),-2), (('F','X'),-5), (('F','Y'),4), + (('F','Z'),-9), (('G','A'),0), (('G','B'),-1), (('G','C'),-6), + (('G','D'),-1), (('G','E'),-2), (('G','F'),-7), (('G','G'),6), + (('G','H'),-6), (('G','I'),-6), (('G','K'),-5), (('G','L'),-7), + (('G','M'),-6), (('G','N'),-1), (('G','P'),-3), (('G','Q'),-4), + (('G','R'),-6), (('G','S'),0), (('G','T'),-3), (('G','V'),-3), + (('G','W'),-10), (('G','X'),-3), (('G','Y'),-9), (('G','Z'),-3), + (('H','A'),-4), (('H','B'),0), (('H','C'),-5), (('H','D'),-1), + (('H','E'),-2), (('H','F'),-4), (('H','G'),-6), (('H','H'),8), + (('H','I'),-6), (('H','K'),-3), (('H','L'),-4), (('H','M'),-6), + (('H','N'),1), (('H','P'),-2), (('H','Q'),2), (('H','R'),0), + (('H','S'),-3), (('H','T'),-4), (('H','V'),-4), (('H','W'),-5), + (('H','X'),-3), (('H','Y'),-1), (('H','Z'),1), (('I','A'),-2), + (('I','B'),-4), (('I','C'),-4), (('I','D'),-5), (('I','E'),-4), + (('I','F'),0), (('I','G'),-6), (('I','H'),-6), (('I','I'),7), + (('I','K'),-4), (('I','L'),1), (('I','M'),1), (('I','N'),-3), + (('I','P'),-5), (('I','Q'),-5), (('I','R'),-3), (('I','S'),-4), + (('I','T'),-1), (('I','V'),3), (('I','W'),-9), (('I','X'),-3), + (('I','Y'),-4), (('I','Z'),-4), (('K','A'),-4), (('K','B'),-1), + (('K','C'),-9), (('K','D'),-2), (('K','E'),-2), (('K','F'),-9), + (('K','G'),-5), (('K','H'),-3), (('K','I'),-4), (('K','K'),6), + (('K','L'),-5), (('K','M'),0), (('K','N'),0), (('K','P'),-4), + (('K','Q'),-1), (('K','R'),2), (('K','S'),-2), (('K','T'),-1), + (('K','V'),-6), (('K','W'),-7), (('K','X'),-3), (('K','Y'),-7), + (('K','Z'),-2), (('L','A'),-4), (('L','B'),-6), (('L','C'),-10), + (('L','D'),-8), (('L','E'),-6), (('L','F'),-1), (('L','G'),-7), + (('L','H'),-4), (('L','I'),1), (('L','K'),-5), (('L','L'),6), + (('L','M'),2), (('L','N'),-5), (('L','P'),-5), (('L','Q'),-3), + (('L','R'),-6), (('L','S'),-6), (('L','T'),-4), (('L','V'),0), + (('L','W'),-4), (('L','X'),-4), (('L','Y'),-4), (('L','Z'),-4), + (('M','A'),-3), (('M','B'),-6), (('M','C'),-9), (('M','D'),-7), + (('M','E'),-4), (('M','F'),-2), (('M','G'),-6), (('M','H'),-6), + (('M','I'),1), (('M','K'),0), (('M','L'),2), (('M','M'),10), + (('M','N'),-5), (('M','P'),-5), (('M','Q'),-2), (('M','R'),-2), + (('M','S'),-3), (('M','T'),-2), (('M','V'),0), (('M','W'),-8), + (('M','X'),-3), (('M','Y'),-7), (('M','Z'),-3), (('N','A'),-2), + (('N','B'),5), (('N','C'),-7), (('N','D'),3), (('N','E'),0), + (('N','F'),-6), (('N','G'),-1), (('N','H'),1), (('N','I'),-3), + (('N','K'),0), (('N','L'),-5), (('N','M'),-5), (('N','N'),6), + (('N','P'),-3), (('N','Q'),-1), (('N','R'),-3), (('N','S'),1), + (('N','T'),0), (('N','V'),-5), (('N','W'),-6), (('N','X'),-2), + (('N','Y'),-3), (('N','Z'),-1), (('P','A'),0), (('P','B'),-4), + (('P','C'),-5), (('P','D'),-4), (('P','E'),-3), (('P','F'),-7), + (('P','G'),-3), (('P','H'),-2), (('P','I'),-5), (('P','K'),-4), + (('P','L'),-5), (('P','M'),-5), (('P','N'),-3), (('P','P'),7), + (('P','Q'),-1), (('P','R'),-2), (('P','S'),0), (('P','T'),-2), + (('P','V'),-3), (('P','W'),-9), (('P','X'),-3), (('P','Y'),-9), + (('P','Z'),-2), (('Q','A'),-2), (('Q','B'),-1), (('Q','C'),-9), + (('Q','D'),0), (('Q','E'),2), (('Q','F'),-9), (('Q','G'),-4), + (('Q','H'),2), (('Q','I'),-5), (('Q','K'),-1), (('Q','L'),-3), + (('Q','M'),-2), (('Q','N'),-1), (('Q','P'),-1), (('Q','Q'),7), + (('Q','R'),0), (('Q','S'),-3), (('Q','T'),-3), (('Q','V'),-4), + (('Q','W'),-8), (('Q','X'),-2), (('Q','Y'),-8), (('Q','Z'),5), + (('R','A'),-4), (('R','B'),-4), (('R','C'),-5), (('R','D'),-6), + (('R','E'),-5), (('R','F'),-7), (('R','G'),-6), (('R','H'),0), + (('R','I'),-3), (('R','K'),2), (('R','L'),-6), (('R','M'),-2), + (('R','N'),-3), (('R','P'),-2), (('R','Q'),0), (('R','R'),8), + (('R','S'),-1), (('R','T'),-4), (('R','V'),-5), (('R','W'),0), + (('R','X'),-3), (('R','Y'),-7), (('R','Z'),-2), (('S','A'),1), + (('S','B'),0), (('S','C'),-1), (('S','D'),-1), (('S','E'),-2), + (('S','F'),-4), (('S','G'),0), (('S','H'),-3), (('S','I'),-4), + (('S','K'),-2), (('S','L'),-6), (('S','M'),-3), (('S','N'),1), + (('S','P'),0), (('S','Q'),-3), (('S','R'),-1), (('S','S'),5), + (('S','T'),2), (('S','V'),-3), (('S','W'),-3), (('S','X'),-1), + (('S','Y'),-5), (('S','Z'),-2), (('T','A'),1), (('T','B'),-1), + (('T','C'),-5), (('T','D'),-2), (('T','E'),-3), (('T','F'),-6), + (('T','G'),-3), (('T','H'),-4), (('T','I'),-1), (('T','K'),-1), + (('T','L'),-4), (('T','M'),-2), (('T','N'),0), (('T','P'),-2), + (('T','Q'),-3), (('T','R'),-4), (('T','S'),2), (('T','T'),6), + (('T','V'),-1), (('T','W'),-8), (('T','X'),-2), (('T','Y'),-4), + (('T','Z'),-3), (('V','A'),-1), (('V','B'),-5), (('V','C'),-4), + (('V','D'),-5), (('V','E'),-4), (('V','F'),-5), (('V','G'),-3), + (('V','H'),-4), (('V','I'),3), (('V','K'),-6), (('V','L'),0), + (('V','M'),0), (('V','N'),-5), (('V','P'),-3), (('V','Q'),-4), + (('V','R'),-5), (('V','S'),-3), (('V','T'),-1), (('V','V'),6), + (('V','W'),-10), (('V','X'),-2), (('V','Y'),-5), (('V','Z'),-4), + (('W','A'),-9), (('W','B'),-7), (('W','C'),-11), (('W','D'),-10), + (('W','E'),-11), (('W','F'),-2), (('W','G'),-10), (('W','H'),-5), + (('W','I'),-9), (('W','K'),-7), (('W','L'),-4), (('W','M'),-8), + (('W','N'),-6), (('W','P'),-9), (('W','Q'),-8), (('W','R'),0), + (('W','S'),-3), (('W','T'),-8), (('W','V'),-10), (('W','W'),13), + (('W','X'),-7), (('W','Y'),-3), (('W','Z'),-10), (('X','A'),-2), + (('X','B'),-2), (('X','C'),-6), (('X','D'),-3), (('X','E'),-3), + (('X','F'),-5), (('X','G'),-3), (('X','H'),-3), (('X','I'),-3), + (('X','K'),-3), (('X','L'),-4), (('X','M'),-3), (('X','N'),-2), + (('X','P'),-3), (('X','Q'),-2), (('X','R'),-3), (('X','S'),-1), + (('X','T'),-2), (('X','V'),-2), (('X','W'),-7), (('X','X'),-3), + (('X','Y'),-5), (('X','Z'),-3), (('Y','A'),-5), (('Y','B'),-4), + (('Y','C'),-2), (('Y','D'),-7), (('Y','E'),-6), (('Y','F'),4), + (('Y','G'),-9), (('Y','H'),-1), (('Y','I'),-4), (('Y','K'),-7), + (('Y','L'),-4), (('Y','M'),-7), (('Y','N'),-3), (('Y','P'),-9), + (('Y','Q'),-8), (('Y','R'),-7), (('Y','S'),-5), (('Y','T'),-4), + (('Y','V'),-5), (('Y','W'),-3), (('Y','X'),-5), (('Y','Y'),9), + (('Y','Z'),-7), (('Z','A'),-1), (('Z','B'),1), (('Z','C'),-9), + (('Z','D'),2), (('Z','E'),5), (('Z','F'),-9), (('Z','G'),-3), + (('Z','H'),1), (('Z','I'),-4), (('Z','K'),-2), (('Z','L'),-4), + (('Z','M'),-3), (('Z','N'),-1), (('Z','P'),-2), (('Z','Q'),5), + (('Z','R'),-2), (('Z','S'),-2), (('Z','T'),-3), (('Z','V'),-4), + (('Z','W'),-10), (('Z','X'),-3), (('Z','Y'),-7), (('Z','Z'),5)] + +-- | Blast defaults, use with gap_open = -5 gap_extend = -3 +-- This should really check for valid nucleotides, and perhaps be more +-- lenient in the case of Ns. Oh well. +blastn_default :: Num a => (Chr,Chr) -> a +blastn_default = simpleMx 1 (-3) + +-- | Construct a simple "matrix" from match score\/mismatch penalty +simpleMx :: (Num a) => a -> a -> (Chr,Chr) -> a +simpleMx match mismatch (x,y) = if x==y || x+32==y || x-32==y + then match else mismatch + +{- +-- Helper function for constructing matrices +readMx s = let mx = filter (\l -> not (null l) && head l /= '#') (lines s) + aas = map head $ words (head mx) + row1 r = zipWith (,) (map (\x -> (head $ head r,x)) aas) (map readInt $ tail r) + readInt :: String -> Int + readInt = read + in map (row1 . words) $ tail mx +-} diff --git a/Bio/Alignment/Multiple.hs b/Bio/Alignment/Multiple.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/Multiple.hs @@ -0,0 +1,37 @@ +{- | + Multiple alignments. + +-} + +module Bio.Alignment.Multiple + +where + +import Bio.Alignment.AlignData +import Bio.Sequence +import Bio.Clustering + +-- | 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,Alignment)) -> [Sequence] -> [String] +progressive = undefined + +-- | Derive alignments indirectly, i.e. calculate A|C using alignments A|B and B|C. +-- This is central for 'Coffee' evaluation of alignments, and T-Coffee construction +-- of alignments. +indirect :: Alignment -> Alignment -> Alignment +indirect (Repl x1 x2:xs) (Repl y1 y2:ys) = Repl x1 y2 : indirect xs ys -- assert x2==y1 +indirect xs@(Repl _ _:_) (Ins y1:ys) = Ins y1 : indirect xs ys +indirect (Repl x1 _:xs) (Del y1:ys) = Del x1 : indirect xs ys + +indirect (Del x1:xs) ys = Del x1 : indirect xs ys -- imply del+ins/=repl + +indirect (Ins x1:xs) (Repl _ y2:ys) = Ins y2 : indirect xs ys -- assert x1==y1 +indirect (Ins x1:xs) (Del y1:ys) = indirect xs ys -- assert x1 == y1 +indirect xs@(Ins _:_) (Ins y1:ys) = Ins y1 : indirect xs ys + +indirect [] ys = ys -- assert: all Ins +indirect xs [] = xs -- assert: all Del + + diff --git a/Bio/Alignment/QAlign.hs b/Bio/Alignment/QAlign.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/QAlign.hs @@ -0,0 +1,166 @@ +{- | + Quality-aware alignments + + Generally, quality data are ignored for alignment\/pattern searching + like Smith-Waterman, Needleman-Wunsch, or BLAST(p|n|x). I believe + that accounting for quality will at the very least affect things like + BLAST statistics, and e.g. is crucial for good EST annotation using Blastx. + + This module performs sequences alignments, takes quality values into + account. +-} +{-# LANGUAGE CPP, ParallelListComp #-} +{-# OPTIONS_GHC -fexcess-precision #-} + +-- #define DEBUG + +module Bio.Alignment.QAlign ( + -- * Smith-Waterman, or locally optimal alignment with affine gaps + local_score, local_align + + -- * Needleman-Wunsch, or globally optimal alignment with affine gaps + , global_score, global_align + + -- * Matrix construction + , qualMx + ) where + +import Data.List (maximumBy,partition,unfoldr,zip4,tails) +import qualified Data.ByteString.Lazy as B +import Data.Array.Unboxed + +import Bio.Sequence.SeqData hiding ((!)) +import Bio.Alignment.AlignData (Chr,Edit(..),SubstMx,Alignment,on,isRepl,showalign) + +-- | The selector must take into account the quality of the sequences +-- on Ins\/Del, the average of qualities surrounding the gap is (should be) used +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 sel z s1 s2 = columns' sel (z,r0,c0) s1' s2' + where (s1',s2') = (zup s1, zup s2) + zup :: Sequence -> [(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) + + -- the first row consists of increasing numbers of deletions + r0 = map (sel . return) + (zip4 (z:r0) (map (Del . fst) s1') (repeat (snd $ head s2')) (map snd s1')) + -- the first column consists of increasing numbers of insertions + c0 = map (sel . return) + (zip4 (z:c0) (map (Ins . fst) s2') (repeat (snd $ head s1')) (map snd s2')) + +columns' :: QSelector a -> (a,[a],[a]) -> [(Chr,Qual)] -> [(Chr,Qual)] -> [[a]] +columns' f (topleft,top,left) s1 s2 = let + c0 = (topleft : left) + -- given the previous column, and the remains of s2, calculate the next column + mkcol (ts, p0:prev, x) = if null x then Nothing else + let (xi,qi) = head x + c = head ts : [f [del,ins,rep] + | del <- zip4 prev (repeat $ Del xi) (repeat qi) (avg2 $ map snd s2) + | ins <- zip4 c (map (Ins . fst) s2) (repeat $ head $ avg2 $ map snd x) (map snd s2) + | rep <- zip4 (p0:prev) (map (Repl xi . fst) s2) (repeat qi) (map snd s2)] + in Just (c, (tail ts, c, tail x)) + in c0 : unfoldr mkcol (top,c0,s1) + +avg2 :: [Qual] -> [Qual] +avg2 = map f . tails + where f (x1:x2:_) = (x1+x2) `div` 2 + f [x] = x + f _ = error "Nasty - incorrect column lenght" + + + +-- | Minus infinity (or an approximation thereof) +minf :: Double +minf = -100000000 + +type QualMx a = Qual -> Qual -> SubstMx a + +qualMx :: Qual -> Qual -> (Chr,Chr) -> Double +qualMx q1 q2 (x,y) = if isN x || isN y then 0.0 else +#ifdef DEBUG + if q1 < 0 || q1 > 99 || q2 < 0 || q2 > 99 + then error ("Qualities out of range: "++show (q1,q2)) + else +#endif + if x==y || x+32==y || x-32==y + then matchtbl!(q1,q2) else mismatchtbl!(q1,q2) + where matchtbl, mismatchtbl :: UArray (Qual,Qual) Double + matchtbl = array ((0,0),(99,99)) [((x,y),adjust True x y) | x <- [0..99], y <- [0..99]] + mismatchtbl = array ((0,0),(99,99)) [((x,y),adjust False x y) | x <- [0..99], y <- [0..99]] + isN c = c `elem` [78,88,110,120] + +adjust :: Bool -> Qual -> Qual -> Double +adjust s q1 q2 = + let fromQual x = 10**(-fromIntegral x/10) + e1 = fromQual q1 + e2 = fromQual q2 + e = (e1+e2-4/3*e1*e2) + in logBase 2 (if s then 4*(1-e) else 4/3*e) + +-- ------------------------------------------------------------ +-- Edit distances + +-- | Calculate global edit distance (Needleman-Wunsch alignment score) +global_score :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> 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 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 :: Double -> QualMx 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] + g = maximum $ minf:[max (sub+go) (gap+ge) | ((sub,gap),_op,_q1,_q2) <- ids] + in (s,g) + +-- ------------------------------------------------------------ +-- Alignments (rip from AAlign) + +-- maximum... +max' :: (Double,Alignment) -> (Double,Alignment) -> (Double,Alignment) +max' (x,ax) (y,yx) = if x>=y then (x,ax) else (y,yx) + +-- ... and addition for compound values +fp :: (Double,Alignment) -> (Double,Edit) -> (Double,Alignment) +fp (x,ax) (s,e) = (x+s,e:ax) + +-- | Calculate global alignment (Needleman-Wunsch) +global_align :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double,Alignment) +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 :: QualMx Double -> (Double,Double) -> Sequence -> Sequence -> (Double, Alignment) +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 + +{- +-- | Calculate the optimal match between a suffix of one sequence and a prefix +-- of the other (useful for e.g. sequence assembly) +overlap_align mx g s1 s2 = + columns (align_select minf mx g) ((0,[]),(fst g,[])) s1 s2 +-} + +-- (maybe better to reverse the inputs for global?) +revsnd (s,a) = (s,reverse a) + +-- | Generic scoring and selection for global and local alignment +align_select :: Double -> QualMx Double -> (Double,Double) -> QSelector ((Double,Alignment),(Double,Alignment)) +align_select minf mx (go,ge) cds = + let (reps,ids) = partition (isRepl.snd') cds + s = maximumBy (compare `on` fst) + $ (minf,[]):[max' sub gap `fp` (mx q1 q2 (x,y),e) | ((sub,gap),e@(Repl x y),q1,q2) <- reps] + g = maximumBy (compare `on` fst) + $ (minf,[]):[max' (sub `fp` (go,e)) (gap `fp` (ge,e)) | ((sub,gap),e,_q1,_q2) <- ids] + in (s,g) + +snd' (_,x,_,_) = x diff --git a/Bio/Alignment/SAlign.hs b/Bio/Alignment/SAlign.hs new file mode 100644 --- /dev/null +++ b/Bio/Alignment/SAlign.hs @@ -0,0 +1,62 @@ +{- | + Simple alignment of sequences + + Standard alignment\/edit distance +-} + +module Bio.Alignment.SAlign + ( + -- * Smith-Waterman, or locally optimal alignment + local_score, local_align + + -- * Needleman-Wunsch, or globally optimal alignment + , global_score, global_align + + ) where + +import Data.List (maximumBy) +import Bio.Sequence.SeqData +import Bio.Alignment.AlignData + +-- ------------------------------------------------------------ +-- Edit distances + +-- | Calculate global edit distance (Needleman-Wunsch alignment score) +global_score :: (Num a, Ord a) => SubstMx a -> a -> Sequence -> Sequence -> 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 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 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 mx g cds = maximum (0:[s+eval mx g e | (s,e) <- cds]) + +-- ------------------------------------------------------------ +-- Alignments + +-- It is of course possible to retain the columns from an alignment score above, and +-- do the usual backtracking, but it is simpler and more efficient on average to +-- store the current alignment (a list; essentially a pointer backwards) along with +-- 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 -> Alignment +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,Alignment) +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 -> Alignment +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,Alignment) +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/Clustering.hs b/Bio/Clustering.hs new file mode 100644 --- /dev/null +++ b/Bio/Clustering.hs @@ -0,0 +1,36 @@ +{- | Implement clustering + +-} + +module Bio.Clustering where + +import qualified Data.Set as S +import Data.List (partition,foldl') + +-- | Data structure for storing hierarchical clusters +data Clustered score datum = Branch score (Clustered score datum) (Clustered score datum) + | Leaf datum deriving Show + +-- | Single linkage agglomerative clustering. +-- Cluster elements by slurping a sorted list of pairs with score (i.e. triples :-) +-- Keeps a set of contained elements at each branch's root, so O(n log n), +-- and requires elements to be in Ord. +-- For this to work, the triples must be sorted on score. Earlier scores in the list will +-- make up the lower nodes, so sort descending for similarity, ascending for distance. +cluster_sl :: (Ord a, Ord s) => [(s,a,a)] -> [Clustered s a] +cluster_sl = map fst . foldl' csl [] + where csl cs (s,a,b) = + -- can be short circuited for more performance + let (acs,tmp) = partition (\(_,objs) -> a `S.member` objs) cs + (bcs,rest) = partition (\(_,objs) -> b `S.member` objs) tmp + in case (acs,bcs) of + ([(ac,ao)],[(bc,bo)]) -> (Branch s ac bc,S.union ao bo):rest + ([(ac,ao)],[]) -> if b `S.member` ao then (ac,ao):rest + else (Branch s ac (Leaf b),S.insert b ao):rest + ([],[(bc,bo)]) -> if a `S.member` bo then (bc,bo):rest + else (Branch s bc (Leaf a),S.insert a bo):rest + ([],[]) -> (Branch s (Leaf a) (Leaf b),S.fromList [a,b]):rest + _ -> error "Grave mistake" + + +-- cluster_gen :: [a] -> (a->a->Bool) -> diff --git a/Bio/Sequence.hs b/Bio/Sequence.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence.hs @@ -0,0 +1,57 @@ +{- | + This is a meta-module importing and re-exporting sequence-related stuff. + + It encompasses the "Bio.Sequence.SeqData", "Bio.Sequence.Fasta", and "Bio.Sequence.TwoBit" modules. +-} + +module Bio.Sequence + ( + -- * Data structures etc ("Bio.Sequence.SeqData") + Sequence(..), Offset, SeqData, Qual, QualData + -- ** Accessor functions + , seqlength, seqlabel, seqheader, seqdata, seqqual, (!) + + -- ** Converting to and from String. + , fromStr, toStr + -- ** Nucleotide functionality. + , compl, revcompl + -- ** Protein sequence functionality + , Amino(..), translate, fromIUPAC, toIUPAC + + -- * File formats + -- ** The Fasta file format ("Bio.Sequence.Fasta") + , readFasta, hReadFasta + , writeFasta, hWriteFasta + -- ** Quality data + -- | Not part of the Fasta format, and treated separately. + , readQual, writeQual, hWriteQual + , readFastaQual + , writeFastaQual, hWriteFastaQual + + -- ** The phd file format ("Bio.Sequence.Phd") + -- | These contain base (nucleotide) calling information, + -- and are generated by @phred@. + , readPhd, hReadPhd + + -- ** TwoBit file format support ("Bio.Seqeunce.TwoBit") + -- | Used by @BLAT@ and related tools. + , decode2Bit, read2Bit, hRead2Bit + -- ,encode2Bit, write2Bit, hWrite2Bit + + -- * Hashing functionality ("Bio.Sequence.HashWord") + -- | Packing words from sequences into integral data types + , HashF (..) + , contigous, rcontig, rcpacked + ) where + +-- basic sequence data structures +import Bio.Sequence.SeqData + +-- file formats +import Bio.Sequence.Fasta +import Bio.Sequence.Phd +import Bio.Sequence.TwoBit + +-- sequence-oriented stuff +import Bio.Sequence.Entropy +import Bio.Sequence.HashWord diff --git a/Bio/Sequence/Entropy.hs b/Bio/Sequence/Entropy.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/Entropy.hs @@ -0,0 +1,22 @@ +module Bio.Sequence.Entropy where + +import Data.List + +class KWords s where + kwords :: Int -> s -> [s] + +instance KWords [a] where + kwords k = dropWhile (( Int -> str -> Double +entropy k s = negate . sum . map nlogn $ probs ls + where ls = map (fromIntegral . length) . group . sort . kwords k $ s + +nlogn :: (Floating a) => a -> a +nlogn x = x*log x/log 2 +probs :: (Fractional a) => [a] -> [a] +probs ls = map (/ sum ls) ls diff --git a/Bio/Sequence/Fasta.hs b/Bio/Sequence/Fasta.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/Fasta.hs @@ -0,0 +1,203 @@ +{- | + Module: Bio.Sequence.Fasta + + This module incorporates functionality for reading and writing + sequence data in the Fasta format. + Each sequence consists of a header (with a '>' prefix) + and a set of lines containing the sequence data. +-} + +module Bio.Sequence.Fasta + ( + -- * Reading and writing plain FASTA files + readFasta, writeFasta, hReadFasta, hWriteFasta + -- * Reading and writing quality files + , readQual, writeQual, hWriteQual + -- * Combining FASTA and quality files + , readFastaQual, hWriteFastaQual, writeFastaQual + -- * Counting sequences in a FASTA file + , countSeqs + -- * Helper function for reading your own sequences + , mkSeqs +) where + + +-- import Data.Char (isSpace) +import Data.List (groupBy,intersperse) +import Data.Int +import Data.Maybe +import System.IO +import System.IO.Unsafe +import Control.Monad + +import qualified Data.ByteString.Char8 as BS +import qualified Data.ByteString.Lazy.Char8 as B +import qualified Data.ByteString.Lazy as BB +import Data.ByteString.Lazy.Char8 (ByteString) + +import Bio.Sequence.SeqData + +splitsAt :: Offset -> ByteString -> [ByteString] +splitsAt n s = let (s1,s2) = B.splitAt n s + in if B.null s2 then [s1] else s1 : splitsAt n s2 + +{- +-- ugly? +class SeqType sd where + toSeq :: sd -> sd -> Sequence + fromSeq :: Sequence -> (sd,sd) + +instance SeqType B.ByteString where + toSeq = Seq + fromSeq (Seq x y) = (x,y) + +instance SeqType BS.ByteString where + toSeq h s = Seq (B.fromChunks [h]) (B.fromChunks [s]) + fromSeq (Seq x y) = (c x, c y) where c = BS.concat . B.toChunks +-} + +-- | Lazily read sequences from a FASTA-formatted file +readFasta :: FilePath -> IO [Sequence] +readFasta f = do + ls <- getLines f + return (mkSeqs ls) + +-- | Write sequences to a FASTA-formatted file. +-- Line length is 60. +writeFasta :: FilePath -> [Sequence] -> 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 f = do + ls <- getLines f + return (mkQual ls) + +-- | Write quality data for sequences to a file. +writeQual :: FilePath -> [Sequence] -> IO () +writeQual f ss = do + h <- openFile f WriteMode + hWriteQual h ss + hClose h + +-- | 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 s q = do + ss <- readFasta s + qs <- readQual q + return ss + -- warning: assumes correct qual file! + let mkseq s1@(Seq x y _) s2@(Seq _ _ (Just z)) + | seqlabel s1 /= seqlabel s2 = error ("mismatching sequence and quality: " + ++show (seqlabel s1,seqlabel s2)) + | B.length y /= B.length z = error ("mismatching sequence and quality lengths:" + ++ show (seqlabel s1,B.length y,B.length z)) + | otherwise = Seq x y (Just z) + return $ zipWith mkseq ss qs + +-- | Write sequence and quality data simulatnously +-- This may be more laziness-friendly. +writeFastaQual :: FilePath -> FilePath -> [Sequence] -> IO () +writeFastaQual f q ss = do + hf <- openFile f WriteMode + hq <- openFile q WriteMode + hWriteFastaQual hf hq ss + hClose hq + hClose hf + +hWriteFastaQual :: Handle -> Handle -> [Sequence] -> 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 h = do + ls <- hGetLines h + return (mkSeqs ls) + +-- | Write sequences in FASTA format to a handle. +hWriteFasta :: Handle -> [Sequence] -> IO () +hWriteFasta h = mapM_ (wFasta h) + +wHead :: Handle -> SeqData -> IO () +wHead h l = do + B.hPut h $ B.pack ">" + B.hPut h l + B.hPut h $ B.pack "\n" + +wFasta :: Handle -> Sequence -> 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 h = mapM_ (wQual h) + +wQual :: Handle -> Sequence -> IO () +wQual h (Seq l _ (Just q)) = do + wHead h l + let qls = splitsAt 20 q + qs = B.pack . unwords . map show . BB.unpack + mapM_ ((\l -> B.hPut h l >> B.hPut h (B.pack "\n")) . qs) qls +wQual _ (Seq _ _ Nothing) = return () + +-- ByteString operations +-- These aren't (or possible weren't) provided by the FPS library. +-- Implement line-based IO (in retrospect it'd be simpler +-- and better to just use 'lines'.) + +-- lazily read lines from file +getLines :: FilePath -> IO [ByteString] +getLines f = do + h <- openFile f ReadMode + hGetLines' (hClose h) h + +-- lazily read lines from handle +hGetLines :: Handle -> IO [ByteString] +hGetLines = hGetLines' (return ()) + +-- add an optional handle-closing parameter +hGetLines' :: IO () -> Handle -> IO [ByteString] +hGetLines' c h = do + e <- hIsEOF h + if e then c >> return [] + else do l' <- BS.hGetLine h + let l = B.fromChunks $ if BS.null l' then [] else [l'] + ls <- unsafeInterleaveIO $ hGetLines' c h + return (l:ls) + +-- | Convert a list of FASTA-formatted lines into a list of sequences. +-- 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 = map mkSeq . blocks + +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 = map f . blocks + where f (l:ls) = Seq (B.drop 1 l) B.empty + (Just $ BB.pack (map int (B.words $ B.unlines ls))) + int = fromIntegral . fst . fromJust' . B.readInt + fromJust' = maybe (error "Error in qual format") id + +-- | Split lines into blocks starting with '>' characters +-- Filter out # comments (but not semicolons?) +blocks :: [ByteString] -> [[ByteString]] +blocks = groupBy (const (('>' /=) . B.head)) . filter ((/='#') . B.head) . filter (not . B.null) + +countSeqs :: FilePath -> IO Int +countSeqs f = do + ss <- B.readFile f + let hdrs = filter (('>'==).B.head) $ filter (not . B.null) $ B.lines ss + return (length hdrs) + diff --git a/Bio/Sequence/GOA.hs b/Bio/Sequence/GOA.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/GOA.hs @@ -0,0 +1,89 @@ +{- | + GOA - parse and index Gene Onthology Annotations + In particular, the file 'gene_association.goa_uniprot' that contains + links between GO terms and UniProt accessions. + + (Where to find the hierarchical relationship between GO terms?) + contains isA relationships + describes the format +-} + +module Bio.Sequence.GOA where +import Data.ByteString.Lazy.Char8 (ByteString,pack,unpack,copy) +import qualified Data.ByteString.Lazy.Char8 as B + +-- | Read the goa_uniprot file (warning: this one is huge!) +readGOA :: FilePath -> IO [Annotation] +readGOA f = B.readFile f >>= + return . map mkAnn . decomment + +-- | Read GO term definitions +readGO :: FilePath -> IO [GoDef] +readGO f = B.readFile f >>= + return . map mkGoDef . decomment + +decomment :: ByteString -> [ByteString] +decomment = filter (\l -> not (B.null l) && B.head l /= '!') . B.lines + +newtype GoTerm = GO Int deriving (Eq,Ord) +type UniProtAcc = ByteString +data GoClass = Func | Proc | Comp + +instance Read GoTerm where + readsPrec n ('G':'O':':':xs) = map (\(i,s)-> (GO i,s)) (readsPrec n xs) + readsPrec n e = error ("couldn't parse GO term: "++show e) +instance Show GoTerm where show (GO x) = "GO:"++show x + +instance Read GoClass where + readsPrec _ ('F':xs) = [(Func,xs)] + readsPrec _ ('P':xs) = [(Proc,xs)] + readsPrec _ ('C':xs) = [(Comp,xs)] + readsPrec _ _ = [] +instance Show GoClass where + show Func = "F" + show Proc = "P" + show Comp = "C" + +-- | GOA Annotation - or multiple annotations? +data Annotation = Ann !UniProtAcc !GoTerm !EvidenceCode deriving (Show) + +mkAnn :: ByteString -> Annotation +mkAnn = pick . B.words + where pick (_db:up:rest) = pick' up $ getGo rest + pick' up' (go:_:ev:_) = Ann (copy up') (read $ unpack go) (read $ unpack ev) + getGo = dropWhile (not . B.isPrefixOf (pack "GO:")) + +-- | GO maps GO terms (GO:xxxx for some number xxxx) to biologically meaningful terms. +-- Defined in +-- The format is "GO:0000000 [tab] text string [tab] F|P|C" +data GoDef = GoDef !GoTerm !ByteString !GoClass deriving (Show) + +mkGoDef :: ByteString -> GoDef +mkGoDef = pick . B.split '\t' + where pick [go,desc,cls] = GoDef (read $ unpack go) (copy desc) (read $ unpack cls) + pick _xs = error ("Couldn't decipher GO definition from: "++show _xs) + +-- | Evidence codes describe the type of support for an annotation +-- +data EvidenceCode = IC -- Inferred by Curator + | IDA -- Inferred from Direct Assay + | IEA -- Inferred from Electronic Annotation + | IEP -- Inferred from Expression Pattern + | IGC -- Inferred from Genomic Context + | IGI -- Inferred from Genetic Interaction + | IMP -- Inferred from Mutant Phenotype + | IPI -- Inferred from Physical Interaction + | ISS -- Inferred from Sequence or Structural Similarity + | NAS -- Non-traceable Author Statement + | ND -- No biological Data available + | RCA -- inferred from Reviewed Computational Analysis + | TAS -- Traceable Author Statement + | NR -- Not Recorded + deriving (Read,Show,Eq) + +-- | The vast majority of GOA data is IEA, while the most reliable information +-- is manually curated. Filtering on this is useful to keep data set sizes +-- manageable, too. +isCurated :: EvidenceCode -> Bool +isCurated = not . (`elem` [ND,IEA]) + diff --git a/Bio/Sequence/HashWord.lhs b/Bio/Sequence/HashWord.lhs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/HashWord.lhs @@ -0,0 +1,217 @@ + | + 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} diff --git a/Bio/Sequence/Phd.hs b/Bio/Sequence/Phd.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/Phd.hs @@ -0,0 +1,51 @@ +{- | + Module: Bio.Sequence.Phd + + Parse phd files (phred base calling output). +-} + +module Bio.Sequence.Phd (readPhd,hReadPhd) where + +import Bio.Sequence.SeqData + +import qualified Data.ByteString.Lazy.Char8 as B +import qualified Data.ByteString.Lazy as BB +import qualified Data.ByteString as BBB -- this is getting silly! + +import System.IO + +-- | Parse a .phd file, extracting the contents as a Sequence +readPhd :: FilePath -> IO Sequence +readPhd f = return . mkPhd =<< B.readFile f + +-- | Parse .phd contents from a handle +hReadPhd :: Handle -> IO Sequence +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 inp = + let (hd:fs) = filter (not . B.null) . B.lines $ inp + (comment,sd) = break (==B.pack "BEGIN_DNA") fs + label = B.drop 15 hd + 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. + +isSubstr :: B.ByteString -> B.ByteString -> Bool +isSubstr s = any (B.isPrefixOf s) . B.tails + +-- | Pack bytestring segments into a single bytestring +-- Allows the (rest of the) file contents to be GC'ed +compact :: B.ByteString -> B.ByteString +compact = B.fromChunks . return . BBB.concat . B.toChunks diff --git a/Bio/Sequence/SeqData.hs b/Bio/Sequence/SeqData.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/SeqData.hs @@ -0,0 +1,202 @@ +{- | + Data structures for manipulating (biological) sequences. + + Generally supports both nucleotide and protein sequences, some functions, + like @revcompl@, only makes sense for nucleotides. +-} + +module Bio.Sequence.SeqData + ( + -- * Data structure + -- | A sequence is a header, sequence data itself, and optional quality data. + -- All items are lazy bytestrings. The Offset type can be used for indexing. + Sequence(..), Offset, SeqData, + + -- | Quality data is normally associated with nucleotide sequences + Qual, QualData + + -- * Accessor functions + , (!), seqlength,seqlabel,seqheader,seqdata + , (?), hasqual, seqqual + + -- * Converting to and from [Char] + , fromStr, toStr + + -- * 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 + + -- * Protein functionality + -- | Proteins are chains of amino acids, represented by the IUPAC alphabet. + , Amino(..), translate, fromIUPAC, toIUPAC -- amino acids + ) where + +import Data.List (unfoldr) +import Data.Char (toUpper) +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 + + +-- | An offset, index, or length of a 'SeqData' +type Offset = Int64 + +-- | The basic data type used in 'Sequence's +type SeqData = B.ByteString + +-- | Basic type for quality data. Range 0..255. Typical Phred output is in +-- the range 6..50, with 20 as the line in the sand separating good from bad. +type Qual = Word8 + +-- | Quality data is a 'Qual' vector, currently implemented as a 'ByteString'. +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) + +-- | Convert a String to 'SeqData' +fromStr :: String -> SeqData +fromStr = B.pack + +-- | Convert a 'SeqData' to a String +toStr :: SeqData -> String +toStr = B.unpack + +-- | Read the character at the specified position in the sequence. +{-# INLINE (!) #-} +(!) :: Sequence -> Offset -> Char +(!) (Seq _ bs _) = B.index bs + +{-# INLINE (?) #-} +(?) :: Sequence -> Offset -> Qual +(?) (Seq _ _ (Just qs)) = BB.index qs + +-- | Return sequence length. +seqlength :: Sequence -> Offset +seqlength (Seq _ bs _) = B.length bs + +-- | Return sequence label (first word of header) +seqlabel :: Sequence -> SeqData +seqlabel (Seq l _ _) = head (B.words l) + +-- | Return full header. +seqheader :: Sequence -> SeqData +seqheader (Seq l _ _) = l + +-- | Return the sequence data. +seqdata :: Sequence -> SeqData +seqdata (Seq _ d _) = d + +-- | Return the quality data, or error if none exist. Use hasqual if in doubt. +seqqual :: Sequence -> QualData +seqqual (Seq _ _ (Just q)) = q + +-- | Check whether the sequence has associated quality data. +hasqual :: Sequence -> Bool +hasqual (Seq _ _ q) = isJust q + +------------------------------------------------------------ +-- Nucleotide (DNA, RNA) specific stuff +------------------------------------------------------------ + +-- | 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) + (maybe Nothing (Just . BB.reverse) qs) + +-- | 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). +compl :: Char -> Char +compl 'A' = 'T' +compl 'T' = 'A' +compl 'C' = 'G' +compl 'G' = 'C' +compl 'a' = 't' +compl 't' = 'a' +compl 'c' = 'g' +compl 'g' = 'c' +compl x = x + +------------------------------------------------------------ +-- Amino acid (protein) stuff +------------------------------------------------------------ + +-- | 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 s' o' = unfoldr triples (s',o') + where triples (s,o) = + if o > seqlength s - 3 then Nothing + else Just (trans1 (map (s!) [o,o+1,o+2]),(s,o+3)) + +-- prop_trans s = tail (translate s 0) == translate s 3 + +trans1 :: String -> Amino +trans1 = maybe Xaa id . flip lookup trans_tbl . map (repUT . toUpper) + where repUT x = if x == 'U' then 'T' else x -- RNA uses U for T + +data Amino = Ala | Arg | Asn | Asp | Cys | Gln | Glu | Gly + | His | Ile | Leu | Lys | Met | Phe | Pro | Ser + | Thr | Tyr | Trp | Val + | STP | Asx | Glx | Xle | Xaa -- unknowns + deriving (Show,Eq) + +trans_tbl :: [(String,Amino)] +trans_tbl = [("GCT",Ala),("GCC",Ala),("GCA",Ala),("GCG",Ala), + ("CGT",Arg),("CGC",Arg),("CGA",Arg),("CGG",Arg),("AGA",Arg),("AGG",Arg), + ("AAT",Asn),("AAC",Asn), + ("GAT",Asp),("GAC",Asp), +-- ("RAT",Asx),("RAC",Asx), -- IUPAC: R is purine (A or G) + ("TGT",Cys),("TGC",Cys), + ("CAA",Gln),("CAG",Gln), + ("GAA",Glu),("GAG",Glu), +-- ("SAA",Glx),("SAG",Glx), -- IUPAC: S is C or G + ("GGT",Gly),("GGC",Gly),("GGA",Gly),("GGG",Gly), + ("CAT",His),("CAC",His), + ("ATT",Ile),("ATC",Ile),("ATA",Ile), + ("TTA",Leu),("TTG",Leu),("CTT",Leu),("CTC",Leu),("CTA",Leu),("CTG",Leu), + ("AAA",Lys),("AAG",Lys), + ("ATG",Met), + ("TTT",Phe),("TTC",Phe), + ("CCT",Pro),("CCC",Pro),("CCA",Pro),("CCG",Pro), + ("TCT",Ser),("TCC",Ser),("TCA",Ser),("TCG",Ser),("AGT",Ser),("AGC",Ser), + ("ACT",Thr),("ACC",Thr),("ACA",Thr),("ACG",Thr), + ("TAT",Tyr),("TAC",Tyr), + ("TGG",Trp), + ("GTT",Val),("GTC",Val),("GTA",Val),("GTG",Val), + ("TAG",STP),("TGA",STP),("TAA",STP) + ] +-- todo: extend with more IUPAC nucleotide wildcards? + +-- | Convert a list of amino acids to a sequence in IUPAC format. +toIUPAC :: [Amino] -> SeqData +toIUPAC = B.pack . map (fromJust . flip lookup iupac) + +-- | Convert a sequence in IUPAC format to a list of amino acids. +fromIUPAC :: SeqData -> [Amino] +fromIUPAC = map (maybe Xaa id . flip lookup iupac') . B.unpack + +iupac :: [(Amino,Char)] +iupac = [(Ala,'A') ,(Arg,'R') ,(Asn,'N') + ,(Asp,'D') ,(Cys,'C') ,(Gln,'Q') + ,(Glu,'E') ,(Gly,'G') ,(His,'H') + ,(Ile,'I') ,(Leu,'L') ,(Lys,'K') + ,(Met,'M') ,(Phe,'F') ,(Pro,'P') + ,(Ser,'S') ,(Thr,'T') ,(Tyr,'Y') + ,(Trp,'W') ,(Val,'V') + ,(Asx,'B') -- Asn or Asp + ,(Glx,'Z') -- Gln or Glu + ,(Xle,'J') -- Ile or Leu + ,(STP,'X') ,(Xaa,'X') + ] +iupac' :: [(Char,Amino)] +iupac' = map (\(a,b)->(b,a)) iupac diff --git a/Bio/Sequence/TwoBit.hs b/Bio/Sequence/TwoBit.hs new file mode 100644 --- /dev/null +++ b/Bio/Sequence/TwoBit.hs @@ -0,0 +1,222 @@ +{- | + This module implements the 2bit format for sequences. + + Based on: + Note! the description is not accurate, it is missing a reserved word + in each sequence record. + + There are also other, completely different ideas of the 2bit format, e.g. + +-} + + +module Bio.Sequence.TwoBit + ( decode2Bit, read2Bit, hRead2Bit + ) where + +import Bio.Sequence.SeqData +import qualified Data.ByteString.Lazy as BB +import qualified Data.ByteString.Lazy.Char8 as B + +import System.IO +import Control.Monad + +import Data.Char +import Data.Binary +import Data.Int +import Data.List +import Data.Bits + +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" + +bswap :: Integral a => Int -> a -> a +bswap n x = let s = bytes x in unbytes . reverse $ s ++ replicate (n-length s) 0 + +bytes :: Integral a => a -> [Word8] +bytes = Data.List.unfoldr (\w -> let (q,r) = quotRem w 256 + in if q == 0 && r == 0 then Nothing + else Just (fromIntegral r,q)) +unbytes :: Integral a => [Word8] -> a +unbytes = Data.List.foldr (\x y -> y*256+x) 0 . map fromIntegral + +instance Arbitrary Word8 where + arbitrary = choose (0,255::Int) >>= return . fromIntegral + +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 +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)) + put (Header m v c r) = do + put default_magic + put default_version + put c + put (0 :: Word32) + +data Entry = Entry { name :: ByteString, offset :: Word32 } + deriving Show + +swapEntry e = e { offset = bswap 4 (offset e) } + +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 + +data Entries = Entries Header [Entry] +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) + +{- + 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 +-} + +data SR = SR { dnaSize, nBlockCount :: Word32 + , nBlockStarts, nBlockSizes :: [Word32] + , maskBlockCount :: Word32 + , maskBlockStarts, 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 + +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) + (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) + +-- used to convert to/from sequence data in the Sequence data structure +fromSR :: SR -> 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 + + 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) + + 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) + +-- | Parse a (lazy) ByteString as sequences in the 2bit format. +decode2Bit :: B.ByteString -> [Sequence] +decode2Bit cs = let (Entries h es) = decode cs :: Entries + ms = map (fromIntegral . offset) es + (c:chunks) = zipWith (-) ms (0:ms) + 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 + +unSRBE (SRBE x) = x +unSRLE (SRLE x) = x + +toSR :: ByteString -> SR +toSR bs = undefined + +-- | Extract sequences from a file in 2bit format. +read2Bit :: FilePath -> IO [Sequence] +read2Bit f = B.readFile f >>= return . decode2Bit + +-- | Extract sequences in the 2bit format from a handle. +hRead2Bit :: Handle -> IO [Sequence] +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 diff --git a/Bio/Util.hs b/Bio/Util.hs new file mode 100644 --- /dev/null +++ b/Bio/Util.hs @@ -0,0 +1,41 @@ +{- | + Utility module, with various useful stuff. +-} + +module Bio.Util where + +import System.IO (stderr, hPutStr, hFlush) +import System.IO.Unsafe (unsafeInterleaveIO) +import Data.List (groupBy) +import qualified Data.ByteString.Lazy.Char8 as B + +-- workaround, current ByteString.Lazy.Char8 contains a bug in 'lines' +lines = case length (B.lines $ B.pack "\n") of 1 -> B.lines + 0 -> mylines +mylines s = case B.elemIndex '\n' s of + Nothing -> if B.null s then [] else [s] + Just i -> B.take i s : mylines (B.drop (i+1) s) + +-- | Break a list of bytestrings on a predicate +splitWhen :: (B.ByteString -> Bool) -> [B.ByteString] -> [[B.ByteString]] +splitWhen p = groupBy (\_ y -> B.null y || not (p y)) + +-- | Output (to stderr) progress while evaluating a lazy list +-- Useful for generating output while (conceptually, at least) in pure code +-- Strictness warning!! This doesn't *quite* work in all cases. Why? +countIO :: String -> String -> Int -> [a] -> IO [a] +countIO msg post step xs = sequence' $ map unsafeInterleaveIO ((blank >> outmsg (0::Int) >> c):cs) + where (c:cs) = ct 0 xs + output = hPutStr stderr + blank = output ('\r':take 70 (repeat ' ')) + outmsg x = output ('\r':msg++show x) >> hFlush stderr + ct s ys = let (a,b) = splitAt (step-1) ys + next = s+step + in case b of [b1] -> map return a ++ [outmsg (s+step) >> hPutStr stderr post >> return b1] + [] -> map return (init a) ++ [outmsg (s+length a) >> hPutStr stderr post >> return (last a)] + _ -> map return a ++ [outmsg s >> return (head b)] ++ ct next (tail b) + +-- lazy sequence +sequence' :: [IO a] -> IO [a] +sequence' ms = foldr k (return []) ms + where k m m' = do { x <- m; xs <- unsafeInterleaveIO m'; return (x:xs) } diff --git a/Bio/Util/Parsex.hs b/Bio/Util/Parsex.hs new file mode 100644 --- /dev/null +++ b/Bio/Util/Parsex.hs @@ -0,0 +1,26 @@ +-- | Lazy "many" combinator for Parsec. +-- Courtesy of Tomasz Zielonka. + +module Bio.Util.Parsex where +import Text.ParserCombinators.Parsec + +lazyMany :: GenParser Char () a -> SourceName -> [Char] -> [a] +lazyMany p file contents = lm state0 + where + Right state0 = parse getParserState file contents + + lm state = case parse p' "" "" of + Left err -> error (show err) + Right x -> x + where + p' = do + setParserState state + choice + [ do + eof + return [] + , do + x <- p + state' <- getParserState + return (x : lm state') + ] diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -0,0 +1,510 @@ + + GNU LESSER GENERAL PUBLIC LICENSE + Version 2.1, February 1999 + + Copyright (C) 1991, 1999 Free Software Foundation, Inc. + 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + Everyone is permitted to copy and distribute verbatim copies + of this license document, but changing it is not allowed. + +[This is the first released version of the Lesser GPL. It also counts + as the successor of the GNU Library Public License, version 2, hence + the version number 2.1.] + + Preamble + + The licenses for most software are designed to take away your +freedom to share and change it. By contrast, the GNU General Public +Licenses are intended to guarantee your freedom to share and change +free software--to make sure the software is free for all its users. + + This license, the Lesser General Public License, applies to some +specially designated software packages--typically libraries--of the +Free Software Foundation and other authors who decide to use it. You +can use it too, but we suggest you first think carefully about whether +this license or the ordinary General Public License is the better +strategy to use in any particular case, based on the explanations +below. + + When we speak of free software, we are referring to freedom of use, +not price. Our General Public Licenses are designed to make sure that +you have the freedom to distribute copies of free software (and charge +for this service if you wish); that you receive source code or can get +it if you want it; that you can change the software and use pieces of +it in new free programs; and that you are informed that you can do +these things. + + To protect your rights, we need to make restrictions that forbid +distributors to deny you these rights or to ask you to surrender these +rights. These restrictions translate to certain responsibilities for +you if you distribute copies of the library or if you modify it. + + For example, if you distribute copies of the library, whether gratis +or for a fee, you must give the recipients all the rights that we gave +you. You must make sure that they, too, receive or can get the source +code. If you link other code with the library, you must provide +complete object files to the recipients, so that they can relink them +with the library after making changes to the library and recompiling +it. And you must show them these terms so they know their rights. + + We protect your rights with a two-step method: (1) we copyright the +library, and (2) we offer you this license, which gives you legal +permission to copy, distribute and/or modify the library. + + To protect each distributor, we want to make it very clear that +there is no warranty for the free library. Also, if the library is +modified by someone else and passed on, the recipients should know +that what they have is not the original version, so that the original +author's reputation will not be affected by problems that might be +introduced by others. + + Finally, software patents pose a constant threat to the existence of +any free program. We wish to make sure that a company cannot +effectively restrict the users of a free program by obtaining a +restrictive license from a patent holder. Therefore, we insist that +any patent license obtained for a version of the library must be +consistent with the full freedom of use specified in this license. + + Most GNU software, including some libraries, is covered by the +ordinary GNU General Public License. This license, the GNU Lesser +General Public License, applies to certain designated libraries, and +is quite different from the ordinary General Public License. We use +this license for certain libraries in order to permit linking those +libraries into non-free programs. + + When a program is linked with a library, whether statically or using +a shared library, the combination of the two is legally speaking a +combined work, a derivative of the original library. The ordinary +General Public License therefore permits such linking only if the +entire combination fits its criteria of freedom. The Lesser General +Public License permits more lax criteria for linking other code with +the library. + + We call this license the "Lesser" General Public License because it +does Less to protect the user's freedom than the ordinary General +Public License. It also provides other free software developers Less +of an advantage over competing non-free programs. These disadvantages +are the reason we use the ordinary General Public License for many +libraries. However, the Lesser license provides advantages in certain +special circumstances. + + For example, on rare occasions, there may be a special need to +encourage the widest possible use of a certain library, so that it +becomes a de-facto standard. To achieve this, non-free programs must +be allowed to use the library. A more frequent case is that a free +library does the same job as widely used non-free libraries. In this +case, there is little to gain by limiting the free library to free +software only, so we use the Lesser General Public License. + + In other cases, permission to use a particular library in non-free +programs enables a greater number of people to use a large body of +free software. For example, permission to use the GNU C Library in +non-free programs enables many more people to use the whole GNU +operating system, as well as its variant, the GNU/Linux operating +system. + + Although the Lesser General Public License is Less protective of the +users' freedom, it does ensure that the user of a program that is +linked with the Library has the freedom and the wherewithal to run +that program using a modified version of the Library. + + The precise terms and conditions for copying, distribution and +modification follow. Pay close attention to the difference between a +"work based on the library" and a "work that uses the library". The +former contains code derived from the library, whereas the latter must +be combined with the library in order to run. + + GNU LESSER GENERAL PUBLIC LICENSE + TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION + + 0. This License Agreement applies to any software library or other +program which contains a notice placed by the copyright holder or +other authorized party saying it may be distributed under the terms of +this Lesser General Public License (also called "this License"). +Each licensee is addressed as "you". + + A "library" means a collection of software functions and/or data +prepared so as to be conveniently linked with application programs +(which use some of those functions and data) to form executables. + + The "Library", below, refers to any such software library or work +which has been distributed under these terms. A "work based on the +Library" means either the Library or any derivative work under +copyright law: that is to say, a work containing the Library or a +portion of it, either verbatim or with modifications and/or translated +straightforwardly into another language. (Hereinafter, translation is +included without limitation in the term "modification".) + + "Source code" for a work means the preferred form of the work for +making modifications to it. For a library, complete source code means +all the source code for all modules it contains, plus any associated +interface definition files, plus the scripts used to control +compilation and installation of the library. + + Activities other than copying, distribution and modification are not +covered by this License; they are outside its scope. The act of +running a program using the Library is not restricted, and output from +such a program is covered only if its contents constitute a work based +on the Library (independent of the use of the Library in a tool for +writing it). Whether that is true depends on what the Library does +and what the program that uses the Library does. + + 1. You may copy and distribute verbatim copies of the Library's +complete source code as you receive it, in any medium, provided that +you conspicuously and appropriately publish on each copy an +appropriate copyright notice and disclaimer of warranty; keep intact +all the notices that refer to this License and to the absence of any +warranty; and distribute a copy of this License along with the +Library. + + You may charge a fee for the physical act of transferring a copy, +and you may at your option offer warranty protection in exchange for a +fee. + + 2. You may modify your copy or copies of the Library or any portion +of it, thus forming a work based on the Library, and copy and +distribute such modifications or work under the terms of Section 1 +above, provided that you also meet all of these conditions: + + a) The modified work must itself be a software library. + + b) You must cause the files modified to carry prominent notices + stating that you changed the files and the date of any change. + + c) You must cause the whole of the work to be licensed at no + charge to all third parties under the terms of this License. + + d) If a facility in the modified Library refers to a function or a + table of data to be supplied by an application program that uses + the facility, other than as an argument passed when the facility + is invoked, then you must make a good faith effort to ensure that, + in the event an application does not supply such function or + table, the facility still operates, and performs whatever part of + its purpose remains meaningful. + + (For example, a function in a library to compute square roots has + a purpose that is entirely well-defined independent of the + application. Therefore, Subsection 2d requires that any + application-supplied function or table used by this function must + be optional: if the application does not supply it, the square + root function must still compute square roots.) + +These requirements apply to the modified work as a whole. If +identifiable sections of that work are not derived from the Library, +and can be reasonably considered independent and separate works in +themselves, then this License, and its terms, do not apply to those +sections when you distribute them as separate works. But when you +distribute the same sections as part of a whole which is a work based +on the Library, the distribution of the whole must be on the terms of +this License, whose permissions for other licensees extend to the +entire whole, and thus to each and every part regardless of who wrote +it. + +Thus, it is not the intent of this section to claim rights or contest +your rights to work written entirely by you; rather, the intent is to +exercise the right to control the distribution of derivative or +collective works based on the Library. + +In addition, mere aggregation of another work not based on the Library +with the Library (or with a work based on the Library) on a volume of +a storage or distribution medium does not bring the other work under +the scope of this License. + + 3. You may opt to apply the terms of the ordinary GNU General Public +License instead of this License to a given copy of the Library. To do +this, you must alter all the notices that refer to this License, so +that they refer to the ordinary GNU General Public License, version 2, +instead of to this License. (If a newer version than version 2 of the +ordinary GNU General Public License has appeared, then you can specify +that version instead if you wish.) Do not make any other change in +these notices. + + Once this change is made in a given copy, it is irreversible for +that copy, so the ordinary GNU General Public License applies to all +subsequent copies and derivative works made from that copy. + + This option is useful when you wish to copy part of the code of +the Library into a program that is not a library. + + 4. You may copy and distribute the Library (or a portion or +derivative of it, under Section 2) in object code or executable form +under the terms of Sections 1 and 2 above provided that you accompany +it with the complete corresponding machine-readable source code, which +must be distributed under the terms of Sections 1 and 2 above on a +medium customarily used for software interchange. + + If distribution of object code is made by offering access to copy +from a designated place, then offering equivalent access to copy the +source code from the same place satisfies the requirement to +distribute the source code, even though third parties are not +compelled to copy the source along with the object code. + + 5. A program that contains no derivative of any portion of the +Library, but is designed to work with the Library by being compiled or +linked with it, is called a "work that uses the Library". Such a +work, in isolation, is not a derivative work of the Library, and +therefore falls outside the scope of this License. + + However, linking a "work that uses the Library" with the Library +creates an executable that is a derivative of the Library (because it +contains portions of the Library), rather than a "work that uses the +library". The executable is therefore covered by this License. +Section 6 states terms for distribution of such executables. + + When a "work that uses the Library" uses material from a header file +that is part of the Library, the object code for the work may be a +derivative work of the Library even though the source code is not. +Whether this is true is especially significant if the work can be +linked without the Library, or if the work is itself a library. The +threshold for this to be true is not precisely defined by law. + + If such an object file uses only numerical parameters, data +structure layouts and accessors, and small macros and small inline +functions (ten lines or less in length), then the use of the object +file is unrestricted, regardless of whether it is legally a derivative +work. (Executables containing this object code plus portions of the +Library will still fall under Section 6.) + + Otherwise, if the work is a derivative of the Library, you may +distribute the object code for the work under the terms of Section 6. +Any executables containing that work also fall under Section 6, +whether or not they are linked directly with the Library itself. + + 6. As an exception to the Sections above, you may also combine or +link a "work that uses the Library" with the Library to produce a +work containing portions of the Library, and distribute that work +under terms of your choice, provided that the terms permit +modification of the work for the customer's own use and reverse +engineering for debugging such modifications. + + You must give prominent notice with each copy of the work that the +Library is used in it and that the Library and its use are covered by +this License. You must supply a copy of this License. If the work +during execution displays copyright notices, you must include the +copyright notice for the Library among them, as well as a reference +directing the user to the copy of this License. Also, you must do one +of these things: + + a) Accompany the work with the complete corresponding + machine-readable source code for the Library including whatever + changes were used in the work (which must be distributed under + Sections 1 and 2 above); and, if the work is an executable linked + with the Library, with the complete machine-readable "work that + uses the Library", as object code and/or source code, so that the + user can modify the Library and then relink to produce a modified + executable containing the modified Library. (It is understood + that the user who changes the contents of definitions files in the + Library will not necessarily be able to recompile the application + to use the modified definitions.) + + b) Use a suitable shared library mechanism for linking with the + Library. A suitable mechanism is one that (1) uses at run time a + copy of the library already present on the user's computer system, + rather than copying library functions into the executable, and (2) + will operate properly with a modified version of the library, if + the user installs one, as long as the modified version is + interface-compatible with the version that the work was made with. + + c) Accompany the work with a written offer, valid for at least + three years, to give the same user the materials specified in + Subsection 6a, above, for a charge no more than the cost of + performing this distribution. + + d) If distribution of the work is made by offering access to copy + from a designated place, offer equivalent access to copy the above + specified materials from the same place. + + e) Verify that the user has already received a copy of these + materials or that you have already sent this user a copy. + + For an executable, the required form of the "work that uses the +Library" must include any data and utility programs needed for +reproducing the executable from it. However, as a special exception, +the materials to be distributed need not include anything that is +normally distributed (in either source or binary form) with the major +components (compiler, kernel, and so on) of the operating system on +which the executable runs, unless that component itself accompanies +the executable. + + It may happen that this requirement contradicts the license +restrictions of other proprietary libraries that do not normally +accompany the operating system. Such a contradiction means you cannot +use both them and the Library together in an executable that you +distribute. + + 7. You may place library facilities that are a work based on the +Library side-by-side in a single library together with other library +facilities not covered by this License, and distribute such a combined +library, provided that the separate distribution of the work based on +the Library and of the other library facilities is otherwise +permitted, and provided that you do these two things: + + a) Accompany the combined library with a copy of the same work + based on the Library, uncombined with any other library + facilities. This must be distributed under the terms of the + Sections above. + + b) Give prominent notice with the combined library of the fact + that part of it is a work based on the Library, and explaining + where to find the accompanying uncombined form of the same work. + + 8. You may not copy, modify, sublicense, link with, or distribute +the Library except as expressly provided under this License. Any +attempt otherwise to copy, modify, sublicense, link with, or +distribute the Library is void, and will automatically terminate your +rights under this License. However, parties who have received copies, +or rights, from you under this License will not have their licenses +terminated so long as such parties remain in full compliance. + + 9. You are not required to accept this License, since you have not +signed it. However, nothing else grants you permission to modify or +distribute the Library or its derivative works. These actions are +prohibited by law if you do not accept this License. Therefore, by +modifying or distributing the Library (or any work based on the +Library), you indicate your acceptance of this License to do so, and +all its terms and conditions for copying, distributing or modifying +the Library or works based on it. + + 10. Each time you redistribute the Library (or any work based on the +Library), the recipient automatically receives a license from the +original licensor to copy, distribute, link with or modify the Library +subject to these terms and conditions. You may not impose any further +restrictions on the recipients' exercise of the rights granted herein. +You are not responsible for enforcing compliance by third parties with +this License. + + 11. If, as a consequence of a court judgment or allegation of patent +infringement or for any other reason (not limited to patent issues), +conditions are imposed on you (whether by court order, agreement or +otherwise) that contradict the conditions of this License, they do not +excuse you from the conditions of this License. If you cannot +distribute so as to satisfy simultaneously your obligations under this +License and any other pertinent obligations, then as a consequence you +may not distribute the Library at all. For example, if a patent +license would not permit royalty-free redistribution of the Library by +all those who receive copies directly or indirectly through you, then +the only way you could satisfy both it and this License would be to +refrain entirely from distribution of the Library. + +If any portion of this section is held invalid or unenforceable under +any particular circumstance, the balance of the section is intended to +apply, and the section as a whole is intended to apply in other +circumstances. + +It is not the purpose of this section to induce you to infringe any +patents or other property right claims or to contest validity of any +such claims; this section has the sole purpose of protecting the +integrity of the free software distribution system which is +implemented by public license practices. Many people have made +generous contributions to the wide range of software distributed +through that system in reliance on consistent application of that +system; it is up to the author/donor to decide if he or she is willing +to distribute software through any other system and a licensee cannot +impose that choice. + +This section is intended to make thoroughly clear what is believed to +be a consequence of the rest of this License. + + 12. If the distribution and/or use of the Library is restricted in +certain countries either by patents or by copyrighted interfaces, the +original copyright holder who places the Library under this License +may add an explicit geographical distribution limitation excluding those +countries, so that distribution is permitted only in or among +countries not thus excluded. In such case, this License incorporates +the limitation as if written in the body of this License. + + 13. The Free Software Foundation may publish revised and/or new +versions of the Lesser General Public License from time to time. +Such new versions will be similar in spirit to the present version, +but may differ in detail to address new problems or concerns. + +Each version is given a distinguishing version number. If the Library +specifies a version number of this License which applies to it and +"any later version", you have the option of following the terms and +conditions either of that version or of any later version published by +the Free Software Foundation. If the Library does not specify a +license version number, you may choose any version ever published by +the Free Software Foundation. + + 14. If you wish to incorporate parts of the Library into other free +programs whose distribution conditions are incompatible with these, +write to the author to ask for permission. For software which is +copyrighted by the Free Software Foundation, write to the Free +Software Foundation; we sometimes make exceptions for this. Our +decision will be guided by the two goals of preserving the free status +of all derivatives of our free software and of promoting the sharing +and reuse of software generally. + + NO WARRANTY + + 15. BECAUSE THE LIBRARY IS LICENSED FREE OF CHARGE, THERE IS NO +WARRANTY FOR THE LIBRARY, TO THE EXTENT PERMITTED BY APPLICABLE LAW. +EXCEPT WHEN OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR +OTHER PARTIES PROVIDE THE LIBRARY "AS IS" WITHOUT WARRANTY OF ANY +KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING, BUT NOT LIMITED TO, THE +IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR +PURPOSE. THE ENTIRE RISK AS TO THE QUALITY AND PERFORMANCE OF THE +LIBRARY IS WITH YOU. SHOULD THE LIBRARY PROVE DEFECTIVE, YOU ASSUME +THE COST OF ALL NECESSARY SERVICING, REPAIR OR CORRECTION. + + 16. IN NO EVENT UNLESS REQUIRED BY APPLICABLE LAW OR AGREED TO IN +WRITING WILL ANY COPYRIGHT HOLDER, OR ANY OTHER PARTY WHO MAY MODIFY +AND/OR REDISTRIBUTE THE LIBRARY AS PERMITTED ABOVE, BE LIABLE TO YOU +FOR DAMAGES, INCLUDING ANY GENERAL, SPECIAL, INCIDENTAL OR +CONSEQUENTIAL DAMAGES ARISING OUT OF THE USE OR INABILITY TO USE THE +LIBRARY (INCLUDING BUT NOT LIMITED TO LOSS OF DATA OR DATA BEING +RENDERED INACCURATE OR LOSSES SUSTAINED BY YOU OR THIRD PARTIES OR A +FAILURE OF THE LIBRARY TO OPERATE WITH ANY OTHER SOFTWARE), EVEN IF +SUCH HOLDER OR OTHER PARTY HAS BEEN ADVISED OF THE POSSIBILITY OF SUCH +DAMAGES. + + END OF TERMS AND CONDITIONS + + How to Apply These Terms to Your New Libraries + + If you develop a new library, and you want it to be of the greatest +possible use to the public, we recommend making it free software that +everyone can redistribute and change. You can do so by permitting +redistribution under these terms (or, alternatively, under the terms +of the ordinary General Public License). + + To apply these terms, attach the following notices to the library. +It is safest to attach them to the start of each source file to most +effectively convey the exclusion of warranty; and each file should +have at least the "copyright" line and a pointer to where the full +notice is found. + + + + Copyright (C) + + This library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + This library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with this library; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +Also add information on how to contact you by electronic and paper mail. + +You should also get your employer (if you work as a programmer) or +your school, if any, to sign a "copyright disclaimer" for the library, +if necessary. Here is a sample; alter the names: + + Yoyodyne, Inc., hereby disclaims all copyright interest in the + library `Frob' (a library for tweaking knobs) written by James + Random Hacker. + + , 1 April 1990 + Ty Coon, President of Vice + +That's all there is to it! + + diff --git a/README b/README new file mode 100644 --- /dev/null +++ b/README @@ -0,0 +1,34 @@ +biolib - a Haskell library for bioinformatics + +This is a collection of data structures and algorithms +I've found useful when building various bioinformatics-related tools +and utilities. + +Current list of features includes: a Sequence data type supporting +protein and nucleotide sequences and conversion between them, 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 +multiple alignment. + +To install, you need to acquire a working GHC (possibly other Haskell +system). You also need the following external libraries: + + QuickCheck - for unit tests + binary - mainly for dealing with the TwoBit sequence format + tagsoup - for parsing XML output from Blast + +You should be able to get what you need from . + +You can then build with 'make', doing either 'make install' if you can sudo, or 'make user_install' if you can not. Of course, the Makefile just proxies for +the regular Cabal routine, which will work just as well: + + runhaskell Setup configure + runhaskell Setup build + sudo runhaskell Setup install + +(Use --prefix=$HOME and remove the sudo, if you don't want to install as root.) + +If that didn't work, mail me at , and we'll try to +work things out. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,3 @@ +#!/usr/bin/env runhaskell +import Distribution.Simple +main = defaultMain diff --git a/bio.cabal b/bio.cabal new file mode 100644 --- /dev/null +++ b/bio.cabal @@ -0,0 +1,48 @@ +Name: bio +Version: 0.3.3 +License: LGPL +License-file: LICENSE +Author: Ketil Malde +Maintainer: ketil@ii.uib.no + +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 + and utilities. + . + Current list of features includes: a Sequence data type supporting + protein and nucleotide sequences and conversion between them, 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 + multiple alignment. + . + The Darcs repository is at: . +Homepage: http://malde.org/~ketil/ + +Tested-With: GHC==6.8.2 +Build-Type: Simple +Build-Depends: base>3, QuickCheck<2, binary, tagsoup>=0.5, bytestring, + containers, array, interlude, parallel, parsec +-- add fps for ghc 6.4.2; change imports in Bio/Sequence/TwoBit.hs if you want QC 2 + +-- We omit the debian/ and Test/ files because those are for development, not installation. +Data-Files: README + +Exposed-modules: Bio.Sequence, + Bio.Sequence.SeqData, + Bio.Sequence.Fasta, Bio.Sequence.TwoBit, Bio.Sequence.Phd, + Bio.Sequence.Entropy, Bio.Sequence.HashWord, + Bio.Sequence.GOA, + 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.Clustering, + Bio.Util, Bio.Util.Parsex + +extensions: CPP, ParallelListComp +ghc-options: -Wall -O2 -fexcess-precision -funbox-strict-fields -auto-all