calcFST() errors, length undefined, with a windowed FST calculation

[bhaller]

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.