Logo Search packages:      
Sourcecode: adun.app version File versions  Download package

PairwiseStrxAlignment.m

/* Copyright 2003  Alexander V. Diemand

    This file is part of MolTalk.

    MolTalk is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    MolTalk is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with MolTalk; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA 
 */
 
/* vim: set filetype=objc: */


#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <time.h>

#include "PairwiseStrxAlignment.oh"
#include "Structure.oh"
#include "Chain.oh"
#include "Residue.oh"
#include "Selection.oh"
#include "Matrix.oh"
#include "Matrix53.oh"
#include "Stream.oh"


#undef DEBUG_COMPUTING_TIME

#define GAPOPENINGPENALTY 4
#define GAPEXTENDPENALTY 1
#define MAXDIST 6.0
#define SIGNIFDISTTHRESHOLD 5.0


@implementation AlPos

-(id)init   //@nodoc
{
      [super init];
      res1 = nil;
      res2 = nil;
      significant = YES;
      distance = -1.0;
      return self;
}


-(void)dealloc    //@nodoc
{
      //printf("AlPos_dealloc\n");
      if (res1)
      {
            RELEASE(res1);
      }     
      if (res2)
      {
            RELEASE(res2);
      }
      [super dealloc];
}


-(NSString*)description
{
      char char1=0;
      char char2=0;
      if (res1 == nil)
      {
            char1 = '-';
      } else {
            char1 = [[res1 oneLetterCode] characterAtIndex:0];
            if (!significant) { char1 += 32; }
      }
      if (res2 == nil)
      {
            char2 = '-';
      } else {
            char2 = [[res2 oneLetterCode] characterAtIndex:0];
            if (!significant) { char2 += 32; }
      }
      if (distance >= 0.0)
      {
            return [NSString stringWithFormat: @"%c %c %1.2f",char1,char2,distance];
      } else {
            return [NSString stringWithFormat: @"%c %c     ",char1,char2];
      }
}


/*
 *   set flag
 */ 
-(id)signify
{
      significant = YES;
      return self;
}


/*
 *   remove flag
 */ 
-(id)designify
{
      significant = NO;
      return self;
}


/*
 *   true if one of the residues is missing, indicating a gap
 */
-(BOOL)isGapped
{
      if (res1==nil || res2==nil)
      {
            return YES;
      }
      return NO;
}


/*
 *   return reference to residue
 */
-(Residue*)res1
{
      return res1;
}


/*
 *   return reference to residue
 */
-(Residue*)res2
{
      return res2;
}


/*
 *   return distance between the two residues, or -1.0 for a gapped position
 */
-(double)distance
{
      if (res1!=nil && res2!=nil)
      {
            //return [res1 distanceCATo:res2];
            return distance;
      }
      return -1.0;
}


/*
 *   create alPos with pair of residues
 */ 
+(AlPos*)alposWithRes1:(Residue*)r1 res2:(Residue*)r2
{
      AlPos *res = [AlPos new];
      res->res1 = r1;
      if (r1)
      {
            RETAIN(res->res1);
      }
      res->res2 = r2;
      if (r2)
      {
            RETAIN(res->res2);
      }
      if (r1 && r2)
      {
            res->distance = [r1 distanceCATo:r2];
      } else {
            res->distance = -1.0;
      }
      return AUTORELEASE(res);
}

@end

@interface AlPos (Private)

/*
 *   setters
 */
-(void)res1:(Residue*)res1;
-(void)res2:(Residue*)res2;
-(void)distance:(double)dist;

@end

@implementation AlPos (Private)

-(void)res1:(Residue*)p_res1
{
      if (p_res1)
      {
            RETAIN(p_res1);
      }
      if (res1)
      {
            RELEASE(res1);
      }
      res1 = p_res1;
}

-(void)res2:(Residue*)p_res2
{
      if (p_res2)
      {
            RETAIN(p_res2);
      }
      if (res2)
      {
            RELEASE(res2);
      }
      res2 = p_res2;
}

-(void)distance:(double)dist
{
      if (res1 && res2)
      {
            distance = dist;
      } else {
            distance = -1.0;
      }
}

@end


/* = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = */

@implementation PairwiseStrxAlignment


-(id)init   //@nodoc
{
      [super init];
      chain1 = nil;
      chain2 = nil;
      positions = nil;
      calculated = NO;
      return self;
}


-(void)dealloc    //@nodoc
{
      if (positions)
      {
            RELEASE(positions);
      }
      if (chain1)
      {
            RELEASE(chain1);
      }     
      if (chain2)
      {
            RELEASE(chain2);
      }
      if (transformation)
      {
            RELEASE(transformation);
      }
      [super dealloc];
}


-(NSString*)description
{
      return @"PairwiseAlignment";
}


-(Chain*)chain1;
{
      return chain1;
}


-(Chain*)chain2
{
      return chain2;
}


-(int)countPairs
{
      if (!calculated)
      {
            NSLog(@"PairwiseStrxAlignment needs to be calculated first.");
            return -1;
      }
      int count=0;
      if (positions)
      {
            count = [positions count];
      }
      return count;
}


-(int)countUngappedPairs
{
      if (!calculated)
      {
            NSLog(@"PairwiseStrxAlignment needs to be calculated first.");
            return -1;
      }
      int count = 0;
      int j;
      int i=[positions count];
      AlPos *alpos;
      for (j=1; j<=i; j++)
      {
            alpos = [positions objectAtIndex: (i-j)];
            if (![alpos isGapped])
            {
                  count++;
            }
      }
      return count;
}


-(Matrix53*)getTransformation
{
      return transformation;
}
            

-(double)calculateRMSD
{
      if (!calculated)
      {
            NSLog(@"PairwiseStrxAlignment needs to be calculated first.");
            return -1.0;
      }
      double deviation = 0.0;
      double tdist;
      int count = 0;
      int i=[positions count];
      int j;
      AlPos *alpos;
      for (j=1; j<=i; j++)
      {
            alpos = [positions objectAtIndex: (i-j)];
            if (![alpos isGapped])
            {
                  tdist = [alpos distance];
                  deviation += (tdist * tdist);
                  count++;
            }
      }
      if (count==0)
      {
            return -1.0;
      }
      return sqrt(deviation / count);
}


/*
 *    after an initial superposition has been calculated, try to improve it by selecting residue pairs
 *    having at most 3.5 A distance and recalculate transformation.
 */
-(void)optimize
{
      if (!calculated)
      {
            NSLog(@"PairwiseStrxAlignment needs to be calculated first.");
            return;
      }
      Selection *matching1;
      Selection *matching2;
      int nowaligned=0;
      int maxaligned=0;
      int counter=1;
      [chain1 selectResiduesCloseTo: chain2 maxDistance: 3.5f sel1:&matching1 sel2:&matching2];
      nowaligned = [matching1 count];
      Matrix53 *rtop;
      while (matching1 && nowaligned>maxaligned)
      {
            maxaligned = nowaligned;
            //printf("####### ##### %d we have %d aligned residues\n",counter,nowaligned);
            rtop = [matching2 alignTo: matching1];
            [chain2 transformBy: rtop];
            [chain1 selectResiduesCloseTo: chain2 maxDistance: 3.5f sel1:&matching1 sel2:&matching2];
            nowaligned = [matching1 count];
            counter++;
      }
      //printf("# at the end we have %d aligned residues after %d runs\n",nowaligned,counter);
}


-(NSArray*)alignmentPositions
{
      if (!calculated)
      {
            NSLog(@"PairwiseStrxAlignment needs to be calculated first.");
            return nil;
      }
      return positions;
}


-(int)countPairsMaxDistance:(double)dist
{
      if (!calculated)
      {
            NSLog(@"PairwiseStrxAlignment needs to be calculated first.");
            return -1;
      }
      int count = 0;
      int i=[positions count];
      int j;
      AlPos *alpos;
      double actdist;
      for (j=1; j<=i; j++)
      {
            alpos = [positions objectAtIndex: (i-j)];
            actdist = [alpos distance];
            if (actdist >= 0.0 && actdist <= dist)
            {
                  count++;
            }
      }
      
      return count;
}


/*
 *   output derived structural alignment from computed superposition as a T-Coffee library.
 */
-(void)toStreamAsTCoffee:(Stream*)stream name1:(NSString*)name1 name2:(NSString*)name2
{
      if (!calculated)
      {
            NSLog(@"PairwiseStrxAlignment needs to be calculated first.");
            return;
      }

      int countA = 1;
      int countB = 1;
      int i=[positions count];
      int j;
      AlPos *alpos;
      double actdist;
      Residue *res1, *res2;
      if (![stream ok])
      {
            return;
      }
      /* write header */
      [stream writeString:@"! T-COFFEE_LIB_FORMAT_01\n"];
      [stream writeString:@"2\n"];
      /* clean up sequence of chain1 */
      NSMutableDictionary *chaindict1 = [NSMutableDictionary new];
      NSEnumerator *enumerator = [chain1 allResidues];
      NSNumber *resnum,*resnum2;
      while (enumerator && (res1 = [enumerator nextObject]))
      {
            if ([res1 isStandardAminoAcid])
            {
                  resnum = [res1 number];
                  [chaindict1 setObject:[NSNumber numberWithInt:countA] forKey:resnum];
                  countA++;
            }
      }
      /* clean up sequence of chain2 */
      NSMutableDictionary *chaindict2 = [NSMutableDictionary new];
      enumerator = [chain2 allResidues];
      while (enumerator && (res2 = [enumerator nextObject]))
      {
            if ([res2 isStandardAminoAcid])
            {
                  resnum2 = [res2 number];
                  [chaindict2 setObject:[NSNumber numberWithInt:countB] forKey:resnum2];
                  countB++;
            }
      }
      NSString *sequence = [chain1 get3DSequence]; // without 'X'
      [stream writeString:[NSString stringWithFormat:@"%@%d %d %@\n",[[chain1 structure] pdbcode],[chain1 code],[sequence length],sequence]];
      sequence = [chain2 get3DSequence]; // without 'X'
      [stream writeString:[NSString stringWithFormat:@"%@%d %d %@\n",[[chain2 structure] pdbcode],[chain2 code],[sequence length],sequence]];
      [stream writeString:@"#1 2\n"]; // as always the case
      
      int score;
      for (j=1; j<=i; j++)
      {
            alpos = [positions objectAtIndex: (i-j)];
            if (![alpos isGapped])
            {
                  res1 = [alpos res1];
                  res2 = [alpos res2];
                  actdist = [alpos distance];
                  if (actdist <= MAXDIST)
                  {
                        score = (int)floor((MAXDIST-actdist)*100.0f/MAXDIST+0.5f);
                  } else {
                        score = 1;
                  }
                  resnum = [chaindict1 objectForKey: [res1 number]];
                  resnum2 = [chaindict2 objectForKey: [res2 number]];
                  [stream writeString:[NSString stringWithFormat:@"%@ %@ %d\n",resnum,resnum2,score]];
            }
      }
      AUTORELEASE(chaindict1);
      AUTORELEASE(chaindict2);
      [stream writeString:@"! SEQ_1_TO_N\n"];
}


/*
 *   input a pairwise sequence alignment from a stream and interpret as a T-Coffee library.<br>
 *   the library indicates the selection of residues pairs which will be included in the calculation of the transformation
 */
-(void)fromStreamAsTCoffee:(Stream*)stream
{
      int nseqs;
      NSString *sname1, *sname2;
      int len1, len2;
      NSString *seq1, *seq2;
      Selection *sel1, *sel2;
      
      CREATE_AUTORELEASE_POOL(pool);
      
      /* get number of sequences */
      NSString *ln = [stream readStringLineLength: 10]; // should do it
      nseqs = [ln intValue];
      if (nseqs!=2)
      {
            NSLog(@"PairwiseStrxAlignment_fromStreamAsTCoffee: requires input file with 2 sequences (%d).",nseqs);
            RELEASE(pool);
            return;
      }
      
      /* read first sequence */
      ln = [stream readStringLineLength: 8192]; // might have a lot of them
      NSScanner *scanner = [NSScanner scannerWithString:ln];
      [scanner scanUpToString:@" " intoString:&sname1];
      [scanner scanInt:&len1];
      [scanner scanUpToString:@"\n" intoString:&seq1];
      if ([sname1 length] != 7)
      {
            sname1 = [[sname1 stringByAppendingString:@"     "] substringToIndex:7];
      }
      
      /* read second sequence */
      ln = [stream readStringLineLength: 8192]; // might have a lot of them
      scanner = [NSScanner scannerWithString:ln];
      [scanner scanUpToString:@" " intoString:&sname2];
      [scanner scanInt:&len2];
      [scanner scanUpToString:@"\n" intoString:&seq2];
      if ([sname2 length] != 7)
      {
            sname2 = [[sname2 stringByAppendingString:@"     "] substringToIndex:7];
      }     
      
      //printf("found: %d sequences\n%@ (%d) ...\n%@ (%d) ...\n",nseqs,sname1,len1,sname2,len2);
      
      /* check whether we are concerned at all */
      /* is sequence1 really our chain1? */
      NSString *t_chainid = [[NSString stringWithFormat:@"%@%d      ",[[chain1 structure] pdbcode],[chain1 code]] substringToIndex:7];
      if (![t_chainid isEqualToString: sname1])
      {
            NSLog(@"sequence1:%@ is not the same as chain1:%@",sname1,t_chainid);
            RELEASE(pool);
            return;
      }
      /* is sequence2 really our chain2? */
      t_chainid = [[NSString stringWithFormat:@"%@%d      ",[[chain2 structure]pdbcode],[chain2 code]] substringToIndex:7];
      if (![t_chainid isEqualToString: sname2])
      {
            NSLog(@"sequence2:%@ is not the same as chain2:%@",sname2,t_chainid);
            RELEASE(pool);
            return;
      }

      /* map sequences to residues */
      int start1, start2;
      NSString *realseq = [chain1 get3DSequence];
      NSRange range = [realseq rangeOfString:seq1];
      //printf("found subsequence from %d to %d\n",range.location,range.location+range.length);
      if (range.length != len1)
      {
            NSLog(@"could not find sequence 1:\n%@\nin chain 1:\n%@\n",seq1,realseq);
            RELEASE(pool);
            return;
      }
      start1 = range.location;
      
      realseq = [chain2 get3DSequence];
      range = [realseq rangeOfString:seq2];
      //printf("found subsequence from %d to %d\n",range.location,range.location+range.length);
      if (range.length != len2)
      {
            NSLog(@"could not find sequence 2:\n%@\nin chain 2:\n%@\n",seq2,realseq);
            RELEASE(pool);
            return;
      }
      start2 = range.location;
      
      /* make selection */
      sel1 = [Selection selectionWithChain: chain1];
      sel2 = [Selection selectionWithChain: chain2];
      NSArray *allresidues1 = [[chain1 allResidues] allObjects];
      NSArray *allresidues2 = [[chain2 allResidues] allObjects];
      Residue *residue;
      int nres1=0,nres2=0;
      int ct1,ct2;
      ct1 = [allresidues1 count]; ct2 = [allresidues2 count];
      while (start1+nres1 < ct1 && start2+nres2 < ct2)
      {
            /* skip comments */
            ln = [stream readStringLineLength: 80];
            if (![stream ok])
            {
                  break;
            }
            if ([ln hasPrefix:@"!"])
            {
                  continue;
            }
            if ([ln hasPrefix:@"#"])
            {
                  /* we assume that everything is correct ;-) !!! */
                  continue;
            }
            scanner = [NSScanner scannerWithString:ln];
            [scanner scanInt: &nres1];
            [scanner scanInt: &nres2];
            
            residue = [allresidues1 objectAtIndex:nres1-1+start1];
            //printf ("res1: %@ at %d\n",residue,(nres1-1+start1));

            while (start1+nres1 < ct1 && residue && ([[residue oneLetterCode] characterAtIndex:0] != [seq1 characterAtIndex:nres1-1]))
            {
                  printf ("gap in sequence 1");
                  start1++;
                  residue = nil;
                  if (start1+nres1 < ct1)
                  {
                        residue = [allresidues1 objectAtIndex:nres1-1+start1];
                  }
            }
            if (residue)
            {
                  [sel1 addResidue:residue];
            } else {
                  NSLog(@"residue %c in sequence 1 at %d not found.",[seq1 characterAtIndex:nres1-1],nres1);
                  RELEASE(pool);
                  return;
            }
            residue = [allresidues2 objectAtIndex:nres2-1+start2];
            //printf ("res2: %@ at %d\n",residue,(nres2-1+start2));
            while (start2+nres2 < ct2 && residue && ([[residue oneLetterCode] characterAtIndex:0] != [seq2 characterAtIndex:nres2-1]))
            {
                  start2++;
                  printf ("gap in sequence 2");
                  residue = nil;
                  if (start2+nres2 < ct2)
                  {
                        residue = [allresidues2 objectAtIndex:nres2-1+start2];
                  }
            }
            if (residue)
            {
                  [sel2 addResidue:residue];
            } else {
                  NSLog(@"residue (%c) in sequence 2 at %d not found.",[seq2 characterAtIndex:nres2-1],nres2);
                  RELEASE(pool);
                  return;
            }
      }

      /* superimpose the two chains given the two selections */
      if (transformation)
      {
            RELEASE(transformation);
      }
      transformation = RETAIN([sel2 alignTo: sel1]);
      
      //printf ("strxal: have RT (%@)\n",transformation);
      
      RELEASE(pool);
}


/*
 *   given a superpositions of the @class(Chain), derive the structural alignment using dynamic programming
 */
-(void)deriveStructuralAlignment
{
      CREATE_AUTORELEASE_POOL(pool);

#ifdef DEBUG_COMPUTING_TIME
      clock_t timebase1,timebase2;
#endif

      int *scorematrix;
      int *hinsert, *vinsert;
      int *tbmatrix;
      NSMutableArray *seq1;
      NSMutableArray *seq2;
      Residue *here, *there;
      float dist;
      int col,row,i,j;
      int h1=0,h2=0,h3=0;
      int ttval,tval;

#ifdef DEBUG_COMPUTING_TIME
      timebase1 = clock ();
#endif      
      if (positions)
      {
            [positions removeAllObjects]; 
            RELEASE(positions); // get rid of old AlPos(itions)
      }
      seq1 = [NSMutableArray arrayWithCapacity:[chain1 countStandardAminoAcids]];
      seq2 = [NSMutableArray arrayWithCapacity:[chain2 countStandardAminoAcids]];   
      
      /* prepare sequence 1 and 2 */
      NSEnumerator *residues1 = [chain1 allResidues];
      while (residues1 && (here = [residues1 nextObject]))
      {
            if ([here isStandardAminoAcid])
            {
                  [seq1 addObject: here];
            }
      }
      NSEnumerator *residues2 = [chain2 allResidues];
      while (residues2 && (there = [residues2 nextObject]))
      {
            if ([there isStandardAminoAcid])
            {
                  [seq2 addObject: there];
            }
      }
      
      int score;
      int len1 = [seq1 count];
      int len2 = [seq2 count];
      
      scorematrix = (int*)calloc((len1+1)*(len2+1),sizeof(int));
      vinsert = (int*)calloc((len1+1)*(len2+1),sizeof(int));
      hinsert = (int*)calloc((len1+1)*(len2+1),sizeof(int));
      tbmatrix = (int*)calloc((len1+1)*(len2+1),sizeof(int));     
      
      if (! (scorematrix && vinsert && hinsert && tbmatrix))
      {
            NSLog(@"failed to allocate scoring matrix.");
            RELEASE(pool);
            return;
      }
      
#ifdef DEBUG_COMPUTING_TIME
      timebase2 = clock ();
      printf("  time spent in allocation: %1.1f ms\n",((timebase2-timebase1)*1000.0f/CLOCKS_PER_SEC));
      timebase1 = timebase2;
#endif
      
      i=0; j=0; // will hold row/col of highest value in matrix
      score=0; // maximum score
      int dir; // direction of transition: -1==down, 1==right, 2==diagonal, 0==end of alignment
      for (col=0; col<len1; col++)
      {
            here = [seq1 objectAtIndex: col];
            for (row=0; row<len2; row++)
            {
                  there = [seq2 objectAtIndex: row];
                  dist = (float)[here distanceCATo: there];

                  h1=0;h2=0;h3=0;
                  if (dist <= MAXDIST)
                  {
                        tval = (int)floor((MAXDIST-dist)*1.0+0.5);
                  } else {
                        tval = (int)floor((MAXDIST-dist)*6.0+0.5);
                  }
                  //printf("%@ - %@ d:%2.2f h1=%2.2f\n",here,there,dist,h1);
                  h1 = scorematrix + tval; // diagonal element
                  h2 = hinsert + tval; // end of gap horizontal
                  h3 = vinsert + tval; // end of gap vertical
                  tval = 0;
                  dir = 0;
                  if (h3>tval)
                  {
                        tval = h3;
                        dir = -1; // vertical
                  }
                  if (h2>tval)
                  {
                        tval = h2;
                        dir = 1; // horizontal
                  }
                  if (h1>=tval)     // overwrite in case of same value
                  {
                        tval = h1;
                        dir = 2;
                  }
                  scorematrix = tval; // save maximum value in matrix
                  tbmatrix = dir; // save traceback direction
                  /* save maximum score */
                  if (tval>score)
                  {
                        i=row; j=col;
                        score=tval;
                  }
                  
                  /* update insertion matrices */
                  tval = scorematrix - GAPOPENINGPENALTY;
                  ttval = hinsert - GAPEXTENDPENALTY;
                  hinsert = (tval>ttval?tval:ttval);
                  tval = scorematrix - GAPOPENINGPENALTY;
                  ttval = vinsert - GAPEXTENDPENALTY;
                  vinsert = (tval>ttval?tval:ttval);
                  
            } /* all rows */
      } /* all columns */
      row=i; col=j; // position on maximum end of alignment
#ifdef DEBUG
      printf ("maximum value:%d in row:%d col:%d\n",score,row,col);
#endif

#ifdef DEBUG_COMPUTING_TIME
      timebase2 = clock ();
      printf("  time spent in matrix calc: %1.1f ms\n",((timebase2-timebase1)*1000.0f/CLOCKS_PER_SEC));
      timebase1 = timebase2;
#endif      
      
      /* write score matrix to file */
/* * * * * * * * * * * * * * * * * *
      //printf ("writing scores\n");

      FILE *outfile = fopen("t_scores.csv","w");
      if (outfile)
      {
            fprintf(outfile,". ");
            for (i=0; i<len1; i++)
            {
                  fprintf(outfile,"   %c  ",[[[seq1 objectAtIndex:i] oneLetterCode] characterAtIndex:0]);
            }
            fprintf(outfile,"\n");
            for (j=0; j<len2; j++)
            {
                  fprintf(outfile,"%c ",[[[seq2 objectAtIndex:j] oneLetterCode] characterAtIndex:0]);
                  for (i=0; i<len1; i++)
                  {
                        fprintf(outfile,"% 3d ",scorematrix[(j+1)*len1+(i+1)]);
                  }
                  fprintf(outfile,"\n");
            }
            fclose(outfile);
      }
      outfile = fopen("t_tb.csv","w");
      if (outfile)
      {
            fprintf(outfile,". ");
            for (i=0; i<len1; i++)
            {
                  fprintf(outfile,"   %c  ",[[[seq1 objectAtIndex:i] oneLetterCode] characterAtIndex:0]);
            }
            fprintf(outfile,"\n");
            for (j=0; j<len2; j++)
            {
                  fprintf(outfile,"%c ",[[[seq2 objectAtIndex:j] oneLetterCode] characterAtIndex:0]);
                  for (i=0; i<len1; i++)
                  {
                        fprintf(outfile,"% 3d ",tbmatrix[(j+1)*len1+(i+1)]);
                  }
                  fprintf(outfile,"\n");
            }
            fclose(outfile);
      }
* * * * * * * * * * * * * * * * * */
            
      /* trace back */
      positions = RETAIN([NSMutableArray arrayWithCapacity:len1]);
      i=0; // counter
      AlPos *alpos;
      tval = scorematrix;

      while (tval > 0 && col >= 0 && row >= 0)
      {
            dir = tbmatrix; // direction of traceback
            //printf("col=%d row=%d dir=%d tval=%d\n",col,row,dir,tval);
            /* follow maximal transition */
            if (dir==-1) {
                  alpos = [AlPos alposWithRes1:nil res2:[seq2 objectAtIndex:row]];
                  [alpos designify];
                  row--;
            } else if (dir==1) {
                  alpos = [AlPos alposWithRes1:[seq1 objectAtIndex:col] res2:nil];
                  [alpos designify];
                  col--;
            } else if (dir==2) {
                  alpos = [AlPos alposWithRes1:[seq1 objectAtIndex:col] res2:[seq2 objectAtIndex:row]];
                  col--;row--;
            } else {
                  alpos = nil;
                  printf("wrong traceback. row=%d col=%d\n",row,col);
            }
            if (alpos)
            {
                  [positions addObject: alpos];
            }
            tval = scorematrix;
      }
      
#ifdef DEBUG
      i=[positions count];    
      for (j=1; j<=i; j++)
      {
            alpos = [positions objectAtIndex: (i-j)];
            if ([alpos distance] > SIGNIFDISTTHRESHOLD)
            {
                  [alpos designify];
            }
            NSLog(@"%03d %@\n",j,alpos);
      }
#endif 
      
#ifdef DEBUG_COMPUTING_TIME
      timebase2 = clock ();
      printf("  time spent in traceback: %1.1f ms\n",((timebase2-timebase1)*1000.0f/CLOCKS_PER_SEC));
      timebase1 = timebase2;
#endif
      
      calculated = YES;
      
      free (hinsert);
      free (vinsert);
      free (scorematrix);
      free (tbmatrix);
      
      RELEASE(pool);

#ifdef DEBUG_COMPUTING_TIME
      timebase2 = clock ();
      printf("  time spent in cleanup: %1.1f ms\n",((timebase2-timebase1)*1000.0f/CLOCKS_PER_SEC));
      timebase1 = timebase2;
#endif      
}


+(PairwiseStrxAlignment*)alignmentBetween:(Chain*)p_chain1 andChain:(Chain*)p_chain2
{
      if (p_chain1 == nil || p_chain2 == nil)
      {
            return nil;
      }
      PairwiseStrxAlignment *res = [PairwiseStrxAlignment new];
      res->chain1 = RETAIN(p_chain1);
      res->chain2 = RETAIN(p_chain2);
      
      return AUTORELEASE(res);
}


@end


Generated by  Doxygen 1.6.0   Back to index