calcFST() errors, length undefined, with a windowed FST calculation #548
Description
Untested by unit tests, apparently, whoops. Just needs to be chromosome.length. Test model:
// set up a simple neutral simulation
initialize() {
initializeMutationRate(1e-7);
// m1 mutation type: neutral initializeMutationType("m1", 0.5, "f", 0.0); // g1 genomic element type: uses m1 for all mutations initializeGenomicElementType("g1", m1, 1.0); // uniform chromosome of length 100 kb with uniform recombination initializeGenomicElement(g1, 0, 99999); initializeRecombinationRate(1e-8); }
// create a population of 500 individuals
1 early() {
sim.addSubpop("p1", 500);
sim.addSubpop("p2", 500);
}
2000 late() {
myCalcFST(p1.haplosomes, p2.haplosomes, NULL, start=100, end=50000);
}
function (float$)myCalcFST(object haplosomes1, object haplosomes2, [No muts = NULL], [Ni$ start = NULL], [Ni$ end = NULL])
{
if ((haplosomes1.length() == 0) | (haplosomes2.length() == 0))
stop("ERROR (calcFST): haplosomes1 and haplosomes2 must both be non-empty.");
species = haplosomes1[0].individual.subpopulation.species; if (community.allSpecies.length() > 1) { if (any(c(haplosomes1, haplosomes2).individual.subpopulation.species != species)) stop("ERROR (calcFST): all haplosomes must belong to the same species."); if (!isNULL(muts)) if (any(muts.mutationType.species != species)) stop("ERROR (calcFST): all mutations must belong to the same species as the haplosomes."); } chromosome = haplosomes1[0].chromosome; if (species.chromosomes.length() > 1) { if (any(c(haplosomes1, haplosomes2).chromosome != chromosome)) stop("ERROR (calcFST): all haplosomes must be associated with the same chromosome."); if (!isNULL(muts)) if (any(muts.chromosome != chromosome)) stop("ERROR (calcFST): all mutations must be associated with the same chromosome as the haplosomes."); } if (isNULL(muts)) muts = species.subsetMutations(chromosome=chromosome); // handle windowing if (!isNULL(start) & !isNULL(end)) { if (start > end) stop("ERROR (calcFST): start must be less than or equal to end."); if ((start < 0) | (end >= chromosome.length)) stop("ERROR (calcFST): start and end must be within the bounds of the focal chromosome"); mpos = muts.position; muts = muts[(mpos >= start) & (mpos <= end)]; } else if (!isNULL(start) | !isNULL(end)) { stop("ERROR (calcFST): start and end must both be NULL or both be non-NULL."); } // do the calculation; if the FST is undefined, return NULL p1_p = haplosomes1.mutationFrequenciesInHaplosomes(muts); p2_p = haplosomes2.mutationFrequenciesInHaplosomes(muts); mean_p = (p1_p + p2_p) / 2.0; H_t = 2.0 * mean_p * (1.0 - mean_p); H_s = p1_p * (1.0 - p1_p) + p2_p * (1.0 - p2_p); mean_H_t = mean(H_t); if (isNULL(mean_H_t)) // occurs if muts is zero-length return NAN; if (mean_H_t == 0) // occurs if muts is not zero-length but all frequencies are zero return NAN; fst = 1.0 - mean(H_s) / mean_H_t; return fst; }
Found by Mingming Zhang in the Helsinki workshop.