src/dist_dna.c

   1 /* dist_dna.c       2010-03-29 */
   2
   3 /* Copyright 2005-2010 Emmanuel Paradis
   4
   5 /* This file is part of the R-package `ape'. */
   6 /* See the file ../COPYING for licensing issues. */
   7
   8 #include <R.h>
   9 #include <R_ext/Lapack.h>
  10
  11 /* from R: print(log(4), d = 22) */
  12 #define LN4 1.386294361119890572454
  13
  14 /* returns 8 if the base is known surely, 0 otherwise */
  15 #define KnownBase(a) (a & 8)
  16
  17 /* returns 1 if the base is adenine surely, 0 otherwise */
  18 #define IsAdenine(a) a == 136
  19
  20 /* returns 1 if the base is guanine surely, 0 otherwise */
  21 #define IsGuanine(a) a == 72
  22
  23 /* returns 1 if the base is cytosine surely, 0 otherwise */
  24 #define IsCytosine(a) a == 40
  25
  26 /* returns 1 if the base is thymine surely, 0 otherwise */
  27 #define IsThymine(a) a == 24
  28
  29 /* returns 1 if the base is a purine surely, 0 otherwise */
  30 #define IsPurine(a) a > 63
  31
  32 /* returns 1 if the base is a pyrimidine surely, 0 otherwise */
  33 #define IsPyrimidine(a) a < 64
  34
  35 /* returns 1 if both bases are different surely, 0 otherwise */
  36 #define DifferentBase(a, b) (a & b) < 16
  37
  38 /* returns 1 if both bases are the same surely, 0 otherwise */
  39 #define SameBase(a, b) KnownBase(a) && a == b
  40
  41 /* computes directly the determinant of a 4x4 matrix */
  42 double detFourByFour(double *x)
  43 {
  44     double det, a33a44, a34a43, a34a42, a32a44, a32a43, a33a42, a34a41, a31a44, a31a43, a33a41, a31a42, a32a41;
  45
  46     a33a44 = x[10]*x[15]; a34a43 = x[14]*x[11];
  47     a34a42 = x[14]*x[7];  a32a44 = x[6]*x[15];
  48     a32a43 = x[6]*x[11];  a33a42 = x[10]*x[7];
  49     a34a41 = x[14]*x[3];  a31a44 = x[2]*x[15];
  50     a31a43 = x[2]*x[11];  a33a41 = x[10]*x[3];
  51     a31a42 = x[2]*x[7];   a32a41 = x[6]*x[3];
  52
  53     det = x[0]*x[5]*(a33a44 - a34a43) + x[0]*x[9]*(a34a42 - a32a44) +
  54       x[0]*x[13]*(a32a43 - a33a42) + x[4]*x[9]*(a31a44 - a34a41) +
  55       x[4]*x[13]*(a33a41 - a31a43) + x[4]*x[1]*(a34a43 - a33a44) +
  56       x[8]*x[13]*(a31a42 - a32a41) + x[8]*x[1]*(a32a44 - a34a42) +
  57       x[8]*x[5]*(a34a41 - a31a44) + x[12]*x[1]*(a33a42 - a32a43) +
  58       x[12]*x[5]*(a31a43 - a33a41) + x[12]*x[9]*(a32a41 - a31a42);
  59
  60     return det;
  61 }
  62
  63 #define CHECK_PAIRWISE_DELETION\
  64     if (KnownBase(x[s1]) && KnownBase(x[s2])) L++;\
  65     else continue;
  66
  67 void distDNA_raw(unsigned char *x, int *n, int *s, double *d, int scaled)
  68 {
  69     int i1, i2, s1, s2, target, Nd;
  70
  71     target = 0;
  72     for (i1 = 1; i1 < *n; i1++) {
  73         for (i2 = i1 + 1; i2 <= *n; i2++) {
  74             Nd = 0;
  75             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n)
  76               if (DifferentBase(x[s1], x[s2])) Nd++;
  77             if (scaled) d[target] = ((double) Nd / *s);
  78             else d[target] = ((double) Nd);
  79             target++;
  80         }
  81     }
  82 }
  83
  84 void distDNA_raw_pairdel(unsigned char *x, int *n, int *s, double *d, int scaled)
  85 {
  86     int i1, i2, s1, s2, target, Nd, L;
  87
  88     target = 0;
  89     for (i1 = 1; i1 < *n; i1++) {
  90         for (i2 = i1 + 1; i2 <= *n; i2++) {
  91             Nd = L = 0;
  92             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
  93                 CHECK_PAIRWISE_DELETION
  94                 if (DifferentBase(x[s1], x[s2])) Nd++;
  95             }
  96             if (scaled) d[target] = ((double) Nd/L);
  97             else d[target] = ((double) Nd);
  98             target++;
  99         }
 100     }
 101 }
 102
 103 #define COMPUTE_DIST_JC69\
 104     p = ((double) Nd/L);\
 105     if (*gamma)\
 106       d[target] = 0.75 * *alpha*(pow(1 - 4*p/3, -1/ *alpha) - 1);\
 107     else d[target] = -0.75*log(1 - 4*p/3);\
 108     if (*variance) {\
 109         if (*gamma) var[target] = p*(1 - p)/(pow(1 - 4*p/3, -2/(*alpha + 1)) * L);\
 110         else var[target] = p*(1 - p)/(pow(1 - 4*p/3, 2)*L);\
 111     }
 112
 113 void distDNA_JC69(unsigned char *x, int *n, int *s, double *d,
 114                   int *variance, double *var, int *gamma, double *alpha)
 115 {
 116     int i1, i2, s1, s2, target, Nd, L;
 117     double p;
 118
 119     L = *s;
 120
 121     target = 0;
 122     for (i1 = 1; i1 < *n; i1++) {
 123         for (i2 = i1 + 1; i2 <= *n; i2++) {
 124             Nd = 0;
 125             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n)
 126               if (DifferentBase(x[s1], x[s2])) Nd++;
 127             COMPUTE_DIST_JC69
 128             target++;
 129         }
 130     }
 131 }
 132
 133 void distDNA_JC69_pairdel(unsigned char *x, int *n, int *s, double *d,
 134                           int *variance, double *var, int *gamma, double *alpha)
 135 {
 136     int i1, i2, s1, s2, target, Nd, L;
 137     double p;
 138
 139     target = 0;
 140     for (i1 = 1; i1 < *n; i1++) {
 141         for (i2 = i1 + 1; i2 <= *n; i2++) {
 142             Nd = L = 0;
 143             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 144                 CHECK_PAIRWISE_DELETION
 145                 if (DifferentBase(x[s1], x[s2])) Nd++;
 146             }
 147             COMPUTE_DIST_JC69
 148             target++;
 149         }
 150     }
 151 }
 152
 153 #define COUNT_TS_TV\
 154     if (SameBase(x[s1], x[s2])) continue;\
 155     Nd++;\
 156     if (IsPurine(x[s1]) && IsPurine(x[s2])) {\
 157         Ns++;\
 158         continue;\
 159     }\
 160     if (IsPyrimidine(x[s1]) && IsPyrimidine(x[s2])) Ns++;
 161
 162 #define COMPUTE_DIST_K80\
 163     P = ((double) Ns/L);\
 164     Q = ((double) (Nd - Ns)/L);\
 165     a1 = 1 - 2*P - Q;\
 166     a2 = 1 - 2*Q;\
 167     if (*gamma) {\
 168         b = -1 / *alpha;\
 169         d[target] = *alpha * (pow(a1, b) + 0.5*pow(a2, b) - 1.5)/2;\
 170     }\
 171     else d[target] = -0.5 * log(a1 * sqrt(a2));\
 172     if (*variance) {\
 173         if (*gamma) {\
 174             b = -(1 / *alpha + 1);\
 175             c1 = pow(a1, b);\
 176             c2 = pow(a2, b);\
 177             c3 = (c1 + c2)/2;\
 178         } else {\
 179           c1 = 1/a1;\
 180           c2 = 1/a2;\
 181           c3 = (c1 + c2)/2;\
 182         }\
 183         var[target] = (c1*c1*P + c3*c3*Q - pow(c1*P + c3*Q, 2))/L;\
 184     }
 185
 186 void distDNA_K80(unsigned char *x, int *n, int *s, double *d,
 187                  int *variance, double *var, int *gamma, double *alpha)
 188 {
 189     int i1, i2, s1, s2, target, Nd, Ns, L;
 190     double P, Q, a1, a2, b, c1, c2, c3;
 191
 192     L = *s;
 193
 194     target = 0;
 195     for (i1 = 1; i1 < *n; i1++) {
 196         for (i2 = i1 + 1; i2 <= *n; i2++) {
 197             Nd = Ns = 0;
 198             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 199                 COUNT_TS_TV
 200             }
 201             COMPUTE_DIST_K80
 202             target++;
 203         }
 204     }
 205 }
 206
 207 void distDNA_K80_pairdel(unsigned char *x, int *n, int *s, double *d,
 208                          int *variance, double *var, int *gamma, double *alpha)
 209 {
 210     int i1, i2, s1, s2, target, Nd, Ns, L;
 211     double P, Q, a1, a2, b, c1, c2, c3;
 212
 213     target = 0;
 214     for (i1 = 1; i1 < *n; i1++) {
 215         for (i2 = i1 + 1; i2 <= *n; i2++) {
 216             Nd = Ns = L = 0;
 217             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 218                 CHECK_PAIRWISE_DELETION
 219                 COUNT_TS_TV
 220             }
 221             COMPUTE_DIST_K80
 222             target++;
 223         }
 224     }
 225 }
 226
 227 #define COMPUTE_DIST_F81\
 228     p = ((double) Nd/L);\
 229     if (*gamma) d[target] = E * *alpha * (pow(1 - p/E, -1/ *alpha) - 1);\
 230     else d[target] = -E*log(1 - p/E);\
 231     if (*variance) {\
 232         if (*gamma) var[target] = p*(1 - p)/(pow(1 - p/E, -2/(*alpha + 1)) * L);\
 233         else var[target] = p*(1 - p)/(pow(1 - p/E, 2)*L);\
 234     }
 235
 236 void distDNA_F81(unsigned char *x, int *n, int *s, double *d, double *BF,
 237                  int *variance, double *var, int *gamma, double *alpha)
 238 {
 239     int i1, i2, s1, s2, target, Nd, L;
 240     double p, E;
 241
 242     L = *s;
 243     E = 1 - BF[0]*BF[0] - BF[1]*BF[1] - BF[2]*BF[2] - BF[3]*BF[3];
 244
 245     target = 0;
 246     for (i1 = 1; i1 < *n; i1++) {
 247         for (i2 = i1 + 1; i2 <= *n; i2++) {
 248             Nd = 0;
 249             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n)
 250               if (DifferentBase(x[s1], x[s2])) Nd++;
 251             COMPUTE_DIST_F81
 252             target++;
 253         }
 254     }
 255 }
 256
 257 void distDNA_F81_pairdel(unsigned char *x, int *n, int *s, double *d, double *BF,
 258                          int *variance, double *var, int *gamma, double *alpha)
 259 {
 260     int i1, i2, s1, s2, target, Nd, L;
 261     double p, E;
 262
 263     E = 1 - BF[0]*BF[0] - BF[1]*BF[1] - BF[2]*BF[2] - BF[3]*BF[3];
 264
 265     target = 0;
 266     for (i1 = 1; i1 < *n; i1++) {
 267         for (i2 = i1 + 1; i2 <= *n; i2++) {
 268             Nd = L = 0;
 269             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 270                 CHECK_PAIRWISE_DELETION
 271                 if (DifferentBase(x[s1], x[s2])) Nd++;
 272             }
 273             COMPUTE_DIST_F81
 274             target++;
 275         }
 276     }
 277 }
 278
 279 #define COUNT_TS_TV1_TV2\
 280     if (SameBase(x[s1], x[s2])) continue;\
 281     Nd++;\
 282     if ((x[s1] | x[s2]) == 152 || (x[s1] | x[s2]) == 104) {\
 283         Nv1++;\
 284         continue;\
 285     }\
 286     if ((x[s1] | x[s2]) == 168 || (x[s1] | x[s2]) == 88) Nv2++;
 287
 288
 289 #define COMPUTE_DIST_K81\
 290     P = ((double) (Nd - Nv1 - Nv2)/L);\
 291     Q = ((double) Nv1/L);\
 292     R = ((double) Nv2/L);\
 293     a1 = 1 - 2*P - 2*Q;\
 294     a2 = 1 - 2*P - 2*R;\
 295     a3 = 1 - 2*Q - 2*R;\
 296     d[target] = -0.25*log(a1*a2*a3);\
 297     if (*variance) {\
 298         a = (1/a1 + 1/a2)/2;\
 299         b = (1/a1 + 1/a3)/2;\
 300         c = (1/a2 + 1/a3)/2;\
 301       var[target] = (a*a*P + b*b*Q + c*c*R - pow(a*P + b*Q + c*R, 2))/2;\
 302     }
 303
 304 void distDNA_K81(unsigned char *x, int *n, int *s, double *d,
 305                  int *variance, double *var)
 306 {
 307     int i1, i2, Nd, Nv1, Nv2, L, s1, s2, target;
 308     double P, Q, R, a1, a2, a3, a, b, c;
 309
 310     L = *s;
 311
 312     target = 0;
 313     for (i1 = 1; i1 < *n; i1++) {
 314         for (i2 = i1 + 1; i2 <= *n; i2++) {
 315             Nd = Nv1 = Nv2 = 0;
 316             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 317                 COUNT_TS_TV1_TV2
 318             }
 319             COMPUTE_DIST_K81
 320             target++;
 321         }
 322     }
 323 }
 324
 325 void distDNA_K81_pairdel(unsigned char *x, int *n, int *s, double *d,
 326                          int *variance, double *var)
 327 {
 328     int i1, i2, Nd, Nv1, Nv2, L, s1, s2, target;
 329     double P, Q, R, a1, a2, a3, a, b, c;
 330
 331     target = 0;
 332     for (i1 = 1; i1 < *n; i1++) {
 333         for (i2 = i1 + 1; i2 <= *n; i2++) {
 334             Nd = Nv1 = Nv2 = L = 0;
 335             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 336                 CHECK_PAIRWISE_DELETION
 337                 COUNT_TS_TV1_TV2
 338             }
 339             COMPUTE_DIST_K81
 340             target++;
 341         }
 342     }
 343 }
 344
 345 #define PREPARE_BF_F84\
 346     A = (BF[0]*BF[2])/(BF[0] + BF[2]) + (BF[1]*BF[3])/(BF[1] + BF[3]);\
 347     B = BF[0]*BF[2] + BF[1]*BF[3];\
 348     C = (BF[0] + BF[2])*(BF[1] + BF[3]);
 349
 350 #define COMPUTE_DIST_F84\
 351    P = ((double) Ns/L);\
 352    Q = ((double) (Nd - Ns)/L);\
 353    d[target] = -2*A*log(1 - P/(2*A) - (A - B)*Q/(2*A*C)) + 2*(A - B - C)*log(1 - Q/(2*C));\
 354    if (*variance) {\
 355        t1 = A*C;\
 356        t2 = C*P/2;\
 357        t3 = (A - B)*Q/2;\
 358        a = t1/(t1 - t2 - t3);\
 359        b = A*(A - B)/(t1 - t2 - t3) - (A - B - C)/(C - Q/2);\
 360        var[target] = (a*a*P + b*b*Q - pow(a*P + b*Q, 2))/L;\
 361    }
 362
 363 void distDNA_F84(unsigned char *x, int *n, int *s, double *d,
 364                  double *BF, int *variance, double *var)
 365 {
 366     int i1, i2, Nd, Ns, L, target, s1, s2;
 367     double P, Q, A, B, C, a, b, t1, t2, t3;
 368
 369     PREPARE_BF_F84
 370     L = *s;
 371
 372     target = 0;
 373     for (i1 = 1; i1 < *n; i1++) {
 374         for (i2 = i1 + 1; i2 <= *n; i2++) {
 375             Nd = Ns = 0;
 376             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 377                 COUNT_TS_TV
 378             }
 379             COMPUTE_DIST_F84
 380             target++;
 381         }
 382     }
 383 }
 384
 385 void distDNA_F84_pairdel(unsigned char *x, int *n, int *s, double *d,
 386                          double *BF, int *variance, double *var)
 387 {
 388     int i1, i2, Nd, Ns, L, target, s1, s2;
 389     double P, Q, A, B, C, a, b, t1, t2, t3;
 390
 391     PREPARE_BF_F84
 392
 393     target = 0;
 394     for (i1 = 1; i1 < *n; i1++) {
 395         for (i2 = i1 + 1; i2 <= *n; i2++) {
 396             Nd = Ns = L = 0;
 397             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 398                 CHECK_PAIRWISE_DELETION
 399                 COUNT_TS_TV
 400             }
 401             COMPUTE_DIST_F84
 402             target++;
 403         }
 404     }
 405 }
 406
 407 #define COMPUTE_DIST_T92\
 408     P = ((double) Ns/L);\
 409     Q = ((double) (Nd - Ns)/L);\
 410     a1 = 1 - P/wg - Q;\
 411     a2 = 1 - 2*Q;\
 412     d[target] = -wg*log(a1) - 0.5*(1 - wg)*log(a2);\
 413     if (*variance) {\
 414         c1 = 1/a1;\
 415         c2 = 1/a2;\
 416         c3 = wg*(c1 - c2) + c2;\
 417         var[target] = (c1*c1*P + c3*c3*Q - pow(c1*P + c3*Q, 2))/L;\
 418     }
 419
 420 void distDNA_T92(unsigned char *x, int *n, int *s, double *d,
 421                  double *BF, int *variance, double *var)
 422 {
 423     int i1, i2, Nd, Ns, L, target, s1, s2;
 424     double P, Q, wg, a1, a2, c1, c2, c3;
 425
 426     L = *s;
 427     wg = 2 * (BF[1] + BF[2]) * (1 - (BF[1] + BF[2]));
 428
 429     target = 0;
 430     for (i1 = 1; i1 < *n; i1++) {
 431         for (i2 = i1 + 1; i2 <= *n; i2++) {
 432             Nd = Ns = 0;
 433             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 434                 COUNT_TS_TV
 435             }
 436             COMPUTE_DIST_T92
 437             target++;
 438         }
 439     }
 440 }
 441
 442 void distDNA_T92_pairdel(unsigned char *x, int *n, int *s, double *d,
 443                          double *BF, int *variance, double *var)
 444 {
 445     int i1, i2, Nd, Ns, L, target, s1, s2;
 446     double P, Q, wg, a1, a2, c1, c2, c3;
 447
 448     wg = 2 * (BF[1] + BF[2]) * (1 - (BF[1] + BF[2]));
 449
 450     target = 0;
 451     for (i1 = 1; i1 < *n; i1++) {
 452         for (i2 = i1 + 1; i2 <= *n; i2++) {
 453             Nd = Ns = L = 0;
 454             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 455                 CHECK_PAIRWISE_DELETION
 456                 COUNT_TS_TV
 457             }
 458             COMPUTE_DIST_T92
 459             target++;
 460         }
 461     }
 462 }
 463
 464 /* returns 1 if one of the base is adenine and
 465    the other one is guanine surely, 0 otherwise */
 466 #define AdenineAndGuanine(a, b) (a | b) == 200
 467
 468 /* returns 1 if one of the base is cytosine and
 469    the other one is thymine surely, 0 otherwise */
 470 #define CytosineAndThymine(a, b) (a | b) == 56
 471
 472 #define PREPARE_BF_TN93\
 473     gR = BF[0] + BF[2];\
 474     gY = BF[1] + BF[3];\
 475     k1 = 2 * BF[0] * BF[2] / gR;\
 476     k2 = 2 * BF[1] * BF[3] / gY;\
 477     k3 = 2 * (gR * gY - BF[0]*BF[2]*gY/gR - BF[1]*BF[3]*gR/gY);
 478
 479 #define COUNT_TS1_TS2_TV\
 480     if (DifferentBase(x[s1], x[s2])) {\
 481         Nd++;\
 482         if (AdenineAndGuanine(x[s1], x[s2])) {\
 483             Ns1++;\
 484         continue;\
 485         }\
 486         if (CytosineAndThymine(x[s1], x[s2])) Ns2++;\
 487     }
 488
 489 #define COMPUTE_DIST_TN93\
 490     P1 = ((double) Ns1/L);\
 491     P2 = ((double) Ns2/L);\
 492     Q = ((double) (Nd - Ns1 - Ns2)/L);\
 493     w1 = 1 - P1/k1 - Q/(2*gR);\
 494     w2 = 1 - P2/k2 - Q/(2*gY);\
 495     w3 = 1 - Q/(2*gR*gY);\
 496     if (*gamma) {\
 497         k4 = 2*(BF[0]*BF[2] + BF[1]*BF[3] + gR*gY);\
 498         b = -1 / *alpha;\
 499         c1 = pow(w1, b);\
 500         c2 = pow(w2, b);\
 501         c3 = pow(w3, b);\
 502         c4 = k1*c1/(2*gR) + k2*c2/(2*gY) + k3*c3/(2*gR*gY);\
 503         d[target] = *alpha * (k1*pow(w1, b) + k2*pow(w2, b) + k3*pow(w3, b) - k4);\
 504     } else {\
 505         k4 = 2*((BF[0]*BF[0] + BF[2]*BF[2])/(2*gR*gR) + (BF[2]*BF[2] + BF[3]*BF[3])/(2*gY*gY));\
 506         c1 = 1/w1;\
 507         c2 = 1/w2;\
 508         c3 = 1/w3;\
 509         c4 = k1 * c1/(2 * gR) + k2 * c2/(2 * gY) + k4 * c3;\
 510         d[target] = -k1*log(w1) - k2*log(w2) - k3*log(w3);\
 511     }\
 512     if (*variance)\
 513       var[target] = (c1*c1*P1 + c2*c2*P2 + c4*c4*Q - pow(c1*P1 + c2*P2 + c4*Q, 2))/L;
 514
 515 void distDNA_TN93(unsigned char *x, int *n, int *s, double *d,
 516                   double *BF, int *variance, double *var,
 517                   int *gamma, double *alpha)
 518 {
 519     int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2;
 520     double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
 521
 522     L = *s;
 523
 524     PREPARE_BF_TN93
 525
 526     target = 0;
 527     for (i1 = 1; i1 < *n; i1++) {
 528         for (i2 = i1 + 1; i2 <= *n; i2++) {
 529             Nd = Ns1 = Ns2 = 0;
 530             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 531                 COUNT_TS1_TS2_TV
 532             }
 533             COMPUTE_DIST_TN93
 534             target++;
 535         }
 536     }
 537 }
 538
 539 void distDNA_TN93_pairdel(unsigned char *x, int *n, int *s, double *d,
 540                           double *BF, int *variance, double *var,
 541                           int *gamma, double *alpha)
 542 {
 543     int i1, i2, k, Nd, Ns1, Ns2, L, target, s1, s2;
 544     double P1, P2, Q, A, B, C, gR, gY, k1, k2, k3, k4, w1, w2, w3, c1, c2, c3, c4, b;
 545
 546     PREPARE_BF_TN93
 547
 548     target = 0;
 549     for (i1 = 1; i1 < *n; i1++) {
 550         for (i2 = i1 + 1; i2 <= *n; i2++) {
 551             Nd = Ns1 = Ns2 = L = 0;
 552             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 553                 CHECK_PAIRWISE_DELETION
 554                 COUNT_TS1_TS2_TV
 555             }
 556             COMPUTE_DIST_TN93
 557             target++;
 558         }
 559     }
 560 }
 561
 562 void distDNA_GG95(unsigned char *x, int *n, int *s, double *d,
 563                   int *variance, double *var)
 564 {
 565     int i1, i2, s1, s2, target, GC, Nd, Ns, tl, npair;
 566     double *theta, gcprop, *P, pp, *Q, qq, *tstvr, svr, A, sum, ma /* mean alpha */, K1, K2;
 567
 568     theta = &gcprop;
 569     P = &pp;
 570     Q = &qq;
 571     tstvr = &svr;
 572
 573     npair = *n * (*n - 1) / 2;
 574
 575     theta = (double*)R_alloc(*n, sizeof(double));
 576     P = (double*)R_alloc(npair, sizeof(double));
 577     Q = (double*)R_alloc(npair, sizeof(double));
 578     tstvr = (double*)R_alloc(npair, sizeof(double));
 579
 580     /* get the proportion of GC (= theta) in each sequence */
 581     for (i1 = 1; i1 <= *n; i1++) {
 582         GC = 0;
 583         for (s1 = i1 - 1; s1 < i1 + *n*(*s - 1); s1 += *n)
 584           if (IsCytosine(x[s1]) || IsGuanine(x[s1])) GC += 1;
 585         theta[i1 - 1] = ((double) GC / *s);
 586     }
 587
 588     /* get the proportions of transitions and transversions,
 589        and the estimates of their ratio for each pair */
 590     target = 0;
 591     for (i1 = 1; i1 < *n; i1++) {
 592         for (i2 = i1 + 1; i2 <= *n; i2++) {
 593             Nd = Ns = 0;
 594             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 595                 COUNT_TS_TV
 596             }
 597             P[target] = ((double) Ns / *s);
 598             Q[target] = ((double) (Nd - Ns) / *s);
 599             A = log(1 - 2*Q[target]);
 600             tstvr[target] = 2*(log(1 - 2*P[target] - Q[target]) - 0.5*A)/A;
 601             target++;
 602         }
 603     }
 604
 605     /* compute the mean alpha (ma) = mean Ts/Tv */
 606     sum = 0;
 607     tl = 0;
 608     for (i1 = 0; i1 < npair; i1++)
 609     /* some values of tstvr are -Inf if there is no
 610        transversions observed: we exclude them */
 611       if (R_FINITE(tstvr[i1])) {
 612           sum += tstvr[i1];
 613           tl += 1;
 614       }
 615     ma = sum/tl;
 616
 617     /* compute the distance for each pair */
 618     target = 0;
 619     for (i1 = 1; i1 < *n; i1++) {
 620         for (i2 = i1 + 1; i2 <= *n; i2++) {
 621             A = 1 - 2*Q[target];
 622             K1 = 1 + ma*(theta[i1 - 1]*(1 - theta[i1 - 1]) + theta[i2 - 1]*(1 - theta[i2 - 1]));
 623             K2 = ma*pow(theta[i1 - 1] - theta[i2 - 1], 2)/(ma + 1);
 624             d[target] = -0.5*K1*log(A) + K2*(1 - pow(A, 0.25*(ma + 1)));
 625             if (*variance)
 626               var[target] = pow(K1 + K2*0.5*(ma + 1)*pow(A, 0.25*(ma + 1)), 2)*Q[target]*(1 - Q[target])/(A*A * *s);
 627             target++;
 628         }
 629     }
 630 }
 631
 632 void distDNA_GG95_pairdel(unsigned char *x, int *n, int *s, double *d,
 633                           int *variance, double *var)
 634 {
 635     int i1, i2, s1, s2, target, *L, length, GC, Nd, Ns, tl, npair;
 636     double *theta, gcprop, *P, pp, *Q, qq, *tstvr, svr, A, sum, ma /* mean alpha */, K1, K2;
 637
 638     theta = &gcprop;
 639     L = &length;
 640     P = &pp;
 641     Q = &qq;
 642     tstvr = &svr;
 643
 644     npair = *n * (*n - 1) / 2;
 645
 646     theta = (double*)R_alloc(*n, sizeof(double));
 647     L = (int*)R_alloc(npair, sizeof(int));
 648     P = (double*)R_alloc(npair, sizeof(double));
 649     Q = (double*)R_alloc(npair, sizeof(double));
 650     tstvr = (double*)R_alloc(npair, sizeof(double));
 651
 652     /* get the proportion of GC (= theta) in each sequence */
 653     for (i1 = 1; i1 <= *n; i1++) {
 654         tl = GC = 0;
 655         for (s1 = i1 - 1; s1 < i1 + *n*(*s - 1); s1 += *n) {
 656             if (KnownBase(x[s1])) tl++;
 657             else continue;
 658             if (IsCytosine(x[s1]) || IsGuanine(x[s1])) GC += 1;
 659         }
 660         theta[i1 - 1] = ((double) GC / tl);
 661     }
 662
 663     /* get the proportions of transitions and transversions,
 664        and the estimates of their ratio for each pair; we
 665        also get the sample size for each pair in L */
 666     target = 0;
 667     for (i1 = 1; i1 < *n; i1++) {
 668         for (i2 = i1 + 1; i2 <= *n; i2++) {
 669             Nd = Ns = L[target] = 0;
 670             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 671                 if (KnownBase(x[s1]) && KnownBase(x[s2])) L[target]++;
 672                 else continue;
 673                 COUNT_TS_TV
 674             }
 675             P[target] = ((double) Ns/L[target]);
 676             Q[target] = ((double) (Nd - Ns)/L[target]);
 677             A = log(1 - 2*Q[target]);
 678             tstvr[target] = 2*(log(1 - 2*P[target] - Q[target]) - 0.5*A)/A;
 679             target++;
 680         }
 681     }
 682
 683     /* compute the mean alpha (ma) = mean Ts/Tv */
 684     sum = 0;
 685     tl = 0;
 686     for (i1 = 0; i1 < npair; i1++)
 687     /* some values of tstvr are -Inf if there is no
 688        transversions observed: we exclude them */
 689       if (R_FINITE(tstvr[i1])) {
 690           sum += tstvr[i1];
 691           tl += 1;
 692       }
 693     ma = sum/tl;
 694
 695     /* compute the distance for each pair */
 696     target = 0;
 697     for (i1 = 1; i1 < *n; i1++) {
 698         for (i2 = i1 + 1; i2 <= *n; i2++) {
 699             A = 1 - 2*Q[target];
 700             K1 = 1 + ma*(theta[i1 - 1]*(1 - theta[i1 - 1]) + theta[i2 - 1]*(1 - theta[i2 - 1]));
 701             K2 = ma*pow(theta[i1 - 1] - theta[i2 - 1], 2)/(ma + 1);
 702             d[target] = -0.5*K1*log(A) + K2*(1 - pow(A, 0.25*(ma + 1)));
 703             if (*variance)
 704               var[target] = pow(K1 + K2*0.5*(ma + 1)*pow(A, 0.25*(ma + 1)), 2)*Q[target]*(1 - Q[target])/(A*A*L[target]);
 705             target++;
 706         }
 707     }
 708 }
 709
 710 #define DO_CONTINGENCY_NUCLEOTIDES\
 711     switch (x[s1]) {\
 712     case 136 : m = 0; break;\
 713     case 72 : m = 1; break;\
 714     case 40 : m = 2; break;\
 715     case 24 : m = 3; break;\
 716     }\
 717     switch (x[s2]) {\
 718     case 72 : m += 4; break;\
 719     case 40 : m += 8; break;\
 720     case 24 : m += 12; break;\
 721     }\
 722     Ntab[m]++;
 723
 724 #define COMPUTE_DIST_LogDet\
 725     for (k = 0; k < 16; k++) Ftab[k] = ((double) Ntab[k]/L);\
 726     d[target] = (-log(detFourByFour(Ftab))/4 - LN4)/4;\
 727     if (*variance) {\
 728         /* For the inversion, we first make U an identity matrix */\
 729         for (k = 1; k < 15; k++) U[k] = 0;\
 730         U[0] = U[5] = U[10] = U[15] = 1;\
 731         /* The matrix is not symmetric, so we use 'dgesv'. */\
 732         /* This subroutine puts the result in U. */\
 733         F77_CALL(dgesv)(&ndim, &ndim, Ftab, &ndim, ipiv, U, &ndim, &info);\
 734         var[target] = (U[0]*U[0]*Ftab[0] + U[1]*U[1]*Ftab[4] +\
 735                        U[2]*U[2]*Ftab[8] + U[3]*U[3]*Ftab[12] +\
 736                        U[4]*U[4]*Ftab[1] + U[5]*U[5]*Ftab[5] +\
 737                        U[6]*U[6]*Ftab[9] + U[7]*U[7]*Ftab[13] +\
 738                        U[8]*U[8]*Ftab[2] + U[9]*U[9]*Ftab[6] +\
 739                        U[10]*U[10]*Ftab[10] + U[11]*U[11]*Ftab[14] +\
 740                        U[12]*U[12]*Ftab[3] + U[13]*U[13]*Ftab[7] +\
 741                        U[14]*U[14]*Ftab[11] + U[15]*U[15]*Ftab[15] - 16)/(L*16);\
 742     }
 743
 744 void distDNA_LogDet(unsigned char *x, int *n, int *s, double *d,
 745                     int *variance, double *var)
 746 {
 747     int i1, i2, k, m, s1, s2, target, L, Ntab[16], ndim = 4, info, ipiv[16];
 748     double Ftab[16], U[16];
 749
 750     L = *s;
 751
 752     target = 0;
 753     for (i1 = 1; i1 < *n; i1++) {
 754         for (i2 = i1 + 1; i2 <= *n; i2++) {
 755             for (k = 0; k < 16; k++) Ntab[k] = 0;
 756             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 757                 DO_CONTINGENCY_NUCLEOTIDES
 758             }
 759             COMPUTE_DIST_LogDet
 760             target++;
 761         }
 762     }
 763 }
 764
 765 void distDNA_LogDet_pairdel(unsigned char *x, int *n, int *s, double *d,
 766                             int *variance, double *var)
 767 {
 768     int i1, i2, k, m, s1, s2, target, L, Ntab[16], ndim = 4, info, ipiv[16];
 769     double Ftab[16], U[16];
 770
 771     target = 0;
 772     for (i1 = 1; i1 < *n; i1++) {
 773         for (i2 = i1 + 1; i2 <= *n; i2++) {
 774             for (k = 0; k < 16; k++) Ntab[k] = 0;
 775             L = 0;
 776             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 777                 CHECK_PAIRWISE_DELETION
 778                 DO_CONTINGENCY_NUCLEOTIDES
 779             }
 780             COMPUTE_DIST_LogDet
 781             target++;
 782         }
 783     }
 784 }
 785
 786 void distDNA_BH87(unsigned char *x, int *n, int *s, double *d,
 787                   int *variance, double *var)
 788 /* <FIXME>
 789    For the moment there is no need to check for pairwise deletions
 790    since DO_CONTINGENCY_NUCLEOTIDES considers only the known nucleotides.
 791    In effect the pairwise deletion has possibly been done before.
 792    The sequence length(s) are used only to compute the variances, which is
 793    currently not available.
 794    </FIXME> */
 795 {
 796     int i1, i2, k, kb, s1, s2, m, Ntab[16], ROWsums[4], ndim = 4, info, ipiv[16];
 797     double P12[16], P21[16], U[16];
 798
 799     for (i1 = 1; i1 < *n; i1++) {
 800         for (i2 = i1 + 1; i2 <= *n; i2++) {
 801             for (k = 0; k < 16; k++) Ntab[k] = 0;
 802             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 803                 DO_CONTINGENCY_NUCLEOTIDES
 804             }
 805
 806             /* get the rowwise sums of Ntab */
 807             ROWsums[0] = Ntab[0] + Ntab[4] + Ntab[8] + Ntab[12];
 808             ROWsums[1] = Ntab[1] + Ntab[5] + Ntab[9] + Ntab[13];
 809             ROWsums[2] = Ntab[2] + Ntab[6] + Ntab[10] + Ntab[14];
 810             ROWsums[3] = Ntab[3] + Ntab[7] + Ntab[11] + Ntab[15];
 811
 812             for (k = 0; k < 16; k++)
 813               P12[k] = ((double) Ntab[k]);
 814
 815             /* scale each element of P12 by its rowwise sum */
 816             for (k = 0; k < 4; k++)
 817               for (kb = 0; kb < 16; kb += 4)
 818                 P12[k + kb] = P12[k + kb]/ROWsums[k];
 819
 820             d[*n*(i2 - 1) + i1 - 1] = -log(detFourByFour(P12))/4;
 821
 822             /* compute the columnwise sums of Ntab: these
 823                are the rowwise sums of its transpose */
 824             ROWsums[0] = Ntab[0] + Ntab[1] + Ntab[2] + Ntab[3];
 825             ROWsums[1] = Ntab[4] + Ntab[5] + Ntab[6] + Ntab[7];
 826             ROWsums[2] = Ntab[8] + Ntab[9] + Ntab[10] + Ntab[11];
 827             ROWsums[3] = Ntab[12] + Ntab[13] + Ntab[14] + Ntab[15];
 828
 829             /* transpose Ntab and store the result in P21 */
 830             for (k = 0; k < 4; k++)
 831                for (kb = 0; kb < 4; kb++)
 832                  P21[kb + 4*k] = Ntab[k + 4*kb];
 833
 834             /* scale as above */
 835             for (k = 0; k < 4; k++)
 836               for (kb = 0; kb < 16; kb += 4)
 837                 P21[k + kb] = P21[k + kb]/ROWsums[k];
 838
 839             d[*n*(i1 - 1) + i2 - 1] = -log(detFourByFour(P21))/4;
 840         }
 841     }
 842 }
 843
 844 #define COMPUTE_DIST_ParaLin\
 845     for (k = 0; k < 16; k++) Ftab[k] = ((double) Ntab[k]/L);\
 846     d[target] = -log(detFourByFour(Ftab)/\
 847                      sqrt(find[0][i1 - 1]*find[1][i1 - 1]*find[2][i1 - 1]*find[3][i1 - 1]*\
 848                           find[0][i2 - 1]*find[1][i2 - 1]*find[2][i2 - 1]*find[3][i2 - 1]))/4;\
 849     if (*variance) {\
 850         /* For the inversion, we first make U an identity matrix */\
 851         for (k = 1; k < 15; k++) U[k] = 0;\
 852         U[0] = U[5] = U[10] = U[15] = 1;\
 853         /* The matrix is not symmetric, so we use 'dgesv'. */\
 854         /* This subroutine puts the result in U. */\
 855         F77_CALL(dgesv)(&ndim, &ndim, Ftab, &ndim, ipiv, U, &ndim, &info);\
 856         var[target] = (U[0]*U[0]*Ftab[0] + U[1]*U[1]*Ftab[4] +\
 857                        U[2]*U[2]*Ftab[8] + U[3]*U[3]*Ftab[12] +\
 858                        U[4]*U[4]*Ftab[1] + U[5]*U[5]*Ftab[5] +\
 859                        U[6]*U[6]*Ftab[9] + U[7]*U[7]*Ftab[13] +\
 860                        U[8]*U[8]*Ftab[2] + U[9]*U[9]*Ftab[6] +\
 861                        U[10]*U[10]*Ftab[10] + U[11]*U[11]*Ftab[14] +\
 862                        U[12]*U[12]*Ftab[3] + U[13]*U[13]*Ftab[7] +\
 863                        U[14]*U[14]*Ftab[11] + U[15]*U[15]*Ftab[15] -\
 864                        4*(1/sqrt(find[0][i1 - 1]*find[0][i2 - 1]) +\
 865                        1/sqrt(find[1][i1 - 1]*find[1][i2 - 1]) +\
 866                        1/sqrt(find[2][i1 - 1]*find[2][i2 - 1]) +\
 867                        1/sqrt(find[3][i1 - 1]*find[3][i2 - 1])))/(L*16);\
 868     }
 869
 870 void distDNA_ParaLin(unsigned char *x, int *n, int *s, double *d,
 871                      int *variance, double *var)
 872 {
 873     int i1, i2, k, s1, s2, m, target, L, Ntab[16], ndim = 4, info, ipiv[16];
 874     double Ftab[16], U[16], *find[4];
 875
 876     L = *s;
 877
 878     for (k = 0; k < 4; k++)
 879       find[k] = (double*)R_alloc(*n, sizeof(double));
 880
 881     for (i1 = 0; i1 < *n; i1++)
 882       for (k = 0; k < 4; k++) find[k][i1] = 0.0;
 883
 884     for (i1 = 0; i1 < *n; i1++) {
 885         for (s1 = i1; s1 < i1 + *n*(*s - 1) + 1; s1+= *n) {
 886             switch (x[s1]) {
 887             case 136 : find[0][i1]++; break;
 888             case 40 : find[1][i1]++; break;
 889             case 72 : find[2][i1]++; break;
 890             case 24 : find[3][i1]++; break;
 891             }
 892         }
 893         for (k = 0; k < 4; k++) find[k][i1] /= L;
 894     }
 895
 896     target = 0;
 897     for (i1 = 1; i1 < *n; i1++) {
 898         for (i2 = i1 + 1; i2 <= *n; i2++) {
 899             for (k = 0; k < 16; k++) Ntab[k] = 0;
 900             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 901                 DO_CONTINGENCY_NUCLEOTIDES
 902             }
 903             COMPUTE_DIST_ParaLin
 904             target++;
 905         }
 906     }
 907 }
 908
 909 void distDNA_ParaLin_pairdel(unsigned char *x, int *n, int *s, double *d,
 910                              int *variance, double *var)
 911 {
 912     int i1, i2, k, s1, s2, m, target, L, Ntab[16], ndim = 4, info, ipiv[16];
 913     double Ftab[16], U[16], *find[4];
 914
 915     L = 0;
 916
 917     for (k = 0; k < 4; k++)
 918       find[k] = (double*)R_alloc(*n, sizeof(double));
 919
 920     for (i1 = 0; i1 < *n; i1++)
 921       for (k = 0; k < 4; k++) find[k][i1] = 0.0;
 922
 923     for (i1 = 0; i1 < *n; i1++) {
 924         L = 0;
 925         for (s1 = i1; s1 < i1 + *n*(*s - 1) + 1; s1+= *n) {
 926             if (KnownBase(x[s1])) {
 927                 L++;
 928                 switch (x[s1]) {
 929                 case 136 : find[0][i1]++; break;
 930                 case 40 : find[1][i1]++; break;
 931                 case 72 : find[2][i1]++; break;
 932                 case 24 : find[3][i1]++; break;
 933                 }
 934             }
 935         }
 936         for (k = 0; k < 4; k++) find[k][i1] /= L;
 937     }
 938
 939     target = 0;
 940     for (i1 = 1; i1 < *n; i1++) {
 941         for (i2 = i1 + 1; i2 <= *n; i2++) {
 942             L = 0;
 943             for (k = 0; k < 16; k++) Ntab[k] = 0;
 944             for (s1 = i1 - 1, s2 = i2 - 1; s1 < i1 + *n*(*s - 1); s1+= *n, s2 += *n) {
 945                 CHECK_PAIRWISE_DELETION
 946                 DO_CONTINGENCY_NUCLEOTIDES
 947             }
 948             COMPUTE_DIST_ParaLin
 949             target++;
 950         }
 951     }
 952 }
 953
 954 void BaseProportion(unsigned char *x, int *n, double *BF, int *freq)
 955 {
 956     int i, m;
 957
 958     m = 0;
 959     for (i = 0; i < *n; i++) {
 960         if (KnownBase(x[i])) {
 961             m++;
 962             switch (x[i]) {
 963             case 136 : BF[0]++; break;
 964             case 40 : BF[1]++; break;
 965             case 72 : BF[2]++; break;
 966             case 24 : BF[3]++; break;
 967             }
 968         }
 969     }
 970     if (! *freq) for (i = 0; i < 4; i++) BF[i] /= m;
 971 }
 972
 973 void SegSites(unsigned char *x, int *n, int *s, int *seg)
 974 {
 975     int i, j;
 976     unsigned char basis;
 977
 978     for (j = 0; j < *s; j++) {
 979         i = *n * j;
 980         while (!KnownBase(x[i])) i++;
 981         basis = x[i];
 982         i++;
 983         while (i < *n * (j + 1)) {
 984             if (!KnownBase(x[i]) || x[i] == basis) i++;
 985             else {
 986                 seg[j] = 1;
 987                 break;
 988             }
 989         }
 990     }
 991 }
 992
 993 void GlobalDeletionDNA(unsigned char *x, int *n, int *s, int *keep)
 994 {
 995     int i, j;
 996
 997     for (j = 0; j < *s; j++) {
 998         i = *n * j;
 999         while (i < *n * (j + 1)) {
1000             if (KnownBase(x[i])) i++;
1001             else {
1002                 keep[j] = 0;
1003                 break;
1004             }
1005         }
1006     }
1007 }
1008
1009 void dist_dna(unsigned char *x, int *n, int *s, int *model, double *d,
1010               double *BF, int *pairdel, int *variance, double *var,
1011               int *gamma, double *alpha)
1012 {
1013     switch (*model) {
1014     case 1 : if (pairdel) distDNA_raw_pairdel(x, n, s, d, 1);
1015              else distDNA_raw(x, n, s, d, 1); break;
1016
1017     case 2 : if (pairdel) distDNA_JC69_pairdel(x, n, s, d, variance, var, gamma, alpha);
1018              else distDNA_JC69(x, n, s, d, variance, var, gamma, alpha); break;
1019
1020     case 3 : if (pairdel) distDNA_K80_pairdel(x, n, s, d, variance, var, gamma, alpha);
1021              else distDNA_K80(x, n, s, d, variance, var, gamma, alpha); break;
1022
1023     case 4 : if (pairdel) distDNA_F81_pairdel(x, n, s, d, BF, variance, var, gamma, alpha);
1024              else distDNA_F81(x, n, s, d, BF, variance, var, gamma, alpha); break;
1025
1026     case 5 : if (pairdel) distDNA_K81_pairdel(x, n, s, d, variance, var);
1027              else distDNA_K81(x, n, s, d, variance, var); break;
1028
1029     case 6 : if (pairdel) distDNA_F84_pairdel(x, n, s, d, BF, variance, var);
1030              else distDNA_F84(x, n, s, d, BF, variance, var); break;
1031
1032     case 7 : if (pairdel) distDNA_T92_pairdel(x, n, s, d, BF, variance, var);
1033              else distDNA_T92(x, n, s, d, BF, variance, var); break;
1034
1035     case 8 : if (pairdel) distDNA_TN93_pairdel(x, n, s, d, BF, variance, var, gamma, alpha);
1036              else distDNA_TN93(x, n, s, d, BF, variance, var, gamma, alpha); break;
1037
1038     case 9 : if (pairdel) distDNA_GG95_pairdel(x, n, s, d, variance, var);
1039              else distDNA_GG95(x, n, s, d, variance, var); break;
1040
1041     case 10 : if (pairdel) distDNA_LogDet_pairdel(x, n, s, d, variance, var);
1042               else distDNA_LogDet(x, n, s, d, variance, var); break;
1043
1044     case 11 : distDNA_BH87(x, n, s, d, variance, var); break;
1045
1046     case 12 : if (pairdel) distDNA_ParaLin_pairdel(x, n, s, d, variance, var);
1047               else distDNA_ParaLin(x, n, s, d, variance, var); break;
1048     case 13 : if (pairdel) distDNA_raw_pairdel(x, n, s, d, 0);
1049              else distDNA_raw(x, n, s, d, 0); break;
1050     }
1051 }