src/dist_dna.c

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