text: /* makeGraph - this module creates a txGraph from a list of 
   * linkedBeds.  The main steps to this process are:
   *   1) Create list of all unique start and end points in
   *      beds.  These are 'hard' if they are at a real splice
   *      site and 'soft' otherwise.
   *   2) Create list of all unique edges. Edges include evidence
   *      if any.
   *   3) Snap soft ends to nearby hard ends of compatible type.
   *   4) For half-soft edges, examine all other edges that share
   *      the hard end, and and that are hard on both ends.  Merge
   *      the half-soft edge with the most similar such double-hard
   *      edge, favoring extensions of edge over truncations. 
   *   5) For edges that are half soft that can't yet be snapped, 
   *      merge all soft ends to a single point that is a consensus favoring 
   *      larger transcripts.
   *    6) Merge together all overlapping edges that are soft on both sides.  Make
   *     the edge boundaries equal to the median value of all the merged boundaries. 
   * Through most of the algorithm the graph is represented as two binary trees, one
   * for vertices and one for edges. At the end this is converted to the txGraph
   * structure. */
  
  #include "common.h"
  #include "hash.h"
  #include "localmem.h"
  #include "rbTree.h"
  #include "rangeTree.h"
  #include "dlist.h"
  #include "dnautil.h"
  #include "bed.h"
  #include "ggTypes.h"
  #include "nibTwo.h"
  #include "txGraph.h"
  #include "txBedToGraph.h"
  
  struct evidence
  /* One piece of evidence. */
      {
      struct evidence *next;
      struct linkedBeds *lb;	/* Associated input data. */
      int start,end;		/* Bounds of evidence. */
      };
  
  struct vertex
  /* A point in our graph */
      {
      int position;		/* Position in genome. */
      enum ggVertexType type;	/* hard/soft? start/end? */
      struct vertex *movedTo;	/* Forwarding address (temporary) */
      struct slRef *waysIn;	/* References to preceeding edges (temporary). */
      struct slRef *waysOut;	/* References to following edges (temporary). */
      int count;			/* Number of times used, or index of vertex (temporary). */
      };
  
  struct edge
  /* An edge in our graph */
      {
      struct edge *next;          /* for forming a singly-linked list of edges */
      struct vertex *start;	/* Starting vertex. */
      struct vertex *end;		/* Ending vertex. */
      struct evidence *evList;	/* List of evidence. */
      };
  
  
  static int vertexCmp(void *va, void *vb)
  /* Return: if a before b, return negative.  If a after b, return positive.  
   * else if a and b have the same position and same type, return 0.
   * else if a and b have same position but different types, return 
   * a nonzero value reflecting the difference in their types. */
  {
  struct vertex *a = va;
  struct vertex *b = vb;
  int diff = a->position - b->position;
  if (diff == 0)
      diff = a->type - b->type;
  return diff;
  }
  
  static int edgeCmp(void *va, void *vb)
  /* Return -1 if a before b,  0 if a and b overlap,
   * and 1 if a after b. */
  {
  struct edge *a = va;
  struct edge *b = vb;
  int diff = vertexCmp(a->start, b->start);
  if (diff == 0)
      diff = vertexCmp(a->end, b->end);
  return diff;
  }
  
  static boolean encloses(struct edge *a, struct edge *b)
  /* Return TRUE if the first edge encloses the second, false otherwise */
  {
  if (a->start->position <= b->start->position
      && a->end->position >= b->end->position)
      {
      return(TRUE);
      } 
  else 
      {
      return(FALSE);
      }
  }
  
  static boolean trustedEdge(struct edge *edge)
  /* Return TRUE if any member of edge evidence is from a trusted source. */
  {
  struct evidence *ev;
  for (ev = edge->evList; ev != NULL; ev = ev->next)
      {
      if (trustedSource(ev->lb->sourceType))
          return TRUE;
      }
  return FALSE;
  }
  
  
  static struct vertex *matchingVertex(struct rbTree *tree, int position, enum ggVertexType type)
  /* Find matching vertex.  Return NULL if none. */
  {
  struct vertex temp;
  temp.position = position;
  temp.type = type;
  return rbTreeFind(tree, &temp);
  }
  
  static struct vertex *addUniqueVertex(struct rbTree *tree, int position, enum ggVertexType type)
  /* Find existing vertex if it exists, otherwise create and return new one. */
  {
  struct vertex *v = matchingVertex(tree, position, type);
  if (v == NULL)
      {
      lmAllocVar(tree->lm, v);
      v->position = position;
      v->type = type;
      rbTreeAdd(tree, v);
      }
  return v;
  }
  
  static struct edge *matchingEdge(struct rbTree *tree, struct vertex *start, struct vertex *end)
  /* Find matching edge. Return NULL if none. */
  {
  struct edge temp;
  temp.start = start;
  temp.end = end;
  return rbTreeFind(tree, &temp);
  }
  
  static struct edge *addUniqueEdge(struct rbTree *tree, struct vertex *start, struct vertex *end,
  	struct linkedBeds *lb)
  /* Find existing edge if it exists.  Otherwise create and return new one. 
   * Regardless add lb as evidence to edge. */
  {
  struct edge *e = matchingEdge(tree, start, end);
  if (e == NULL)
      {
      lmAllocVar(tree->lm, e);
      e->start = start;
      e->end = end;
      e->next = NULL;
      rbTreeAdd(tree, e);
      }
  struct evidence *ev;
  lmAllocVar(tree->lm, ev);
  ev->lb = lb;
  ev->start = start->position;
  ev->end = end->position;
  slAddHead(&e->evList, ev);
  return e;
  }
  
  static struct rbTree *makeVertexTree(struct linkedBeds *lbList)
  /* Make tree of unique vertices. */
  {
  struct rbTree *vertexTree = rbTreeNew(vertexCmp);
  struct linkedBeds *lb;
  for (lb = lbList; lb != NULL; lb = lb->next)
      {
      struct bed *bed;
      for (bed = lb->bedList; bed != NULL; bed = bed->next)
          {
  	/* Add very beginning and end, they'll be soft. */
  	addUniqueVertex(vertexTree, bed->chromStart, ggSoftStart);
  	addUniqueVertex(vertexTree, bed->chromEnd, ggSoftEnd);
  
  	/* Add internal hard ends. */
  	int i, lastBlock = bed->blockCount-1;
  	for (i=0; i<lastBlock; ++i)
  	    {
  	    addUniqueVertex(vertexTree, 
  	    	bed->chromStart + bed->chromStarts[i] + bed->blockSizes[i], ggHardEnd);
  	    addUniqueVertex(vertexTree, 
  	    	bed->chromStart + bed->chromStarts[i+1], ggHardStart);
  	    }
  	}
      }
  return vertexTree;
  }
  
  static struct rbTree *makeEdgeTree(struct linkedBeds *lbList, struct rbTree *vertexTree)
  /* Make tree of unique edges. */
  {
  struct rbTree *edgeTree = rbTreeNew(edgeCmp);
  struct linkedBeds *lb;
  for (lb = lbList; lb != NULL; lb = lb->next)
      {
      struct bed *bed, *nextBed;
      for (bed = lb->bedList; bed != NULL; bed = nextBed)
          {
  	nextBed = bed->next;
  
  	/* Loop to add all introns and all but last exon. */
  	struct vertex *start = matchingVertex(vertexTree, bed->chromStart, ggSoftStart);
  	int i, lastBlock = bed->blockCount-1;
  	for (i=0; i<lastBlock; ++i)
  	    {
  	    /* Add exon */
  	    struct vertex *end = matchingVertex(vertexTree,
  		start->position + bed->blockSizes[i], ggHardEnd);
  	    addUniqueEdge(edgeTree, start, end, lb);
  
  	    /* Add intron */
  	    start = matchingVertex(vertexTree, 
  		    bed->chromStart + bed->chromStarts[i+1], ggHardStart);
  	    addUniqueEdge(edgeTree, end, start, lb);
  	    }
  
  	/* Add final exon */
  	struct vertex *end = matchingVertex(vertexTree, bed->chromEnd, ggSoftEnd);
  	addUniqueEdge(edgeTree, start, end, lb);
  
  	/* If there's another bed to go, add a soft intron connecting it. */
  	if (nextBed != NULL)
  	    {
  	    start = matchingVertex(vertexTree, nextBed->chromStart, ggSoftStart);
  	    addUniqueEdge(edgeTree, end, start, lb);
  	    }
  	}
      }
  return edgeTree;
  }
  
  static void incVertexUses(void *item)
  /* Callback to clear edge->start->useCount and edge->end->useCount. */
  {
  struct edge *edge = item;
  edge->start->count += 1; 
  edge->end->count += 1;
  }
  
  static void removeUnusedVertices(struct rbTree *vertexTree, struct rbTree *edgeTree)
  /* Remove vertices not connected to any edges. */
  {
  /* Get vertex list and clear counts. */
  struct slRef *vRef, *vRefList = rbTreeItems(vertexTree);
  for (vRef = vRefList; vRef != NULL; vRef = vRef->next)
      {
      struct vertex *v = vRef->val;
      v->count = 0;
      }
  
  /* Inc counts of vertices connected to edges. */
  rbTreeTraverse(edgeTree, incVertexUses);
  
  /* Remove unused vertices. */
  for (vRef = vRefList; vRef != NULL; vRef = vRef->next)
      {
      struct vertex *v = vRef->val;
      if (v->count == 0)
          rbTreeRemove(vertexTree, v);
      }
  
  slFreeList(&vRefList);
  }
  
  static void removeEmptyEdges(struct rbTree *vertexTree, struct rbTree *edgeTree)
  /* Remove edges that are zero or negative in length. */
  {
  int removeCount = 0;
  struct slRef *edgeRef, *edgeRefList = rbTreeItems(edgeTree);
  for (edgeRef = edgeRefList; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
      if (edge->start->position >= edge->end->position)
          {
  	removeCount += 1;
  	rbTreeRemove(edgeTree, edge);
  	}
      }
  if (removeCount)
      removeUnusedVertices(vertexTree, edgeTree);
  slFreeList(&edgeRefList);
  }
  
  
  static struct dlList *sortedListFromTree(struct rbTree *tree)
  /* Create a double-linked list from tree. List will be sorted.  */
  {
  struct slRef *ref, *refList = rbTreeItems(tree);
  struct dlList *list = dlListNew();
  for (ref = refList; ref != NULL; ref = ref->next)
      dlAddValTail(list, ref->val);
  slFreeList(&refList);
  return list;
  }
  
  static void addWaysInAndOut(struct rbTree *vertexTree, struct rbTree *edgeTree,
  	struct lm *lm)
  /* Put a bunch of ways in and out of the graph onto the edges. */
  {
  struct slRef *edgeRef, *edgeRefList = rbTreeItems(edgeTree);
  for (edgeRef = edgeRefList; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
      struct slRef *inRef, *outRef;
      lmAllocVar(lm, inRef);
      lmAllocVar(lm, outRef);
      inRef->val = outRef->val = edge;
      slAddHead(&edge->start->waysOut, outRef);
      slAddHead(&edge->end->waysIn, inRef);
      }
  slFreeList(&edgeRefList);
  }
  
  static boolean checkSeqSimilar(struct dnaSeq *a, struct dnaSeq *b, int minScore)
  /* Check that two sequences are similar.  Assumes sequences are same size. */
  {
  int size = a->size;
  int size1 = size-1;
  return dnaScoreMatch(a->dna, b->dna, size) >= minScore 
  	|| dnaScoreMatch(a->dna+1, b->dna, size1) >= minScore
  	|| dnaScoreMatch(a->dna, b->dna+1, size1) >= minScore;
  }
  
  static boolean uncuttableEdge(struct edge *e)
  /* Return TRUE if any evidence for edge is uncuttable (ccds) */
  {
  struct evidence *ev;
  for (ev = e->evList; ev != NULL; ev = ev->next)
      if (noCutSource(ev->lb->sourceType))
          return TRUE;
  return FALSE;
  }
  
  static boolean vertexPartOfUncuttableEdge(struct vertex *v, boolean isStart)
  /* See if vertex is part of an uncuttable edge. */
  {
  struct slRef *eRef, *eRefList = (isStart ? v->waysOut : v->waysIn);
  assert(eRefList != NULL);
  for (eRef = eRefList; eRef != NULL; eRef = eRef->next)
      if (uncuttableEdge(eRef->val))
  	return TRUE;
  return FALSE;
  }
  
  static boolean checkSnapOk(struct vertex *vOld, struct vertex *vNew, boolean isRev, 
  	int bleedSize, int maxUncheckedSize, struct nibTwoCache *seqCache, char *chromName)
  /* Load sequence that corresponds to bleed-over, and  make sure that sequence of next
   * exon is similar. */
  {
  if (vertexPartOfUncuttableEdge(vOld, isRev))
      return FALSE;
  if (bleedSize <= maxUncheckedSize)
      return TRUE;
  int minScore = bleedSize-2;
  boolean similar = FALSE;
  if (isRev)
      {
      int oldStart = vOld->position;
      struct slRef *eRef;
      struct dnaSeq *oldSeq = nibTwoCacheSeqPartExt(seqCache, chromName, oldStart, bleedSize, FALSE, NULL);
      for (eRef = vNew->waysIn; eRef != NULL; eRef = eRef->next)
          {
  	struct edge *edge = eRef->val;
  	struct vertex *vRest = edge->start;
  	int newStart = vRest->position - bleedSize;
  	struct dnaSeq *newSeq = nibTwoCacheSeqPartExt(seqCache, chromName, newStart, bleedSize, FALSE, NULL);
  	similar = checkSeqSimilar(oldSeq, newSeq, minScore);
  	dnaSeqFree(&newSeq);
  	if (similar)
  	    break;
  	}
      dnaSeqFree(&oldSeq);
      }
  else
      {
      int oldStart = vOld->position - bleedSize;
      struct slRef *eRef;
      struct dnaSeq *oldSeq = nibTwoCacheSeqPartExt(seqCache, chromName, oldStart, bleedSize, FALSE, NULL);
      for (eRef = vNew->waysOut; eRef != NULL; eRef = eRef->next)
          {
  	struct edge *edge = eRef->val;
  	struct vertex *vRest = edge->end;
  	int newStart = vRest->position;
  	struct dnaSeq *newSeq = nibTwoCacheSeqPartExt(seqCache, chromName, newStart, bleedSize, FALSE, NULL);
  	similar = checkSeqSimilar(oldSeq, newSeq, minScore);
  	dnaSeqFree(&newSeq);
  	if (similar)
  	    break;
  	}
      }
  return similar;
  }
  
  static boolean snapVertex(struct dlNode *oldNode, int maxSnapSize, int maxUncheckedSnapSize,
  	struct nibTwoCache *seqCache, char *chromName)
  /* Snap vertex to nearby previous hard vertex. */
  {
  /* Figure out hard type we want to snap to, and return if not
   * soft to begin with. */
  struct vertex *vOld = oldNode->val;
  int oldPos = vOld->position;
  enum ggVertexType currentType = vOld->type;
  
  /* If no seqCache, then set up things so not bothering to check it. */
  if (seqCache == NULL)
      maxUncheckedSnapSize = maxSnapSize;
  
  /* Search backwards */
  if (currentType == ggSoftEnd)
      {
      struct dlNode *node;
      for (node = oldNode->prev; !dlStart(node); node = node->prev)
  	{
  	struct vertex *v = node->val;
  	int newPos = v->position;
  	int dif = oldPos - newPos;
  	if (dif > maxSnapSize)
  	    break;
  	if (v->type == ggHardEnd)
  	    {
  	    if (checkSnapOk(vOld, v, FALSE, dif, maxUncheckedSnapSize, seqCache, chromName))
  		{
  		vOld->movedTo = v;
  		verbose(3, "snapping vertex %d to %d\n", oldPos, newPos);
  		return TRUE;
  		}
  	    }
  	}
      }
  
  /* Search forward */
  else if (currentType == ggSoftStart)
      {
      struct dlNode *node;
      for (node = oldNode->next; !dlEnd(node); node = node->next)
  	{
  	struct vertex *v = node->val;
  	int newPos = v->position;
  	int dif = newPos - oldPos;
  	if (dif > maxSnapSize)
  	    break;
  	if (v->type == ggHardStart)
  	    {
  	    if (checkSnapOk(vOld, v, TRUE, dif, maxUncheckedSnapSize, seqCache, chromName))
  		{
  		vOld->movedTo = v;
  		return TRUE;
  		}
  	    }
  	}
      }
  return FALSE;
  }
  
  static void mergeOrAddEdge(struct rbTree *edgeTree, struct edge *edge)
  /* Add edge back if it is still unique, otherwise move evidence from
   * edge into existing edge. */
  {
  struct edge *existing = rbTreeFind(edgeTree, edge);
  if (existing)
      {
      existing->evList = slCat(existing->evList, edge->evList);
      edge->evList = NULL;
      }
  else
      rbTreeAdd(edgeTree, edge);
  }
  
  static void updateForwardedEdges(struct rbTree *edgeTree)
  /* Go through edges, following the movedTo's in the
   * vertices if need be. */
  {
  struct slRef *ref, *refList = rbTreeItems(edgeTree);
  int forwardCount = 0;
  for (ref = refList; ref != NULL; ref = ref->next)
      {
      struct edge *edge = ref->val;
      struct vertex *start = edge->start, *end = edge->end;
      if (start->movedTo || end->movedTo)
          {
  	++forwardCount;
  	rbTreeRemove(edgeTree, edge);
  	if (start->movedTo)
  	    edge->start = start->movedTo;
  	if (end->movedTo)
  	    edge->end = end->movedTo;
  	mergeOrAddEdge(edgeTree, edge);
  	}
      }
  if (forwardCount > 0)
      verbose(3, "Forwarded %d edges.\n", forwardCount);
  slFreeList(&refList);
  }
  
  static void snapSoftToCloseHard(struct rbTree *vertexTree, struct rbTree *edgeTree, int maxSnapSize,
  	int maxUncheckedSnapSize, struct nibTwoCache *seqCache, char *chromName)
  /* Snap hard vertices to nearby soft vertices of same type. */
  {
  struct lm *lm = lmInit(0);
  addWaysInAndOut(vertexTree, edgeTree, lm);
  struct dlList *vList = sortedListFromTree(vertexTree);
  struct dlNode *node;
  int snapCount = 0;
  for (node = vList->head; !dlEnd(node); node = node->next)
      {
      if (snapVertex(node, maxSnapSize, maxUncheckedSnapSize, seqCache, chromName))
  	{
  	rbTreeRemove(vertexTree, node->val);
          ++snapCount;
  	}
      }
  /* Clean up ways in and out since have removed some nodes. */
  for (node = vList->head; !dlEnd(node); node = node->next)
      {
      struct vertex *v = node->val;
      v->waysIn = v->waysOut = NULL;
      }
  if (snapCount > 0)
      {
      verbose(3, "Snapped %d close edges, now have %d vertices\n", snapCount, vertexTree->n);
      updateForwardedEdges(edgeTree);
      }
  dlListFree(&vList);
  lmCleanup(&lm);
  }
  
  int edgeRefCmpEnd(const void *va, const void *vb)
  /* Comparison to sort a slRef list of edges
   * by the edge end location. */
  {
  const struct slRef *a = *((struct slRef **)va);
  const struct slRef *b = *((struct slRef **)vb);
  struct edge *aEdge = a->val;
  struct edge *bEdge = b->val;
  return aEdge->end->position - bEdge->end->position;
  }
  
  int edgeRefCmpStartRev(const void *va, const void *vb)
  /* Comparison to sort a slRef list of edges
   * by the edge start location, with biggest biggest
   * starts first. */
  {
  const struct slRef *a = *((struct slRef **)va);
  const struct slRef *b = *((struct slRef **)vb);
  struct edge *aEdge = a->val;
  struct edge *bEdge = b->val;
  return bEdge->end->position - aEdge->end->position;
  }
  
  
  
  int snapHalfHardForward(struct vertex *v, struct rbTree *edgeTree,
  	enum ggVertexType softType, enum ggVertexType hardType)
  /* V is a hard vertex. Try to snap soft end vertices connected to v
   * to nearest hard vertex connected to v. */
  {
  int snapCount = 0;
  // enum ggVertex otherHardType = (hardType == ggHardStart ? ggHardEnd : ggHardStart);
  slSort(&v->waysOut, edgeRefCmpEnd); 
  struct slRef *hardRef = v->waysOut, *softRef = v->waysOut;
  for (;;)
      {
      /* Advance softRef to next soft ended edge. */
      for (;softRef != NULL; softRef = softRef->next)
          {
  	struct edge *edge = softRef->val;
  	if (edge->end->type == softType)
  	    break;
  	}
      if (softRef == NULL)
          break;
      struct edge *softEdge = softRef->val;
  
      /* If hardRef is before softRef (or it's not hard) advance it */
      for (;hardRef != NULL; hardRef = hardRef->next)
          {
  	struct edge *edge = hardRef->val;
  	if (edge->end->type == hardType && edge->end->position >= softEdge->end->position)
  	   break;
  	}
      if (hardRef == NULL)
          break;
      rbTreeRemove(edgeTree, softEdge);
      struct edge *hardEdge = hardRef->val;
      verbose(3, "Snapping half-hard edge ending at %d to %d\n", softEdge->end->position, 
  	    hardEdge->end->position); 
      softEdge->end = hardEdge->end;
      ++snapCount;
      mergeOrAddEdge(edgeTree, softEdge);
      softRef = softRef->next;
      }
  if (snapCount > 0)
      verbose(3, "Snapped %d forward\n", snapCount);
  return snapCount;
  }
  
  int snapHalfHardBackward(struct vertex *v, struct rbTree *edgeTree,
  	enum ggVertexType softType, enum ggVertexType hardType)
  /* V is a hard vertex. Try to snap soft start vertices connected to v
   * to nearest hard vertex connected to v. */
  {
  int snapCount = 0;
  // enum ggVertex otherHardType = (hardType == ggHardStart ? ggHardEnd : ggHardStart);
  slSort(&v->waysIn, edgeRefCmpStartRev); 
  struct slRef *hardRef = v->waysIn, *softRef = v->waysIn;
  for (;;)
      {
      /* Advance softRef to next soft ended edge. */
      for (;softRef != NULL; softRef = softRef->next)
          {
  	struct edge *edge = softRef->val;
  	if (edge->start->type == softType)
  	    break;
  	}
      if (softRef == NULL)
          break;
      struct edge *softEdge = softRef->val;
  
      /* If hardRef is before softRef (or it's not hard) advance it */
      for (;hardRef != NULL; hardRef = hardRef->next)
          {
  	struct edge *edge = hardRef->val;
  	if (edge->start->type == hardType && edge->start->position <= softEdge->start->position)
  	   break;
  	}
      if (hardRef == NULL)
          break;
      rbTreeRemove(edgeTree, softEdge);
      struct edge *hardEdge = hardRef->val;
      verbose(3, "Snapping half-hard edge starting at %d to %d\n", softEdge->start->position, 
  	    hardEdge->start->position); 
      softEdge->start = hardEdge->start;
      ++snapCount;
      mergeOrAddEdge(edgeTree, softEdge);
      softRef = softRef->next;
      }
  if (snapCount > 0)
      verbose(3, "Snapped %d reverse\n", snapCount);
  return snapCount;
  }
  
  static int snapHalfHard(struct vertex *v, struct rbTree *edgeTree)
  /* Snap soft ends of half-hard edges connected to v to 
   * compatible hard end.  Returns number of snaps it's made. */
  {
  enum ggVertexType currentType = v->type;
  enum ggVertexType softType, hardType;
  
  /* Figure out what types look for at the other end.  If we're
   * a soft vertex we can't do anything. */
  if (currentType == ggHardStart)
      {
      softType = ggSoftEnd;
      hardType = ggHardEnd;
      }
  else if (currentType == ggHardEnd)
      {
      softType = ggSoftStart;
      hardType = ggHardStart;
      }
  else
      return 0;
  
  int snapCount = 0;
  snapCount += snapHalfHardForward(v, edgeTree, softType, hardType);
  snapCount += snapHalfHardBackward(v, edgeTree, softType, hardType);
  return snapCount;
  }
  
  static void snapHalfHards(struct rbTree *vertexTree, struct rbTree *edgeTree)
  /* Snap edges that are half hard to edges that are fully hard and that
   * share the original hard end */
  {
  struct lm *lm = lmInit(0);
  addWaysInAndOut(vertexTree, edgeTree, lm);
  struct slRef *vRef, *vRefList = rbTreeItems(vertexTree);
  for (vRef = vRefList; vRef != NULL; vRef = vRef->next)
      {
      struct vertex *v = vRef->val;
      snapHalfHard(v, edgeTree);
      v->waysIn = v->waysOut = NULL;	/* These are trashed so remove them. */
      }
  slFreeList(&vRefList);
  removeUnusedVertices(vertexTree, edgeTree);
  lmCleanup(&lm);
  }
  
  struct sourceAndPos
  /* A vertex source and position. */
      {
      struct sourceAndPos *next;
      int position;		/* Position in genome. */
      bool trustedSource;		/* If true this is trusted source (refSeq) */
      };
  
  int sourceAndPosCmp(const void *va, const void *vb)
  /* Compare two sourceAndPos by position. */
  {
  const struct sourceAndPos *a = *((struct sourceAndPos **)va);
  const struct sourceAndPos *b = *((struct sourceAndPos **)vb);
  return a->position - b->position;
  }
  
  int sourceAndPosCmpRev(const void *va, const void *vb)
  /* Compare two sourceAndPos in reverse direction. */
  {
  const struct sourceAndPos *a = *((struct sourceAndPos **)va);
  const struct sourceAndPos *b = *((struct sourceAndPos **)vb);
  return b->position - a->position;
  }
  
  static struct vertex *consensusVertex(struct rbTree *vertexTree, struct sourceAndPos *list, 
  	int listSize, enum ggVertexType softType)
  /* Return first trusted (refSeq) vertex, or if no trusted one, then 
   * a soft vertex 1/4 of way (rounded down) through list. */
  {
  struct sourceAndPos *el;
  for (el = list; el != NULL; el = el->next)
      {
      if (el->trustedSource)
          break;
      }
  if (el == NULL)
      el = slElementFromIx(list, listSize/4);
  return matchingVertex(vertexTree, el->position, softType);
  }
  
  static void halfConsensusForward(struct vertex *v, 
  	struct rbTree *vertexTree, struct rbTree *edgeTree,
  	enum ggVertexType softType, struct lm *lm)
  /* Figure out consensus end of all edges beginning at v that have soft end. */
  {
  /* Collect a list of all attached softies. */
  struct sourceAndPos *list = NULL, *el;
  struct slRef *edgeRef;
  int softCount = 0;
  for (edgeRef = v->waysOut; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
      struct vertex *v = edge->end;
      if (v->type == softType)
          {
  	struct evidence *ev;
  	for (ev = edge->evList; ev != NULL; ev = ev->next)
  	    {
  	    lmAllocVar(lm, el);
  	    el->position = ev->end;
  	    el->trustedSource = trustedSource(ev->lb->sourceType);
  	    slAddHead(&list, el);
  	    ++softCount;
  	    }
  	}
      }
  
  /* See if have enough elements to make consensus forming
   * worthwhile. */
  if (softCount > 1)
      {
      slSort(&list, sourceAndPosCmpRev);
      struct vertex *end = consensusVertex(vertexTree, list, softCount, softType);
      for (edgeRef = v->waysOut; edgeRef != NULL; edgeRef = edgeRef->next)
  	{
  	struct edge *edge = edgeRef->val;
  	struct vertex *v = edge->end;
  	if (v != end && v->type == softType)
  	    {
  	    rbTreeRemove(edgeTree, edge);
  	    verbose(3, "Performing half-hard consensus: moving edge end from %d to %d\n",
  		    edge->end->position, end->position);
  	    edge->end = end;
  	    mergeOrAddEdge(edgeTree, edge);  // Will always merge. 
  	    }
  	}
      }
  }
  
  
  static void halfConsensusBackward(struct vertex *v, 
  	struct rbTree *vertexTree, struct rbTree *edgeTree,
  	enum ggVertexType softType, struct lm *lm)
  /* Figure out consensus start of all edges end at v that have soft start. */
  {
  /* Collect a list of all attached softies. */
  struct sourceAndPos *list = NULL, *el;
  struct slRef *edgeRef;
  int softCount = 0;
  for (edgeRef = v->waysIn; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
      struct vertex *v = edge->start;
      if (v->type == softType)
          {
  	struct evidence *ev;
  	for (ev = edge->evList; ev != NULL; ev = ev->next)
  	    {
  	    lmAllocVar(lm, el);
  	    el->position = ev->start;
  	    el->trustedSource = trustedSource(ev->lb->sourceType);
  	    slAddHead(&list, el);
  	    ++softCount;
  	    }
  	}
      }
  
  /* See if have enough elements to make consensus forming
   * worthwhile. */
  if (softCount > 1)
      {
      slSort(&list, sourceAndPosCmp);
      struct vertex *start = consensusVertex(vertexTree, list, softCount, softType);
      for (edgeRef = v->waysIn; edgeRef != NULL; edgeRef = edgeRef->next)
  	{
  	struct edge *edge = edgeRef->val;
  	struct vertex *v = edge->start;
  	if (v != start && v->type == softType)
  	    {
  	    rbTreeRemove(edgeTree, edge);
  	    verbose(3, "Performing half-hard consensus: moving edge start from %d to %d\n",
  		    edge->start->position, start->position);
  	    edge->start = start;
  	    mergeOrAddEdge(edgeTree, edge);  // Will always merge. 
  	    }
  	}
      }
  }
  
  static void halfHardConsensus(struct vertex *v, 
  	struct rbTree *vertexTree, struct rbTree *edgeTree,
  	struct lm *lm)
  /* Form consensus of soft vertices connected to hard vertex v. 
   * Consensus is:
   *    1) The largest refSeq.
   *    2) If no refSeq, something 3/4 of the way to the largest. */
  {
  enum ggVertexType currentType = v->type;
  enum ggVertexType softType;
  
  /* Figure out what types look for at the other end.  If we're
   * a soft vertex we can't do anything. */
  if (currentType == ggHardStart)
      softType = ggSoftEnd;
  else if (currentType == ggHardEnd)
      softType = ggSoftStart;
  else
      return;
  
  halfConsensusForward(v, vertexTree, edgeTree, softType, lm);
  halfConsensusBackward(v, vertexTree, edgeTree, softType, lm);
  }
  
  static void halfHardConsensuses(struct rbTree *vertexTree, struct rbTree *edgeTree)
  /* Collect soft edges that share a hard edge. Move all of them to a consensus
   * vertex, which is either:
   *    1) The largest refSeq.
   *    2) If no refSeq, something 3/4 of the way to the largest. */
  {
  struct lm *lm = lmInit(0);
  addWaysInAndOut(vertexTree, edgeTree, lm);
  struct slRef *vRef, *vRefList = rbTreeItems(vertexTree);
  for (vRef = vRefList; vRef != NULL; vRef = vRef->next)
      {
      struct vertex *v = vRef->val;
      halfHardConsensus(v, vertexTree, edgeTree, lm);
      v->waysIn = v->waysOut = NULL;	/* These are trashed so remove them. */
      }
  slFreeList(&vRefList);
  removeUnusedVertices(vertexTree, edgeTree);
  lmCleanup(&lm);
  }
  
  struct edge *edgeFromConsensusOfEvidence(struct rbTree *vertexTree, struct evidence *evList,
  	struct lm *lm)
  /* Attempt to create a single edge from a list of overlapping evidence ranges.
   * The start will be the consensus of all evidence starts.  Likewise
   * the end will be the consensus of all evidence ends.  The evidence that
   * overlaps this edge will be included in the edge. */
  {
  /* Gather up lists of starts and ends. */
  struct sourceAndPos *startList = NULL, *endList = NULL;
  struct evidence *ev, *nextEv;
  int listSize = 0;
  for (ev = evList; ev != NULL; ev = ev->next)
      {
      struct sourceAndPos *x;
      boolean trusted = trustedSource(ev->lb->sourceType);
      lmAllocVar(lm, x);
      x->position = ev->start;
      x->trustedSource = trusted;
      slAddHead(&startList, x);
      lmAllocVar(lm, x);
      x->position = ev->end;
      x->trustedSource = trusted;
      slAddHead(&endList, x);
      ++listSize;
      }
  
  /* Get consensus starts and ends. */
  slSort(&startList, sourceAndPosCmp);
  struct vertex *start = consensusVertex(vertexTree, startList, listSize, ggSoftStart);
  slSort(&endList, sourceAndPosCmpRev);
  struct vertex *end = consensusVertex(vertexTree, endList, listSize, ggSoftEnd);
  
  /* Make edge */
  struct edge *edge;
  AllocVar(edge);
  edge->start = start;
  edge->end = end;
  edge->next = NULL;
  
  /* Add overlapping evidence to edge. */
  for (ev = evList; ev != NULL; ev = nextEv)
      {
      nextEv = ev->next;
      if (rangeIntersection(ev->start, ev->end, start->position, end->position) > 0)
          slAddHead(&edge->evList, ev);
      }
  
  return edge;
  }
  
  static void removeEnclosedDoubleSofts(struct rbTree *vertexTree, struct rbTree *edgeTree, 
  	int maxBleedOver, double singleExonMaxOverlap)
  /* Move double-softs that overlap spliced things to a very great extent into
   * the spliced things. Also remove tiny double-softs (no more than 2*maxBleedOver). */
  {
  /* Traverse graph and build up range tree covering spliced exons.  For each 
   * range of overlapping exons, assemble a singly-linked list of all exons in 
   * the range */
  struct rbTree *rangeTree = rangeTreeNew(0);
  struct slRef *edgeRef, *edgeRefList = rbTreeItems(edgeTree);
  int removedCount = 0;
  for (edgeRef = edgeRefList; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
      struct vertex *start = edge->start;
      struct vertex *end = edge->end;
      if (start->type == ggHardStart || end->type == ggHardEnd)
  	{
  	rangeTreeAddValList(rangeTree, start->position, end->position, edge);
  	}
      }
  
  /* Traverse graph yet one more time looking for doubly-soft exons
   * that are overlapping the spliced exons. */
  for (edgeRef = edgeRefList; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
      struct vertex *start = edge->start;
      struct vertex *end = edge->end;
      if (start->type == ggSoftStart && end->type == ggSoftEnd)
          {
  	int s = start->position;
  	int e = end->position;
  	int size = e - s;
  	if (size <= maxBleedOver+maxBleedOver)
  	     {
  	     /* Tiny case, just remove edge and forget it. */
  	     verbose(3, "Removing tiny double-soft edge from %d to %d\n", s, e);
  	     rbTreeRemove(edgeTree, edge);
  	     ++removedCount;
  	     }
  	else
  	     {
  	     /* Normal case, look for exon list that encloses us, and
  	      * if any single exon in that list encloses us, merge into it. */
  	     int splicedOverlap = rangeTreeOverlapSize(rangeTree, s, e);
  	     if (splicedOverlap > 0 && splicedOverlap > singleExonMaxOverlap*size)
  	         {
  		 if (!trustedEdge(edge))
  		     {
  		     /* Once we find a range that overlaps the doubly-soft edge, find 
  		      * (half-hard or better) edge from that range that encloses the 
  		      * doubly soft edge. */
  		     struct range *r = rangeTreeMaxOverlapping(rangeTree, s, e);
  		     struct edge *nextEdge, *edgeList = r->val;
  		     struct edge *enclosingEdge = NULL;
  		     for (nextEdge = edgeList; edgeList != NULL; edgeList = edgeList->next)
  			 {
  			 if (encloses(nextEdge, edge))
  			     {
  			     enclosingEdge = nextEdge;
  			     }
  			 }
  		     if (enclosingEdge != NULL) 
  			 {
  			 enclosingEdge->evList = slCat(enclosingEdge->evList, edge->evList);
  			 edge->evList = NULL;
  			 verbose(3, "Removing doubly-soft edge %d-%d, reassigning to %d-%d\n",
  				 s, e, enclosingEdge->start->position, 
  				 enclosingEdge->end->position);
  			 rbTreeRemove(edgeTree, edge);
  			 ++removedCount;
  			 }
  		     }
  		 }
  	     }
  	}
      }
  
  /* Clean up and go home. */
  if (removedCount > 0)
      removeUnusedVertices(vertexTree, edgeTree);
  for (edgeRef = edgeRefList; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *nextEdge, *edge = edgeRef->val;
      while (edge != NULL) 
  	{
  	nextEdge = edge->next;
  	edge->next = NULL;
  	edge = nextEdge;
  	}
      }
  slFreeList(&edgeRefList);
  rbTreeFree(&rangeTree);
  }
  
  static void mergeDoubleSofts(struct rbTree *vertexTree, struct rbTree *edgeTree)
  /* Merge together overlapping edges with soft ends. */
  {
  struct mergedEdge
  /* Hold together info on a merged edge. */
      {
      struct evidence *evidence;
      };
  
  /* Traverse graph and build up range tree.  Each node in the range tree
   * will represent the bounds of coordinates of overlapping double softs */
  struct rbTree *rangeTree = rangeTreeNew(0);
  struct slRef *edgeRef, *edgeRefList = rbTreeItems(edgeTree);
  for (edgeRef = edgeRefList; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
      struct vertex *start = edge->start;
      struct vertex *end = edge->end;
      if (start->type == ggSoftStart && end->type == ggSoftEnd)
          rangeTreeAdd(rangeTree, start->position, end->position);
      }
  
  /* Traverse graph again merging edges */
  for (edgeRef = edgeRefList; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
      struct vertex *start= edge->start;
      struct vertex *end = edge->end;
      if (start->type == ggSoftStart && end->type == ggSoftEnd)
          {
  	struct range *r = rangeTreeFindEnclosing(rangeTree,
  		start->position, end->position);
  	assert(r != NULL);
  	/* At this point, r represents the bounds of a double-soft
  	 * region that encompasses this edge.  Collect the set of
  	 * evidence of edges overlapping this range */
          struct mergedEdge *mergeEdge = r->val;
          if (mergeEdge == NULL)
              {
              lmAllocVar(rangeTree->lm, mergeEdge);
              r->val = mergeEdge;
              }
          mergeEdge->evidence = slCat(edge->evList, mergeEdge->evidence);
  	verbose(3, "Merging doubly-soft edge (%d,%d) into range (%d,%d)\n", 
  		start->position, end->position, r->start, r->end);
          edge->evList = NULL;
          rbTreeRemove(edgeTree, edge);
  	}
      }
  
  /* Traverse merged edge list, making a single edge from each range. At this point,
   * each range will have some evidence attached to it, from each of the double softs
   * that fall within the range.  From all of this evidence, make a single consensus edge */
  struct range *r;
  struct lm *lm = lmInit(0);
  for (r = rangeTreeList(rangeTree); r != NULL; r = r->next)
      {
      struct mergedEdge *mergedEdge = r->val;
      struct edge *edge = edgeFromConsensusOfEvidence(vertexTree, mergedEdge->evidence, lm);
      if (edge != NULL)
          rbTreeAdd(edgeTree, edge);
      verbose(3, "Deriving edge (%d,%d) from all the double softs in range (%d,%d)\n", 
  	    edge->start->position, edge->end->position, r->start, r->end);
      }
  
  
  /* Clean up and go home. */
  lmCleanup(&lm);
  removeUnusedVertices(vertexTree, edgeTree);
  slFreeList(&edgeRefList);
  rbTreeFree(&rangeTree);
  }
  
  #ifdef DEBUG
  static void dumpVertices(struct rbTree *vertexTree)
  {
  struct slRef *vRef, *vRefList = rbTreeItems(vertexTree);
  for (vRef = vRefList; vRef != NULL; vRef = vRef->next)
      {
      struct vertex *v = vRef->val;
      printf(" %d(%d)", v->position, v->type);
      }
  printf("\n");
  slFreeList(&vRefList);
  }
  #endif /* DEBUG */
  
  struct txGraph *treeTreeToTxg(struct rbTree *vertexTree, struct rbTree *edgeTree, char *name,
  	struct linkedBeds *lbList)
  /* Convert from vertexTree/edgeTree representation to txGraph */
  {
  /* Allocate txGraph, and fill in some of the relatively easy fields. */
  struct txGraph *txg;
  AllocVar(txg);
  struct bed *bed = lbList->bedList;
  txg->name = cloneString(name);
  txg->tName = cloneString(bed->chrom);
  txg->strand[0] = bed->strand[0];
  txg->vertexCount = vertexTree->n;
  txg->edgeCount = edgeTree->n;
  
  /* Get vertex list and number sequentially. Fill in vertex array */
  int i;
  struct slRef *vRef, *vRefList = rbTreeItems(vertexTree);
  struct txVertex *tv = AllocArray(txg->vertices, vertexTree->n);
  int tStart = BIGNUM, tEnd = -BIGNUM;
  for (vRef = vRefList, i=0; vRef != NULL; vRef = vRef->next, i++)
      {
      struct vertex *v = vRef->val;
      int position = v->position;
      v->count = i;
      tv->position = position;
      tv->type = v->type;
      tStart = min(tStart, position);
      tEnd = max(tEnd, position);
      ++tv;
      }
  slFreeList(&vRefList);
  
  /* Fill in overall bounds. */
  txg->tStart = tStart;
  txg->tEnd = tEnd;
  
  /* Deal with sources. */
  int sourceCount = 0;
  struct linkedBeds *lb;
  for (lb = lbList; lb != NULL; lb = lb->next)
      lb->id = sourceCount++;
  struct txSource *source = AllocArray(txg->sources, sourceCount);
  for (lb = lbList; lb != NULL; lb = lb->next)
      {
      source->accession = cloneString(lb->bedList->name);
      source->type = cloneString(lb->sourceType);
      ++source;
      }
  txg->sourceCount = sourceCount;
  
  /* Convert edges */
  struct slRef *edgeRef, *edgeRefList = rbTreeItems(edgeTree);
  for (edgeRef = edgeRefList; edgeRef != NULL; edgeRef = edgeRef->next)
      {
      struct edge *edge = edgeRef->val;
  
      /* Allocate edge and fill in start and end. */
      struct txEdge *te;
      AllocVar(te);
      struct vertex *start = edge->start;
      struct vertex *end = edge->end;
      te->startIx = start->count;
      te->endIx = end->count;
  
      /* Figure out whether it's an intron or exon. */
      enum ggVertexType startType = edge->start->type;
      if (startType == ggHardStart || startType == ggSoftStart)
          te->type = ggExon;
      else
          te->type = ggIntron;
  
      /* Convert the evidence. */
      struct evidence *ev;
      for (ev = edge->evList; ev != NULL; ev = ev->next)
          {
  	struct txEvidence *tev;
  	AllocVar(tev);
  	tev->sourceId = ev->lb->id;
  	tev->start = ev->start;
  	tev->end = ev->end;
  	slAddHead(&te->evList, tev);
  	te->evCount += 1;
  	}
      slReverse(&te->evList);
  
      slAddHead(&txg->edgeList, te);
      }
  slReverse(&txg->edgeList);
  return txg;
  }
  
  struct txGraph *makeGraph(struct linkedBeds *lbList, int maxBleedOver, 
  	int maxUncheckedBleed, struct nibTwoCache *seqCache,
  	double singleExonMaxOverlap, char *name)
  /* Create a graph corresponding to linkedBedsList.
   * The maxBleedOver parameter controls how much of a soft edge that
   * can be cut off when snapping to a hard edge.  The singleExonMaxOverlap
   * controls what ratio of a single exon transcript can overlap spliced 
   * transcripts */
  {
  char *chromName = lbList->bedList->chrom;
  
  /* Create tree of all unique vertices. */
  struct rbTree *vertexTree = makeVertexTree(lbList);
  verbose(2, "%d unique vertices\n", vertexTree->n);
  
  /* Create tree of all unique edges */
  struct rbTree *edgeTree = makeEdgeTree(lbList, vertexTree);
  verbose(2, "%d unique edges\n", edgeTree->n);
  
  snapSoftToCloseHard(vertexTree, edgeTree, maxBleedOver, maxUncheckedBleed, seqCache, chromName);
  verbose(2, "%d edges, %d vertices after snapSoftToCloseHard\n", 
  	edgeTree->n, vertexTree->n);
  
  removeEmptyEdges(vertexTree, edgeTree);
  verbose(2, "%d edges, %d vertices after removeEmptyEdges\n",
  	edgeTree->n, vertexTree->n);
  
  snapHalfHards(vertexTree, edgeTree);
  verbose(2, "%d edges, %d vertices after snapHalfHards\n", 
  	edgeTree->n, vertexTree->n);
  
  halfHardConsensuses(vertexTree, edgeTree);
  verbose(2, "%d edges, %d vertices after medianHalfHards\n", 
  	edgeTree->n, vertexTree->n);
  
  removeEnclosedDoubleSofts(vertexTree, edgeTree, maxBleedOver, singleExonMaxOverlap);
  verbose(2, "%d edges, %d vertices after mergeEnclosedDoubleSofts\n",
  	edgeTree->n, vertexTree->n);
  
  mergeDoubleSofts(vertexTree, edgeTree);
  verbose(2, "%d edges, %d vertices after mergeDoubleSofts\n",
  	edgeTree->n, vertexTree->n);
  
  struct txGraph *txg = treeTreeToTxg(vertexTree, edgeTree, name, lbList);
  
  /* Clean up and go home. */
  rbTreeFree(&vertexTree);
  rbTreeFree(&edgeTree);
  return txg;
  }
  