risib.py

#!/usr/bin/env python
#
# init and step probabilities for RI by sib-mating
# values returned are log of probabilities
#
# True genotypes BB or DD
# forward direction means BxD for cross

import math # for log


def init(True_gen, is_x_chr=False, forward_direction=True):
    '''log initial probability for RIL by sib-mating'''
    if is_x_chr:
        if forward_direction:
            if True_gen == 'BB':
                return math.log(2.0)-math.log(3.0)
            if True_gen == 'DD':
                return -math.log(3.0)
        else: # backward direction
            if True_gen == BB:
                return math.log(2.0)-math.log(3.0)
            if True_gen == AA:
                return -math.log(3.0)
    else: # autosome
        return -math.log(2.0)

    return null # can't get here


def step(gen_left, gen_right, rec_frac, is_x_chr=False, forward_direction=True):
    '''log transition probability for RIL by sib-mating'''
    if is_x_chr:
        if forward_direction:
            if gen_left == 'BB':
                if gen_right == 'BB':
                    return math.log(1.0 + 2.0*rec_frac) - math.log(1.0 + 4.0*rec_frac)
                if gen_right == 'DD':
                    return math.log(2.0*rec_frac) - math.log(1.0 + 4.0*rec_frac)

            if gen_left == 'DD':
                if gen_right == 'DD':
                    return -math.log(1.0 + 4.0*rec_frac)
                if gen_right == 'BB':
                    return math.log(4.0*rec_frac) - math.log(1.0 + 4.0*rec_frac)
        else: # backward direction
            if gen_left == 'BB':
                if gen_right == 'BB':
                    return -math.log(1.0 + 4.0*rec_frac)
                if gen_right == 'DD':
                    return math.log(4.0*rec_frac) - math.log(1.0 + 4.0*rec_frac)

            if gen_left == 'DD':
                if gen_right == 'DD':
                    return math.log(1.0 + 2.0*rec_frac) - math.log(1.0 + 4.0*rec_frac)
                if gen_right == 'BB':
                    return math.log(2.0*rec_frac) - math.log(1.0 + 4.0*rec_frac)

    else: # autosome
        R = 4.0*rec_frac/(1+6.0*rec_frac)

        if gen_left == gen_right:
            return math.log(1.0-R)
        else:
            return math.log(R)

    return null # can't get here