40ed3e546ef868bffe4c338d93d6528138ecfc44
angie
  Wed Sep 30 09:25:12 2015 -0700
Add test cases (and bug fixes!) to make sure that we get the desired
behavior when an insertion is adjacent to a non-insertion, and when
an insertion falls at the start or end of the search region (if applicable):

- Include insertions that fall at the start or end of the search region

- If the primary row is an insertion, keep secondary non-insertion rows
to the left and right.

- If the primary row is a non-insertion, keep secondary insertion rows
at its start and end.

(i.e. keep insertions at boundaries -- but don't let any non-insertions
slip through)

The VCF logic is more complicated because VCF indels always include an
extra base to the left, so they appear to start before they actually do,
and can be interspersed with non-indels that start there.

diff --git src/lib/annoGrator.c src/lib/annoGrator.c
index c395b0d..413e484 100644
--- src/lib/annoGrator.c
+++ src/lib/annoGrator.c
@@ -64,33 +64,36 @@
 	errAbort("annoGrator %s: Unsorted input from internal source %s (%s < %s)",
 		 self->streamer.name, self->mySource->name, newRow->chrom, self->qTail->chrom);
     else if (cDifNewTail == 0 && newRow->start < self->qTail->start)
 	errAbort("annoGrator %s: Unsorted input from internal source %s (%s, %u < %u)",
 		 self->streamer.name, self->mySource->name,
 		 newRow->chrom, newRow->start, self->qTail->start);
     }
 }
 
 INLINE void agFetchToEnd(struct annoGrator *self, char *chrom, uint start, uint end)
 /* Fetch rows until we are sure we have all items that start to the left of end,
  * i.e. we have an item that starts at/after end or we hit eof. */
 {
 while (!self->eof &&
        (self->qTail == NULL || strcmp(self->qTail->chrom, chrom) < 0 ||
-	(sameString(self->qTail->chrom, chrom) && self->qTail->start < end)))
+	(sameString(self->qTail->chrom, chrom) && self->qTail->start <= end)))
     {
-    struct annoRow *newRow = self->mySource->nextRow(self->mySource, chrom, start, self->qLm);
+    // If primary row is a zero-length insertion, we will want to pick up items adjacent
+    // to it on the left, so move minEnd one base to the left:
+    uint minEnd = (start == end) ? start - 1 : start;
+    struct annoRow *newRow = self->mySource->nextRow(self->mySource, chrom, minEnd, self->qLm);
     if (newRow == NULL)
 	self->eof = TRUE;
     else
 	{
 	agCheckInternalSorting(self, newRow);
 	int cDifNewP = strcmp(newRow->chrom, chrom);
 	if (cDifNewP >= 0)
 	    {
 	    // Add newRow to qTail
 	    if (self->qTail == NULL)
 		{
 		if (self->qHead != NULL)
 		    errAbort("annoGrator %s: qTail is NULL but qHead is non-NULL",
 			     self->streamer.name);
 		self->qHead = self->qTail = newRow;
@@ -107,68 +110,91 @@
 	// If we're skipping past large regions, keep qLm size under control:
 	else
 	    {
 	    self->qSkippedCount++;
 	    if (self->qSkippedCount > 1024 && self->qHead == NULL && self->qTail == NULL)
 		{
 		lmCleanup(&(self->qLm));
 		self->qLm = lmInit(0);
 		self->qSkippedCount = 0;
 		}
 	    }
 	}
     }
 }
 
+static boolean handleRJFail(boolean rjFail, boolean *retRJFilterFailed)
+/* If retRJFilterFailed is non-NULL and rjFail is TRUE then set *retRJFilterFailed to TRUE
+ * and return TRUE.  Otherwise do nothing and return FALSE. */
+{
+if (retRJFilterFailed != NULL && rjFail)
+    {
+    *retRJFilterFailed = TRUE;
+    return TRUE;
+    }
+return FALSE;
+}
+
 struct annoRow *annoGratorIntegrate(struct annoGrator *self, struct annoStreamRows *primaryData,
 				    boolean *retRJFilterFailed, struct lm *callerLm)
 /* Given a single row from the primary source, get all overlapping rows from internal
  * source, and produce joined output rows.
  * If retRJFilterFailed is non-NULL:
  * - any overlapping row has a rightJoin filter failure (see annoFilter.h), or
  * - overlap rule is agoMustOverlap and no rows overlap, or
  * - overlap rule is agoMustNotOverlap and any overlapping row is found,
  * then set retRJFilterFailed and stop. */
 {
 struct annoRow *primaryRow = primaryData->rowList;
 struct annoRow *rowList = NULL;
 agCheckPrimarySorting(self, primaryRow);
-// In order to catch the intersection of two 0-length elements (i.e. two insertions),
-// we have to broaden our search a little:
 int pStart = primaryRow->start, pEnd = primaryRow->end;
-if (pStart == pEnd)
-    {
-    pStart--;
-    pEnd++;
-    }
 char *pChrom = primaryRow->chrom;
 agTrimToStart(self, pChrom, pStart);
 agFetchToEnd(self, pChrom, pStart, pEnd);
 boolean rjFailHard = (retRJFilterFailed != NULL);
 if (rjFailHard)
     *retRJFilterFailed = FALSE;
+int numCols = self->mySource->numCols;
+enum annoRowType rowType = self->mySource->rowType;
 struct annoRow *qRow;
+if (pStart == pEnd)
+    {
+    // If the primary feature is an insertion, catch any feature that is adjacent to it,
+    // not only those that overlap.
     for (qRow = self->qHead;  qRow != NULL;  qRow = qRow->next)
         {
-    if (qRow->start < pEnd && qRow->end > pStart && sameString(qRow->chrom, pChrom))
+        if (qRow->start <= pEnd && qRow->end >= pStart)
             {
-	int numCols = self->mySource->numCols;
-	enum annoRowType rowType = self->mySource->rowType;
             slAddHead(&rowList, annoRowClone(qRow, rowType, numCols, callerLm));
-	if (rjFailHard && qRow->rightJoinFail)
+            if (handleRJFail(qRow->rightJoinFail, retRJFilterFailed))
+                break;
+            }
+        }
+    }
+else
     {
-	    *retRJFilterFailed = TRUE;
+    for (qRow = self->qHead;  qRow != NULL;  qRow = qRow->next)
+        {
+        if (((qRow->start < pEnd && qRow->end > pStart) ||
+             // Make sure to include q insertions at pStart or pEnd:
+             (qRow->start == qRow->end &&
+              (qRow->start == pEnd || qRow->end == pStart))) &&
+            sameString(qRow->chrom, pChrom))
+            {
+            slAddHead(&rowList, annoRowClone(qRow, rowType, numCols, callerLm));
+            if (handleRJFail(qRow->rightJoinFail, retRJFilterFailed))
                 break;
             }
         }
     }
 slReverse(&rowList);
 // If no rows overlapped primary, and there is a right-join, !isExclude (i.e. isInclude) filter,
 // then we need to set retRJFilterFailed because the condition was not met to include
 // the primary item.
 if (retRJFilterFailed &&
     ((rowList == NULL && (self->haveRJIncludeFilter || self->overlapRule == agoMustOverlap)) ||
      (rowList != NULL && self->overlapRule == agoMustNotOverlap)))
     *retRJFilterFailed = TRUE;
 return rowList;
 }