diff --git a/doc/rarecoal.tex b/doc/rarecoal.tex index 25a1726..4b15079 100644 --- a/doc/rarecoal.tex +++ b/doc/rarecoal.tex @@ -8,6 +8,8 @@ \date{} \maketitle +\section{Introduction} +Rarecoal is a method to estimate a population history model, parameterized by population splits, population sizes and admixture edges, from the joint rare site frequency spectrum. \section{The rarecoal coalescent framework} This model describes a coalescent framework for rare alleles. We define rare alleles roughly by requiring i) the allele count of the derived mutation to be small, typically not larger than 10, and ii) the total number of samples to be much larger, say 100 or more. The idea is to provide a general approach of computing the joint allele frequency spectrum for rare alleles under an arbitrary demographic model under population splits and population size changes. Migration and admixture will be incorporated in the future. diff --git a/src-rarecoal/Main.hs b/src-rarecoal/Main.hs index ca589ee..07d7b64 100644 --- a/src-rarecoal/Main.hs +++ b/src-rarecoal/Main.hs @@ -67,7 +67,7 @@ parseProb :: OP.Parser Command parseProb = CmdProb <$> parseProbOpt parseProbOpt :: OP.Parser ProbOpt -parseProbOpt = ProbOpt <$> parseTheta <*> parseModelDesc <*> parseLinGen <*> parseNVec <*> parseKVec +parseProbOpt = ProbOpt <$> parseTheta <*> parseModelDesc <*> parseLinGen <*> parseBranchnames <*> parseNVec <*> parseKVec where parseNVec = OP.argument OP.auto (OP.metavar "[N1,N2,...]" <> OP.help "number of samples, \ \comma-separated, without spaces, and surrounded by square brackets, e.g. \ @@ -76,6 +76,7 @@ parseProbOpt = ProbOpt <$> parseTheta <*> parseModelDesc <*> parseLinGen <*> par \alleles in each population, same format as for NVec, e.g. [1,2] for allele \ \count 3 shared with one sample from the first and two from the second \ \population.") + parseBranchnames = OP.option (splitOn "," <$> OP.str) (OP.help "string of branch names" <> OP.long "branchnames" <> OP.metavar "Pop1,Pop2,...") parseTheta :: OP.Parser Double parseTheta = OP.option OP.auto $ OP.short 't' <> OP.long "theta" <> OP.hidden <> OP.metavar "FLOAT" diff --git a/src-rarecoal/Prob.hs b/src-rarecoal/Prob.hs index 45e60ed..94b3348 100644 --- a/src-rarecoal/Prob.hs +++ b/src-rarecoal/Prob.hs @@ -9,12 +9,13 @@ data ProbOpt = ProbOpt { prTheta :: Double, prModelDesc :: ModelDesc, prLinGen :: Int, + prBranchnames :: [String], prNvec :: [Int], prKvec :: [Int] } runProb :: ProbOpt -> Script () runProb opts = do - modelSpec <- getModelSpec (prModelDesc opts) [] (prTheta opts) (prLinGen opts) + modelSpec <- getModelSpec (prModelDesc opts) (prBranchnames opts) (prTheta opts) (prLinGen opts) val <- tryRight $ getProb modelSpec (prNvec opts) False (prKvec opts) scriptIO $ print val