view release on metacpan or  search on metacpan

GH/Sim4/sim4.2002-03-03/sim4b1.c  view on Meta::CPAN

  int diag_dist, diff;

  exon_list = NULL; 
  if (last_msp<0) return;

  /* Note: in new_exon, the 'flag' and 'length' fields need not be computed */
  mp = msp[last_msp];
  exon_list = new_exon (mp->pos1, mp->pos2, mp->pos1+mp->len-1, 
                        mp->pos2+mp->len-1, -1, 
                        (mp->len*MATCH_SCORE-mp->score)/(MATCH_SCORE-MISMATCH_SCORE), 0, exon_list);

  last_msp = mp->prev;

  while (last_msp>=0) {
    mp = msp[last_msp]; 
    if (((diag_dist=abs((exon_list->from2-exon_list->from1)-(mp->pos2-mp->pos1)))<=L) 
        && (exon_list->from2-(mp->pos2+mp->len-1))<MAX_INTERNAL_GAP) {
          /* merge with previous exon */
          exon_list->edist += diag_dist;
          exon_list->edist += (mp->len*MATCH_SCORE-mp->score)/(MATCH_SCORE-MISMATCH_SCORE);
          if ((diff=mp->pos2+mp->len-exon_list->from2)>0) {   /* overlap */
             int dist1, dist2;
             dist1 = get_edist(exon_list->from1,mp->pos2+mp->len-diff,
                               exon_list->from1+diff-1,mp->pos2+mp->len-1,s1,s2);
             dist2 = get_edist(mp->pos1+mp->len-diff,mp->pos2+mp->len-diff,
                               mp->pos1+mp->len-1,mp->pos2+mp->len-1,s1,s2);
             exon_list->edist -= max(dist1,dist2);
          } else if (diff<0) {  /* gap */
             exon_list->edist += 0.5*P*(-1)*diff; 
          }
          exon_list->to1 = max(exon_list->to1,mp->pos1+mp->len-1);
          exon_list->to2 = max(exon_list->to2,mp->pos2+mp->len-1);
          exon_list->from1 = min(exon_list->from1,mp->pos1);
          exon_list->from2 = min(exon_list->from2,mp->pos2);
    } else {
          /* new exon */
          exon_list = new_exon (mp->pos1, mp->pos2, mp->pos1+mp->len-1,
                                mp->pos2+mp->len-1, -1,
                                (mp->len*MATCH_SCORE-mp->score)/(MATCH_SCORE-MISMATCH_SCORE),
                                0, exon_list);
    }
    last_msp = mp->prev;
  } 
}

static int get_edist(int f1, int f2, int t1, int t2, uchar *seq1, uchar *seq2)
{
    uchar *s1, *s2, *q1, *q2;
    int dist=0;

    s1 = seq1+f1+1;   /* bc at this stage, the msp pos do not have added +1 */
    s2 = seq2+f2+1;
    q1 = seq1+t1+1;
    q2 = seq2+t2+1;

    while (s1<=q1 && s2<=q2) { dist += (*s1!=*s2); s1++; s2++; } 
    
    return dist;
}

/* ----------------------  print endpoints of exons  --------------------*/

#ifdef AUXUTILS 
static void find_introns(Exon *eleft, Intron **Ilist)
{
  Exon   *tmp_exon, *tmp_exon1;
  Intron *new, *tail;
  int     GTAG_score, CTAC_score;

  *Ilist = tail = NULL;
  if (!eleft) fatal("sim4b1.c: Something wrong in the exon list.\n");

  tmp_exon = eleft->next_exon;
  while ((tmp_exon!=NULL) && (tmp_exon1=tmp_exon->next_exon) &&
         tmp_exon1->to1) {
     
     new = (Intron *)ckalloc(sizeof(Intron)); 
     new->from1 = tmp_exon->to1+1;
     new->to1 = tmp_exon1->from1-1;
     new->from2 = tmp_exon->to2;
     new->to2 = tmp_exon1->from2;
     new->length = new->to1-new->from1+1;
     new->next_intron = NULL;
     if (!tail) *Ilist = new;
     else tail->next_intron = new; 
     tail = new;

     /* find orientation */
     GTAG_score = CTAC_score = 0;
     if (*(seq1+new->from1-1)=='G') GTAG_score++;
     else 
       if (*(seq1+new->from1-1)=='C') CTAC_score++;

     if (*(seq1+new->from1-1)=='T')   
          { GTAG_score++;  CTAC_score++; }

     if (*(seq1+new->to1-1)=='A') 
          { GTAG_score++;  CTAC_score++; }

     if (*(seq1+new->to1-1)=='G')  GTAG_score++;
     else 
       if (*(seq1+new->to1-1)=='C') CTAC_score++;

     if  (GTAG_score>=CTAC_score) 
          new->orientation = '+';
     else new->orientation = 'c';

     tmp_exon = tmp_exon1;
  }
 
}
#endif 


/* should only be called when (file_type==EST_GEN) && (match_ori==BWD) */
void complement_exons(Exon **left, int M, int N)
{ 
  Exon  *tmp_block, *right;
  char   prev, ch;

  prev = 'U'; /* unknown; should trigger error */