intervalFind.c

Go to the documentation of this file.

#include "log.h"
#include "format.h"
#include "linestream.h"
#include "numUtil.h"
#include "intervalFind.h"



typedef struct {
  char* chromosome;
  int start;
  int end;
  Array sublist;  // of type Interval*
} SuperInterval;



static Array intervals = NULL;
static Array superIntervals = NULL;
static int superIntervalAssigned = 0;



int intervalFind_getNumberOfIntervals (void)
{
  return arrayMax (intervals);
}




Array intervalFind_getAllIntervals (void)
{
  return arrayCopy (intervals);
}



static void processCommaSeparatedList (Array results, char *str) 
{
  WordIter w;
  char *tok;

  w = wordIterCreate (str,",",0);
  while (tok = wordNext (w)) {
    if (tok[0] == '\0') {
      continue;
    }
    array (results,arrayMax (results),int) = atoi (tok);
  }
  wordIterDestroy (w);
}



void intervalFind_addIntervalsToSearchSpace (char* fileName, int source)
{
  LineStream ls;
  WordIter w;
  char *line;
  Array subIntervalStarts;
  Array subIntervalEnds;
  SubInterval *currSubInterval;
  Interval *currInterval;
  int i;

  if (intervals == NULL) {
    intervals = arrayCreate (100000,Interval);
  }
  ls = ls_createFromFile (fileName);
  while (line = ls_nextLine (ls)) {
    if (line[0] == '\0') {
      continue;
    }
    w = wordIterCreate (line,"\t",0);
    currInterval = arrayp (intervals,arrayMax (intervals),Interval);
    currInterval->source = source;
    currInterval->name = hlr_strdup (wordNext (w));
    currInterval->chromosome = hlr_strdup (wordNext (w));
    currInterval->strand = wordNext (w)[0];
    currInterval->start = atoi (wordNext (w));
    currInterval->end = atoi (wordNext (w));
    currInterval->subIntervalCount = atoi (wordNext (w));
    subIntervalStarts = arrayCreate (currInterval->subIntervalCount,int);
    subIntervalEnds = arrayCreate (currInterval->subIntervalCount,int);
    processCommaSeparatedList (subIntervalStarts,wordNext (w));
    processCommaSeparatedList (subIntervalEnds,wordNext (w));
    if (arrayMax (subIntervalStarts) != arrayMax (subIntervalEnds)) {
      die ("Unequal number of subIntervalStarts and subIntervalEnds");
    }
    currInterval->subIntervals = arrayCreate (currInterval->subIntervalCount,SubInterval);
    for (i = 0; i < currInterval->subIntervalCount; i++) {
      currSubInterval = arrayp (currInterval->subIntervals,arrayMax (currInterval->subIntervals),SubInterval);
      currSubInterval->start = arru (subIntervalStarts,i,int);
      currSubInterval->end = arru (subIntervalEnds,i,int);
    }
    arrayDestroy (subIntervalStarts);
    arrayDestroy (subIntervalEnds);
    wordIterDestroy (w);
  }
  ls_destroy (ls);
}



static int sortIntervalsByChromosomeAndStartAndEnd (Interval *a, Interval *b)
{
  int diff;

  diff = strcmp (a->chromosome,b->chromosome);
  if (diff != 0) {
    return diff;
  } 
  diff = a->start - b->start;
  if (diff != 0) {
    return diff;
  }
  return b->end - a->end;
}



static int sortSuperIntervalsByChromosomeAndStartAndEnd (SuperInterval *a, SuperInterval *b)
{
  int diff;

  diff = strcmp (a->chromosome,b->chromosome);
  if (diff != 0) {
    return diff;
  } 
  diff = a->start - b->start;
  if (diff != 0) {
    return diff;
  }
  return b->end - a->end;
}



static void assignSuperIntervals (void)
{
  int i,j;
  Interval *currInterval,*nextInterval;
  SuperInterval *currSuperInterval;

  superIntervals = arrayCreate (100000,SuperInterval);
  arraySort (intervals,(ARRAYORDERF)sortIntervalsByChromosomeAndStartAndEnd);
  i = 0;
  while (i < arrayMax (intervals)) {
    currInterval = arrp (intervals,i,Interval);
    currSuperInterval = arrayp (superIntervals,arrayMax (superIntervals),SuperInterval);
    currSuperInterval->chromosome = hlr_strdup (currInterval->chromosome);
    currSuperInterval->start = currInterval->start; 
    currSuperInterval->end = currInterval->end;
    currSuperInterval->sublist = arrayCreate (10,Interval*);
    array (currSuperInterval->sublist,arrayMax (currSuperInterval->sublist),Interval*) = currInterval;
    j = i + 1;
    while (j < arrayMax (intervals)) {
      nextInterval = arrp (intervals,j,Interval);
      if (strEqual (currInterval->chromosome,nextInterval->chromosome) &&
          currInterval->start <= nextInterval->start && 
          currInterval->end >= nextInterval->end) { 
        array (currSuperInterval->sublist,arrayMax (currSuperInterval->sublist),Interval*) = nextInterval;
      }
      else {
        break;
      }
      j++;
    }
    i = j;
  }
  arraySort (superIntervals,(ARRAYORDERF)sortSuperIntervalsByChromosomeAndStartAndEnd);
}



static void addIntervals (Array matchingIntervals, Array sublist, int start, int end) 
{
  int i;
  Interval *currInterval;

  for (i = 0; i < arrayMax (sublist); i++) {
    currInterval = arru (sublist,i,Interval*);
    if (rangeIntersection (currInterval->start,currInterval->end,start,end) >= 0) {
      array (matchingIntervals,arrayMax (matchingIntervals),Interval*) = currInterval;
    }
  }
}



Array intervalFind_getOverlappingIntervals (char* chromosome, int start, int end)
{
  SuperInterval testSuperInterval;
  SuperInterval *currSuperInterval;
  int i,index;
  static Array matchingIntervals = NULL;

  if (superIntervalAssigned == 0) {
     assignSuperIntervals ();
     superIntervalAssigned = 1;
  }
  if (matchingIntervals == NULL) {
    matchingIntervals = arrayCreate (20,Interval*);
  }
  else {
    arrayClear (matchingIntervals);
  }
  testSuperInterval.chromosome = chromosome;
  testSuperInterval.start = start;
  testSuperInterval.end = end;
  arrayFind (superIntervals,&testSuperInterval,&index,(ARRAYORDERF)sortSuperIntervalsByChromosomeAndStartAndEnd);
  // Index points to the location where testSuperInterval would be inserted
  i = index;
  while (i >= 0) {
    currSuperInterval = arrp (superIntervals,i,SuperInterval);
    if (!strEqual (currSuperInterval->chromosome,chromosome) || currSuperInterval->end < start) {
      break;
    }
    addIntervals (matchingIntervals,currSuperInterval->sublist,start,end); 
    i--;
  }
  i = index + 1;
  while (i < arrayMax (superIntervals)) {
    currSuperInterval = arrp (superIntervals,i,SuperInterval);
    if (!strEqual (currSuperInterval->chromosome,chromosome) ||  currSuperInterval->start > end) {
      break;
    }
    addIntervals (matchingIntervals,currSuperInterval->sublist,start,end); 
    i++;
  }
  return matchingIntervals;
}

---