Jonatan Blais
Guest

Hi Sergei,
Thank you so much for the bootstrapping program.
We managed to insert it to replace the old permutation routine in the F_ST.bf and it seems that HYPHY reads it through (after correcting a little mistake by replacing “basePartition1,2” by “basePartition2,1” in “resampleP2 = Random (basePartition2,1)”).
However, the output doesn’t seem to make sense. All Fst estimates from bootstrapped samples are lower than the observed one. Normally of course, some should be lower and some higher so a confidence interval could be obtained. I’ve pasted the last part of the F_ST.bf containing the newly inserted bootstrapping code. I suspect there is something wrong in the way the bootstrapping is done but can’t manage to figure out what it could be.
Cheers,
Jon
F_ST_OBS = {{pi_D/(resMx[1]+pi_D),pi_D/(resMx[1]+resMx[2]),1resMx[1]/resMx[0]}};
fprintf (stdout, "\n\nPopulation characterisitcs:", "\nMetapopulation diversity (pi_T) = ", resMx[0], "\nMean subpopulation diversity (pi_S) = ", resMx[1], "\nMean interpopulation diversity (pi_B) = ", resMx[2], "\n\nF_ST\n", "\nHudson, Slatkin and Madison (Genetics 132:583589): ",F_ST_OBS[0], "\nSlatkin (Evolution 47:264279) : ",F_ST_OBS[1], "\nHudson, Boos and Kaplan (Mol Bio Evol 9: 138151) : ",F_ST_OBS[2], "\n"); F_ST_OBS = F_ST_OBS; ChoiceList (resample,"Permutation Test",1,SKIP_NONE, "Skip","Do not perform a permutation test.", "But of course","Randomly allocate sequences into subpopulations and tabulate the distribution of various F_ST statistics.");
if (resample < 1) { return 0; }
sampleCount = 0; while (sampleCount < 1) { fprintf (stdout, "\nHow many permutations would you like to perform?"); fscanf (stdin,"Number",sampleCount); }
F_ST_1 = {sampleCount,2}; F_ST_2 = {sampleCount,2}; F_ST_3 = {sampleCount,2}; step = sampleCount$20;
fprintf (stdout, "Running the simulations...\n");
/* New bootstrapping within groups replacing old permutation code */
saveDM = distanceMatrix;
p1Size = Abs (p1vector); p2Size = Abs (p2vector); basePartition1 = {1,p1Size}; basePartition2 = {1,p2Size};
for (k=0; k<p1Size;k=k+1) { basePartition1[k] = p1vector[k]; } for (k=0; k<p2Size;k=k+1) { basePartition2[k] = p2vector[k]; }
for (sampleID = 0; sampleID < sampleCount; sampleID = sampleID + 1) {
/* repeat the following code for each replicate */ resampleP1 = Random (basePartition1,1); /* resample vector of indices with replacement */ resampleP2 = Random (basePartition2,1);
for (k=0; k<p1Size;k=k+1) /* repopulate distances for population 1 */ { ki = resampleP1[k]; for (k2 = k+1; k2 < p1Size; k2=k2+1) { k2i = resampleP1[k2]; distanceMatrix [k][k2] = saveDM[ki][ki2]; distanceMatrix [k2][k] = saveDM[ki2][ki]; } }
for (k=0; k<p2Size;k=k+1) /* repopulate distances for population 2 */ { ki = resampleP2[k]; for (k2 = k+1; k2 < p2Size; k2=k2+1) { k2i = resampleP2[k2]; distanceMatrix [k][k2] = saveDM[ki][ki2]; distanceMatrix [k2][k] = saveDM[ki2][ki]; } }
for (k=0; k<p1Size;k=k+1) /* repopulate interpopulation distances */ { ki = resampleP1[k]; for (k2 = 0; k2 < p2Size; k2=k2+1) { k2i = resampleP2[k2]; distanceMatrix [k][k2] = saveDM[ki][ki2]; distanceMatrix [k2][k] = saveDM[ki2][ki]; } } /*End of new bootstrapping within groups code*/
resMx = computeCompartentValues (p1_1,p2_1);*/ resMx = computeCompartentValues (p1vector,p2vector); pi_D = resMx[2]resMx[1]; F_ST_1[sampleID][1] = pi_D/(resMx[1]+pi_D); F_ST_2[sampleID][1] = pi_D/(resMx[1]+resMx[2]); F_ST_3[sampleID][1] = 1resMx[1]/resMx[0]; if ((sampleID+1)%step == 0) { fprintf (stdout, Format ((sampleID+1)*100/sampleCount, 6, 2), "% done\n"); } }
F_ST_1 = F_ST_1%1; F_ST_2 = F_ST_2%1; F_ST_3 = F_ST_3%1;
pv = {{sampleCount,sampleCount,sampleCount}};
for (sampleID = 0; sampleID < sampleCount; sampleID = sampleID + 1) { step = sampleID/sampleCount; F_ST_1[sampleID][0] = step; F_ST_2[sampleID][0] = step; F_ST_3[sampleID][0] = step; if (pv[0]==sampleCount) { if (F_ST_1[sampleID][1] > F_ST_OBS[0]) { pv[0] = sampleID; } } if (pv[1]==sampleCount) { if (F_ST_2[sampleID][1] > F_ST_OBS[1]) { pv[1] = sampleID; } } if (pv[2]==sampleCount) { if (F_ST_3[sampleID][1] > F_ST_OBS[2]) { pv[2] = sampleID; } } }
fprintf (stdout, "\n\nProb {Random F_ST > Observed F_ST}\n", "\nHudson, Slatkin and Madison : ",(sampleCountpv[0])/sampleCount, "\nSlatkin : ",(sampleCountpv[1])/sampleCount, "\nHudson, Boos and Kaplan : ",(sampleCountpv[2])/sampleCount, "\n");
labels = {{"Cumulative Weight","F_ST"}};
OpenWindow (CHARTWINDOW,{{"Hudson, Slatkin and...
