cistematic/core/motifextension.c

   1 /*##########################################################################*/
   2 /*#                                                                         #*/
   3 /*# C O P Y R I G H T   N O T I C E                                         #*/
   4 /*#  Copyright (c) 2003-10 by:                                              #*/
   5 /*#    * California Institute of Technology                                 #*/
   6 /*#                                                                         #*/
   7 /*#    All Rights Reserved.                                                 #*/
   8 /*#                                                                         #*/
   9 /*# Permission is hereby granted, free of charge, to any person             #*/
  10 /*# obtaining a copy of this software and associated documentation files    #*/
  11 /*# (the "Software"), to deal in the Software without restriction,          #*/
  12 /*# including without limitation the rights to use, copy, modify, merge,    #*/
  13 /*# publish, distribute, sublicense, and/or sell copies of the Software,    #*/
  14 /*# and to permit persons to whom the Software is furnished to do so,       #*/
  15 /*# subject to the following conditions:                                    #*/
  16 /*#                                                                         #*/
  17 /*# The above copyright notice and this permission notice shall be          #*/
  18 /*# included in all copies or substantial portions of the Software.         #*/
  19 /*#                                                                         #*/
  20 /*# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,         #*/
  21 /*# EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF      #*/
  22 /*# MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND                   #*/
  23 /*# NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS     #*/
  24 /*# BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN      #*/
  25 /*# ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN       #*/
  26 /*# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE        #*/
  27 /*# SOFTWARE.                                                               #*/
  28 /*###########################################################################*/
  29 /*# C extension for commonly used motif methods.                             */
  30
  31 #include <Python.h>
  32 #include <stdio.h>
  33 #include <stdlib.h>
  34 #include <math.h>
  35
  36 static char module_doc[] =
  37 "This module only implements the locateMotif and correlateMotifs in C for now.";
  38
  39 double probMatchPWM_func(double *PWM, int index, char *aSeq, long seqPos) {
  40     double aval, cval, gval, tval;
  41
  42     aval = PWM[4 * index];
  43     cval = PWM[4 * index + 1];
  44     gval = PWM[4 * index + 2];
  45     tval = PWM[4 * index + 3];
  46
  47     switch (aSeq[seqPos]) {
  48         case 'A': case 'a':
  49             return aval;
  50             break;
  51         case 'C': case 'c':
  52             return cval;
  53             break;
  54         case 'G': case 'g':
  55             return gval;
  56             break;
  57         case 'T': case 't':
  58             return tval;
  59             break;
  60         case 'N': case 'n':
  61             return 1.0;
  62             break;
  63         case 'S': case 's':
  64             return cval + gval;
  65             break;
  66         case 'W': case 'w':
  67             return aval + tval;
  68             break;
  69         case 'R': case 'r':
  70             return aval + gval;
  71             break;
  72         case 'Y': case 'y':
  73             return cval + tval;
  74             break;
  75         case 'M': case 'm':
  76             return aval + cval;
  77             break;
  78         case 'K': case 'k':
  79             return gval + tval;
  80             break;
  81         case 'B': case 'b':
  82             return cval + gval + tval;
  83             break;
  84         case 'D': case 'd':
  85             return aval + gval + tval;
  86             break;
  87         case 'H': case 'h':
  88             return aval + cval + tval;
  89             break;
  90         case 'V': case 'v':
  91             return aval + cval + gval;
  92             break;
  93     }
  94     return 0.0;
  95 }
  96
  97 double probMatchDPWM_func(double *PWM, int index, char *aSeq, long seqPos) {
  98     double aval, cval, gval, tval;
  99     int prevIndex;
 100
 101     aval = 0.0;
 102     cval = 0.0;
 103     gval = 0.0;
 104     tval = 0.0;
 105
 106     if (index == 0) {
 107         for (prevIndex = 0; prevIndex < 4; prevIndex++) {
 108             aval += PWM[prevIndex * 4];
 109             cval += PWM[prevIndex * 4 + 1];
 110             gval += PWM[prevIndex * 4 + 2];
 111             tval += PWM[prevIndex * 4 + 3];
 112         }
 113     } else {
 114         if (seqPos > 0) {
 115             switch(aSeq[seqPos - 1]) {
 116                 case 'A': case 'a':
 117                     aval = PWM[16 * index];
 118                     cval = PWM[16 * index + 1];
 119                     gval = PWM[16 * index + 2];
 120                     tval = PWM[16 * index + 3];
 121                     break;
 122                 case 'C': case 'c':
 123                     aval = PWM[16 * index + 4];
 124                     cval = PWM[16 * index + 5];
 125                     gval = PWM[16 * index + 6];
 126                     tval = PWM[16 * index + 7];
 127                     break;
 128                 case 'G': case 'g':
 129                     aval = PWM[16 * index + 8];
 130                     cval = PWM[16 * index + 9];
 131                     gval = PWM[16 * index + 10];
 132                     tval = PWM[16 * index + 11];
 133                     break;
 134                 case 'T': case 't':
 135                     aval = PWM[16 * index + 12];
 136                     cval = PWM[16 * index + 13];
 137                     gval = PWM[16 * index + 14];
 138                     tval = PWM[16 * index + 15];
 139                     break;
 140             }
 141         }
 142     }
 143
 144     switch (aSeq[seqPos]) {
 145         case 'A': case 'a':
 146             return aval;
 147             break;
 148         case 'C': case 'c':
 149             return cval;
 150             break;
 151         case 'G': case 'g':
 152             return gval;
 153             break;
 154         case 'T': case 't':
 155             return tval;
 156             break;
 157     }
 158
 159     return 0.0001;
 160 }
 161
 162 void scoreMot_func(char *theSeq, double *mPWM, double *rPWM, long pos, long motLen, float *bestCons, float *maxDiff, float *score, char *sense)
 163 {
 164     float forScore, revScore;
 165     long index, currentPos;
 166
 167     forScore = 0.0;
 168     revScore = 0.0;
 169     *sense = 'F';
 170
 171     for (index = 0; index < motLen; index++) {
 172         currentPos = pos + index;
 173         forScore += probMatchPWM_func(mPWM, index, theSeq, currentPos);
 174         revScore += probMatchPWM_func(rPWM, index, theSeq, currentPos);
 175     }
 176
 177     if (((forScore + *maxDiff) < *bestCons) && ((revScore + *maxDiff) < *bestCons)) {
 178         *score = -1.0;
 179         return;
 180     }
 181
 182     if (forScore > revScore) {
 183         *score = forScore;
 184     } else {
 185         *score = revScore;
 186         *sense = 'R';
 187     }
 188 }
 189
 190 void scoreDPWM_func(char *theSeq, double *mDPWM, double *rDPWM, long pos, long motLen, float *bestCons, float *score, char *sense)
 191 {
 192     float forScore, revScore;
 193     long index, currentPos;
 194
 195     forScore = 0.0;
 196     revScore = 0.0;
 197     *sense = 'F';
 198
 199     for (index = 0; index < motLen; index++) {
 200         currentPos = pos + index;
 201         forScore -= log(probMatchDPWM_func(mDPWM, index, theSeq, currentPos)) / log(2.0);
 202         revScore -= log(probMatchDPWM_func(rDPWM, index, theSeq, currentPos)) / log(2.0);
 203     }
 204
 205     if ((forScore > *bestCons) && (revScore > *bestCons)) {
 206         *score = -1.0;
 207         return;
 208     }
 209
 210     if (forScore < revScore) {
 211         *score = forScore;
 212     } else {
 213         *score = revScore;
 214         *sense = 'R';
 215     }
 216 }
 217
 218 static PyObject*
 219 the_func(PyObject *self, PyObject *args)
 220 {
 221     PyObject *pySeq, *motifPWM, *revPWM, *mScore, *mDiff, *results;
 222     double *mPWM, *rPWM;
 223     float maxScore, maxDiff, seqScore;
 224     long pos, maxPos;
 225     int ok, index, indexMax, ntIndex, motLen, seqLen, skipping;
 226     char *seq, sense;
 227
 228     results = Py_BuildValue("[]");
 229
 230     ok = PyArg_UnpackTuple(args, "ref", 5, 5, &pySeq, &motifPWM, &revPWM, &mScore, &mDiff);
 231
 232     seq = PyString_AsString(pySeq);
 233     motLen = PyList_Size(motifPWM);
 234     seqLen = PyString_Size(pySeq);
 235     maxScore = PyFloat_AsDouble(mScore);
 236     maxDiff = PyFloat_AsDouble(mDiff);
 237
 238     indexMax = 4 * motLen;
 239     mPWM = malloc(indexMax * sizeof(double));
 240     rPWM = malloc(indexMax * sizeof(double));
 241
 242     for (index = 0; index < motLen; index++) {
 243         for (ntIndex = 0; ntIndex < 4; ntIndex++) {
 244             mPWM[4 * index + ntIndex] = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(motifPWM, index), ntIndex));
 245             rPWM[4 * index + ntIndex] = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(revPWM, index), ntIndex));
 246         }
 247     }
 248
 249     pos = 0;
 250     maxPos = seqLen - motLen;
 251
 252     while (pos <= maxPos) {
 253         skipping = 0;
 254         for (index = 0; index < motLen; index++) {
 255             if ((seq[pos + index] == 'N') || (seq[pos + index] == 'n')) {
 256                 skipping = index + 1;
 257             }
 258         }
 259
 260         if (skipping) {
 261             pos += skipping;
 262             continue;
 263         }
 264
 265         scoreMot_func(seq, mPWM, rPWM, pos, motLen, &maxScore, &maxDiff, &seqScore, &sense);
 266
 267         if (seqScore > 1.0) {
 268             PyList_Append(results, Py_BuildValue("(i, c)", pos, sense));
 269         }
 270         pos++;
 271     }
 272
 273     free(mPWM);
 274     free(rPWM);
 275     return results;
 276 }
 277
 278 static PyObject*
 279 dpwm_func(PyObject *self, PyObject *args)
 280 {
 281     PyObject *pySeq, *motifPWM, *revPWM, *mScore, *results;
 282     double *mDPWM, *rDPWM;
 283     float maxScore, seqScore;
 284     long pos, maxPos;
 285     int ok, index, indexMax, prevIndex, currentIndex, motLen, seqLen, skipping;
 286     char *seq, sense;
 287
 288     results = Py_BuildValue("[]");
 289
 290     ok = PyArg_UnpackTuple(args, "ref", 4, 4, &pySeq, &motifPWM, &revPWM, &mScore);
 291
 292     seq = PyString_AsString(pySeq);
 293     motLen = PyList_Size(motifPWM);
 294     seqLen = PyString_Size(pySeq);
 295     maxScore = PyFloat_AsDouble(mScore);
 296
 297     indexMax = 16 * motLen;
 298     mDPWM = malloc(indexMax * sizeof(double));
 299     rDPWM = malloc(indexMax * sizeof(double));
 300
 301     for (index = 0; index < motLen; index++) {
 302         for (prevIndex = 0; prevIndex < 4; prevIndex++) {
 303             for (currentIndex = 0; currentIndex < 4; currentIndex++) {
 304                 mDPWM[16 * index + 4 * prevIndex + currentIndex] = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(PyList_GetItem(motifPWM, index), prevIndex), currentIndex));
 305                 rDPWM[16 * index + 4 * prevIndex + currentIndex] = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(PyList_GetItem(revPWM, index), prevIndex), currentIndex));
 306             }
 307         }
 308     }
 309
 310     pos = 0;
 311     maxPos = seqLen - motLen;
 312
 313     while (pos <= maxPos) {
 314         skipping = 0;
 315         for (index = 0; index < motLen; index++) {
 316             if ((seq[pos + index] == 'N') || (seq[pos + index] == 'n')) {
 317                 skipping = index + 1;
 318             }
 319         }
 320
 321         if (skipping) {
 322             pos += skipping;
 323             continue;
 324         }
 325
 326         scoreDPWM_func(seq, mDPWM, rDPWM, pos, motLen, &maxScore, &seqScore, &sense);
 327
 328         if (seqScore > 1.0) {
 329             PyList_Append(results, Py_BuildValue("(i, c)", pos, sense));
 330         }
 331         pos++;
 332     }
 333
 334     free(mDPWM);
 335     free(rDPWM);
 336     return results;
 337 }
 338
 339 static PyObject*
 340 mer_func(PyObject *self, PyObject *args)
 341 {
 342     PyObject *pySeq, *forwardMer, *revcompMer, *maxMismatches, *results;
 343     long pos, maxPos;
 344     int ok, index, motLen, seqLen, skipping, mismatches, fmis, rmis;
 345     char *fMer, *rMer;
 346     char *seq, sense;
 347
 348     results = Py_BuildValue("[]");
 349
 350     ok = PyArg_UnpackTuple(args, "ref", 4, 4, &pySeq, &forwardMer, &revcompMer, &maxMismatches);
 351
 352     seq = PyString_AsString(pySeq);
 353     fMer = PyString_AsString(forwardMer);
 354     rMer = PyString_AsString(revcompMer);
 355     motLen = PyString_Size(forwardMer);
 356     seqLen = PyString_Size(pySeq);
 357     mismatches = PyInt_AsLong(maxMismatches);
 358
 359     pos = 0;
 360     maxPos = seqLen - motLen;
 361
 362     while (pos <= maxPos) {
 363         skipping = 0;
 364         for (index = 0; index < motLen; index++) {
 365             if ((seq[pos + index] == 'N') || (seq[pos + index] == 'n')) {
 366                 skipping = index + 1;
 367             }
 368         }
 369
 370         if (skipping) {
 371             pos += skipping;
 372             continue;
 373         }
 374
 375         fmis = 0;
 376         rmis = 0;
 377         for (index = 0; index < motLen; index++) {
 378             if (seq[pos + index] != fMer[index]) {
 379                 fmis++;
 380             }
 381             if (seq[pos + index] != rMer[index]) {
 382                 rmis++;
 383             }
 384             if ((fmis > mismatches) && (rmis > mismatches)) {
 385                 break;
 386             }
 387         }
 388
 389         if ((fmis <= mismatches) || (rmis <= mismatches)) {
 390             sense = 'F';
 391             if (rmis < fmis) {
 392                 sense = 'R';
 393             }
 394             PyList_Append(results, Py_BuildValue("(i, c)", pos, sense));
 395         }
 396         pos++;
 397     }
 398
 399     return results;
 400 }
 401
 402 double pearson_func(double *colA, double *colB, int pos)
 403 {
 404     double c, numerator;
 405     double meanA, denominatorA;
 406     double meanB, denominatorB;
 407     long index;
 408
 409     meanA = 0.0;
 410     meanB = 0.0;
 411
 412     for (index = 0; index < 4; index++) {
 413         meanA += colA[pos + index];
 414         meanB += colB[pos + index];
 415     }
 416
 417     meanA /= 4;
 418     meanB /= 4;
 419
 420     denominatorA = 0.0;
 421     denominatorB = 0.0;
 422     numerator = 0.0;
 423
 424     for (index = 0; index < 4; index++) {
 425         numerator += (colA[pos + index] - meanA) * (colB[pos + index] - meanB);
 426         denominatorA += (colA[pos + index] - meanA) * (colA[pos + index] - meanA);
 427         denominatorB += (colB[pos + index] - meanB) * (colB[pos + index] - meanB);
 428     }
 429
 430     if (denominatorA == 0.0 || denominatorB == 0.0) {
 431         c = 0.0;
 432     } else {
 433         c = numerator / sqrt(denominatorA * denominatorB);
 434     }
 435
 436     return c;
 437 }
 438
 439 static PyObject*
 440 corr_func(PyObject *self, PyObject *args)
 441 {
 442     PyObject *PyaPWM, *PybPWM, *PycPWM, *PyMaxSlide;
 443     double *aPWM, *bPWM, *cPWM, *tempA, *tempB, *tempC;
 444     float fscore, rscore, bestScore;
 445     long maxSlide;
 446     int ok, index, indexMax, bIndexMax, ntIndex, motLen, padLen, slide, adjustedPadLen, adjustedSlide, tempMax, tempSize;
 447
 448     bestScore = 0.0;
 449
 450     ok = PyArg_UnpackTuple(args, "corr", 4, 4, &PyaPWM, &PybPWM, &PycPWM, &PyMaxSlide);
 451
 452     motLen = PyList_Size(PyaPWM);
 453     padLen = motLen - PyList_Size(PybPWM);
 454     maxSlide = PyInt_AsLong(PyMaxSlide);
 455
 456     if (maxSlide > motLen) {
 457         maxSlide = motLen - 1;
 458     }
 459
 460     indexMax = 4 * motLen;
 461     bIndexMax = 4 * (motLen - padLen);
 462     aPWM = malloc(indexMax * sizeof(double));
 463     bPWM = malloc(bIndexMax * sizeof(double));
 464     cPWM = malloc(bIndexMax * sizeof(double));
 465
 466     for (index = 0; index < (motLen - padLen); index++) {
 467         for (ntIndex = 0; ntIndex < 4; ntIndex++) {
 468             aPWM[4 * index + ntIndex] = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(PyaPWM, index), ntIndex));
 469             bPWM[4 * index + ntIndex] = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(PybPWM, index), ntIndex));
 470             cPWM[4 * index + ntIndex] = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(PycPWM, index), ntIndex));
 471         }
 472     }
 473
 474     for (; index < motLen; index++) {
 475         for (ntIndex = 0; ntIndex < 4; ntIndex++) {
 476             aPWM[4 * index + ntIndex] = PyFloat_AsDouble(PyList_GetItem(PyList_GetItem(PyaPWM, index), ntIndex));
 477         }
 478     }
 479
 480     for (slide = -1 * maxSlide; slide < (maxSlide + padLen + 1); slide++ ) {
 481         tempA = malloc(3 * indexMax * sizeof(double));
 482         tempB = malloc(3 * indexMax * sizeof(double));
 483         tempC = malloc(3 * indexMax * sizeof(double));
 484         if (slide < 0) {
 485             tempSize = abs(slide) + motLen;
 486             tempMax = 4 * tempSize;
 487             for (index = 0; index < 4 * abs(slide); index++) {
 488                 tempA[index] = 0.25;
 489             }
 490             for (; index < tempMax; index++) {
 491                 tempA[index] = aPWM[index - 4 * slide];
 492             }
 493             for (index = 0; index < bIndexMax; index++) {
 494                 tempB[index] = bPWM[index];
 495                 tempC[index] = cPWM[index];
 496             }
 497             for (; index < tempMax; index++) {
 498                 tempB[index] = 0.25;
 499                 tempC[index] = 0.25;
 500             }
 501         } else if ((slide > 0) && (slide <= maxSlide)) {
 502             if (padLen > 0) {
 503                 if (padLen >= slide) {
 504                     adjustedPadLen = padLen - slide;
 505                     adjustedSlide = 0;
 506                 } else {
 507                     adjustedPadLen = 0;
 508                     adjustedSlide = slide - padLen;
 509                 }
 510                 tempSize = motLen + adjustedSlide;
 511                 tempMax = indexMax + 4 * adjustedSlide;
 512                 for (index = 0; index < indexMax; index++) {
 513                     tempA[index] = aPWM[index];
 514                 }
 515                 for (; index < tempMax; index++) {
 516                     tempA[index] = 0.25;
 517                 }
 518                 tempMax = 4 * slide;
 519                 for (index = 0; index < tempMax; index++) {
 520                     tempB[index] = 0.25;
 521                     tempC[index] = 0.25;
 522                 }
 523                 tempMax += bIndexMax;
 524                 for (; index < tempMax; index++) {
 525                     tempB[index] = bPWM[index - 4 * slide];
 526                     tempC[index] = cPWM[index - 4 * slide];
 527                 }
 528                 tempMax = indexMax + 4 * adjustedSlide;
 529                 for (; index < tempMax; index++) {
 530                     tempB[index] = 0.25;
 531                     tempC[index] = 0.25;
 532                 }
 533             } else {
 534                 tempSize = motLen + slide;
 535                 tempMax = indexMax + 4 * slide;
 536                 for (index = 0; index < indexMax; index++) {
 537                     tempA[index] = aPWM[index];
 538                 }
 539                 for (; index < tempMax; index++) {
 540                     tempA[index] = 0.25;
 541                 }
 542                 tempMax = 4 * slide;
 543                 for (index = 0; index < tempMax; index++) {
 544                     tempB[index] = 0.25;
 545                     tempC[index] = 0.25;
 546                 }
 547                 tempMax += bIndexMax;
 548                 for (; index < tempMax; index++) {
 549                     tempB[index] = bPWM[index - 4 * slide];
 550                     tempC[index] = cPWM[index - 4 * slide];
 551                 }
 552             }
 553         } else if (slide > maxSlide) {
 554             tempSize = motLen + maxSlide;
 555             tempMax = indexMax + 4 * maxSlide;
 556             for (index = 0; index < indexMax; index++) {
 557                 tempA[index] = aPWM[index];
 558             }
 559             for (; index < tempMax; index++) {
 560                 tempA[index] = 0.25;
 561             }
 562             tempMax = 4 * slide;
 563             for (index = 0; index < tempMax; index++) {
 564                 tempB[index] = 0.25;
 565                 tempC[index] = 0.25;
 566             }
 567             tempMax += bIndexMax;
 568             for (; index < tempMax; index++) {
 569                 tempB[index] = bPWM[index - 4 * slide];
 570                 tempC[index] = cPWM[index - 4 * slide];
 571             }
 572             tempMax = indexMax + 4 * maxSlide;
 573             for (; index < tempMax; index++) {
 574                 tempB[index] = 0.25;
 575                 tempC[index] = 0.25;
 576             }
 577         } else {
 578             tempSize = motLen;
 579
 580             for (index = 0; index < indexMax; index++) {
 581                 tempA[index] = aPWM[index];
 582             }
 583             for (index = 0; index < bIndexMax; index++) {
 584                 tempB[index] = bPWM[index];
 585                 tempC[index] = cPWM[index];
 586             }
 587             for (; index < indexMax; index++) {
 588                 tempB[index] = 0.25;
 589                 tempC[index] = 0.25;
 590             }
 591         }
 592         fscore = 0.0;
 593         rscore = 0.0;
 594
 595         for (index = 0; index <tempSize; index++) {
 596             fscore += pearson_func(tempA, tempB, 4 * index);
 597             rscore += pearson_func(tempA, tempC, 4 * index);
 598             /* fprintf(stdout,"A: %f %f %f %f\n", tempA[4 * index], tempA[4 * index + 1], tempA[4 * index + 2], tempA[4 * index + 3]); */
 599             /* fprintf(stdout,"B: %f %f %f %f\n", tempB[4 * index], tempB[4 * index + 1], tempB[4 * index + 2], tempB[4 * index + 3]); */
 600         }
 601         fscore = fscore / tempSize;
 602         rscore = rscore / tempSize;
 603         /* fprintf(stdout,"slide %d %f %f\n", slide, fscore, rscore); */
 604         if ((fscore < rscore) && (rscore > bestScore)) {
 605             bestScore = rscore;
 606         }  else if (fscore > bestScore) {
 607             bestScore = fscore;
 608         }
 609
 610         free(tempA);
 611         free(tempB);
 612         free(tempC);
 613     }
 614
 615     free(aPWM);
 616     free(bPWM);
 617     free(cPWM);
 618
 619     return PyFloat_FromDouble(bestScore);
 620 }
 621
 622 static char the_func_doc[] =
 623 "returns a list of positions on aSeq that match the PWM within a Threshold, given as a percentage of the optimal consensus score.";
 624
 625 static char dpwm_func_doc[] =
 626 "returns a list of positions on aSeq that match the DPWM within a given Fold of the optimal consensus score.";
 627
 628 static char mer_func_doc[] =
 629 "returns a list of positions on aSeq that match an N-mer within M mismatches. Assumes Mers and Seq are in the same case.";
 630
 631 static char corr_func_doc[] =
 632 "returns a pearson-correlation coefficient based similarity value between -1 (anti-correlated) and +1 (identical) for two motifs.";
 633
 634 static PyMethodDef module_methods[] = {
 635     {"locateMotif", the_func, METH_VARARGS, the_func_doc},
 636     {"locateMarkov1", dpwm_func, METH_VARARGS, dpwm_func_doc},
 637     {"locateMer", mer_func, METH_VARARGS, mer_func_doc},
 638     {"correlateMotifs", corr_func, METH_VARARGS, corr_func_doc},
 639     {NULL, NULL}
 640 };
 641
 642 PyMODINIT_FUNC
 643 init_motif(void)
 644 {
 645     Py_InitModule3("_motif", module_methods, module_doc);
 646 }