1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
|
import Probability.Distribution.PairwiseAlignment
import Bio.Alignment
import Bio.Alphabet
import Effect
import IModel
import MCMC
import Probability
import SModel
import Tree
import Tree.Newick
import Data.IntMap (IntMap)
import qualified Data.IntMap as IntMap
import qualified Data.JSON as J
import Probability.Logger
import qualified Data.Text.IO as T
import Data.Array
import System.Environment
import System.IO
-- Read the tree from a file
getTree args = do
let (filename:_) = args
rtree <- dropInternalLabels <$> readBranchLengthTree filename
return rtree
getStartLength args = read startLength' :: Int
where (_:startLength':_) = args
-- Sample substitution model parameters and define the substitution model
getSmodel = do
kappaPur <- sample $ logNormal (log 2) 0.25
kappaPyr <- sample $ logNormal (log 2) 0.25
pi <- sample $ symmetric_dirichlet_on (letters dna) 2
return $ tn93' dna kappaPur kappaPyr pi
model rootedTree startLength = do
smodel <- getSmodel
-- Sample ancestral sequence STATES
sequences <- sample $ phyloCTMC rootedTree startLength smodel 1
-- Return the AlignedCharacterData
return sequences
main = do
-- 1. Read the tree and get the starting sequence length
args <- getArgs
rootedTree <- getTree args
let startLength = getStartLength args
alignedSequences <- runRandomLazy $ model rootedTree startLength
T.putStr $ toFasta $ alignedSequences
|
