ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gffread/gff_utils.cpp
(Generate patch)
# Line 3 | Line 3
3   extern bool verbose;
4   extern bool debugMode;
5  
6 + //bool debugState=false;
7 +
8   void printFasta(FILE* f, GStr& defline, char* seq, int seqlen) {
9   if (seq==NULL) return;
10   int len=(seqlen>0)?seqlen:strlen(seq);
# Line 315 | Line 317
317   }
318  
319  
320 < void printLocus(GffLocus* loc, const char* pre=NULL) {
320 > void printLocus(GffLocus* loc, const char* pre) {
321    if (pre!=NULL) fprintf(stderr, "%s", pre);
322    GMessage(" [%d-%d] : ", loc->start, loc->end);
323    GMessage("%s",loc->rnas[0]->getID());
# Line 339 | Line 341
341     }
342   }
343  
344 < void placeGf(GffObj* t, GenomicSeqData* gdata, bool doCluster, bool collapseRedundant, bool matchAllIntrons, bool fuzzSpan) {
345 <  //GMessage(">>Placing transcript %s\n", t->getID());
344 > void placeGf(GffObj* t, GenomicSeqData* gdata, bool doCluster, bool collapseRedundant,
345 >                                               bool matchAllIntrons, bool fuzzSpan) {
346    GTData* tdata=new GTData(t);
347    gdata->tdata.Add(tdata);
348    int tidx=-1;
349 +  /*
350 +  if (debug) {
351 +     GMessage(">>Placing transcript %s\n", t->getID());
352 +     debugState=true;
353 +     }
354 +    else debugState=false;
355 +   */
356    if (t->exons.Count()>0)
357                tidx=gdata->rnas.Add(t); //added it in sorted order
358              else {
# Line 356 | Line 365
365         //GMessage("  <<make it first locus %d-%d \n",t->start, t->end);
366         return;
367         }
368 +   /*    
369    //DEBUG: show available loci:
370 <  // GMessage("  [%d loci already:\n", gdata->loci.Count());
371 <  //for (int l=0;l<gdata->loci.Count();l++) {
372 <  //    printLocus(gdata->loci[l]);
373 <  //    }
370 >   if (debug) {
371 >    GMessage("  [%d loci already:\n", gdata->loci.Count());
372 >    for (int l=0;l<gdata->loci.Count();l++) {
373 >       printLocus(gdata->loci[l]);
374 >       }
375 >    }
376 >  */
377    int nidx=qsearch_gloci(t->end, gdata->loci); //get index of nearest locus starting just ABOVE t->end
378    //GMessage("\tlooking up end coord %d in gdata->loci.. (qsearch got nidx=%d)\n", t->end, nidx);
379    if (nidx==0) {
380       //cannot have any overlapping loci
381 <     //GMessage("  <<no ovls possible, create locus %d-%d \n",t->start, t->end);
381 >     //if (debug) GMessage("  <<no ovls possible, create locus %d-%d \n",t->start, t->end);
382       gdata->loci.Add(new GffLocus(t));
383       return;
384       }
# Line 373 | Line 386
386    int lfound=0; //count of parent loci
387    GArray<int> mrgloci(false);
388    GList<GffLocus> tloci(true); //candidate parent loci to adopt this
389 <  //GMessage("\tchecking all loci from %d to 0\n",nidx-1);
389 >  //if (debug) GMessage("\tchecking all loci from %d to 0\n",nidx-1);
390    for (int l=nidx-1;l>=0;l--) {
391        GffLocus& loc=*(gdata->loci[l]);
392        if (loc.strand!='.' && t->strand!='.'&& loc.strand!=t->strand) continue;
# Line 381 | Line 394
394             if (t->start-loc.start>GFF_MAX_LOCUS) break; //give up already
395             continue;
396             }
397 <      if (loc.start>t->end) continue;
398 <          //this should never be the case if nidx was found correctly
399 <      //GMessage(" !range overlap found with locus ");
400 <      //printLocus(&loc);
397 >      if (loc.start>t->end) {
398 >               //this should never be the case if nidx was found correctly
399 >               GMessage("Warning: qsearch_gloci found loc.start>t.end!(t=%s)\n", t->getID());
400 >               continue;
401 >               }
402 >      /*
403 >      if (debug) {
404 >          GMessage(" !range overlap found with locus ");
405 >          printLocus(&loc);
406 >          }
407 >      */
408        if (loc.add_RNA(t)) {
409           //will add this transcript to loc
410           lfound++;
# Line 410 | Line 430
430                }//for each transcript in the exon-overlapping locus
431            } //if doCollapseRedundant
432           } //overlapping locus
433 <      }
433 >      } //for each existing locus
434    if (lfound==0) {
435        //overlapping loci not found, create a locus with only this mRNA
436 <      //GMessage("  overlapping locus not found, create locus %d-%d \n",t->start, t->end);
436 >      /* if (debug) {
437 >        GMessage("  overlapping locus not found, create locus %d-%d \n",t->start, t->end);
438 >        }
439 >      */
440        int addidx=gdata->loci.Add(new GffLocus(t));
441        if (addidx<0) {
442 +         //should never be the case!
443           GMessage("  WARNING: new GffLocus(%s:%d-%d) not added!\n",t->getID(), t->start, t->end);
444           }
445        }
446 <   else if (lfound>1) {
447 <      //more than one loci found parenting this mRNA, merge loci
448 <      //if (lfound>2) GMessage(" merging %d loci \n",lfound);
449 <       lfound--;
446 >   else { //found at least one overlapping locus
447 >     lfound--;
448 >     int locidx=mrgloci[lfound];
449 >     GffLocus& loc=*(gdata->loci[locidx]);
450 >     //last locus index found is also the smallest index
451 >     if (lfound>0) {
452 >       //more than one loci found parenting this mRNA, merge loci
453 >       /* if (debug)
454 >          GMessage(" merging %d loci \n",lfound);
455 >       */
456         for (int l=0;l<lfound;l++) {
457 <          int mlidx=mrgloci[l]; //largest indices first, so it's safe to remove
458 <          gdata->loci[mrgloci[lfound]]->addMerge(*(gdata->loci[mlidx]), t);
459 <          gdata->loci.Delete(mlidx);
457 >          int mlidx=mrgloci[l];
458 >          loc.addMerge(*(gdata->loci[mlidx]), t);
459 >          gdata->loci.Delete(mlidx); //highest indices first, so it's safe to remove
460            }
461 <      }
461 >       }
462 >     int i=locidx;  
463 >     while (i>0 && loc.start<gdata->loci[i-1]->start) {
464 >       //bubble down until it's in the proper order
465 >       i--;
466 >       gdata->loci.Swap(i,i+1);
467 >       }
468 >     }//found at least one overlapping locus
469   }
470  
471   void collectLocusData(GList<GenomicSeqData>& ref_data) {

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines