Skip to content

Commit

Permalink
Merge pull request #3 from kern-lab/rngMutationIssues
Browse files Browse the repository at this point in the history 
Rng mutation issues
  • Loading branch information
Kern Lab authored Nov 24, 2017
2 parents 86585be + edd8742 commit e4a8647
Show file tree
Hide file tree
Showing 4 changed files with 108 additions and 34 deletions.
2 changes: 1 addition & 1 deletion discoal.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ typedef struct rootedNode
{
struct rootedNode *leftParent, *rightParent, *leftChild, *rightChild;
double time, branchLength, blProb;
float muts[MAXMUTS];
double muts[MAXMUTS];
#ifdef BIG
uint16_t ancSites[MAXSITES];
#else
Expand Down
123 changes: 93 additions & 30 deletions discoalFunctions.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -439,7 +439,7 @@ int recombineAtTimePopn(double cTime, int popn){
int xOver;

aNode = pickNodePopn(popn);
// printf("picked: %p\n",aNode);
// printf("recombine\n");
xOver = ignuin(0,nSites-1);
// printf("xo: %d test1: %d test2: %d\n",xOver,siteBetweenChunks(aNode, xOver),isActive(xOver) );
// printNode(aNode);
Expand Down Expand Up @@ -508,7 +508,7 @@ void geneConversionAtTimePopn(double cTime, int popn){
int tractL;

aNode = pickNodePopn(popn);
// printf("picked: %p\n",aNode);
// printf("GC\n",aNode);
xOver = ignuin(0,nSites);
tractL = (int) ceil(log(genunf(0,1))/log(1.0-(1.0/gcMean)));
// printf("xo: %d test1: %d tractL: %d gcMean: %d\n",xOver,siteBetweenChunks(aNode, xOver),tractL,gcMean );
Expand Down Expand Up @@ -804,7 +804,7 @@ double neutralPhaseGeneralPopNumber(int *bpArray,double startTime, double endTim
for(i=0;i<npops;i++){
mRate[i]=0.0;
}
// printf("currPopSize[0]: %d currPopSize[1]: %d alleleNumber: %d totNodeNumber: %d activeSites: %d cTime: %f\n",popnSizes[0],popnSizes[1],alleleNumber,totNodeNumber, activeSites, cTime);
//printf("currPopSize[0]: %d currPopSize[1]: %d alleleNumber: %d totNodeNumber: %d activeSites: %d cTime: %f\n",popnSizes[0],popnSizes[1],alleleNumber,totNodeNumber, activeSites, cTime);
for(i=0;i<npops;i++){
cRate[i] = popnSizes[i] * (popnSizes[i] - 1) * 0.5 / sizeRatio[i];
rRate[i] = rho * popnSizes[i] * 0.5;// * ((float)activeSites/nSites);
Expand All @@ -818,6 +818,7 @@ double neutralPhaseGeneralPopNumber(int *bpArray,double startTime, double endTim
totGCRate += gcRate[i];
totRate += cRate[i] + rRate[i] + mRate[i] + gcRate[i];
}
//printf("totRate: %f totCRate: %f totRRate: %f \n",totRate,totCRate,totRRate);

//find time of next event
waitTime = genexp(1.0) * (1.0/ totRate);
Expand Down Expand Up @@ -1849,7 +1850,7 @@ void coalesceAtTimePopnSweep(double cTime, int popn, int sp){
void dropMutations(){
int i, j, m;
double p;
float mutSite, error;
double mutSite, error;

//get time and set probs
coaltime = totalTimeInTree();
Expand All @@ -1858,26 +1859,26 @@ void dropMutations(){
//add Mutations
p = allNodes[i]->blProb * theta * 0.5;
if(p>0.0){
m = ignpoi(p);
while(m>0){
mutSite = genunf((float)allNodes[i]->lLim / nSites, (float) allNodes[i]->rLim / nSites);
while(isAncestralHere(allNodes[i],mutSite) != 1){
if (allNodes[i]->lLim == allNodes[i]->rLim){
if (mutSite*nSites < (float)allNodes[i]->lLim)
error = (float)allNodes[i]->lLim-(mutSite*nSites);
else
error = 0.0;
p = allNodes[i]->rLim + (1.0/nSites);
mutSite = genunf(((double)allNodes[i]->lLim + error) / nSites, (double)(p + error) / nSites);

}
else
mutSite = genunf((double)allNodes[i]->lLim / nSites, (double) allNodes[i]->rLim / nSites);
}
// printf("mut: %f\n",mutSite);
addMutation(allNodes[i],mutSite);
m--;
}
m = ignpoi(p);
while(m>0){
mutSite = genunf((float)allNodes[i]->lLim / nSites, (float) allNodes[i]->rLim / nSites);
while(isAncestralHere(allNodes[i],mutSite) != 1){
// printf("in here\n");
if (allNodes[i]->lLim == allNodes[i]->rLim){
if (mutSite*nSites < (float)allNodes[i]->lLim)
error = (float)allNodes[i]->lLim-(mutSite*nSites);
else
error = 0.0;
p = allNodes[i]->rLim + (1.0/nSites);
mutSite = genunf(((double)allNodes[i]->lLim + error) / nSites, (double)(p + error) / nSites);
}
else
mutSite = genunf((double)allNodes[i]->lLim / nSites, (double) allNodes[i]->rLim / nSites);
}
//printf("new mut allele: %d mut: %f\n",i,mutSite);
addMutation(allNodes[i],mutSite);
m--;
}
}
}
// printf("coaltime: %f totM: %d\n",coaltime,tm);
Expand All @@ -1894,6 +1895,49 @@ void dropMutations(){
}
}

void dropMutationsRecurse(){
int i, j, m;
double p;
float mutSite, error;

//get time and set probs
coaltime = totalTimeInTree();
//printf("%f\n",coaltime);
for(i=0;i<totNodeNumber;i++){
//add Mutations
p = allNodes[i]->blProb * theta * 0.5;
if(p>0.0){
m = ignpoi(p);
while(m>0){
mutSite = genunf((float)allNodes[i]->lLim / nSites, (float) allNodes[i]->rLim / nSites);
while(isAncestralHere(allNodes[i],mutSite) != 1){
// printf("in here\n");
if (allNodes[i]->lLim == allNodes[i]->rLim){
if (mutSite*nSites < (float)allNodes[i]->lLim)
error = (float)allNodes[i]->lLim-(mutSite*nSites);
else
error = 0.0;
p = allNodes[i]->rLim + (1.0/nSites);
mutSite = genunf(((double)allNodes[i]->lLim + error) / nSites, (double)(p + error) / nSites);
}
else
mutSite = genunf((double)allNodes[i]->lLim / nSites, (double) allNodes[i]->rLim / nSites);
}
// printf("mut: %f\n",mutSite);
addMutation(allNodes[i],mutSite);
m--;
}
}
}
// printf("coaltime: %f totM: %d\n",coaltime,tm);
//push mutations down tree
for(i=totNodeNumber-1;i>=0;i--){
for(j=0;j<allNodes[i]->mutationNumber;j++){
mutSite = allNodes[i]->muts[j];
recurseTreePushMutation(allNodes[i],mutSite);
}
}
}
//calculates the total time in the tree, then set blProbs for each node
double totalTimeInTree(){
double tTime, siteLength;
Expand All @@ -1910,14 +1954,22 @@ double totalTimeInTree(){
return tTime;
}


void addMutation(rootedNode *aNode, float site){
void recurseTreePushMutation(rootedNode *aNode, float site){
if(aNode->leftChild != NULL && isAncestralHere(aNode->leftChild,site))
recurseTreePushMutation(aNode->leftChild,site);
if(aNode->rightChild != NULL && isAncestralHere(aNode->rightChild,site))
recurseTreePushMutation(aNode->rightChild,site);
if( isLeaf(aNode) && isAncestralHere(aNode,site)){
if(hasMutation(aNode,site)==0) addMutation(aNode, site);
}
}
void addMutation(rootedNode *aNode, double site){
aNode->muts[aNode->mutationNumber] = site;
aNode->mutationNumber += 1;

}

int hasMutation(rootedNode *aNode, float site){
int hasMutation(rootedNode *aNode, double site){
int i;
for (i = 0; i < aNode->mutationNumber; i++){
if (aNode->muts[i] == site) return 1;
Expand Down Expand Up @@ -1979,14 +2031,15 @@ void newickRecurse(rootedNode *aNode, float site, float tempTime){
/*makeGametesMS-- MS style sample output */
void makeGametesMS(int argc,const char *argv[]){
int i,j, size, count, k, mutNumber;
float allMuts[MAXMUTS];
double allMuts[MAXMUTS];

/* get unique list of muts */
size = 0;
for (i = 0; i < sampleSize; i++){
for (j = 0; j < allNodes[i]->mutationNumber; j++){
count = 0;
for( k = 0; k < size; k++){
// printf("%f %f %d\n",allMuts[k],allNodes[i]->muts[j],allMuts[k] == allNodes[i]->muts[j]);
if (allMuts[k] == allNodes[i]->muts[j]){
count += 1;
}
Expand All @@ -2000,11 +2053,11 @@ void makeGametesMS(int argc,const char *argv[]){
}

mutNumber = size;
qsort(allMuts, size, sizeof(allMuts[0]), compare_floats);
qsort(allMuts, size, sizeof(allMuts[0]), compare_doubles);
printf("\n//\nsegsites: %d",mutNumber);
if(mutNumber > 0) printf("\npositions: ");
for(i = 0; i < mutNumber; i++)
fprintf(stdout,"%7.5lf ",allMuts[i] );
fprintf(stdout,"%6.6lf ",allMuts[i] );
fprintf(stdout,"\n");

/* go through sites, if there is a mut is allMuts print out segSite */
Expand All @@ -2026,6 +2079,16 @@ void makeGametesMS(int argc,const char *argv[]){
}
}

void errorCheckMutations(){
int i,j;

for (i = 0; i < sampleSize; i++){
printf("allNodes[%d]:\n",i);
for (j = 0; j < allNodes[i]->mutationNumber; j++){
printf("muts[%d]=%lf\n",j,allNodes[i]->muts[j]);
}
}
}

//dropMutationsUntilTime-- places mutationson the tree iff they occur before time t
//this will not return the "correct" number of mutations conditional on theta
Expand Down
7 changes: 5 additions & 2 deletions discoalFunctions.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -21,9 +21,12 @@ int nAncestorsHere(rootedNode *aNode, float site);

int siteBetweenChunks(rootedNode *aNode, int xOverSite);
void dropMutations();
void addMutation(rootedNode *aNode, float site);
int hasMutation(rootedNode *aNode, float site);
void addMutation(rootedNode *aNode, double site);
int hasMutation(rootedNode *aNode, double site);
void makeGametesMS(int argc,const char *argv[]);
void dropMutationsRecurse();
void recurseTreePushMutation(rootedNode *aNode, float site);
void errorCheckMutations();

void mergePopns(int popnSrc, int popnDest);
void admixPopns(int popnSrc, int popnDest1, int popnDest2, double admixProp);
Expand Down
10 changes: 9 additions & 1 deletion discoal_multipop.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -253,8 +253,10 @@ int main(int argc, const char * argv[]){
}
else{
//Hudson style output
//errorCheckMutations();
makeGametesMS(argc,argv);
}
//printf("rep: %d\n",i);
i++;
}
else{
Expand Down Expand Up @@ -289,8 +291,14 @@ void getParameters(int argc,const char **argv){
if( argc < 3){
usage();
}
//locusNumber = atoi(argv[1]);

sampleSize = atoi(argv[1]);
if(sampleSize > 254){
#ifndef BIG
printf("Error: sampleSize > 254. recompile discoal and set the -DBIG flag\n");
exit(666);
#endif
}
sampleNumber = atoi(argv[2]);
nSites = atoi(argv[3]);
if(nSites>MAXSITES){
Expand Down

0 comments on commit e4a8647

Please sign in to comment.