simuMigration.py

Given subpopulation size, number of subpopulations, the migration rate m, generations to evolve, this script shows changes of allele frequencies among subpopulations at single locus, which results from the effect of migration. All subpopulations have equal number of individuals. Migration is depicted in the form of a N*N matrix, which is consisted with identical non-zero values at non-diagonal positions and zero on the diagonal line, where N is the number of subpopulations. Sum of values of elements on row i (i = 1, 2, , N) of the migration matrix equals to the migration rate, thus it ensures that every individual in subpopulation i has probability m to become a migrant, and such a migrant has equal possibility of entering any other subpopulation. Initialization of allele frequencies among subpopulations is to use a list of numbers, which have counts equal to number of subpopulations and been set evenly distributed between 0 and 1, including 0 and 1. Therefore, the average allele frequency named theoretical value that all subpopulations share after a long term of evolution will be 0.5. The main part of this function is shown below. The full Python code is available for download.

def simuMigration(subPopSize, numOfSubPops, m, generations):
    '''Simulate the change of allele frequencies among subpopulations as a result of migration.'''
    # diploid population, one chromosome with 1 locus
    # random mating with sex
    pop = Population(size=[subPopSize]*numOfSubPops, loci=[1], infoFields=['migrate_to'])
    # set initial allele frequencies to each subpopulation
    # Rule of initialization is to use a list of numbers, which have counts equal to number of
    # subpopulations and been set evenly distributed between 0 and 1.
    # Therefore, the average allele frequency named theoretical value among all subpopulations will be 0.5
    for i in range(numOfSubPops):
        initGenotype(pop, freq=[i*1./(numOfSubPops-1), 1 - i*1./(numOfSubPops-1)], subPops=[i])

    # check if plot
    if useRPy:
        plotter = VarPlotter('[0.5] + [subPop[i]["alleleFreq"][0][0] for i in range(%d)]' % numOfSubPops,
            ylim=[0, 1], update=generations-1, legend=['Theoretical'] + ['Allele freq at subpop %d' % x for x in range(numOfSubPops)],
	    xlab="generation", ylab="alleleFrequency", saveAs='migration.png')
    else:
        plotter = NoneOp()

    a=[]
    # set migration rate matrix a[]. Within each row i of matrix a[], all elements have the same value which is
    # m divided by number of subpopulations minus one, except the diagonal element, whose value is set to be 0.
    # This setting ensures that every individual in subpopulation i has probability m to become the Migrator,
    # and such a migrator.has equal possibility of entering any other subpopulation. 
    for i in range(numOfSubPops):
        b=[];
        for j in range(numOfSubPops):
            if j==i:
                b.append(0)
            else:
                b.append(m/(numOfSubPops-1))
        a.append(b)
    # if number of generations is smaller than 200, step will be set to 10 generations, otherwise it will be 20.
    if generations <= 200:
        s = 10
    else:
        s = 20
    pop.evolve(
        initOps = InitSex(),
        preOps = [
            Stat(alleleFreq=[0], vars=['alleleFreq_sp']), 
            PyEval(r'"Frequency at generation %d: %s\n" % (gen, ", ".join(["%.2f" % x for x in freq]))',
                stmts = 'freq = [subPop[x]["alleleFreq"][0][0] for x in range(%i)]' % numOfSubPops, step = s),
        ],
        matingScheme = RandomMating(),
        postOps = [
            Migrator(rate = a),
            plotter,
            ],
        gen = generations,
        )

This script will produce figures such as

Leave a Comment

Comment
Author
Enter code 977