Probabilistic programming in BAli-Phy

Use model expressions or Haskell programs to describe Bayesian models and run inference with BAli-Phy. You can combine existing evolutionary models or define new models with random trees, changing dependencies, and unknown numbers of parameters.

BAli-Phy implements a universal probabilistic programming language: a model defines a distribution over program execution histories. Random choices can determine which branches of a program are taken and which further random variables are created. MCMC samples these histories conditional on the observed data.

  • Modular models: combine smaller model components in new analyses.
  • Changing structure: let random choices determine dependencies, such as which node is a parent in a random tree.
  • Variable dimensions: infer an unknown number of components without writing reversible-jump proposals yourself.
  • Random data structures: work with trees, lists, and user-defined types as well as numeric parameters.

For the conceptual explanation, see Programs and inference: execution histories, conditioning on data, and how BAli-Phy updates a model efficiently.

From model expressions to programs

Command-line model expressions provide a compact way to construct models using BAli-Phy’s probabilistic programming capabilities. For example, this command combines the HKY85 substitution model with gamma-distributed variation in rates across sites:

bali-phy sequences.fasta -S 'HKY85 +> ASRV.Gamma'

The command estimates an alignment and tree from the sequences in sequences.fasta. You can also evaluate a model expression and display its value:

bali-phy print HKY85 -A DNA

Here, the alphabet is supplied because HKY85 depends on it. See the command reference and user's guide for more about using model expressions.

An unknown number of rate categories

Model expressions can also describe models whose number of parameters is inferred from the data. For example, ASRV.Free estimates the rates and weights of rate categories. Giving its category count n a prior distribution lets inference explore different numbers of categories:

bali-phy sequences.fasta -S 'HKY85 +> ASRV.Free(n=1 + ~Geometric(0.5))' --log-format=json

The ~ means “sample from”: ~Geometric(0.5) samples a nonnegative integer, and adding one makes the category count positive. During MCMC, the number of categories can change along with their rates and weights.

Use JSON logging for this model: unlike TSV, it can record the changing number of category-specific values.

Writing a program directly

BAli-Phy translates model expressions into Haskell, which it uses to represent probabilistic models and perform inference. Writing a Haskell program directly gives you more control over how a model is constructed, what data it observes, and what results it records. The examples below illustrate this approach.

Example programs

Download a program and its data into the same directory, then run the commands there. These are short trial runs, not sufficient to establish convergence. For program options, use bali-phy run PROGRAM.hs --help.

Linear regression

This model describes observations around a line f(x) = b*x + a. The slope b, intercept a, and residual standard deviation sigma have prior distributions. The program reads the x and y columns of xy.csv.

model xs ys = do

    b     <- prior $ normal 0 1

    a     <- prior $ normal 0 1

    sigma <- prior $ exponential 1

    let f x = b * x + a

    observe ys $ independent [ normal (f x) sigma | x <- xs ]

    return ["b" %=% b, "a" %=% a, "sigma" %=% sigma]
  • prior $ normal 0 1 samples a parameter and contributes its prior density.
  • observe ys $ independent [...] contributes the likelihood of the data.
  • Each sampled slope and intercept defines a line. Observations vary around that line according to normal (f x) sigma.
  • The returned named quantities are passed to the program's logging setup, which records them during MCMC.
Complete program: LinearRegression.hs 
module LinearRegression where

import           BAliPhy.Run
import           MCMC (runMCMC)
import           Options.Applicative
import           Probability
import           Data.Frame

model xs ys = do

    b     <- prior $ normal 0 1

    a     <- prior $ normal 0 1

    sigma <- prior $ exponential 1

    let f x = b * x + a

    observe ys $ independent [ normal (f x) sigma | x <- xs ]

    return ["b" %=% b, "a" %=% a, "sigma" %=% sigma]

main = do
  options <- execParser $ modelRunParser "LinearRegression" 200000

  runInfo <- initializeModelRun (runMode options)

  xy_data <- readTable "xy.csv"

  let xs = xy_data $$ "x" :: [Double]
      ys = xy_data $$ "y" :: [Double]

  mcmcState <- makeLoggedMCMCState runInfo (logFormats options) $ model xs ys

  case runInfo of
    TestRun -> printInitialModel (logFormats options) mcmcState
    MCMCRun directory -> do
      reportModelRun (iterations options) (logFormats options) directory
      runMCMC (iterations options) mcmcState

Download LinearRegression.hs and xy.csv. Keep the name xy.csv, which is used by the program.

bali-phy run LinearRegression.hs --iterations=1000 --name=Regression
statreport Regression-1/C1.log.json

Inspect the posterior summaries for a, b, and sigma, together with their sampling diagnostics.

Source in the BAli-Phy repository.

Tree and alignment inference

This program samples a tree and alignment for DNA sequences. It combines a tree prior, an indel process, and a TN93 substitution model.

branch_length_dist topology branch = gamma (1/2) (2/fromIntegral n) where n = numBranches topology

model seq_data = do
    let taxa            = getTaxa seq_data
        tip_seq_lengths = getSequenceLengths seq_data

    -- Tree
    scale <- prior $ gamma (1/2) 2
    tree  <- prior $ uniformLabelledTree'' taxa branch_length_dist

    -- Indel model
    indel_rate   <- prior $ logLaplace (-4) 0.707
    mean_length <- (1 +) <$> sample (exponential 10)
    let imodel = rs07 indel_rate mean_length tree

    -- Substitution model
    freqs  <- prior $ symmetricDirichletOn (letterSet dna) 1
    kappa1 <- prior $ logNormal 0 1
    kappa2 <- prior $ logNormal 0 1
    let tn93_model = tn93' dna kappa1 kappa2 freqs

    -- Alignment
    alignment <- prior $ phyloAlignment tree imodel scale tip_seq_lengths

    -- Observation
    observe seq_data $ phyloCTMC tree alignment tn93_model scale

    return
        [ "tree" %=% writeNewick tree
        , "log(indel_rate)" %=% log indel_rate
        , "mean_length" %=% mean_length
        , "kappa1" %=% kappa1
        , "kappa2" %=% kappa2
        , "frequencies" %=% freqs
        , "scale" %=% scale
        , "|T|" %=% treeLength tree
        , "scale*|T|" %=% treeLength tree * scale
        , "|A|" %=% alignmentLength alignment
        ]
  • The tree prior supplies topology and branch lengths; scale multiplies those lengths for the sequence model.
  • The indel model and tree define the distribution of the alignment through phyloAlignment.
  • phyloCTMC connects the tree, alignment, and substitution model to the observed sequences.
Complete program: InferTreeAlignment.hs 
module Model where

import           BAliPhy.Run
import           MCMC (runMCMC)
import           Options.Applicative
import           Probability
import           Bio.Alignment
import           Bio.Alphabet
import           Bio.Sequence
import           Tree
import           Tree.Newick
import           SModel
import           IModel

branch_length_dist topology branch = gamma (1/2) (2/fromIntegral n) where n = numBranches topology

model seq_data = do
    let taxa            = getTaxa seq_data
        tip_seq_lengths = getSequenceLengths seq_data

    -- Tree
    scale <- prior $ gamma (1/2) 2
    tree  <- prior $ uniformLabelledTree'' taxa branch_length_dist

    -- Indel model
    indel_rate   <- prior $ logLaplace (-4) 0.707
    mean_length <- (1 +) <$> sample (exponential 10)
    let imodel = rs07 indel_rate mean_length tree

    -- Substitution model
    freqs  <- prior $ symmetricDirichletOn (letterSet dna) 1
    kappa1 <- prior $ logNormal 0 1
    kappa2 <- prior $ logNormal 0 1
    let tn93_model = tn93' dna kappa1 kappa2 freqs

    -- Alignment
    alignment <- prior $ phyloAlignment tree imodel scale tip_seq_lengths

    -- Observation
    observe seq_data $ phyloCTMC tree alignment tn93_model scale

    return
        [ "tree" %=% writeNewick tree
        , "log(indel_rate)" %=% log indel_rate
        , "mean_length" %=% mean_length
        , "kappa1" %=% kappa1
        , "kappa2" %=% kappa2
        , "frequencies" %=% freqs
        , "scale" %=% scale
        , "|T|" %=% treeLength tree
        , "scale*|T|" %=% treeLength tree * scale
        , "|A|" %=% alignmentLength alignment
        ]

main = do
    (options, filename) <- execParser $
      modelRunParserWith "Model" 200000 $
        strArgument (metavar "SEQUENCES" <> help "Unaligned DNA sequences")

    runInfo <- initializeModelRun (runMode options)

    seq_data <- mkUnalignedCharacterData dna <$> loadSequences filename

    mcmcState <- makeLoggedMCMCState runInfo (logFormats options) $ model seq_data

    case runInfo of
      TestRun -> printInitialModel (logFormats options) mcmcState
      MCMCRun directory -> do
        reportModelRun (iterations options) (logFormats options) directory
        runMCMC (iterations options) mcmcState

Download InferTreeAlignment.hs and 5d-muscle.fasta.

bali-phy run InferTreeAlignment.hs 5d-muscle.fasta --iterations=1000 --name=TreeAlignment
statreport --ignore=tree TreeAlignment-1/C1.log.json

The JSON log includes sampled trees as text. Exclude that field when summarizing numeric parameters such as tree length, alignment length, and indel rate.

Source in the BAli-Phy repository.

A custom codon model

The function gtr_m7_model combines GTR nucleotide exchangeabilities with a beta distribution of dN/dS (ω) across sites. M7 describes variation in selective constraint: ω lies between 0 and 1.

gtr_m7_model codons = do
    let nucs = getNucleotides codons

    -- GTR model parameters
    sym <- sample $ symmetricDirichletOn (Set.fromList $ letter_pair_names nucs) 1
    pi <- sample $ symmetricDirichletOn (letterSet nucs) 1

    let posSelModel w = gtr' sym pi nucs +> x3 codons +> dNdS w

    -- Independent uniform priors on the beta mean and normalized variance.
    mu <- sample $ uniform 0 1
    v <- sample $ uniform 0 1
    let m7Model = posSelModel +> m7 mu v 4

    let loggers =
            [ "gtr:sym" %=% sym
            , "gtr:pi" %=% pi
            , "m7:mu" %=% mu
            , "m7:v" %=% v
            ]

    return (m7Model, loggers)

The function defines the substitution model and its priors, returning both the model and quantities to log. The complete program uses this model to infer a tree and alignment, with branch-specific indel rates.

Complete program: M7.hs 
module Model where

import BAliPhy.Run
import Bio.Alignment
import Bio.Alphabet
import qualified Data.Set as Set
import IModel
import MCMC
import Options.Applicative
import Probability
import SModel
import SModel.Parsimony
import Tree
import Tree.Newick

gtr_m7_model codons = do
    let nucs = getNucleotides codons

    -- GTR model parameters
    sym <- sample $ symmetricDirichletOn (Set.fromList $ letter_pair_names nucs) 1
    pi <- sample $ symmetricDirichletOn (letterSet nucs) 1

    let posSelModel w = gtr' sym pi nucs +> x3 codons +> dNdS w

    -- Independent uniform priors on the beta mean and normalized variance.
    mu <- sample $ uniform 0 1
    v <- sample $ uniform 0 1
    let m7Model = posSelModel +> m7 mu v 4

    let loggers =
            [ "gtr:sym" %=% sym
            , "gtr:pi" %=% pi
            , "m7:mu" %=% mu
            , "m7:v" %=% v
            ]

    return (m7Model, loggers)

model sequenceData = do
    let taxa = getTaxa sequenceData

    tree <- sample $ uniformLabelledTree taxa (gamma 0.5 (1 / fromIntegral (length taxa)))
    let tlength = treeLength tree

    sigma <- sample $ logLaplace (-3) 1
    indelRates <- fmap (** sigma) <$> sample (iidMap (getUEdgesSet tree) (logNormal 0 1))
    let indelTree = addBranchRates indelRates tree

    scale <- sample $ gamma 0.5 2
    addMove 2 (scaleGroupsSlice [scale] (branchLengths tree))
    addMove 1 (scaleGroupsMH [scale] (branchLengths tree))

    let codons = mkCodons dna (geneticCode "standard")
    (m7_model, log_m7_model) <- gtr_m7_model codons

    rate <- sample $ logLaplace (-4) 0.707
    meanLength <- sample $ shifted_exponential 10 1
    let imodel = IModel.rs07 rate meanLength tree

    let sequenceLengths = getSequenceLengths sequenceData
    (alignment, propertiesA) <- sampleWithProps (phyloAlignment indelTree imodel scale sequenceLengths)
    properties <- observe sequenceData (phyloCTMC tree alignment m7_model scale)

    let alignment_length = alignmentLength alignment
    let num_indels = totalNumIndels alignment
    let total_length_indels = totalLengthIndels alignment
    let prior_A = ln (probability propertiesA)
    let ancStates = prop_anc_cat_states properties
    let ancAlignment = toFasta $ ancestralAlignment tree alignment (getSMap m7_model) codons ancStates
    let substs = parsimony tree (unitCostMatrix codons) (sequenceData, alignment)

    let loggers =
            [ "indelRates:sigma" %=% sigma
            , "S1" %>% log_m7_model
            , "rs07:rate" %=% rate
            , "rs07:mean_length" %=% meanLength
            , "scale" %=% scale
            , "scale*|T|" %=% (scale * tlength)
            , "|A|" %=% alignment_length
            , "#indels" %=% num_indels
            , "|indels|" %=% total_length_indels
            , "#substs" %=% substs
            , "prior_A" %=% prior_A
            ]

    return loggers

main = do
    (options, filename) <-
        execParser $
            modelRunParserWith "Model" 200000 $
                strArgument (metavar "SEQUENCES" <> help "Unaligned coding sequences")

    runInfo <- initializeModelRun (runMode options)

    sequenceData <-
        mkUnalignedCharacterData (mkCodons dna standard_code)
            <$> loadSequences filename

    mcmcState <- makeLoggedMCMCState runInfo (logFormats options) $ model sequenceData

    case runInfo of
        TestRun -> printInitialModel (logFormats options) mcmcState
        MCMCRun directory -> do
            reportModelRun (iterations options) (logFormats options) directory
            runMCMC (iterations options) mcmcState

Download M7.hs and bglobin.fasta.

bali-phy run M7.hs bglobin.fasta --iterations=100 --name=CodonModel
statreport CodonModel-1/C1.log.json

Inspect the beta-distribution parameters and other logged quantities. With this short run, some credible intervals may be unavailable. Directory suffixes increase on subsequent runs; use the directory reported by the program.

Source in the BAli-Phy repository.

Further reading

Programs and inference explains traces, changing model structure, and incremental computation. The user's guide covers built-in evolutionary models and ordinary analyses, and the command reference describes run and print.