/* This file contains new versions of the routines segup, segdown, */ /* and some auxiliary routines for use with */ /* LODSCORE, ILINK, LINKMAP, and MLINK */ /* The modifications are described in the papers: */ /* R. W. Cotiingham, Jr., R. M. Idury, and A. A. Schaffer, */ /* Faster Sequential Genetic Linkage Computations */ /* American Journal of Human Genetics, 53(1993), pp. 252-263*/ /* and A. A. Schaffer, S. K. Gupta, K. Shriram, and R. W. Cottingham, Jr. */ /* Avoiding Recomputation in Linkage Analysis, */ /* Human Heredity 44(1994), pp. 225-237. */ /* The versions in this file do not incorporate some of the changes */ /* described in the paper that require lots of memory */ /* Faster versions that use more memory are in automodified.c */ /*Most of the code in this file was written by R. M. Idury*/ #include "commondefs.h" /*segsum2 is used in segup to compute some common subexpressions in the probability updates. first and second are the joint genotypes of the two parents. fslength is the number of probabilities needed for the combined isozygote class of both parental genotypes; fs stands for the product of first and second. LINK is used to pass the genetic information for parents and children */ Local void segsum2(first,second,fslength) unsigned first,second, fslength; { int g1, g2, g3, g4; /*indices to store gene numbers*/ int f1, f2, s1, s2; /*indices to store haplotype numbers*/ int index; /*counter for combined isozygote equivalence class, will range from 0 to fslength-1*/ int i, j, k, l; /*loop indices*/ int FORLIM; /*limit on a for loop*/ gennurec *WITH; /*temporary pointer for gene records*/ int FORLIM1, FORLIM2; /*lower and upper limits on a for loop*/ int *TEMPGENE1, *TEMPGENE2; /*temporary pointers into genenumber table*/ double *tempwith3; /*temporarily stores genarray for current child*/ WITH = gennustruct; /*get genenumber records*/ FORLIM = fence[first]; /*end of the isozygote class for first*/ FORLIM1 = base[second]; /*start of the isozygote class for second*/ FORLIM2 = fence[second]; /*end of the isozygote class for second*/ index = 0; /*start with first isozygote in the class*/ /*iterate over isozygotes of i*/ for (i = base[first]; i < FORLIM; i++) { f1 = haps1[i]; /*retrieve first haplotype of joint genotype i*/ f2 = haps2[i]; /*retrieve second haplotype*/ TEMPGENE1 = WITH->genenumber[f1 - 1]; /*lookup partial information in genenumber table*/ TEMPGENE2 = WITH->genenumber[f2 - 1]; /*iterate over isozygotes of j*/ for (j = FORLIM1; j < FORLIM2; j++) { s1 = haps1[j]; /*retrieve first haplotype of joint genotype j*/ s2 = haps2[j]; /*retrieve second haplotype*/ /*lookup the four ways to combine one haplotype from i and one from j*/ g1 = TEMPGENE1[s1 - 1]; g2 = TEMPGENE1[s2 - 1]; g3 = TEMPGENE2[s1 - 1]; g4 = TEMPGENE2[s2 - 1]; /*iterate over children; update partial computations of new probabilities and store in array tempseg note that tempseg has nchild*fslength relevant entries The fslength entries for child 1 come first, then the fslength entries, for child 2, etc. This is why the increment on l is fslength for each change of child*/ for (l = index, k = 0 ; k < nchild; l += fslength, k++) { tempwith3 = thischild[k]->genarray; /*retrieve genarray*/ /* sum the probabilities for the four joint genotypes that the child can get from the current parental joint genotypes*/ tempseg[l] = (tempwith3[g1 - 1] + tempwith3[g2 - 1] + tempwith3[g3 - 1] + tempwith3[g4 - 1]); } index++; /*increment isozygote class counter*/ } } } /*segsumdown2 is used in segdown to compute some common subexpressions in the probability updates. first and second are the joint genotypes of the two parents. fslength is the number of probabilities needed for the combined isozygote class for both parental genotypes; fs stands for the product of first and second. LINK is used to pass the genetic information for parents and children */ Local void segsumdown2(first,second,fslength) unsigned first,second, fslength; { int g1, g2, g3, g4; /*indices to store gene numbers*/ int f1, f2, s1, s2; /*indices to store haplotype numbers*/ int index; /*counter for combined isozygote equivalence class, will range from 0 to fslength-1*/ int index2; /*used as index for isozygote class*/ int i, j, k, l; /*loop indices*/ gennurec *WITH; /*temporary pointer for gene records*/ int FORLIM; /*limit on a for loop*/ int FORLIM1, FORLIM2; /*lower and upper limits on a for loop*/ int *TEMPGENE1, *TEMPGENE2; /*temporary pointers into genenumber table*/ double *tempwith3; /*temporarily stores genarray for current child*/ WITH = gennustruct; /*get genenumber records*/ FORLIM = fence[first]; /*end of the isozygote class for first*/ FORLIM1 = base[second]; /*start of the isozygote class for second*/ FORLIM2 = fence[second]; /*end of the isozygote class for second*/ index = 0; /*start with first isozygote in the class*/ index2 = 0; /*iterate over isozygotes of i*/ for (i = base[first]; i < FORLIM; i++) { f1 = haps1[i]; /*retrieve first haplotype of joint genotype i*/ f2 = haps2[i]; /*retrieve second haplotype*/ TEMPGENE1 = WITH->genenumber[f1 - 1]; /*lookup partial information in genenumber table*/ TEMPGENE2 = WITH->genenumber[f2 - 1]; /*iterate over isozygotes of j*/ for (j = FORLIM1; j < FORLIM2; j++) { s1 = haps1[j]; /*retrieve first haplotype of joint genotype j*/ s2 = haps2[j]; /*retrieve second haplotype*/ /*lookup the four ways to combine one haplotype from i and one from j*/ g1 = TEMPGENE1[s1 - 1]; g2 = TEMPGENE1[s2 - 1]; g3 = TEMPGENE2[s1 - 1]; g4 = TEMPGENE2[s2 - 1]; /*store these gene numbers for later use; the key point is that we will want these numbers consecutively later, so we store them consecutively in segindex*/ segindex[index2++] = g1; segindex[index2++] = g2; segindex[index2++] = g3; segindex[index2++] = g4; /*iterate over children; update partial computations of new probabilities and store in array tempseg note that tempseg has nchild*fslength relevant entries The fslength entries for child 1 come first, then the fslength entries, for child 2, etc. This is why the increment on l is fslength for each change of child*/ for (l = index, k = 0 ; k < nchild; l += fslength, k++) { tempwith3 = thischild[k]->genarray; /*retrieve genarray*/ /* sum the probabilities for the four joint genotypes that the child can get from the current parental joint genotypes*/ tempseg[l] = (tempwith3[g1 - 1] + tempwith3[g2 - 1] + tempwith3[g3 - 1] + tempwith3[g4 - 1]); } index++; /*increment isozygote class counter*/ } } } /*lsegfun2 does a logical test similar to the computation done in the original segfun to determine whether segfun would return 0.0. If segfun would return 0.0, then lsegfun2 return 0 (FALSE), while if segfun would not return 0.0., lsegfun2 returns 1 (TRUE). Given, a combined isozygote class, we want to know whether any elements of that isozygote class are possible joint genotypes for the parents. This will be the case if and only if each child has a nonzero probability for at least one of the joint genotypes that the parents can produce given their isozygote class. first and second are the joint genotypes for the parents.*/ unsigned lsegfun2(first,second) unsigned first, second; { int g1, g2, g3, g4; /*four gene numbers*/ int i, j, k; /* loop indices*/ int f1, f2, s1, s2; /*haplotype numbers*/ int FORLIM; /*loop bound*/ gennurec *WITH; /*pointer to gene number records*/ int FORLIM1, FORLIM2; /*loop bounds*/ int *TEMPGENE1, *TEMPGENE2; /* store pointers into genenumber*/ unsigned char *tempflag3; /*stores sparsity pattern for child's genarray*/ WITH = gennustruct; /*fetch gene records*/ FORLIM = fence[first]; /*find end of isozygote class for first*/ FORLIM1 = base[second];/*find beginning and end of isozygote class for second*/ FORLIM2 = fence[second]; /*try to find a non-zero value for each child*/ for (k = 0; k < nchild; k++) { /*code for the non-boolean version is shown in comments tempwith3 = thischild[k]->genarray;*/ /*retrieve sparsity pattern for child k*/ tempflag3 = thischild[k]->sparseflag; /*iterate over all recombined isozygotes of first*/ for (i = base[first]; i < FORLIM; i++) { f1 = haps1[i]; /*retrieve the haplotypes of this genotype*/ f2 = haps2[i]; TEMPGENE1 = WITH->genenumber[f1 - 1];/*get pointer into genenumber for this haplotype*/ TEMPGENE2 = WITH->genenumber[f2 - 1]; /*iterate over all the recombined isozygotes of second*/ for (j = FORLIM1; j < FORLIM2; j++) { s1 = haps1[j]; /*get haplotypes of this genotype*/ s2 = haps2[j]; /*retrieve the four genes that this combination of joint haplotypes can produce*/ g1 = TEMPGENE1[s1 - 1]; g2 = TEMPGENE1[s2 - 1]; g3 = TEMPGENE2[s1 - 1]; g4 = TEMPGENE2[s2 - 1]; /*if(tempwith3[g1 - 1] != 0.0 || tempwith3[g2 - 1] != 0.0 || tempwith3[g3 - 1] != 0.0 || tempwith3[g4 - 1] != 0.0) goto notzero;*/ /*if any of the flags is TRUE, then this child can have this genotype, and we move immediately to testing the next child*/ if(tempflag3[g1 - 1] != 0 || tempflag3[g2 - 1] != 0 || tempflag3[g3 - 1] != 0 || tempflag3[g4 - 1] != 0) goto notzero; } } return 0; /*this child failed to have any of the genotypes for the isozygote class, so this isozygote class is not possible*/ notzero: continue; } return 1; /*all children passed the test and have a possible joint genotype*/ } /*segup updates the genotype probabilities for a parent (p) based on the probabilities of children. The parameter LINK includes all the genetic information about p, p's spouse q, and their children*/ Void segup(LINK) struct LOC_seg *LINK; { int findex, sindex; int FORLIM; /*loop bound*/ int nonzindex, nonzcount; /*loop index and counter for nonzero values*/ unsigned int nonzgens[maxfem]; /*list of genotypes having nonzero prob. */ int step1, step2; /*size of isozygote classes*/ double val, temp1; /*temporaries to store intermediate values*/ unsigned int i, j, first, second; /*genotype indices*/ unsigned int fslength; /*size of product isozygote class for p and q*/ unsigned k; /*index to loop over children*/ unsigned char skip; /*used to skip iterations of a loop*/ thisarray *WITH2; /*stores genetic information about p*/ thisarray *WITH3; /*stores gnetic information about q*/ double *newwith2, *newwith3, *newwithr, *newsegr; /*store genarrays for p,q, and children*/ thetavalues *WITH4; /*stores theta values for p*/ thetavalues *WITH5; /*store theta values for q*/ boolean depend; /*used to handle loops*/ unsigned char *newflag2, *newflag3; /*store sparsity patterns for p and q genarrays*/ /* newsegprob, newsegprob1, and newsegprob2 are used to hold segprob arrays, which contain the probabilities for various patterns of recombination events*/ double *newsegprob1, *newsegprob2; double *tempprob; /*temporary holder for probability array*/ initseg(LINK); /*get data about this p,q,children triple*/ /*get sparsity patterns for p and q genarrays*/ newflag2 = (*LINK->p)->gen->sparseflag; newflag3 = (*LINK->q)->gen->sparseflag; WITH2 = (*LINK->p)->gen; /*get genetic data for p*/ WITH3 = (*LINK->q)->gen; /*get genetic data for q*/ newwith2=WITH2->genarray; /*get genarray for p*/ newwith3=WITH3->genarray; /*get genarray for q*/ WITH4 = LINK->firstsex; /*get recomb. probabilities for p*/ WITH5 = LINK->secondsex; /*get recomb. probabilities for q*/ /*The case of 1 child is handled specially because the subcomputations are much simpler. In particular we do not need to multiply the probabilities among the different children. In a typical pedigree many pairs will have only one child about which information is known (in essence this is the child that keeps the pedigree connected), so this is an important and common special case. */ { /*nchild is bigger than 1*/ /*find nonzero entries in q's genarray and make a list of them stored in nonzgens; just get one per isozygote class*/ nonzcount = 0; /*iterate over genotypes for q*/ for(i = 0; i < fgeno; i += step2) { /*number of distinct probabilities needed for i's isoz. class*/ step2 = probend[i] - probstart[i] + 1; for(j = i; j < i+step2; j++) if(newflag3[j] != 0) { nonzgens[nonzcount++] = i; /*store index of nonzero value*/ break; /*go to next isozygote class*/ } } newsegprob1 = WITH4->segprob; /*get p recomb. probabilities*/ newsegprob2 = WITH5->segprob; /*get q recomb. probabilties*/ /*iterate over genotypes for p*/ for (first = 0; first < fgeno; first += step1) { /*number of distinct probabilties needed for first's isoz. class*/ step1 = probend[first] - probstart[first] + 1; skip = 1; /*work only on those isozygotes that are possible*/ for(i = first; i < first+step1; i++) if(newflag2[i] != 0) { skip = 0; break; /*go to next isozygote in class*/ } if(skip) continue; /*initialize update multiple for each isozygote in class*/ for(i = 0; i < step1; i++) segval[i] = 0.0; /*iterate over the genotypes representing isozygote classes that q may have*/ for (nonzindex = 0; nonzindex < nonzcount; nonzindex++) { second = nonzgens[nonzindex]; /*check if this first, second pair yield a nonzero value among children's probabilities*/ if(lsegfun2(first,second) == 0) continue; /*number of distinct probabilties needed for second's isoz. class*/ step2 = probend[second] - probstart[second] + 1; fslength = step1 * step2; /*number of probs. for combined isozygote class*/ /*call segsum2 to compute the part of the conditional probability update that is common to all the members of the combined isozygote class defined by first and second */ segsum2(first,second,fslength); newsegr = tempseg2; tempprob = tempseg; for(k = 0; k < nchild; k++) { for(i = 0; i < step1; i++) { newwithr = tempprob; findex = probstart[first+i]-1; temp1 = newsegprob1[findex]; for(sindex = 0; sindex < step2; sindex++) newsegr[sindex] = temp1*(*newwithr++); FORLIM = probend[first+i]; for(findex++; findex < FORLIM; findex++) { temp1 = newsegprob1[findex]; for(sindex = 0; sindex < step2; sindex++) newsegr[sindex] += temp1*(*newwithr++); } newsegr += step2; } tempprob += fslength; } /*now specialize update for each member of first's class*/ for(i = 0; i < step1; i++) { if(newflag2[first+i] == 0) { continue; /*skip if this isozygote is not possible*/ } /*further specialize update for each member of second's isozygote class*/ for(j = 0; j < step2; j++) { /*skip if this isozygote not possible*/ if(newflag3[second+j] == 0) continue; /*get offset into probtable; the offset depends on the isozygote class size and index of each parent note that fisozygoteindex gets incremented by the size of the joint class each time, so it represents a sum of all the numbers of distinct probabilities needed for all joint iso. classes considered before the current p isozygote class*/ /*combine for all children*/ /*due to the arrangement in segsum all the probability contributions for a given child are contiguous in the tempseg array. the number of contributions is fslength which is the number of probabilities needed for the joint isozygote class of the parents. We get the contribution of the first child (index 0) and then loop over the rest*/ val = 1.0; newwithr = tempseg2 + i*step2; tempprob = newsegprob2 + probstart[second+j] - 1; for(k = 0; k < nchild; k++) { temp1 = 0.0; for(sindex = 0; sindex < step2; sindex++) temp1 += tempprob[sindex] * newwithr[sindex]; val *= temp1; newwithr += fslength; } /*update segval entry for this isozygote of first*/ segval[i] += newwith3[second+j] * val; } } } /*update p's genarray for each isozygote in this class*/ for(i = 0; i < step1; i++) newwith2[first+i] *= segval[i] * segscale; } } /* If any of the nonzero entries in p's genarray became 0, we want to set them to zero to avoid computations on subsequent calls*/ for(i = 0; i < fgeno; i++) if((newflag2[i] != 0) && (newwith2[i] == 0.0)) newflag2[i] = 0; /*Added by Alex to handle loops */ if (loopfirstgen && (!(*LINK->p)->loopdepend)) { depend = false; for(i=0; iloopdepend) { depend = true; break; } depend = (depend || ((*LINK->q)->loopdepend)); if (depend) { (*LINK->p)->loopdepend = depend; (*LINK->p)->loopneeded = false; } } if ((*LINK->p)->loopdepend) { if (!((*LINK->q)->loopdepend)) (*LINK->q)->loopneeded = true; for(i=0; i< nchild; i++) if (!(childarray[i]->loopdepend)) childarray[i]->loopneeded = true; } cleanup(LINK->q, LINK->LINK); exitseg(LINK); } /*segup*/ /* segdown updates the genotype probabilities for a child based on the probabilities of parents and other children. The parameter LINK includes all the genetic information about p, p's spouse q, and their children */ Void segdown(LINK) struct LOC_seg *LINK; { int findex, sindex; int FORLIM; /*loop bound*/ int nonzindex, nonzcount; /*loop index and counter for nonzero values*/ unsigned int nonzgens[maxfem]; /*list of genotypes having nonzero prob. */ int step1, step2; /*size of isozygote classes*/ genotype gene; /*stores newly forming genarray for parent p*/ double valtemp; /*intermediate value in probability updates*/ double val, temp1; /*temporaries to store intermediate values*/ unsigned int f1,f2; unsigned int here, i, j, first, second; /*genotype indices*/ unsigned int fslength; /*number of probs. for product isoz. class of p,q*/ unsigned jointisoindex; /*index to work within isozygote classes*/ unsigned currentindex; /*index to update genarray within isozygote class*/ unsigned k; /*index to loop over children*/ unsigned char skip; /*used to skip iterations of a loop*/ thisarray *WITH2; /*stores genetic information about p*/ thisarray *WITH3; /*stores gnetic information about q*/ double *newwith2, *newwith3, *newwithr, *newsegr, *newwithc; /*store genarrays for p,q, and children*/ thetavalues *WITH4; /*stores theta values for p*/ thetavalues *WITH5; /*store theta values for q*/ char flag[maxfem]; /*flag is used to weed out uninteresting haplotypes more precisely, flag[i] will be true if and only if some child has a non-zero probability of having a joint genotype, one of whose haplotypes is i*/ unsigned int c1, c2; /*haplotypes*/ unsigned char *newflag2, *newflag3, *newflagr; /*store sparsity patterns for p and q genarrays*/ /* newsegprob, newsegprob1, and newsegprob2 are used to hold segprob arrays, which contain the probabilities for various patterns of recombination events*/ double *newsegprob, *newsegprob1, *newsegprob2; double *tempprob; /*temporary holder for probability array*/ boolean depend; /*used to handle loops*/ /* The arrays psumcache and qsumcache store conditional probabilities of different haplotypes being passed on from p and q respectively */ double psumcache[maxhap], qsumcache[maxhap]; initseg(LINK); /*get data about this p,q,children triple*/ /*get sparsity patterns for p, q, and child genarrays*/ newflag2 = (*LINK->p)->gen->sparseflag; newflag3 = (*LINK->q)->gen->sparseflag; newflagr = (*LINK->r)->gen->sparseflag; WITH2 = (*LINK->p)->gen; /*get genetic data for p*/ WITH3 = (*LINK->q)->gen; /*get genetic data for q*/ newwith2 = WITH2->genarray; /*get genarray for p*/ newwith3 = WITH3->genarray; /*get genarray for q*/ newwithr = (*LINK->r)->gen->genarray; /*get genarray for first child*/ WITH4 = LINK->firstsex; /*get p recomb. probabilities*/ WITH5 = LINK->secondsex; /*get q recomb. probabilities*/ /*The case of 1 child (nchild==0) is handled specially because the subcomputations are much simpler. In particular we do not need to multiply the probabilities among the different children. In a typical pedigree many pairs will have only one child about which information is known (in essence this is the child that keeps the pedigree connected), so this is an important and common special case. */ if(nchild == 0) { /*initialize cache data structures*/ for(i = 0;i < maxhap;i++) { flag[i] = 0; psumcache[i] = 0.0; qsumcache[i] = 0.0; } /*initialize gene array and set up flag array*/ for(i = 0; i < fgeno; i++) { gene[i] = 0.0; flag[invgenenum1[i]-1] = 1; flag[invgenenum2[i]-1] = 1; } /*This section of the code precomputes for each haplotype the the probability that the child will inherit this haplotype from p. Each genotype has two haplotypes, but can produce others by recombination. Therefore, for each genotype we must sum over the different haplotypes that can be produced by its isozygote class. The contributions for each haplotype are stored in psumcache. Afterwards a similar computation is done for the inheritance from q with the results stored in qsumcache.*/ newsegprob = WITH4->segprob; /*get probabilities for recomb. patterns for p*/ for (first = 0; first < fgeno; first++) { if(newflag2[first] == 0) continue; /*use only possible genotypes*/ FORLIM = fence[first]; /*find end of isozygote class of first*/ valtemp = newwith2[first]; /*probability of getting this genotype*/ /*iterate over all members of first's isozygote calss*/ for (i = base[first]; i < FORLIM; i++) { f1 = haps1[i]; /*get haplotypes*/ f2 = haps2[i]; if ((flag[f1-1] !=0) || (flag[f2-1] != 0)) { /*get probability offset for i*/ /*multiply probability of getting genotype times probability of this recombination pattern and haplo. choice*/ val = valtemp * newsegprob[hind[i]]; /*add to psumcache*/ if(flag[f1-1] != 0) { psumcache[f1-1] += val; } if(flag[f2-1] != 0) { psumcache[f2-1] += val; } } } } newsegprob = WITH5->segprob; /*get probabilities for q*/ /*iterate over all joint genotypes*/ for (first = 0; first < fgeno; first++ ) { if(newflag3[first] == 0) continue; /*check if joint genotype is possible*/ valtemp = newwith3[first]; /*get cond. prob. that q has genotype first*/ FORLIM = fence[first]; /*find bounds for the isozygote class of first*/ /*iterate over all the recombined genotypes of the isozygote class*/ for (i = base[first]; i < FORLIM; i++) { f1 = haps1[i]; f2 = haps2[i]; if((flag[f1-1] !=0) || (flag[f2-1] !=0)) { /*condition probability of first as a genotype multiplied by probability of passing on f1 (alternatively f2) as a haplotype*/ /*get probability offset for i*/ val = valtemp*newsegprob[hind[i]]; /*store in qsumcache*/ if(flag[f1-1] != 0) { qsumcache[f1 - 1] += val; } if(flag[f2-1] != 0) { qsumcache[f2 - 1] += val; } } } } /*In this section of the code we update the probabilities for the child based on the probabilities for the parents*/ /*Iterate over all joint genotypes of the child*/ for(here = 0; here < fgeno; here++) { if(newflagr[here] == 0) continue; c1 = invgenenum1[here]; c2 = invgenenum2[here]; /*probability of child getting genotype here as a result of p passing on c1 and q passing on c2 is summed to gene[here] */ gene[here] += psumcache[c1-1] * qsumcache[c2-1]; /*if c1 is distinct from c2 we need to do the same computation reversing the roles of the two haplotypes*/ if(c1 != c2) { gene[here] += psumcache[c2-1] * qsumcache[c1-1]; } } /*set up new genarray for r; it is gene scaled by segscale*/ for (first = 0; first < fgeno; first++) { if(newflagr[first] == 0) continue; if(gene[first] == 0.0) newflagr[first] = 0; newwithr[first] *= segscale*gene[first]; } } else { /*nchild is bigger than 0*/ /*initialize genarray entries for child to 0*/ for(i = 0; i < fgeno; i++) { gene[i] = 0.0; } /*find nonzero entries in q's genarray and make a list of them stored in nonzgens; just get one per isozygote class*/ nonzcount = 0; newwith3=(*LINK->q)->gen->genarray; /*iterate over genotypes for q*/ for(i = 0; i < fgeno; i += step2) { /*number of distinct probabilties needed for i's isoz. class*/ step2 = probend[i] - probstart[i] + 1; for(j = i; j < i+step2; j++) if(newflag3[j] != 0) { nonzgens[nonzcount++] = i; /*store index of nonzero value*/ break; /*go to next isozygote class*/ } } newsegprob1 = WITH4->segprob; /*get p probabilities*/ newsegprob2 = WITH5->segprob; /*get q probabilties*/ /*iterate over genotypes for p*/ for (first = 0; first < fgeno; first += step1) { /*number of distinct probabilties needed for first's isoz. class*/ step1 = probend[first] - probstart[first] + 1; skip = 1; /*work only on those isozygotes that are possible*/ for(i = first; i < first+step1; i++) if(newflag2[i] != 0) { skip = 0; break; /*go to next isozygote in class*/ } if(skip) continue; /*iterate over the genotypes representing isozygote classes that q may have*/ for (nonzindex = 0; nonzindex < nonzcount; nonzindex++) { second = nonzgens[nonzindex]; /*check if this first, second pair yield a nonzero value among children's probabilities*/ if(lsegfun2(first,second) == 0) continue; /*number of distinct probabilties needed for second's isoz. class*/ step2 = probend[second] - probstart[second] + 1; fslength = step1 * step2; /*number of probs. of combined isozygote class*/ /*call segsumdown2 to compute the part of the conditional probability update that is common to all the members of the combined isozygote class defined by first and second */ segsumdown2(first,second,fslength); for(jointisoindex = 0; jointisoindex < fslength; jointisoindex++) segval[jointisoindex] = 0.0; newsegr = tempseg2; tempprob = tempseg; for(k = 0; k < nchild; k++) { for(i = 0; i < step1; i++) { newwithr = tempprob; findex = probstart[first+i]-1; temp1 = newsegprob1[findex]; for(sindex = 0; sindex < step2; sindex++) newsegr[sindex] = temp1*(*newwithr++); for(findex++; findex < probend[first+i]; findex++) { temp1 = newsegprob1[findex]; for(sindex = 0; sindex < step2; sindex++) newsegr[sindex] += temp1*(*newwithr++); } newsegr += step2; } tempprob += fslength; } /*now specialize update for each member of first's class*/ for(i = 0; i < step1; i++) { if(newflag2[first+i] == 0) { continue; /*skip if this isozygote is not possible*/ } /*further specialize update for each member of second's isozygote class*/ for(j = 0; j < step2; j++) { /*skip if this isozygote not possible*/ if(newflag3[second+j] == 0) { continue; } /*get offset into probtable; the offset depends on the isozygote class size and index of each parent note that fisozygoteindex gets incremented by the size of the joint class each time, so it represents a sum of all the numbers of distinct probabilities needed for all joint iso. classes considered before the current p isozygote class*/ /*combine for all children*/ val = 1.0; newwithc = tempseg2 + i * step2; tempprob = newsegprob2 + probstart[second+j] - 1; for(k = 0; k < nchild; k++) { temp1 = 0.0; for(sindex = 0; sindex < step2; sindex++) temp1 += tempprob[sindex] * newwithc[sindex]; val *= temp1; newwithc += fslength; } /*probability of this combination of parent genotypes*/ val *= newwith2[first+i] * newwith3[second+j]; jointisoindex = 0; for(findex = probstart[first+i]-1; findex < probend[first+i]; findex++) for(sindex = probstart[second+j]-1; sindex < probend[second+j]; sindex++) /*probability of this recombination pattern (based on other children)*/ segval[jointisoindex++] += newsegprob1[findex] * newsegprob2[sindex] * val; } } /*update the probabilities of four joint genotypes the child might get; each different choice of recombination pattern will lead to a different set of four genotypes*/ currentindex = 0; for(jointisoindex = 0; jointisoindex < fslength; jointisoindex++) { temp1 = segval[jointisoindex]; gene[segindex[currentindex++]-1] += temp1; gene[segindex[currentindex++]-1] += temp1; gene[segindex[currentindex++]-1] += temp1; gene[segindex[currentindex++]-1] += temp1; } } } /*finally update child's real genarray by coppy gene multiplied by scale factor segscale*/ newwithr = (*LINK->r)->gen->genarray; /*get genarray for first child*/ for(i = 0; i < fgeno; i++) { if(gene[i] == 0.0) newflagr[i] = 0; /*if probability changes from nonzero to 0.0 change flag to 0*/ newwithr[i] *= segscale * gene[i]; } } /*Added by A. A. Schaffer to handle loops */ if (loopfirstgen && (!((*LINK->r)->loopdepend))) { depend = false; for(i=0; iloopdepend) { depend = true; break; } depend = (depend || ((*LINK->p)->loopdepend) || ((*LINK->q)->loopdepend)); if (depend) { (*LINK->r)->loopdepend = depend; (*LINK->r)->loopneeded = false; } } if ((*LINK->r)->loopdepend) { if (!((*LINK->p)->loopdepend)) (*LINK->p)->loopneeded = true; if (!((*LINK->q)->loopdepend)) (*LINK->q)->loopneeded = true; for(i=0; i< nchild; i++) if (!(childarray[i]->loopdepend)) childarray[i]->loopneeded = true; } cleanup(LINK->p, LINK->LINK); cleanup(LINK->q, LINK->LINK); exitseg(LINK); } /*segdown*/