4a9f5caa2ea145d8ef8d86b6d4abb9a5d70bc345
hiram
  Tue Jul 11 14:54:48 2017 -0700
adding score calculation and second PAM sequence refs #18969

diff --git src/utils/crisprKmers/crisprKmers.c src/utils/crisprKmers/crisprKmers.c
index d4ea4f0..4947dc4 100644
--- src/utils/crisprKmers/crisprKmers.c
+++ src/utils/crisprKmers/crisprKmers.c
@@ -47,39 +47,41 @@
 static char *ranges = NULL;	/* use ranges <file> to limit scanning */
 static struct hash *rangesHash = NULL;	/* ranges into hash + binkeeper */
 static char *dumpKmers = NULL;	/* file name to write out kmers from scan */
 static char *loadKmers = NULL;	/* file name to read in kmers from previous scan */
 
 /* Command line validation table. */
 static struct optionSpec options[] = {
    {"bed", OPTION_STRING},
    {"offTargets", OPTION_STRING},
    {"ranges", OPTION_STRING},
    {"dumpKmers", OPTION_STRING},
    {"loadKmers", OPTION_STRING},
    {NULL, 0},
 };
 
+// sizeof(struct crispr): 32
 struct crispr
 /* one chromosome set of crisprs */
     {
     struct crispr *next;		/* Next in list. */
     long long sequence;	/* sequence value in 2bit format */
     long long start;		/* chromosome start 0-relative */
     char strand;			/* strand: + or - */
     };
 
+// sizeof(struct crisprList): 72
 struct crisprList
 /* all chromosome sets of crisprs */
     {
     struct crisprList *next;	/* Next in list. */
     char *chrom;		/* chrom name */
     int size;			/* chrom size */
     struct crispr *chromCrisprs;	/* all the crisprs on this chrom */
     long long crisprCount;	/* number of crisprs on this chrom */
     long long *sequence;	/* array of the sequences */
     long long *start;		/* array of the starts */
     char *strand;		/* array of the strand characters */
     int **offBy;		/* offBy[5][n] */
     };
 
 #define A_BASE	0
@@ -110,37 +112,41 @@
 int i;
 for (i=0; i<ArraySize(orderedNtVal); i++)
     orderedNtVal[i] = -1;
 orderedNtVal['A'] = orderedNtVal['a'] = A_BASE;
 orderedNtVal['C'] = orderedNtVal['c'] = C_BASE;
 orderedNtVal['G'] = orderedNtVal['g'] = G_BASE;
 orderedNtVal['T'] = orderedNtVal['t'] = T_BASE;
 orderedNtVal['U'] = orderedNtVal['u'] = T_BASE;
 bases[A_BASE] = 'A';
 bases[C_BASE] = 'C';
 bases[G_BASE] = 'G';
 bases[T_BASE] = 'T';
 }	//	static void initOrderedNtVal()
 
 static void timingMessage(char *prefix, long long count, char *message,
-    long ms, char *units)
+    long ms, char *units, char *invUnits)
 {
 double perSecond = 0.0;
-if (ms > 0)
+double inverse = 0.0;
+if ((ms > 0) && (count > 0))
+    {
     perSecond = 1000.0 * count / ms;
+    inverse = 1.0 / perSecond;
+    }
 
-verbose(1, "# %s: %lld %s @ %ld ms -> %.2f %s\n", prefix, count, message, ms, perSecond, units);
+verbose(1, "# %s: %lld %s @ %ld ms -> %.2f %s == %g %s\n", prefix, count, message, ms, perSecond, units, inverse, invUnits);
 }
 
 static struct hash *readRanges(char *bedFile)
 /* Read bed and return it as a hash keyed by chromName
  * with binKeeper values.  (from: src/hg/bedIntersect/bedIntersec.c) */
 {
 struct hash *hash = newHash(0);	/* key is chromName, value is binkeeper data */
 struct lineFile *lf = lineFileOpen(bedFile, TRUE);
 char *row[3];
 
 while (lineFileNextRow(lf, row, 3))
     {
     struct binKeeper *bk;
     struct bed3 *item;
     struct hashEl *hel = hashLookup(hash, row[0]);
@@ -260,75 +266,79 @@
     int r = 0;
     for (r = 0; r < 5; ++r)
         cl->offBy[r] = (int *)needLargeZeroedMem(memSize);
 
     long long i = 0;
     struct crispr *c = NULL;
     for (c = cl->chromCrisprs; c; c = c->next)
         {
 	++itemsCopied;
         cl->sequence[i] = c->sequence;
         cl->start[i] = c->start;
         cl->strand[i++] = c->strand;
         }
     }
 long elapsedMs = clock1000() - startTime;
-timingMessage("copyToArray:", itemsCopied, "items copied", elapsedMs, "items/sec");
+timingMessage("copyToArray:", itemsCopied, "items copied", elapsedMs, "items/sec", "seconds/item");
 }	//	static void copyToArray(struct crisprList *list)	*/
 
 static struct crisprList *generateKmers(struct dnaSeq *seq)
 {
 struct crispr *crisprSet = NULL;
 struct crisprList *returnList = NULL;
 AllocVar(returnList);
 returnList->chrom = cloneString(seq->name);
 returnList->size = seq->size;
 int i;
 DNA *dna;
 long long chromPosition = 0;
 long long startGap = 0;
 long long gapCount = 0;
 int kmerLength = 0;
 long long kmerVal = 0;
+long long endsAG = (A_BASE << 2) | G_BASE;
 long long endsGG = (G_BASE << 2) | G_BASE;
+long long beginsCT = (long long)((C_BASE << 2) | T_BASE) << 42;
 long long beginsCC = (long long)((C_BASE << 2) | C_BASE) << 42;
 long long reverseMask = (long long)0xf << 42;
+verbose(4, "#   endsAG: %032llx\n", endsAG);
 verbose(4, "#   endsGG: %032llx\n", endsGG);
+verbose(4, "# beginsCT: %032llx\n", beginsCT);
 verbose(4, "# beginsCC: %032llx\n", beginsCC);
 verbose(4, "#  46 bits: %032llx\n", (long long) fortySixBits);
 
 dna=seq->dna;
 for (i=0; i < seq->size; ++i)
     {
     int val = orderedNtVal[(int)dna[i]];
     if (val >= 0)
 	{
         kmerVal = fortySixBits & ((kmerVal << 2) | val);
         ++kmerLength;
 	if (kmerLength > 22)
 	    {
-	    if (endsGG == (kmerVal & 0xf))	/* check positive strand */
-                {
+	    if ((endsAG == (kmerVal & 0xf)) || (endsGG == (kmerVal & 0xf)))
+                {	/* have match for positive strand */
                 struct crispr *oneCrispr = NULL;
                 AllocVar(oneCrispr);
                 oneCrispr->start = chromPosition - 22;
                 oneCrispr->strand = '+';
                 oneCrispr->sequence = kmerVal;
                 slAddHead(&crisprSet, oneCrispr);
                 }
-	    if (beginsCC == (kmerVal & reverseMask))	/* check neg strand */
-                {
+	    if ((beginsCT == (kmerVal & reverseMask)) || (beginsCC == (kmerVal & reverseMask)))
+                {	/* have match for negative strand */
                 struct crispr *oneCrispr = NULL;
                 AllocVar(oneCrispr);
                 oneCrispr->start = chromPosition - 22;
                 oneCrispr->strand = '-';
                 oneCrispr->sequence = revComp(kmerVal);
                 slAddHead(&crisprSet, oneCrispr);
                 }
 	    }
         }	// if (val >= 0)
         else
             {
             ++gapCount;
             startGap = chromPosition;
             /* skip all N's == any value not = 0,1,2,3 */
             while ( ((i+1) < seq->size) && (orderedNtVal[(int)dna[i+1]] < 0))
@@ -366,120 +376,111 @@
     return notUnique;
 else
     {
     int offSum = offBy[1][index] + offBy[2][index] + offBy[3][index] +
 		offBy[4][index];
     if (offSum < 100)
        return highScore;
     else if (offSum < 150)
        return mediumScore;
     else if (offSum < 250)
        return lowScore;
     else
        return twoMany;
     }
 }
-// sample of offBy sums from the chrM vs. hg38 run:
-// Q1 91.000000
-// median 150.000000
-// Q3 251.000000
-// average 238.378477
-// min 4.000000
-// max 9023.000000
-// count 2193
-// total 522764.000000
-// standard deviation 424.309155
 
 static void countsOutput(struct crisprList *all)
 /* everything has been scanned and counted, print out all the data from arrays*/
 {
 long startTime = clock1000();
 long long itemsOutput = 0;
 
 FILE *bedFH = NULL;
 if (bedFileOut)
     bedFH = mustOpen(bedFileOut, "w");
 
 struct crisprList *list;
 long long totalOut = 0;
 for (list = all; list; list = list->next)
     {
     long long c = 0;
     for (c = 0; c < list->crisprCount; ++c)
         {
 	++itemsOutput;
 	int negativeOffset = 0;
         if (list->strand[c] == '-')
 	    negativeOffset = 3;
 	long long txStart = list->start[c] + negativeOffset;
 	long long txEnd = txStart + 20;
+
         int totalOffs = list->offBy[0][c] + list->offBy[1][c] +
 	    list->offBy[2][c] + list->offBy[3][c] + list->offBy[4][c];
 
         char *color = itemColor(list->offBy, c);
 
         if (0 == totalOffs)
 verbose(1, "# PERFECT score %s:%lld %c\t%s\n", list->chrom, list->start[c], list->strand[c], kmerValToString(list->sequence[c], 3));
 
 	if (bedFH)
 	    fprintf(bedFH, "%s\t%lld\t%lld\t%d,%d,%d,%d,%d\t%d\t%c\t%lld\t%lld\t%s\t%s\t%s\t%s\t%s\t%d\t%d\t%d\t%d\n", list->chrom, list->start[c], list->start[c]+23, list->offBy[0][c], list->offBy[1][c], list->offBy[2][c], list->offBy[3][c], list->offBy[4][c],  list->offBy[0][c], list->strand[c], txStart, txEnd, color, color, color, kmerValToString(list->sequence[c], 3), kmerPAMString(list->sequence[c]), list->offBy[1][c], list->offBy[2][c], list->offBy[3][c], list->offBy[4][c]);
 
         if (list->offBy[0][c])
            verbose(3, "# array identical: %d %s:%lld %c\t%s\n", list->offBy[0][c], list->chrom, list->start[c], list->strand[c], kmerValToString(list->sequence[c], 3));
 	++totalOut;
 	}
     }
 if (bedFH)
     carefulClose(&bedFH);
 
 long elapsedMs = clock1000() - startTime;
-timingMessage("copyToArray:", itemsOutput, "items output", elapsedMs, "items/sec");
+timingMessage("copyToArray:", itemsOutput, "items output", elapsedMs, "items/sec", "seconds/item");
 }	// static void countsOutput(struct crisprList *all)
 
 static struct crisprList *scanSequence(char *inFile)
 /* scan the given file, return list of crisprs */
 {
 verbose(1, "#\tscanning sequence file: %s\n", inFile);
 dnaUtilOpen();
 struct dnaLoad *dl = dnaLoadOpen(inFile);
 struct dnaSeq *seq;
 struct crisprList *listReturn = NULL;
 
 long elapsedMs = 0;
 long scanStart = clock1000();
 long long totalCrisprs = 0;
 
 /* scanning all sequences, setting up crisprs on the listReturn */
 while ((seq = dnaLoadNext(dl)) != NULL)
     {
     if (startsWithNoCase("chrUn", seq->name) ||
          rStringIn("hap", seq->name) || rStringIn("_alt", seq->name) )
 	{
 	verbose(1, "# skip chrom: %s\n", seq->name);
 	continue;
 	}
 
     long startTime = clock1000();
     struct crisprList *oneList = generateKmers(seq);
     slAddHead(&listReturn, oneList);
     totalCrisprs += oneList->crisprCount;
     elapsedMs = clock1000() - startTime;
-    timingMessage(seq->name, oneList->crisprCount, "crisprs", elapsedMs, "crisprs/sec");
+    timingMessage(seq->name, oneList->crisprCount, "crisprs", elapsedMs, "crisprs/sec", "seconds/crispr");
     }
 
 elapsedMs = clock1000() - scanStart;
-timingMessage("scanSequence", totalCrisprs, "total crisprs", elapsedMs, "crisprs/sec");
+timingMessage("scanSequence", totalCrisprs, "total crisprs", elapsedMs, "crisprs/sec", "seconds/crispr");
 return listReturn;
 }	/*	static crisprList *scanSequence(char *inFile)	*/
 
 static struct crisprList *rangeExtraction(struct crisprList **allReference)
 /* given ranges in global rangesHash, construct new list of crisprs that
  *     have any type of overlap with the ranges, also extract those items from
  *         the all list.  Returns new list.
  */
 {
 struct crisprList *all = *allReference;
 struct crisprList *listReturn = NULL;
 struct crisprList *list = NULL;
 int inputChromCount = slCount(all);
 long long returnListCrisprCount = 0;
 long long examinedCrisprCount = 0;
@@ -544,60 +545,117 @@
 	else
 	    {
 	    list->crisprCount = slCount(list->chromCrisprs);
             long long removedCrisprs = beforeCrisprCount - list->crisprCount;
 	    verbose(1, "# range scan chrom %s had %lld crisprs, removed %lld leaving %lld target\n",
 		list->chrom, beforeCrisprCount, removedCrisprs, list->crisprCount);
 	    }
 	}
     prevChromList = list;
     }	//	for (list = all; list; list = list->next)
 
 elapsedMs = clock1000() - scanStart;
 verbose(1, "# range scanning %d chroms, return %d selected chroms, leaving %d chroms\n",
     inputChromCount, slCount(listReturn), slCount(all));
 long long targetCrisprCount = examinedCrisprCount - returnListCrisprCount;
-timingMessage("range scan", examinedCrisprCount, "examined crisprs", elapsedMs, "crisprs/sec");
-timingMessage("range scan", targetCrisprCount, "remaining target crisprs", elapsedMs, "crisprs/sec");
-timingMessage("range scan", returnListCrisprCount, "returned query crisprs", elapsedMs, "crisprs/sec");
+timingMessage("range scan", examinedCrisprCount, "examined crisprs", elapsedMs, "crisprs/sec", "seconds/crispr");
+timingMessage("range scan", targetCrisprCount, "remaining target crisprs", elapsedMs, "crisprs/sec", "seconds/crispr");
+timingMessage("range scan", returnListCrisprCount, "returned query crisprs", elapsedMs, "crisprs/sec", "seconds/crispr");
 if (NULL == all)
     {
     allReference = NULL;	// they have all been removed
     verbose(1, "# range scan, every single chrom has been removed\n");
     }
 else if (*allReference != all)
     {
     verbose(1, "# range scan, first chrom list has been moved from %p to %p\n", (void *)*allReference, (void *)all);
     *allReference = all;
     }
 return listReturn;
 }	//	static crisprList *rangeExtraction(crisprList *all)
 
+static int hitScoreM[20] = {0,0,0.014,0,0,0.395,0.317,0,0.389,0.079,0.445,0.508,0.613,0.851,0.732,0.828,0.615,0.804,0.685,0.583};
+
+static int calcHitScore(long long sequence1, long long sequence2)
+{
+double score1 = 1.0;
+int mmCount = 0;
+int lastMmPos = -1;
+/* the XOR will determine differences in two sequences, the shift
+ # right 6 removes the PAM sequence */
+long long misMatch = (sequence1 ^ sequence2) >> 6;
+int distCount = 0;
+int distSum = 0;
+int pos = 0;
+for (pos = 0; pos < 20; ++pos)
+    {
+    int diff = misMatch & 0x3;
+    if (diff)
+	{
+	++mmCount;
+        if (lastMmPos != -1)
+	    {
+	    distSum += pos-lastMmPos;
+	    ++distCount;
+	    }
+	score1 *= 1 - hitScoreM[pos];
+	lastMmPos = pos;
+	}
+    misMatch >>= 2;
+    }
+
+double score2 = 1.0;
+if (distCount > 1)
+    {
+    double avgDist = (double)distSum / distCount;
+    score2 = 1.0 / (((19-avgDist)/19.0) * 4 + 1);
+    }
+double score3 = 1.0;
+if (mmCount > 0)
+    score3 = 1.0 / (mmCount * mmCount);
+
+return (int)(score1 * score2 * score3 * 100);
+}
+
 static void recordOffTargets(struct crisprList *query,
     struct crisprList *target, int bitsOn, long long qIndex,
 	long long tIndex)
-/* bitsOn is from 1 to 4, record this match when less than 10 count */
+/* bitsOn is from 1 to 4, record this match when less than 1000 total */
 {
+if (query->offBy[0][qIndex] ) // no need to accumulate if 0-mismatch > 0
+    return;
+
 if (offFile)
     {
-    if (query->offBy[bitsOn][qIndex] < 10)
-	{/* needs to be two prints as the kmerValToString return is confused */
-	fprintf(offFile, "%s:%lld %c %s\t", query->chrom, query->start[qIndex],
-	    query->strand[qIndex], kmerValToString(query->sequence[qIndex], 3));
-        fprintf(offFile, "%s:%lld %c %s\t%d\n", target->chrom,
+    int i = 0;
+    int bitsOnSum = 0;
+    for (i = 1; i < 5; ++i)
+        bitsOnSum += query->offBy[i][qIndex];
+
+    if (bitsOnSum < 1000)	// could be command line option limit
+        { /* needs to be two fprintfs, as the kmer*String() function
+           * returns are confused due to static returns */
+        int hitScore =
+		calcHitScore(query->sequence[qIndex], target->sequence[tIndex]);
+        fprintf(offFile, "%s:%lld %c %s %s %d\t", query->chrom,
+            query->start[qIndex], query->strand[qIndex],
+                kmerValToString(query->sequence[qIndex], 3),
+		    kmerPAMString(query->sequence[qIndex]), hitScore);
+        fprintf(offFile, "%s:%lld %c %s %s\t%d\n", target->chrom,
             target->start[tIndex], target->strand[tIndex],
-		kmerValToString(target->sequence[tIndex], 3), bitsOn);
+                kmerValToString(target->sequence[tIndex], 3),
+		    kmerPAMString(target->sequence[tIndex]), bitsOn);
         }
     }
 }	//	static void recordOffTargets(struct crisprList *query,
 	//	    struct crisprList *target, int bitsOn, long long qIndex,
 	//		long long tIndex)
 
 static void queryVsAll(struct crisprList *query, struct crisprList *target)
 /* run the query crisprs list against the target list in the array structures */
 {
 struct crisprList *qList = NULL;
 long long totalCrisprsQuery = 0;
 long long totalCrisprsTarget = 0;
 long long totalCompares = 0;
 
 long processStart = clock1000();
@@ -626,47 +684,47 @@
                  *  the shift right 6 removes the PAM sequence
                  */
                 long long misMatch =
                     (qList->sequence[qCount] ^ tList->sequence[tCount]) >> 6;
                 if (misMatch)
                     {
                     /* possible misMatch bit values: 01 10 11
                      *  turn those three values into just: 01
                      */
                     misMatch = (misMatch | (misMatch >> 1)) & 0x5555555555;
                     int bitsOn = _mm_popcnt_u64(misMatch);
 		    if (bitsOn < 5)
 			{
 			recordOffTargets(qList, tList, bitsOn, qCount, tCount);
 			qList->offBy[bitsOn][qCount] += 1;
-			tList->offBy[bitsOn][tCount] += 1;
+//			tList->offBy[bitsOn][tCount] += 1; not needed
 			}
                     }
                 else
                     { 	/* no misMatch, identical crisprs */
                     qList->offBy[0][qCount] += 1;
-                    tList->offBy[0][tCount] += 1;
+//                  tList->offBy[0][tCount] += 1;	not needed
                     }
                 } // for (tCount = 0; tCount < tList->crisprCount; ++tCount)
             }	//	for (tList = target; tList; tList = tList->next)
 	}	//	for (qCount = 0; qCount < qList->crisprCount; ++qCount)
     }	//	for (qList = query; qList; qList = qList->next)
 verbose(1, "# done with scanning, check timing\n");
 elapsedMs = clock1000() - processStart;
-timingMessage("queryVsAll", totalCrisprsQuery, "query crisprs processed", elapsedMs, "crisprs/sec");
-timingMessage("queryVsAll", totalCrisprsTarget, "vs target crisprs", elapsedMs, "crisprs/sec");
-timingMessage("queryVsAll", totalCompares, "total comparisons", elapsedMs, "compares/sec");
+timingMessage("queryVsAll", totalCrisprsQuery, "query crisprs processed", elapsedMs, "crisprs/sec", "seconds/crispr");
+timingMessage("queryVsAll", totalCrisprsTarget, "vs target crisprs", elapsedMs, "crisprs/sec", "seconds/crispr");
+timingMessage("queryVsAll", totalCompares, "total comparisons", elapsedMs, "compares/sec", "seconds/compare");
 }	/* static struct crisprList *queryVsAll(struct crisprList *query,
 	    struct crisprList *target) */
 
 static void allVsAll(struct crisprList *all)
 /* run this 'all' list vs. itself avoiding self to self comparisons */
 {
 struct crisprList *qList = NULL;
 long long totalCrisprsQuery = 0;
 long long totalCrisprsCompare = 0;
 
 long processStart = clock1000();
 long elapsedMs = 0;
 
 /* query runs through all chroms */
 for (qList = all; qList; qList = qList->next)
@@ -705,32 +763,32 @@
 			qList->offBy[bitsOn][qCount] += 1;
 			tList->offBy[bitsOn][tCount] += 1;
 			}
                     }
                 else
                     { 	/* no misMatch, identical crisprs */
                     qList->offBy[0][qCount] += 1;
                     tList->offBy[0][tCount] += 1;
                     }
                 } // for (tCount = 0; tCount < tList->crisprCount; ++tCount)
                 tStart = 0;	/* following chroms run through all */
             }	//	for (tList = target; tList; tList = tList->next)
 	}	//	for (qCount = 0; qCount < qList->crisprCount; ++qCount)
     }	//	for (qList = query; qList; qList = qList->next)
 elapsedMs = clock1000() - processStart;
-timingMessage("allVsAll", totalCrisprsQuery, "crisprs processed", elapsedMs, "crisprs/sec");
-timingMessage("allVsAll", totalCrisprsCompare, "total comparisons", elapsedMs, "compares/sec");
+timingMessage("allVsAll", totalCrisprsQuery, "crisprs processed", elapsedMs, "crisprs/sec", "seconds/crispr");
+timingMessage("allVsAll", totalCrisprsCompare, "total comparisons", elapsedMs, "compares/sec", "seconds/compare");
 }	/* static struct crisprList *allVsAll(struct crisprList *query,
 	    struct crisprList *target) */
 
 static struct crisprList *readKmers(char *fileIn)
 /* read in kmer list from 'fileIn', return list structure */
 {
 verbose(1, "# reading crisprs from: %s\n", fileIn);
 struct crisprList *listReturn = NULL;
 struct lineFile *lf = lineFileOpen(fileIn, TRUE);
 char *row[10];
 int wordCount = 0;
 long long crisprsInput = 0;
 
 long startMs = clock1000();
 
@@ -757,61 +815,61 @@
         AllocVar(oneCrispr);
         oneCrispr->sequence = sqlLongLong(row[0]);
         oneCrispr->start = sqlLongLong(row[1]);
         oneCrispr->strand = row[2][0];
         slAddHead(&newItem->chromCrisprs, oneCrispr);
 	}
     if (verifyCount != newItem->crisprCount)
 	errAbort("expecting %lld kmer items at line %d in file %s, found: %lld",
 	    newItem->crisprCount, lf->lineIx, fileIn, verifyCount);
     crisprsInput += verifyCount;
     }
 
 lineFileClose(&lf);
 
 long elapsedMs = clock1000() - startMs;
-timingMessage("readKmers", crisprsInput, "crisprs read in", elapsedMs, "crisprs/sec");
+timingMessage("readKmers", crisprsInput, "crisprs read in", elapsedMs, "crisprs/sec", "seconds/crispr");
 
 return listReturn;
 }	//	static struct crisprList *readKmers(char *fileIn)
 
 static void writeKmers(struct crisprList *all, char *fileOut)
 /* write kmer list 'all' to 'fileOut' */
 {
 FILE *fh = mustOpen(fileOut, "w");
 struct crisprList *list = NULL;
 long long crisprsWritten = 0;
 
 long startMs = clock1000();
 
 slReverse(&all);
 for (list = all; list; list = list->next)
     {
     fprintf(fh, "%s\t%lld\t%d\n", list->chrom, list->crisprCount, list->size);
     struct crispr *c = NULL;
     slReverse(&list->chromCrisprs);
     for (c = list->chromCrisprs; c; c = c->next)
 	{
 	fprintf(fh, "%lld\t%lld\t%c\n", c->sequence,
 	    c->start, c->strand);
 	++crisprsWritten;
 	}
     }
 carefulClose(&fh);
 
 long elapsedMs = clock1000() - startMs;
-timingMessage("writeKmers", crisprsWritten, "crisprs written", elapsedMs, "crisprs/sec");
+timingMessage("writeKmers", crisprsWritten, "crisprs written", elapsedMs, "crisprs/sec", "seconds/crispr");
 }	//	static void writeKmers(struct crisprList *all, char *fileOut)
 
 static void crisprKmers(char *sequence)
 /* crisprKmers - find and annotate crispr sequences. */
 {
 struct crisprList *queryCrisprs = NULL;
 // struct crisprList *countedCrisprs = NULL;
 struct crisprList *allCrisprs = NULL;
 
 if (loadKmers)
     allCrisprs = readKmers(loadKmers);
 else
     allCrisprs = scanSequence(sequence);
 
 if (dumpKmers)