Haskell syntax

It probably helps to know that

• f x and f $ x both mean the function f applied to x
• let x = y is a simple assignment (no side-effects)
• do...x <- y... performs some kind of action and assigns the result to x.
  The action could be random sampling, an IO operation, etc, depending on the context.
• map f [x1,x2,...] applies the function f to every element of the list.
  The result looks like [f x1, f x2, ...].
  It is used instead of for-loops.

For quick introductions to Haskell syntax you might want to take a look at the short interactive tutorial at tryhaskell.org or the quick tour of Haskell syntax.

Linear regression

Here is a short program that performs linear regression. Here the goal is to find a line f(x)=a*x+b that best predicts y[i] from x[i]. The data ys gives y[i] at each location x[i] in xs.

module LinearRegression where

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
  xy_data <- readTable "xy.csv"

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

  return $ model xs ys

• sampling from a distribution looks like b <- prior $ normal 0 1.
  (This specifies a prior term.)
• observing data from a distribution looks like observe data $ distribution.
  (This specifies a likelihood term.)
• defining a function looks like let f x = b*x + a.
  (This defines the best fit line.)
• logging parameters is done by the code return ["b" %=% b, "a" %=% a, ...].
  (This writes a corresponding JSON object each MCMC iteration.)

Note that for each x, the distribution of y(x) is normal (f x) sigma. The term f x is the location predicted by the best-fit line, but there is a distribution because the observed point may not fall exactly on the line.

You can find this file in bali-phy/tests/prob_prog/regression/ and run it as bali-phy -m LinearRegression.hs --iter=1000.

Tree and alignment inference

Here is a short program that infers the tree and alignment from a FASTA file given on the command line.

module Model where

import           Probability
import           Bio.Alignment
import           Bio.Alphabet
import           Bio.Sequence
import           Tree
import           Tree.Newick
import           SModel
import           IModel
import           System.Environment  -- for getArgs

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 = get_sequence_lengths seq_data

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

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

    -- Substitution model
    freqs  <- prior $ symmetric_dirichlet_on ["A", "C", "G", "T"] 1
    kappa1 <- prior $ log_normal 0 1
    kappa2 <- prior $ log_normal 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" %=% write_newick tree
        , "log(indel_rate)" %=% log indel_rate
        , "mean_length" %=% mean_length
        , "kappa1" %=% kappa1
        , "kappa2" %=% kappa2
        , "frequencies" %=% freqs
        , "scale" %=% scale
        , "|T|" %=% tree_length tree
        , "scale*|T|" %=% tree_length tree * scale
        , "|A|" %=% alignmentLength alignment
        ]

main = do
    [filename] <- getArgs

    seq_data <- mkUnalignedCharacterData dna <$> load_sequences filename

    return $ model seq_data

You can find this file in bali-phy/tests/prob_prog/infer_tree/1/ and run it as bali-phy -m Model.hs 5d-muscle.fasta --iter=1000.

Tree and alignment inference under a dN/dS-across-sites model

Some of the value of this framework can be seen by showing how we can do alignment inference under a more complicated substitution model. Here we have factored the substitution model and its prior out into a function gtr_m7_model.

This uses a model of positive selection that (i) is based on GTR, not HKY and (ii) allows dNdS to have a Beta distribution across sites. This is based on the Haskell generated by the command

bali-phy bglobin.fasta --test --smodel 'function(w: gtr +> x3 +> dNdS(w)) +> m7'

module M7 where

import Bio.Alignment
import Bio.Alphabet
import IModel
import Probability
import SModel
import System.Environment -- for getArgs
import Tree
import Tree.Newick

gtr_m7_model codons = do
    -- GTR model parameters
    let nucs = getNucleotides codons
    sym <- prior $ symmetric_dirichlet_on (letter_pair_names nucs) 1
    pi <- prior $ symmetric_dirichlet_on (letters nucs) 1
    let gtr_model = gtr' sym pi nucs

    -- Positive selection model based on GTR
    let pos_sel_model w = gtr_model +> x3 codons +> dNdS w

    -- M7 model parameters
    mu <- prior $ uniform 0 1 -- mean of dN/dS
    gamma <- prior $ beta 1 10 -- spread of dN/dS
    let m7_model = pos_sel_model +> m7 mu gamma 4

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

    return (m7_model, loggers)

branch_length_dist topology b = gamma 0.5 (2 / fromIntegral n)
  where
    n = numBranches topology

model sequenceData = do
    let the_codons = codons dna standard_code
        taxa = getTaxa sequenceData

    -- Tree
    scale <- prior $ gamma 0.5 2
    tree <- prior $ uniform_labelled_tree taxa branch_length_dist

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

    -- Substitution model
    (m7_model, log_m7_smodel) <- gtr_m7_model the_codons

    -- Alignment
    let sequence_lengths = get_sequence_lengths sequenceData
    alignment <- prior $ phyloAlignment tree imodel scale sequence_lengths

    -- Observation
    observe sequenceData $ phyloCTMC tree alignment m7_model scale

    let alignment_length = alignmentLength alignment
        num_indels = totalNumIndels alignment
        total_length_indels = totalLengthIndels alignment

    return
        [ "tree" %=% write_newick tree
        , "scale" %=% scale
        , "S1" %>% log_m7_smodel
        , "|T|" %=% tree_length tree
        , "scale*|T|" %=% scale * tree_length tree
        , "|A|" %=% alignment_length
        , "#indels" %=% num_indels
        , "|indels|" %=% total_length_indels
        ]

main = do
    [filename] <- getArgs

    sequenceData <- mkUnalignedCharacterData (codons dna standard_code) <$> load_sequences filename

    return $ model sequenceData

You can run it as bali-phy --iter=1000 -m M7 bglobin.fasta.