simseq (empty) → 0.0
raw patch · 10 files changed
+865/−0 lines, 10 filesdep +basedep +biodep +bytestringsetup-changed
Dependencies added: base, bio, bytestring, random
Files
- LICENSE +340/−0
- README +44/−0
- Setup.hs +5/−0
- idea.txt +35/−0
- simseq.cabal +28/−0
- src/Options.hs +151/−0
- src/R454.hs +38/−0
- src/Sanger.hs +33/−0
- src/Simseq.hs +31/−0
- src/UnfoldMut.hs +160/−0
+ LICENSE view
@@ -0,0 +1,340 @@+ GNU GENERAL PUBLIC LICENSE+ Version 2, June 1991++ Copyright (C) 1989, 1991 Free Software Foundation, Inc.+ 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA+ Everyone is permitted to copy and distribute verbatim copies+ of this license document, but changing it is not allowed.++ Preamble++ The licenses for most software are designed to take away your+freedom to share and change it. By contrast, the GNU General Public+License is intended to guarantee your freedom to share and change free+software--to make sure the software is free for all its users. This+General Public License applies to most of the Free Software+Foundation's software and to any other program whose authors commit to+using it. (Some other Free Software Foundation software is covered by+the GNU Library General Public License instead.) You can apply it to+your programs, too.++ When we speak of free software, we are referring to freedom, 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 or use pieces of it+in new free programs; and that you know you can do these things.++ To protect your rights, we need to make restrictions that forbid+anyone to deny you these rights or to ask you to surrender the rights.+These restrictions translate to certain responsibilities for you if you+distribute copies of the software, or if you modify it.++ For example, if you distribute copies of such a program, whether+gratis or for a fee, you must give the recipients all the rights that+you have. You must make sure that they, too, receive or can get the+source code. And you must show them these terms so they know their+rights.++ We protect your rights with two steps: (1) copyright the software, and+(2) offer you this license which gives you legal permission to copy,+distribute and/or modify the software.++ Also, for each author's protection and ours, we want to make certain+that everyone understands that there is no warranty for this free+software. If the software is modified by someone else and passed on, we+want its recipients to know that what they have is not the original, so+that any problems introduced by others will not reflect on the original+authors' reputations.++ Finally, any free program is threatened constantly by software+patents. We wish to avoid the danger that redistributors of a free+program will individually obtain patent licenses, in effect making the+program proprietary. To prevent this, we have made it clear that any+patent must be licensed for everyone's free use or not licensed at all.++ The precise terms and conditions for copying, distribution and+modification follow.++ GNU GENERAL PUBLIC LICENSE+ TERMS AND CONDITIONS FOR COPYING, DISTRIBUTION AND MODIFICATION++ 0. This License applies to any program or other work which contains+a notice placed by the copyright holder saying it may be distributed+under the terms of this General Public License. The "Program", below,+refers to any such program or work, and a "work based on the Program"+means either the Program or any derivative work under copyright law:+that is to say, a work containing the Program or a portion of it,+either verbatim or with modifications and/or translated into another+language. (Hereinafter, translation is included without limitation in+the term "modification".) Each licensee is addressed as "you".++Activities other than copying, distribution and modification are not+covered by this License; they are outside its scope. The act of+running the Program is not restricted, and the output from the Program+is covered only if its contents constitute a work based on the+Program (independent of having been made by running the Program).+Whether that is true depends on what the Program does.++ 1. You may copy and distribute verbatim copies of the Program's+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 give any other recipients of the Program a copy of this License+along with the Program.++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 Program or any portion+of it, thus forming a work based on the Program, 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) You must cause the modified files to carry prominent notices+ stating that you changed the files and the date of any change.++ b) You must cause any work that you distribute or publish, that in+ whole or in part contains or is derived from the Program or any+ part thereof, to be licensed as a whole at no charge to all third+ parties under the terms of this License.++ c) If the modified program normally reads commands interactively+ when run, you must cause it, when started running for such+ interactive use in the most ordinary way, to print or display an+ announcement including an appropriate copyright notice and a+ notice that there is no warranty (or else, saying that you provide+ a warranty) and that users may redistribute the program under+ these conditions, and telling the user how to view a copy of this+ License. (Exception: if the Program itself is interactive but+ does not normally print such an announcement, your work based on+ the Program is not required to print an announcement.)++These requirements apply to the modified work as a whole. If+identifiable sections of that work are not derived from the Program,+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 Program, 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 Program.++In addition, mere aggregation of another work not based on the Program+with the Program (or with a work based on the Program) on a volume of+a storage or distribution medium does not bring the other work under+the scope of this License.++ 3. You may copy and distribute the Program (or a work based on it,+under Section 2) in object code or executable form under the terms of+Sections 1 and 2 above provided that you also do one of the following:++ a) 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; or,++ b) Accompany it with a written offer, valid for at least three+ years, to give any third party, for a charge no more than your+ cost of physically performing source distribution, a complete+ machine-readable copy of the corresponding source code, to be+ distributed under the terms of Sections 1 and 2 above on a medium+ customarily used for software interchange; or,++ c) Accompany it with the information you received as to the offer+ to distribute corresponding source code. (This alternative is+ allowed only for noncommercial distribution and only if you+ received the program in object code or executable form with such+ an offer, in accord with Subsection b above.)++The source code for a work means the preferred form of the work for+making modifications to it. For an executable work, 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 executable. However, as a+special exception, the source code 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.++If distribution of executable or 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 counts as+distribution of the source code, even though third parties are not+compelled to copy the source along with the object code.++ 4. You may not copy, modify, sublicense, or distribute the Program+except as expressly provided under this License. Any attempt+otherwise to copy, modify, sublicense or distribute the Program 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.++ 5. 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 Program or its derivative works. These actions are+prohibited by law if you do not accept this License. Therefore, by+modifying or distributing the Program (or any work based on the+Program), you indicate your acceptance of this License to do so, and+all its terms and conditions for copying, distributing or modifying+the Program or works based on it.++ 6. Each time you redistribute the Program (or any work based on the+Program), the recipient automatically receives a license from the+original licensor to copy, distribute or modify the Program 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 to+this License.++ 7. 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 Program at all. For example, if a patent+license would not permit royalty-free redistribution of the Program 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 Program.++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.++ 8. If the distribution and/or use of the Program is restricted in+certain countries either by patents or by copyrighted interfaces, the+original copyright holder who places the Program 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.++ 9. The Free Software Foundation may publish revised and/or new versions+of the 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 Program+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 Program does not specify a version number of+this License, you may choose any version ever published by the Free Software+Foundation.++ 10. If you wish to incorporate parts of the Program into other free+programs whose distribution conditions are different, 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++ 11. BECAUSE THE PROGRAM IS LICENSED FREE OF CHARGE, THERE IS NO WARRANTY+FOR THE PROGRAM, TO THE EXTENT PERMITTED BY APPLICABLE LAW. EXCEPT WHEN+OTHERWISE STATED IN WRITING THE COPYRIGHT HOLDERS AND/OR OTHER PARTIES+PROVIDE THE PROGRAM "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 PROGRAM IS WITH YOU. SHOULD THE+PROGRAM PROVE DEFECTIVE, YOU ASSUME THE COST OF ALL NECESSARY SERVICING,+REPAIR OR CORRECTION.++ 12. 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 PROGRAM 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 PROGRAM (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 PROGRAM TO OPERATE WITH ANY OTHER+PROGRAMS), 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 Programs++ If you develop a new program, and you want it to be of the greatest+possible use to the public, the best way to achieve this is to make it+free software which everyone can redistribute and change under these terms.++ To do so, attach the following notices to the program. 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.++ <one line to give the program's name and a brief idea of what it does.>+ Copyright (C) <year> <name of author>++ This program is free software; you can redistribute it and/or modify+ it under the terms of the GNU General Public License as published by+ the Free Software Foundation; either version 2 of the License, or+ (at your option) any later version.++ This program 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 General Public License for more details.++ You should have received a copy of the GNU General Public License+ along with this program; if not, write to the Free Software+ Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA+++Also add information on how to contact you by electronic and paper mail.++If the program is interactive, make it output a short notice like this+when it starts in an interactive mode:++ Gnomovision version 69, Copyright (C) year name of author+ Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.+ This is free software, and you are welcome to redistribute it+ under certain conditions; type `show c' for details.++The hypothetical commands `show w' and `show c' should show the appropriate+parts of the General Public License. Of course, the commands you use may+be called something other than `show w' and `show c'; they could even be+mouse-clicks or menu items--whatever suits your program.++You should also get your employer (if you work as a programmer) or your+school, if any, to sign a "copyright disclaimer" for the program, if+necessary. Here is a sample; alter the names:++ Yoyodyne, Inc., hereby disclaims all copyright interest in the program+ `Gnomovision' (which makes passes at compilers) written by James Hacker.++ <signature of Ty Coon>, 1 April 1989+ Ty Coon, President of Vice++This General Public License does not permit incorporating your program into+proprietary programs. If your program is a subroutine library, you may+consider it more useful to permit linking proprietary applications with the+library. If this is what you want to do, use the GNU Library General+Public License instead of this License.
+ README view
@@ -0,0 +1,44 @@++Simseq - SIMulate SEQuences. Yep, that's real creative.++ Synopsis+ --------++Generates a bunch of sequences from a set of reference sequences.+For ESTs, NCBI's refseq transcripts are probably good choices.++The generated sequences are generated using a model that specifies+priming conditions and error generation.++Currently, this is not very refined, you can try++ simseq --model=sanger:n,d reference.fasta++Where n indicates the number of sequences to generate, starting points+drawn from a uniform distribution, and d probability of being in the +forward direction. Or, even more experimentally:++ simseq --model=454:n,d++Which implemets a completely unfounded and baseless model of 454/Roche+pyrosequencing. (Okay, actually based on a paper by Marguiles et al, but+more data is definitely a requirement).++Solexa will be installed as soon as anybody says something definitive+about the error modes.++In any case, running out of sequence results in X's, indicating vector,+which I hope makes sense for Sanger, at least.++ Install+ -------++The usual Cabal routine. Get a working GHC compiler, install+my 'bio' library, and do:++ chmod +x Setup.hs+ ./Setup.hs configure+ ./Setup.hs build+ sudo ./Setup.hs install++Mail me if it didn't work - <ketil at malde.org>.
+ Setup.hs view
@@ -0,0 +1,5 @@+#!/usr/bin/env runhaskell+ +import Distribution.Simple+main = defaultMain+
+ idea.txt view
@@ -0,0 +1,35 @@++Simple tool to generate simulated ESTs.++Must support (minimum, generalize?):++ sanger: initial medium/high error rate,+ stable low error rate (mostly subst?)+ slowly increasing to one from pos 750-900++ 454: medium error rate, mostly ins/del for monomers++ Solexa: any idea of error rate?++Generalized:+ general error models, training on data sets?+ HMM-like?++** PLAN **++ primer: determines starting point (initial state) from sequences+ transcribe: unfold new sequence, while modifying state+ error model: depends on state+ terminator: ends transcript ++Specialized: homopolymers for 454 (must be part of state)+++ Terminator = Distribution+ Error Model = Distribution + Mutator+ Primer = ...?++Mutators can be 'ins(c,p)', 'del(p)', 'subst(c,p) etc.+Or should the distribution multiplier apply already to distrib?++-k
+ simseq.cabal view
@@ -0,0 +1,28 @@+Name: simseq+Version: 0.0+License: GPL+License-File: LICENSE++Author: Ketil Malde+Maintainer: Ketil Malde <ketil@malde.org>++Stability: Alpha+Category: Bioinformatics+Synopsis: Simulate sequencing with different models for priming and errors+Description: This is a simulator that can generate simulated sequences -- primarily EST type sequences,+ but quite possibly other types as well. Mail me for further information on usage etc.+Homepage: http://malde.org/~ketil/+++Build-Depends: base, bio >= 0.2, random, bytestring+ -- add fps for ghc-6.4.2+Build-Type: Simple++data-files: README, idea.txt++Executable: simseq+Main-Is: Simseq.hs+Hs-Source-Dirs: src+Other-modules: Options, R454, UnfoldMut, Sanger++Ghc-Options: -Wall -funbox-strict-fields
+ src/Options.hs view
@@ -0,0 +1,151 @@+{- Handle Options++ Options:++ --prime model for starting points [uniform]+ --error error model [uniform,gradient,..]+ --terminator termination condition [sigma..]+ --model shortcut for useful combinations [sanger|454|solexa]++ -}++module Options where++import System.Environment (getArgs)+import System.Exit (exitWith,ExitCode(..))+import System.Console.GetOpt+import System.IO+import Data.List (unfoldr)+import Control.Monad (when)+import Bio.Sequence++import UnfoldMut+import Sanger+import R454++version, usagemsg :: String+version = "0.0"+usagemsg = "Usage: simseq [options] FILE.."++data Opts = Opts { primer :: Primer+ , errmod :: [Mutator]+ , term :: Terminator+ }++getOpts :: IO (Opts,[Sequence])+getOpts = do+ (opt1,fs,err) <- getArgs >>= (return . getOpt Permute options)+ opts <- parseargs opt1+ when (not $ null err) (error $ usage err)+ when (forceOpts opts) (error "Impossible!")+ ss <- if null fs then hReadFasta stdin+ else return . concat =<< mapM readFasta fs+ return (opts, ss)+ where+ forceOpts (Opts x y z) = x `seq` y `seq` z `seq` False++parseargs :: [Opts -> IO Opts] -> IO Opts+parseargs args = foldl (>>=) (return defaultopts) args++usage :: [String] -> String+usage errs = usageInfo (concat errs ++ usagemsg) options++------------------------------------------------------------++options :: [OptDescr (Opts -> IO Opts)]+options = [ Option [] ["primer"] (ReqArg (setPrimer . processArgs) "p")+ ("Model for initiating sequencing "++show (map fst primers))+ , Option [] ["error"] (ReqArg (setMutator . processArgs) "e")+ ("Error model "++show (map fst mutators))+ , Option [] ["terminator"] (ReqArg (setTerm . processArgs) "t")+ ("Termination condition "++show (map fst terminators))+ , Option [] ["model"] (ReqArg (setModel . processArgs) "m")+ ("Shortcut for useful combinations "++show (map fst models))+ , Option ['h'] ["help"] (ReqArg help "option") ("Help with 'option'")+ ]++defaultopts :: Opts+defaultopts = Opts (error "no primer specified")+ (error "no error model specified")+ (error "no terminator specified")++help :: String -> Opts -> IO Opts+help arg _ = do putStrLn usagemsg+ case arg of+ "primer" -> list "specifies when sequence generation starts" primers+ "error" -> list "adds an error model for sequence generation" mutators+ "terminator" -> list "specifies when sequence generation terminates" terminators+ "model" -> list "specifies a complete model including primer, error model, and terminator" models+ _ -> putStrLn (" "++arg++": no such option")+ exitWith ExitSuccess++ where list str tbl = putStrLn (" "++arg ++ ": "++str++"\n Valid alternatives are: "++unwords (map fst tbl))+++------------------------------------------------------------++-- | Break a string argument representing a function and arguments+-- into components. Syntax is func:arg1,arg2,..+processArgs :: String -> (String,[String])+processArgs ss = let (fn,as) = break (==':') ss+ in (fn,splitStr ',' as)++splitStr :: Char -> String -> [String]+splitStr c = unfoldr (\s -> if null s then Nothing+ else let (x,y) = break (==c) (drop 1 s)+ in Just (x,y))++------------------------------------------------------------+-- Lookup tables for selecting the desired behavior+------------------------------------------------------------++notfound :: String -> [(String,a)] -> String -> String+notfound what whats which = "'"++which++"': no such "++what++".\n"+++ "Available "++what++"s are: "++ (unwords $ map fst whats)++setPrimer :: (String,[String]) -> Opts -> IO Opts+setPrimer (prim,args) (Opts _ m t) = return (Opts p m t)+ where p = maybe (error $ notfound "primer" primers prim)+ (\f -> f args) (lookup prim primers)++primers :: [(String,[String]->Primer)]+primers = [("uniform",p_uniform)]+-- help text? "uniform:n,p","select n random starting points, with probability p of forward direction, 1-p or reverse"++------------------------------------------------------------++setTerm :: (String,[String]) -> Opts -> IO Opts+setTerm (trm,args) (Opts p m _) = return (Opts p m t)+ where t = maybe (error $ notfound "terminator" terminators trm)+ (\f -> f args) (lookup trm terminators)++terminators :: [(String,[String]->Terminator)]+terminators =+ [("sigma",\args -> case args of+ [s,w] -> \ (MS _ _ l) -> sigma (read s) (read w) $ fromIntegral l+ _ -> error "terminator 'sigma' needs two arguments")+ ]++------------------------------------------------------------++setMutator :: (String,[String]) -> Opts -> IO Opts+setMutator (mut,args) (Opts p m t) = return (Opts p m' t)+ where m' = maybe (error $ notfound "mutator" mutators mut)+ (\f -> f args : m) (lookup mut mutators)+++mutators :: [(String,[String]->Mutator)]+mutators = [] -- ("uniform",...)]++------------------------------------------------------------++setModel :: (String,[String]) -> Opts -> IO Opts+setModel (md,args) _ = return (maybe (error $ notfound "model" models md)+ (\f -> case f args of (p,m,t) -> Opts {primer=p,errmod=m,term=t})+ (lookup md models))++models :: [(String,[String]->Model)]+models = [("sanger",sanger)+ ,("stest",\[] -> sanger ["1","1"])+ ,("454",r454)+ ]
+ src/R454.hs view
@@ -0,0 +1,38 @@+{- Modelling 454 sequencing+ (or Roche-something, as it is now called)++ From Marguiles M. et al:+ insertion error rate 1.67%+ deletion error rate 1.60%+ substitution error rate 0.68%++ No quantification of homopolymer length errors.++ Primer is uniform, this is even more so for 454 than for Sanger+ (but the paper has a more elaborate stats to be incorporated: TODO)+-}++module R454 where++import UnfoldMut+import Bio.Sequence.SeqData (Offset, Sequence)++r454_dist :: Offset -> Prob+r454_dist = (*0.02) . uniform -- 2% base rate of error++r454_hp :: MState -> Prob+r454_hp = (*0.2) . homopolymer -- 0-40% chance for 0..20+ equal bases++basic, hp :: [Mutator]+basic = mkmut r454_dist (subst "ACGTN" ++ ins "ACGTN" ++ (concat $ replicate 5 del))++-- | Depending on homopolymer length, duplicate or delete+-- hmm..ugly+hp = [\m -> let (s,m') = (head dup) m in (r454_hp m,s,m')+ ,\m -> let (s,m') = (head del) m in (r454_hp m,s,m')]++terminator :: MState -> Prob+terminator (MS _ _ l) = gradient 90 150 $ fromIntegral l++r454 :: [String] -> ([Sequence] -> IO [MState], [Mutator], MState -> Prob)+r454 = \[n,d] -> (p_uniform [n,d],basic++hp,terminator)
+ src/Sanger.hs view
@@ -0,0 +1,33 @@+{- Implement Sanger type sequencing+ + Random primer position+ Terminate after ...850±50 nucleotides(?)+ initial errors, uniform errors, terminal errors+-} ++module Sanger (sanger) where+import UnfoldMut++-- Sanger error distribution. High initial rate,+-- long stretch of high quality, then degradation at the end+sanger_dist :: Distribution+sanger_dist = combine [ (*0.05) . uniform -- 0.5% error rate in HQ region+ , (*0.30) . gradient 25 0 -- high initial error rate+ , (*0.40) . sigma 700 25] -- terminal degradation++mods :: [Char]+mods = "ACGTN"++-- | Relationship of 3:1 between substitutions and indels+s_mutate :: [Mutator]+s_mutate = mkmut ((*0.60) . sanger_dist) (subst mods) -- 3/4 subst+ ++ mkmut ((*0.10) . sanger_dist) (ins mods) -- 1/8 insert + ++ mkmut ((*0.10) . sanger_dist) del -- 1/8 delete+ ++ mkmut ((*0.30) . sigma 700 25) dup -- terminal duplications+ ++ mkmut ((*0.20) . sigma 700 25) del -- ...and deletions++terminator :: Terminator+terminator (MS _ _ l) = sigma 850 25 $ fromIntegral l++sanger :: [String] -> Model+sanger = \[n,d] -> (p_uniform [n,d],s_mutate,terminator)
+ src/Simseq.hs view
@@ -0,0 +1,31 @@+{-+ Main function for simseq++ Read a set of sequences, apply an error model to generate simulated+ sequences.++ Usage: simseq [options] fasta-file..++ See Options.hs for option details.+-}++module Main where++import Bio.Sequence++import System.IO+import Data.ByteString.Lazy.Char8 (pack)++import UnfoldMut+import Options++main :: IO ()+main = do+ (Opts pri mut trm,ss) <- getOpts+ ms <- pri ss+ let ts = mkseqs $ map (transcribe mut trm) ms+ hWriteFasta stdout ts++mkseqs :: [String] -> [Sequence]+mkseqs = map ($Nothing) . zipWith Seq [pack (">tr_"++show x) | x <- ([1..]::[Int])] . map pack+
+ src/UnfoldMut.hs view
@@ -0,0 +1,160 @@+{- Transcribe from a base sequence using unfoldr++ Probably more complicated than necessary.++ 'prime' selects starting points (in the form of initial state)+ 'transcribe' takes a list of mutators, a terminator, and an+ initial state and generates the transcript.+-}++module UnfoldMut where++import Bio.Sequence hiding ((!))+import Data.Char (toUpper) -- for testing+import Data.Int+import Data.List (sort)+import Data.List (unfoldr)+import Data.Maybe (isJust)+import System.Random+import qualified Bio.Sequence as S++at :: Cursor -> Char+at (s,o,d) | o < 0 || o >= seqlength s = 'X'+ | d == Fwd = (S.!) s o+ | otherwise = compl $ (S.!) s o++type Prob = Double++-- | The cursor tracks position in the sequence (simulating polymerase)+data Direction = Fwd | Rev deriving (Eq,Show)+type Cursor = (Sequence,Offset,Direction)++-- | A primer determines when to initiate transcripts+type Primer = [Sequence] -> IO [MState]++------------------------------------------------------------+-- | A mutator inspects the state, returns the probability of an action,+-- and the action, in the form of a string and a modified state+type Mutator = MState -> (Prob,[Char], MState)+type MutFunc = MState -> ([Char],MState)++combine :: [Distribution] -> Distribution+combine = foldr (\f g -> \x -> f x + g x) (const 0)++mkmut :: Distribution -> [MutFunc] -> [Mutator]+mkmut dist fs =+ map (\f -> \m -> let (s,m') = f m in (t m, s, m')) fs+ where t = (/fromIntegral (length fs)) . dist . tlen+ tlen (MS _ _ l) = fromIntegral l++-- | Uniform probs for substitution and insertion+subst, ins :: String -> [MutFunc]+subst = map (\x -> \m -> ([x],fwd m))+ins = map (\x -> \m -> ([x],m))++-- | Deletion of current base+del :: [MutFunc]+del = [\m -> ([],fwd m)]++-- | Duplicate current base without forwarding+dup :: [MutFunc]+dup = [\m@(MS _ c _) -> ([toUpper $ at c],m)]++------------------------------------------------------------+-- | The terminator examines the state, and the probablity of+-- transcript termination+type Terminator = MState -> Prob++-- | A Model is shorthand for a complete parameter set+type Model = (Primer,[Mutator],Terminator)++-- | The state of the current transcription+data MState = MS StdGen Cursor Int+instance Show MState where show = const "xxx"++------------------------------------------------------------+-- | A distribution is a probability as a function of length+type Distribution = Offset -> Prob++-- | The sigmoid function provides a soft threshold+-- center on s, scale by w+sigma :: Double -> Double -> Distribution+sigma s w = \t -> 1/(1+exp ((s-fromIntegral t)/w))++uniform :: Distribution+uniform = const 1++gradient :: Double -> Double -> Distribution+gradient s e = \t -> max 0 . min 1 $ (fromIntegral t-s)/(e-s)++-- | For 454 (and other?) sequencing, we need a distribution+-- that depends on sequence content+homopolymer :: MState -> Prob+homopolymer (MS _ c@(s,o,_) _) =+ (/20) . fromIntegral . length . takeWhile (== at c) .+ map ((S.!) s) $ [o..min (o+19) (seqlength s-1)]++------------------------------------------------------------+-- | Increments the cursor in the MState+fwd :: MState -> MState+fwd (MS sg (s,o,d) l) = (MS sg (s,if d == Fwd then o+1 else o-1,d) l)++transcribe :: [Mutator] -> Terminator -> MState -> [Char]+transcribe ms t is = concat . unfoldr (sequencer ms t) $ is++-- | walk a list of mutator funcitions, pick one based on p+-- whoa! this is almost :: [Mutator] -> Mutator!+sequencer :: [Mutator] -> Terminator -> MState -> Maybe ([Char],MState)+sequencer ms t =+ \m@(MS g c _) ->+ let (p,g') = randomR (0,1) g+ tp = t m+ new = mutate m (p-tp) ms+ in if p < tp then Nothing+ else update_ms g' $ if isJust new then new+ else Just ([at c], fwd m)++-- | update rndgen and length+update_ms :: StdGen -> Maybe (String,MState) -> Maybe (String,MState)+update_ms g (Just (xs, MS _ c l)) = Just (xs, MS g c (l+length xs))+update_ms _ Nothing = error "this is impossible"++-- | try all mutators in a list+mutate :: MState -> Prob -> [Mutator] -> Maybe (String,MState)+mutate st p (m:ms) = let (t,ch,st') = m st+ p' = p-t+ in if p' < 0 then Just (ch,st')+ else mutate st p' ms+mutate _st _p [] = Nothing++------------------------------------------------------------+-- | A primer with a uniform (poisson) distribution+p_uniform :: [String] -> Primer+p_uniform [n,d] = \ss -> do+ (g,g') <- return . split =<< newStdGen+ let l = sum $ map seqlength ss+ os = sort . take (read n) . randomRs (0,l-1) $ g+ return (gen_ms (read d) g' 0 ss os)+p_uniform _ = error "primer 'uniform' needs two arguments"++-- | Generate MStates for transcript initiation+gen_ms :: Double -> StdGen -> Offset -> [Sequence] -> [Offset] -> [MState]+gen_ms d g c (s:ss) (p:ps)+ | p-c >= seqlength s = gen_ms d g (c+seqlength s) ss (p:ps)+ | otherwise = let (g1,g') = split g+ (z,g2) = randomR (0,1) g1+ in MS g2 (s,p-c,if z<d then Fwd else Rev) 0+ : gen_ms d g' c (s:ss) ps+gen_ms _ _ _ _ [] = []+gen_ms _ _ _ [] _ = error "Impossible: gen_ms ran out of sequences?!"++-- ripped from http://darcs.haskell.org/binary/tests/QuickCheckUtils.hs++instance Random Int64 where+ randomR = integralRandomR+ random = randomR (minBound,maxBound)++integralRandomR :: (Integral a, RandomGen g) => (a,a) -> g -> (a,g)+integralRandomR (a,b) g = case randomR (fromIntegral a :: Integer,+ fromIntegral b :: Integer) g of+ (x,g') -> (fromIntegral x, g')