Skip Navigation
  Print Page

Biostatistics and Bioinformatics Branch (BBB)

Skip sharing on social media links
Share this:
Skip Internal Navigation

Extended Homozygosity Score Test (EHST)

EGHH.cc

#include "EGHH.h"

int main(int argc, char **argv)
{
EGHH eghh;

eghh.load_info(argv[1]);
eghh.load_data(argv[2]);
for (int i = 0; i < eghh.snp_num; i++)
eghh.compute_allele_freq(i);

int len = strlen(argv[3]);
char *str_egh = new char[len + 8 + 1];
str_egh = strcpy(str_egh, argv[3]);
char *str_hmm = new char[len + 8 + 1];
str_hmm = strcpy(str_hmm, argv[3]);
char *str_ehh = new char[len + 8 + 1];
str_ehh = strcpy(str_ehh, argv[3]);
//char *str_FY = new char[len + 7 + 1];
//str_FY = strcpy(str_FY, argv[3]);

eghh.out_egh = strcat(str_egh, "_egh.out");
FILE * output_egh = fopen(eghh.out_egh, "wt");
eghh.output_egh = output_egh;
fprintf(output_egh, "Marker location freq T mu sigma_2 S");
fprintf(output_egh, " p-value\n");

/*eghh.out_hmm = strcat(argv[3], "_hmm.out");
FILE * output_hmm = fopen(eghh.out_hmm, "wt");
eghh.output_hmm = output_hmm;
fprintf(output_hmm, "Marker location freq T mu sigma_2 S");
fprintf(output_hmm, " p-value\n");

eghh.out_ehh = strcat(argv[3], "_ehh.out");
FILE * output_ehh = fopen(eghh.out_ehh, "wt");
eghh.output_ehh = output_ehh;
fprintf(output_ehh, "Marker location freq T mu sigma_2 S");
fprintf(output_ehh, " p-value\n");*/

eghh.out_hmm = strcat(str_hmm, "_hmm.out");
FILE * output_hmm = fopen(eghh.out_hmm, "wt");
eghh.output_hmm = output_hmm;
fprintf(output_hmm, "Marker location freq T mu sigma_2 S");
fprintf(output_hmm, " p-value\n");

eghh.out_ehh = strcat(str_ehh, "_ehh.out");
FILE * output_ehh = fopen(eghh.out_ehh, "wt");
eghh.output_ehh = output_ehh;
fprintf(output_ehh, "Marker location freq T mu sigma_2 S");
fprintf(output_ehh, " p-value\n");

/*eghh.out_FY = strcat(str_FY, "_FY.out");
FILE * output_FY = fopen(eghh.out_FY, "wt");
eghh.output_FY = output_FY;
fprintf(output_FY, "Marker location freq T mu sigma_2 S");
fprintf(output_FY, " p-value\n");*/

for (int i = 0; i < eghh.snp_num; i++)
       {
              eghh.compute_eghst(i);
       eghh.compute_hmmst(i);
       eghh.compute_ehhst(i);
       //eghh.compute_FYst(i);
       }

fclose(output_egh);
fclose(output_hmm);
fclose(output_ehh);
//fclose(output_FY);

delete [] str_egh;
delete [] str_hmm;
delete [] str_ehh;
//delete [] str_FY;

return 0;
}

void EGHH::load_info(char *infofile)
          {
          FILE *fp = fopen(infofile, "r");
          if (fp == 0)
             return;
          char line[128];

          fgets(line, 128, fp);//skip the first line, revised by Ming
          while (fgets(line, 128, fp))
             {
             char *str = strtok(line, "\t"); // the space in "\t" is a tab
             string marker(str);

             snp_list.push_back(marker);

             str = strtok(NULL, "\t");
             phy_map.push_back(atoi(str));

             str = strtok(NULL, "\t");
             allele_coded_0.push_back(str);

             str = strtok(NULL, "\t");
             allele_coded_1.push_back(str);
             }
          fclose(fp);
          }

int EGHH::load_data(char *datafile)
          {
          FILE *fp = fopen(datafile, "r");
          if (fp == 0)
             return -1;

          vector<char> line;

          while (1)
             {
             char in = fgetc(fp);
             if (isalnum(in))
                      {
             line.push_back(in);
             continue;
             }
          if (in == '\n' || in == '\r' || in == EOF)
             {
             // convert line to single char array
             if (line.size() > 0)
                      {
                      char *pl = new char[line.size()+1];
                      memset(pl,0,line.size()+1);
                      for (int i = 0; i < line.size(); i++)
                                        {
                               pl[i] = line[i];
                               }
                      data.push_back(pl);
                               }
                               line.clear();
                      }

             if (in == EOF)
                      break;
             }

          fclose(fp); // added by Andrew Redd on March 6 2009//

          if (data.size() < 1)
             return -2;

          snp_num = strlen(data[0]);
          //printf("snp_num = %d, data.size() = %d \n", snp_num, data.size());
          return 1;
          }

void EGHH::compute_allele_freq(int index)
          {
          double freq = 0.0;
          for (int i = 0; i < data.size(); i++)
             {
             if (data[i][index] == '0')
                freq++;
             }

          p.push_back( freq / data.size() );
          }

void EGHH::compute_hap_freq(int index)
          {
          hap00 = hap01 = hap10 = hap11 = 0.0;
          for (int i = 0; i < data.size(); i++)
             {
             if (data[i][index] == '0' && data[i][index + 1] == '0')
                hap00++;
                   else if (data[i][index] == '0' && data[i][index + 1] == '1')
                      hap01++;
             else if (data[i][index] == '1' && data[i][index + 1] == '0')
                      hap10++;
                   else
                hap11++;
             }
          hap00 = hap00 / data.size();
          hap01 = hap01 / data.size();
          hap10 = hap10 / data.size();
          hap11 = hap11 / data.size();
          }

void EGHH::compute_FYst(int index)
          {
          }

void EGHH::compute_ehhst(int index)
          {
          printf("index = %d \r", index);
          double T = count_homozygotes( index );

          double mu, EM;
          mu = EM = mean_M ( index );

          //calculate ER and ER^2
          mean_R = mean_Rsq = 0.0;
          mean_R_Rsq_hap ( index );
          if (mean_R * mean_R > mean_Rsq + epsilon)
                    printf("ERROR: mean_R * mean_R = %f > mean_Rsq = %f\n", mean_R * mean_R, mean_Rsq);

          // calculate EL and EL^2
          mean_L = mean_Lsq = 0.0;
          mean_L_Lsq_hap ( index );
          if (mean_L * mean_L > mean_Lsq + epsilon)
                    printf("ERROR: mean_L * mean_L = %f > mean_Lsq = %f\n", mean_L * mean_L, mean_Lsq);

          // calculate ELR
          /*double mean_LR = expectation_LR ( index );
          printf("%11s, ", snp_list[index].c_str());
          printf("mean_R = %f, mean_Rsq = %f, ", mean_R, mean_Rsq);
          printf("mean_L = %f, mean_Lsq = %f, ", mean_L, mean_Lsq);
          printf("mean_LR = %f\n\n", mean_LR);*/

          mu += mean_R + mean_L;

          double sigma = EM - EM * EM;
          sigma += (mean_Rsq - mean_R * mean_R) + (mean_Lsq - mean_L * mean_L);
          sigma += 2.0 * (mean_R + mean_L) * (1.0 - EM);
          //sigma += 2.0 * (mean_LR - mean_R * mean_L);
          sigma += 2.0 * mean_R * mean_L * ( 1.0/EM -1.0 );

          double S = T - (data.size()/2.0) * mu;
          S = S / sqrt( data.size() * sigma/2.0 );
          double pvalue = 2 * ndist(S, true);

          fprintf(output_ehh, "%11s %10d %.3f %7.0f %7.3f %7.3f %7.3f ",
                    snp_list[index].c_str(), phy_map[index], p[index], T, mu, sigma, S);
          if (pvalue > 0.05)
                    fprintf(output_ehh, "\n");
          else if (pvalue > 0.0001)
              fprintf(output_ehh, " %7.5f \n", pvalue);
          else
              fprintf(output_ehh, " %7.1e \n", pvalue);
          }

void EGHH::compute_eghst(int index)
          {
          double T = count_homozygotes( index ); //count all homozygotes around locus "index"//

          double mu, EM;
          mu = EM = mean_M (index);

          //calculate ER and ER^2
          mean_R = mean_Rsq = 0.0;
          mean_R_Rsq_gen ( index );
          if (mean_R * mean_R > mean_Rsq + epsilon)
                    printf("ERROR: mean_R * mean_R = %f > mean_Rsq = %f\n", mean_R * mean_R, mean_Rsq);

          // calculate EL and EL^2
          mean_L = mean_Lsq = 0.0;
          mean_L_Lsq_gen ( index );
          if (mean_L * mean_L > mean_Lsq + epsilon)
                    printf("ERROR: mean_L * mean_L = %f > mean_Lsq = %f\n", mean_L * mean_L, mean_Lsq);

          //calculate the mean
          mu += mean_R + mean_L;

          //calculate the variance
          double sigma = EM - EM * EM;
          sigma += (mean_Rsq - mean_R * mean_R) + (mean_Lsq - mean_L * mean_L);
          sigma += 2 * (mean_R + mean_L) * (1.0 - EM);
          sigma += 2 * mean_R * mean_L * (1.0/EM - 1.0);

          double S = T - (data.size()/2.0) * mu;
          S = S / sqrt( data.size() * sigma/2.0);
          double pvalue = 2 * ndist(S, true);

          fprintf(output_egh, "%11s %10d %.3f %7.0f %7.3f %7.3f %7.3f ",
                    snp_list[index].c_str(), phy_map[index], p[index], T, mu, sigma, S);

          if (pvalue > 0.05)
                    fprintf(output_egh, "\n");
          else if (pvalue > 0.0001)
              fprintf(output_egh, " %7.5f \n", pvalue);
          else
              fprintf(output_egh, " %7.1e \n", pvalue);
          }

void EGHH::compute_hmmst(int index)
          {
          printf("index = %d \r", index);
          double T = count_homozygotes( index );

          double mu, EM;
          mu = EM = mean_M ( index );

          //calculate ER and ER^2
          mean_R = mean_Rsq = 0.0;
          mean_R_Rsq_hmm( index );
          if (mean_R * mean_R > mean_Rsq + epsilon)
                    printf("ERROR: mean_R * mean_R = %f > mean_Rsq = %f\n", mean_R * mean_R, mean_Rsq);

          // calculate EL and EL^2
          mean_L = mean_Lsq = 0.0;
          mean_L_Lsq_hmm( index );
          if (mean_L * mean_L > mean_Lsq + epsilon)
                    printf("ERROR: mean_L * mean_L = %f > mean_Lsq = %f\n", mean_L * mean_L, mean_Lsq);

          mu += mean_R + mean_L;

          double sigma = EM - EM * EM;
          sigma += (mean_Rsq - mean_R * mean_R) + (mean_Lsq - mean_L * mean_L);
          sigma += 2.0 * (mean_R + mean_L) * (1.0 - EM);
          sigma += 2.0 * mean_R * mean_L * ( 1.0/EM -1.0 );

          double S = T - (data.size()/2.0) * mu;
          S = S / sqrt( data.size() * sigma/2.0 );
          double pvalue = 2 * ndist(S, true);

          fprintf(output_hmm, "%11s %10d %.3f %7.0f %7.3f %7.3f %7.3f ",
                    snp_list[index].c_str(), phy_map[index], p[index], T, mu, sigma, S);
          if (pvalue > 0.05)
                    fprintf(output_hmm, "\n");
          else if (pvalue > 0.0001)
              fprintf(output_hmm, " %7.5f \n", pvalue);
          else
              fprintf(output_hmm, " %7.1e \n", pvalue);
          }

//calculate the mean of M, under HWE
double EGHH::mean_M (int index)
          {
                 return p[index] * p[index] + (1 - p[index]) * (1 - p[index]);
          }

double EGHH::count_homozygotes(int index)
          {
          double T =0;
          int num_ind = data.size()/2;

          for (int i = 0; i < num_ind; i++)
                   {
                   if (data[2*i][index] == data[2*i + 1][index])
                           {
                                   T++;
                      for (int j = index + 1; j <= snp_num; j++) //count right-hand side homozygotes
                                          {//if (index == 0 && 2*i < 8) printf("T = %f, 2*i = %d, j = %d, data[2*i][j] = %c, data[2*i + 1][j] = %c \n", T, 2* i, j,data[2*i][j], data[2*i + 1][j]);
                           if (data[2*i][j] == data[2*i + 1][j])
                                                   T++;
                                              else
                                                   break;
                                              }
                                for (int j = index - 1; j >= 0; j--) //count left-hand side homozygotes
                                      {
                           if (data[2*i][j] == data[2*i + 1][j])
                                                   T++;
                                              else
                                                   break;
                                              }
                           }
                  }

          return T;
          }

// functions for genotype data
void EGHH::mean_R_Rsq_gen (int index)
          {
          double base = p[index] * p[index] + (1 - p[index]) * (1 - p[index]);
          for (int r = 1; r <= snp_num - index - 1; r++)
                  {
                  base *= p[index + r] * p[index + r] + (1 - p[index + r]) * (1 - p[index + r]);
                  double PR = 2.0 * p[index + r +1] * (1 - p[index + r +1]) * base;

                           mean_R += r * PR;
                           mean_Rsq += r * r * PR;

                           if (r * r * PR < epsilon && p[index + r + 1] > epsilon
                                 && p[index + r + 1] < 1.0 - epsilon)
                                break;
                           }
          }

void EGHH::mean_L_Lsq_gen (int index)
          {
          double base = p[index] * p[index] + (1 - p[index]) * (1 - p[index]);
          for (int ell = 1; ell < index; ell++)
                  {
                  base *= p[index - ell] * p[index - ell]
                                    + (1 - p[index - ell]) * (1 - p[index - ell]);
                  double PL = 2.0 * p[index - ell - 1] * (1 - p[index - ell - 1]) * base;

                           mean_L += ell * PL;
                           mean_Lsq += ell * ell *PL;

                           if (ell * ell * PL < epsilon && p[index - ell -1] > epsilon
                                    && p[index - ell - 1] < 1.0 - epsilon)
                                    break;
                           }
          }

//functions for hidden Markov model
void EGHH::mean_R_Rsq_hmm (int index)
          {
          double a1, a2, b1, b2, c, d;
          a1 = b2 = 1.0;
          a2 = b1 = 0.0;

          for (int r = 1; r <= snp_num - index - 1; r++)
                  {
                           compute_hap_freq (index + r - 1);

                  if (p[index + r - 1] == 1.0)
                                    {
                     c = a1 * hap00 * hap00;
                                    a2 = a1 * hap01 * hap01;
                     a1 = c;

                     d = b1 * hap00 * hap00;
                                    b2 = b1 * hap01 * hap01;
                     b1 = d;
                                    }
                           else if (p[index + r - 1] == 0.0)
                                    {
                     c = a2 * hap10 * hap10;
                                    a2 = a2 * hap11 * hap11;
                     a1 = c;

                     d = b2 * hap10 * hap10;
                                    b2 = b2 * hap11 * hap11;
                     b1 = d;
                                    }
                           else
                                    {
                     c = a1 * pow( hap00 / p[index + r - 1], 2.0)
                                                   + a2 * pow( hap10 / ( 1.0 - p[index + r - 1]), 2.0);
                                    a2 = a1 * pow( hap01 / p[index + r - 1], 2.0)
                                          + a2 * pow( hap11 / (1.0 - p[index + r - 1]), 2.0) ;
                     a1 = c;

                     d = b1 * pow( hap00 / p[index + r - 1], 2.0)
                                          + b2 * pow( hap10 / (1.0 - p[index + r - 1]), 2.0);
                                    b2 = b1 * pow( hap01 / p[index + r - 1], 2.0)
                                          + b2 * pow( hap11 / (1.0 - p[index + r - 1]), 2.0);
                     b1 = d;
                                    }

                  double PR = pow( p[ index ], 2.0 ) * (a1 + a2)
                                                   + pow( 1.0 - p[ index ], 2.0) * (b1 + b2);

                     mean_R += PR;
                     mean_Rsq += (2.0 * r - 1.0) * PR;

                     if ( (2.0 * r - 1.0) * PR < epsilon )
                                    break;
                     }
          }

void EGHH::mean_L_Lsq_hmm (int index)
          {
          double a1, a2, b1, b2, c, d;
          a1 = b2 = 1.0;
          a2 = b1 = 0.0;

          for (int ell = 1; ell <= index; ell++)
               {
                     compute_hap_freq (index - ell);

                     if (p[index - ell + 1] == 1.0)
                           {
                     c = a1 * hap00 * hap00;
                           a2 = a1 * hap10 * hap10;
                           a1 = c;

                           d = b1 * hap00 * hap00;
                           b2 = b1 * hap10 * hap10;
                     b1 = d;
                           }
                       else if (p[index - ell + 1] == 0.0)
                           {
                     c = a2 * hap01 * hap01;
                           a2 = a2 * hap11 * hap11;
                           a1 = c;

                           d = b2 * hap01 * hap01;
                           b2 = b2 * hap11 * hap11;
                     b1 = d;
                           }
                       else
                           {
                     c = a1 * pow( hap00 / p[index - ell + 1], 2.0)
                                          + a2 * pow( hap01 / (1.0 - p[index - ell + 1]), 2.0);
                           a2 = a1 * pow( hap10 / p[index - ell + 1], 2.0)
                                          + a2 * pow( hap11 / (1.0 - p[index - ell + 1]), 2.0);
                           a1 = c;

                           d = b1 * pow( hap00 / p[index - ell + 1], 2.0)
                                          + b2 * pow( hap01 / (1.0 - p[index - ell + 1]), 2.0);
                           b2 = b1 * pow( hap10 / p[index - ell + 1], 2.0)
                                          + b2 * pow( hap11 / (1.0 - p[index - ell + 1]), 2.0);
                     b1 = d;
                           }

          double PL = pow( p[ index ], 2.0 ) * (a1 + a2)
                                                   + pow( 1.0 - p[ index ], 2.0) * (b1 + b2);

                        mean_L += PL;
                        mean_Lsq += (2.0 * ell - 1.0) * PL;

                        if ( (2.0 * ell - 1.0) * PL < epsilon )
                                          break;
                        }
          }

// functions for haplotype data
void EGHH::mean_R_Rsq_hap (int index)
          {
          int test = 0;

          for (int r = 3; r <= snp_num - index + 1; r++)
               {
               map<const char*, int, strCmp> hset;
               map<const char*, int, strCmp>::iterator xx;
               map<const char*, int, strCmp>::iterator yy;

                    for (int i = 0; i < data.size(); i++)
                    {
                              char *str = new char[r + 1];
                         memset(str, 0, r+1);
                    for (int j = 0; j < r; j++)
                              str[j] = data[i][index + j];

                    //insert a key when necessary
                              xx = hset.find(str);
                    if (xx == hset.end())
                              hset[str] = 1;
                         else
                              {
                       (*xx).second += 1;
                                 delete [] str;
                                 }
                    }

               xx = hset.begin();
                    double PR = 0.0;
                    while (xx != hset.end())
                    {
                         for (yy = hset.begin(); yy != hset.end(); yy++)
                              {
                                     if ((*xx).first != (*yy).first && memcmp((*xx).first, (*yy).first, r - 1) == 0)
                                                PR += 1.0 * (*xx).second * (*yy).second / data.size() / data.size();
                                     }
                            xx++;
                    }

                    //delete the keys to save memory and to avoid leaking problems, method by Xiong Hao
                    for (map<const char*, int, strCmp>::const_iterator xx = hset.begin(); xx != hset.end(); ++xx)
                              delete [] xx->first;

              hset.clear();

              mean_R += (r - 2.0) * PR;
              mean_Rsq += (r - 2.0) * (r - 2.0) * PR;

              if ( (r - 2.0) * (r - 2.0) * PR < epsilon )
                              test +=1;
                     else
                          test = 0;

                     if (test > TEST)
                              break;
              }
          }

void EGHH::mean_L_Lsq_hap (int index)
          {
          int test = 0;

          for (int ell = 3; ell <= index + 1; ell++)
              {
              map<const char*, int, strCmp> hset;
              map<const char*, int, strCmp>::iterator xx;
              map<const char*, int, strCmp>::iterator yy;

                     for (int i = 0; i < data.size(); i++)
                     {
                     char *str = new char[ell + 1];
                     memset(str, 0, ell + 1);
                     for (int j = 0; j < ell; j++)
                                     str[j] = data[i][index - j];

                                xx = hset.find(str);
                     if (xx == hset.end())
                           hset[str] = 1;
                                else
                                     {
                                          (*xx).second += 1;
                                     delete[](str);
                                }
                     }

              xx = hset.begin();
                           double PL = 0.0;
                           while (xx != hset.end())
                           {
                                      for (yy = hset.begin(); yy != hset.end(); yy++)
                                               {
                                                        if ((*xx).first != (*yy).first && memcmp((*xx).first, (*yy).first, ell - 1) == 0)
                                                        PL += 1.0 * (*xx).second * (*yy).second / data.size() / data.size();
                                               }
                                 xx++;
                           }

                           //delete the keys to save memory and to avoidleaking problems, method by Xiong Hao
                           for (map<const char*, int, strCmp>::const_iterator xx = hset.begin(); xx != hset.end(); ++xx)
                                 delete [] xx->first;

              hset.clear();

              mean_L += (ell - 2.0) * PL;
                       mean_Lsq += (ell - 2.0) * (ell - 2.0) * PL;

              if ( (ell - 2.0) * (ell - 2.0) * PL < epsilon )
                                test +=1;
                       else
                           test = 0;

                       if (test > TEST)
                           break;
              }
          }

// The normal distribution function
//#define LOWER_TAIL_ONE 7.5
//#define UPPER_TAIL_ZERO 500

double EGHH::ndist(double z, bool upper)
          {
          // C version of ID Hill, "The Normal Integral"
          // Applied Statistics, Vol 22, pp. 424-427

          // If 7 digit accuracy is enough, alternative is
          // return erfcc(x / M_SQRT2) * 0.5;

          if (z < 0)
              {
              upper = !upper;
              z = -z;
              }

          //if (z > LOWER_TAIL_ONE && !upper || z > UPPER_TAIL_ZERO)
              // return (upper) ? 0.0 : 1.0;

          double p, y = 0.5 * z * z;

          if (z < 1.28)
              {
              p = 0.5 - z * (0.398942280444 - 0.399903438504 * y /
                                  (y + 5.75885480458 - 29.8213557808 /
                                  (y + 2.62433121679 + 48.6959930692 /
                                  (y + 5.92885724438))));
              }
          else
              {
              p = 0.398942270385 * exp (-y) /
                    (z - 2.8052e-8 + 1.00000615302 /
                    (z + 3.98064794e-4 + 1.98615381364 /
                    (z - 0.151679116635 + 5.29330324926 /
                    (z + 4.8385912808 - 15.1508972451 /
                    (z + 0.742380924027 + 30.789933034 /
                    (z + 3.99019417011))))));
              }

          return p;
          //return (upper) ? p : 1 - p;
          }

/*double EGHH::expectation_LR (int index)
          {
          int test = 0;
          int test1 = 0;
          double mean_LR = 0.0;

          for (int r = 3; r <= snp_num - index; r++)
              {
              double mean_LR_r = 0.0;
                    int ell_max = index + 1;
                    for (int ell = 3; ell <= index + 1 && ell <= ell_max; ell++)
                    {
                    map<const char*, int, strCmp> hset;
                    map<const char*, int, strCmp>::iterator xx;
                    map<const char*, int, strCmp>::iterator yy;

                          for (int i = 0; i < data.size(); i++)
                          {
                          char *str = new char[ell + r];
                          memset(str, 0, ell + r);
                          for (int j = 0; j < ell + r - 1; j++)
                                            str[j] = data[i][index - ell + j +1];

                                    xx = hset.find(str);
                          if (xx == hset.end())
                              hset[str] = 1;
                                    else
                                             {
                              (*xx).second += 1;
                                             delete[](str);
                              }
                          }

                    xx = hset.begin();
                          double PLR = 0.0;
                          while (xx != hset.end())
                          {if (r < 10 and ell < 10) fprintf(output,"r = %d, ell = %d, (*xx).first = %s, (*xx).second = %d, PLR = %f \n", r, ell, (*xx).first, (*xx).second,PLR);
                                            for (yy = hset.begin(); yy != hset.end(); yy++)
                                                    {if (r < 10 and ell < 10) fprintf(output,"r = %d, ell = %d, (*yy).first = %s, (*yy).second = %d, PLR = %f \n", r, ell,(*yy).first, (*yy).second, PLR);
                                                    if (Memcmp((*xx).first, (*yy).first, 1, ell + r - 3) == 0)
                                                            {
                                    if ((*xx).first[0] == '0' && (*xx).first[ell + r - 2] == '0' && (*yy).first[0] == '1' && (*yy).first[ell + r - 2] == '1')
                                                                    PLR += 1.0 * (*xx).second * (*yy).second / data.size() / data.size();
                                                                    else if ((*xx).first[0] == '1' && (*xx).first[ell + r - 2] == '1' && (*yy).first[0] == '0' && (*yy).first[ell + r - 2] == '0')
                                                                             PLR += 1.0 * (*xx).second * (*yy).second / data.size() / data.size();
                                                                    else if ((*xx).first[0] == '0' && (*xx).first[ell + r - 2] == '1' && (*yy).first[0] == '1' && (*yy).first[ell + r - 2] == '0')
                                                                             PLR += 1.0 * (*xx).second * (*yy).second / data.size() / data.size();
                                                                    else if ((*xx).first[0] == '1' && (*xx).first[ell + r - 2] == '0' && (*yy).first[0] == '0' && (*yy).first[ell + r - 2] == '1')
                                                                             PLR += 1.0 * (*xx).second * (*yy).second / data.size() / data.size();
                                                    }
                                                         }
                                    xx++;
                          }

                    //delete the keys to save memory and to avoidleaking problems
                          for (map<const char*, int, strCmp>::const_iterator xx = hset.begin(); xx != hset.end(); ++xx)
                                             delete [] xx->first;
                    hset.clear();

                                mean_LR_r += (ell - 2.0) * PLR;
fprintf(output, "PLR = %f, mean_LR_r = %f \n", PLR, mean_LR_r);

                    if ( (ell - 2.0) * PLR < epsilon )
                                     test1 +=1;
                          else
                                test1 = 0;

                          if (test1 > TEST) 
                                     {
                          ell_max = ell;
                                               break;
                                               }
                          }

                    mean_LR += (r - 2.0) * mean_LR_r;
              if ( (r - 2.0) * mean_LR_r < epsilon )
                                     test +=1;
                          else
                                 test = 0;

                   if (test > TEST)
                                 break;
              }

          return mean_LR;
          }
*/

Last Updated Date: 04/12/2013
Last Reviewed Date: 04/12/2013

Contact Information

Name: Dr Paul Albert
Chief and Senior Investigator
Biostatistics and Bioinformatics Branch
Phone: 301-496-5582
E-mail: albertp@mail.nih.gov

Staff Directory
Vision National Institutes of Health Home BOND National Institues of Health Home Home Storz Lab: Section on Environmental Gene Regulation Home Machner Lab: Unit on Microbial Pathogenesis Home Division of Intramural Population Health Research Home Bonifacino Lab: Section on Intracellular Protein Trafficking Home Lilly Lab: Section on Gamete Development Home Lippincott-Schwartz Lab: Section on Organelle Biology