Probability

12 problems from Rosalind — Bioinformatics Stronghold. Press Run on any block to execute it in your browser.

# Rosalind: IPRB — Mendel's First Law
# https://rosalind.info/problems/iprb/
#
# Given: Positive integers k (homozygous dominant), m (heterozygous) and
# n (homozygous recessive) organisms.
# Return: The probability that two randomly chosen mating organisms produce an
# individual with a dominant allele.

let k = 2.0
let m = 2.0
let n = 2.0
let total = k + m + n

# Work out the recessive probability and subtract: only mm, mn and nn pairings
# can produce a recessive offspring.
# A continued expression keeps its operator at the end of the line.
let mm = (m / total) * ((m - 1.0) / (total - 1.0)) * 0.25
let mn = (m / total) * (n / (total - 1.0)) * 0.5
let nm = (n / total) * (m / (total - 1.0)) * 0.5
let nn = (n / total) * ((n - 1.0) / (total - 1.0))
let p_recessive = mm + mn + nm + nn

let result = 1.0 - p_recessive

println("Result:   " + str(round(result, 5)))
println("Expected: 0.78333")

fn test_iprb_dominant_probability() {
    assert round(result, 5) == 0.78333, "IPRB: got " + str(result)
}

# Rosalind: IEV — Calculating Expected Offspring
# https://rosalind.info/problems/iev/
#
# Given: Six integers, the number of couples with each genotype pairing.
# Return: The expected number of offspring displaying the dominant phenotype,
# assuming every couple has exactly two offspring.

let couples = [1, 0, 0, 1, 0, 1]

# Probability a single offspring shows the dominant phenotype, per pairing:
# AA-AA, AA-Aa, AA-aa are certain; Aa-Aa is 3/4; Aa-aa is 1/2; aa-aa is 0.
let dominant_probability = [1.0, 1.0, 1.0, 0.75, 0.5, 0.0]
let offspring_per_couple = 2.0

let result = range(0, 6)
    |> map(|i| float(couples[i]) * dominant_probability[i] * offspring_per_couple)
    |> sum()

println("Result:   " + str(result))
println("Expected: 3.5")

fn test_iev_expected_offspring() {
    assert result == 3.5, "IEV: got " + str(result)
}

# Rosalind: LIA — Independent Alleles
# https://rosalind.info/problems/lia/
#
# Given: Integers k <= 7 and N <= 2^k. Tom is AaBb; every organism mates with
# an AaBb partner and has two children.
# Return: The probability that at least N organisms in generation k are AaBb.

let k = 2
let n = 1

# Mendel's second law makes the two loci independent, so each child is AaBb
# with probability 1/4 regardless of generation. Generation k holds 2^k
# organisms, so the count is Binomial(2^k, 1/4).
let population = int(pow(2.0, float(k)))
let p = 0.25

fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) }
fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) }

# P(at least N) = 1 - P(fewer than N).
let below = range(0, n)
    |> map(|i| float(choose(population, i)) * pow(p, float(i)) * pow(1.0 - p, float(population - i)))
    |> sum()

let result = 1.0 - below

println("Result:   " + str(round(result, 3)))
println("Expected: 0.684")

fn test_lia_independent_alleles() {
    assert round(result, 3) == 0.684, "LIA: got " + str(result)
}

# Rosalind: SEXL — Sex-Linked Inheritance
# https://rosalind.info/problems/sexl/
#
# Given: An array A of allele frequencies for a gene on the X chromosome.
# Return: For each frequency, the probability that a randomly selected female
# is a carrier of the recessive allele.

let freqs = [0.1, 0.5, 0.8]

# Males carry one X, so the recessive allele frequency is the given value. A
# female is heterozygous with probability 2q(1-q).
let results = freqs |> map(|q| 2.0 * q * (1.0 - q))
let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: 0.18 0.5 0.32")

fn test_sexl_carrier_probabilities() {
    assert round(results[0], 3) == 0.18, "SEXL[0]: got " + str(results[0])
    assert round(results[1], 3) == 0.5, "SEXL[1]: got " + str(results[1])
    assert round(results[2], 3) == 0.32, "SEXL[2]: got " + str(results[2])
}

# Rosalind: AFRQ — Counting Disease Carriers
# https://rosalind.info/problems/afrq/
#
# Given: An array A giving, for each factor, the proportion of the population
# that is homozygous recessive.
# Return: The proportion carrying at least one copy of the recessive allele,
# assuming Hardy-Weinberg equilibrium.

let homozygous_recessive = [0.1, 0.25, 0.5]

# Under Hardy-Weinberg, q^2 is given, so q = sqrt(q^2) and the carriers plus
# the affected are everyone who is not homozygous dominant: 1 - (1 - q)^2.
let results = homozygous_recessive |> map(|q2| {
    let q = sqrt(q2)
    1.0 - (1.0 - q) * (1.0 - q)
})
let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: 0.532 0.75 0.914")

fn test_afrq_carrier_proportions() {
    assert round(results[0], 3) == 0.532, "AFRQ[0]: got " + str(results[0])
    assert round(results[1], 3) == 0.75, "AFRQ[1]: got " + str(results[1])
    assert round(results[2], 3) == 0.914, "AFRQ[2]: got " + str(results[2])
}

# Rosalind: PROB — Introduction to Random Strings
# https://rosalind.info/problems/prob/
#
# Given: A DNA string s and an array A of GC contents.
# Return: For each GC content, the common logarithm of the probability that a
# random string with that GC content matches s exactly.

let s = "ACGATACAA"
let gc_contents = [0.129, 0.287, 0.423, 0.476, 0.641, 0.742, 0.783]

# At GC content x, each of G and C appears with probability x/2 and each of
# A and T with probability (1-x)/2. Sum the logs rather than multiplying the
# probabilities, which would underflow on a long string.
let results = gc_contents |> map(|x| {
    let gc_probability = x / 2.0
    let at_probability = (1.0 - x) / 2.0
    range(0, len(s))
        |> map(|i| {
            let base = substr(s, i, 1)
            let is_gc = base == "G" or base == "C"
            log10(if is_gc then gc_probability else at_probability)
        })
        |> sum()
})

let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: -5.737 -5.217 -5.263 -5.36 -5.958 -6.628 -7.009")

fn test_prob_random_string_logs() {
    assert round(results[0], 3) == -5.737, "PROB[0]: got " + str(results[0])
    assert round(results[3], 3) == -5.36, "PROB[3]: got " + str(results[3])
    assert round(results[6], 3) == -7.009, "PROB[6]: got " + str(results[6])
}

# Rosalind: EVAL — Expected Number of Restriction Sites
# https://rosalind.info/problems/eval/
#
# Given: A positive integer n, a DNA string s, and an array A of GC contents.
# Return: For each GC content, the expected number of times s appears as a
# substring of a random string of length n.

let n = 10
let s = "AG"
let gc_contents = [0.25, 0.5, 0.75]

# Each of the (n - |s| + 1) starting positions carries the same probability of
# spelling s, and expectation is additive whether or not those events overlap.
let positions = n - len(s) + 1

let results = gc_contents |> map(|x| {
    # No product() builtin, so fold the multiplication.
    let probability = range(0, len(s)) |> reduce(|acc, i| {
        let base = substr(s, i, 1)
        let is_gc = base == "G" or base == "C"
        acc * (if is_gc then x / 2.0 else (1.0 - x) / 2.0)
    }, 1.0)
    float(positions) * probability
})

let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: 0.422 0.563 0.422")

fn test_eval_expected_sites() {
    assert round(results[0], 3) == 0.422, "EVAL[0]: got " + str(results[0])
    assert round(results[1], 3) == 0.563, "EVAL[1]: got " + str(results[1])
    assert round(results[2], 3) == 0.422, "EVAL[2]: got " + str(results[2])
}

# Rosalind: EBIN — Wright-Fisher's Expected Behavior
# https://rosalind.info/problems/ebin/
#
# Given: A positive integer n and an array P of probabilities.
# Return: For each probability, the expected number of successes in n
# independent Bernoulli trials.

let n = 17
let probabilities = [0.1, 0.2, 0.3]

# The expectation of a binomial is n*p, whether or not the trials interact —
# expectation is linear, so no summation over outcomes is needed.
let results = probabilities |> map(|p| float(n) * p)
let formatted = results |> map(|r| str(round(r, 3))) |> join(" ")

println("Result:   " + formatted)
println("Expected: 1.7 3.4 5.1")

fn test_ebin_expected_successes() {
    assert round(results[0], 3) == 1.7, "EBIN[0]: got " + str(results[0])
    assert round(results[1], 3) == 3.4, "EBIN[1]: got " + str(results[1])
    assert round(results[2], 3) == 5.1, "EBIN[2]: got " + str(results[2])
}

# Rosalind: INDC — Independent Segregation of Chromosomes
# https://rosalind.info/problems/indc/
#
# Given: A positive integer n (n <= 50).
# Return: An array A of length 2n where A[k] is the common logarithm of the
# probability that at least k of Tom's 2n chromosomes came from his father.

let n = 5
let trials = 2 * n

# Each chromosome is inherited from either parent with probability 1/2 and
# independently of the others, so the count is Binomial(2n, 1/2).
fn binomial_row(width) {
    range(0, width) |> reduce(
        |row, _| {
            let extended = concat([0], row)
            range(0, len(row) + 1) |> map(|i| {
                let left = extended[i]
                let right = if i < len(row) then row[i] else 0
                left + right
            })
        },
        [1]
    )
}

let coefficients = binomial_row(trials)
let total = float(sum(coefficients))

# P(at least k) sums the tail from k upwards.
let results = range(1, trials + 1) |> map(|k| {
    let tail = range(k, trials + 1) |> map(|i| float(coefficients[i])) |> sum()
    log10(tail / total)
})

println("Result:   " + (results |> map(|r| str(round(r, 3))) |> join(" ")))
println("Expected: the binomial tail of 10 fair trials, ending log10(1/1024) = -3.010")

fn test_indc_independent_segregation() {
    assert len(results) == trials, "INDC: got " + str(len(results)) + " entries"
    # P(at least 1) = 1 - 1/1024, and P(at least 2n) = 1/1024 exactly.
    assert round(results[0], 4) == -0.0004, "INDC: first was " + str(results[0])
    assert round(results[trials - 1], 3) == -3.01, "INDC: last was " + str(results[trials - 1])
    # The tail must be non-increasing: a larger threshold cannot be likelier.
    let i = 1
    while i < trials {
        assert results[i] <= results[i - 1], "INDC: tail increased at " + str(i)
        i = i + 1
    }
}

# Rosalind: WFMD — The Wright-Fisher Model of Genetic Drift
# https://rosalind.info/problems/wfmd/
#
# Given: Positive integers N, m, g and k.
# Return: The probability that in a population of N diploid individuals
# initially carrying m copies of a dominant allele, at least k copies of the
# recessive allele remain after g generations.

let n = 4
let m = 6
let g = 2
let k = 1

let alleles = 2 * n

fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) }
fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) }

# One generation: given i recessive copies now, the next count is
# Binomial(2N, i / 2N) — every allele is redrawn independently.
fn step(distribution) {
    range(0, alleles + 1) |> map(|j| {
        range(0, alleles + 1) |> map(|i| {
            let p = float(i) / float(alleles)
            distribution[i] * float(choose(alleles, j)) * pow(p, float(j)) * pow(1.0 - p, float(alleles - j))
        }) |> sum()
    })
}

# The population starts with m dominant copies, so 2N - m recessive ones.
let start = range(0, alleles + 1) |> map(|i| if i == alleles - m then 1.0 else 0.0)
let after = range(0, g) |> reduce(|d, _| step(d), start)

let result = range(k, alleles + 1) |> map(|i| after[i]) |> sum()

println("Result:   " + str(round(result, 3)))
println("Expected: 0.772")

fn test_wfmd_genetic_drift() {
    assert round(result, 3) == 0.772, "WFMD: got " + str(result)
}

# Rosalind: RSTR — Matching Random Motifs
# https://rosalind.info/problems/rstr/
#
# Given: A positive integer N, a number x between 0 and 1, and a DNA string s.
# Return: The probability that at least one of N random strings constructed
# with GC content x matches s exactly.

let n = 90
let gc_fraction = 0.6
let s = "ATAGCCGA"

# At GC content x each of G and C appears with probability x/2, and each of A
# and T with (1-x)/2.
let single = range(0, len(s)) |> reduce(|acc, i| {
    let base = substr(s, i, 1)
    let is_gc = base == "G" or base == "C"
    acc * (if is_gc then gc_fraction / 2.0 else (1.0 - gc_fraction) / 2.0)
}, 1.0)

# The N strings are independent, so the complement is the clean way in: the
# chance that none of them matches is (1-p)^N.
let result = 1.0 - pow(1.0 - single, float(n))

println("P(one string matches): " + str(single))
println("Result:   " + str(round(result, 5)))
println("Expected: 0.00117 — s has 4 weak bases and 4 strong ones, so")
println("          p = 0.2^4 x 0.3^4 = 1.296e-05, and 1-(1-p)^90 = 0.001166")

fn test_rstr_random_motif_match() {
    # p is the product of four A/T probabilities and four G/C ones.
    assert abs(single - 0.00001296) < 0.0000001, "RSTR: p = " + str(single)
    assert round(result, 5) == 0.00117, "RSTR: got " + str(result)
    # More strings can only raise the chance of a match.
    let fewer = 1.0 - pow(1.0 - single, 10.0)
    assert fewer < result, "RSTR: fewer trials did not lower the probability"
}

# Rosalind: FOUN — The Founder Effect and Genetic Drift
# https://rosalind.info/problems/foun/
#
# Given: Positive integers N and m, and an array A of m allele counts.
# Return: An N x m matrix whose (i, j) entry is the common logarithm of the
# probability that the recessive allele is lost entirely after i+1 generations,
# starting from A[j] copies in a population of N diploid individuals.

let n = 4
let counts = [0, 1, 2]
let alleles = 2 * n

fn factorial(x) { range(1, x + 1) |> reduce(|acc, i| acc * i, 1) }
fn choose(total, taken) { factorial(total) / (factorial(taken) * factorial(total - taken)) }

# One Wright-Fisher generation: from i copies, the next count is
# Binomial(2N, i / 2N) — every allele is drawn afresh.
fn step(distribution) {
    range(0, alleles + 1) |> map(|j| {
        range(0, alleles + 1) |> map(|i| {
            let p = float(i) / float(alleles)
            distribution[i] * float(choose(alleles, j)) * pow(p, float(j)) * pow(1.0 - p, float(alleles - j))
        }) |> sum()
    })
}

# Each row is one further generation; each column one starting count.
let rows = counts |> map(|start| {
    let distribution = range(0, alleles + 1) |> map(|i| if i == start then 1.0 else 0.0)
    range(0, n) |> map(|_| {
        distribution = step(distribution)
        log10(distribution[0])
    })
})

println("Result (rows are starting counts, columns are generations):")
range(0, len(counts)) |> each(|j| {
    println("  from " + str(counts[j]) + ": " + (rows[j] |> map(|v| str(round(v, 6))) |> join(" ")))
})
println("")
println("Expected, generation 1: 0.0  -0.463935  -0.999510")
println("  from 1 copy:  (7/8)^8 = 0.343609, log10 = -0.463936")
println("  from 2 copies: (6/8)^8 = 0.100113, log10 = -0.999510")

fn test_foun_founder_effect() {
    # An allele already absent stays absent, so the probability is 1 and its
    # logarithm 0 in every generation.
    # Bound first: `x |> count_if(f) == n` reads the comparison as count_if's
    # second argument rather than as a test of its result.
    let certain = rows[0] |> count_if(|v| v == 0.0)
    assert certain == n, "FOUN: the zero-copy column is not certain"
    assert round(rows[1][0], 6) == -0.463936, "FOUN: one copy, one generation = " + str(rows[1][0])
    assert round(rows[2][0], 6) == -0.99951, "FOUN: two copies, one generation = " + str(rows[2][0])
    # Loss can only become more likely with time, so each row is non-decreasing.
    let i = 1
    while i < n {
        assert rows[1][i] >= rows[1][i - 1], "FOUN: loss became less likely at generation " + str(i)
        i = i + 1
    }
}