simuLDDecay.py

This script demonstrates the decay of linkage disequilibrium due to recombination. Given population size, generations to evolve, recombinatin rate, number of replicates and method to measure linkage disequilibrium, the core part of this recipe is shown below. The full Python code is available for download.

def simuLDDecay(popSize, gen, recRate, numRep, method, saveFigure, useRPy):
    '''Simulate the decay of linkage disequilibrium as a result
    of recombination.
    '''
    # diploid population, one chromosome with 2 loci
    # random mating with sex
    simu = Simulator(
        Population(size=popSize, ploidy=2, loci=[2]),
        rep = numRep)

    # get method value used to plot and evolve
    if method=="D'":
        methodplot = "LD_prime[0][1]"
        upperlim = 1
        methodeval = r"'%.4f\t' % LD_prime[0][1]"
    elif method=='R2':
        methodplot = "R2[0][1]"
        upperlim = 1
        methodeval = r"'%.4f\t' % R2[0][1]"
    else:
        methodplot = "LD[0][1]"
        upperlim = 0.25
        methodeval = r"'%.4f\t' % LD[0][1]"

    if useRPy:
        print saveFigure
        plotter = VarPlotter(methodplot, 
            ylim = [0, upperlim], saveAs=saveFigure,
            update = gen - 1, ylab=method,
            main="Decay of Linkage Disequilibrium r=%f" % recRate)
    else:
        plotter = NoneOp()

    simu.evolve(
        # everyone will have the same genotype: 01/10
        initOps = [
            InitSex(),
            InitGenotype(genotype=[0,1,1,0])
        ],
        matingScheme = RandomMating(ops=Recombinator(rates = recRate)), 
        postOps = [
            Stat(alleleFreq=[0], LD=[0, 1]),
            PyEval(methodeval),
            PyOutput('\n', reps=-1),
            plotter
        ],
        gen = gen
    )

Using option --saveFigure=decay.png, this script will produce figures such as

Leave a Comment

Comment
Author
Enter code 187