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

Chain.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 "Chain.oh"
#include "privateChain.oh"
#include "Structure.oh"
#include "Residue.oh"
#include "privateResidue.oh"
#include "Coordinates.oh"
#include "Matrix.oh"
#include "Matrix53.oh"
#include "Selection.oh"


/* max range of coordinates that are mapped to the hash */
#define MAX_RANGE 640.0


@implementation Chain


-(id)init    // @nodoc
{
      [super init];
      residues = RETAIN([NSMutableArray arrayWithCapacity: 100]);
      solvent = [NSMutableDictionary new] ;
      heterogens = [NSMutableDictionary new] ;
      residueKeys = [NSMutableDictionary new] ;
      residuehash = nil;
      return self;
}


-(void)dealloc    // @nodoc
{
      //printf("Chain_dealloc '%s'\n",[[self description]cString]);
      if (seqres)
      {     
            RELEASE(seqres);
      }
      if (residuehash)
      {
            [residuehash removeAllObjects];
            RELEASE(residuehash);
      }
      if (residueKeys)
      {
            [residueKeys removeAllObjects];
            RELEASE(residueKeys);
      }
      if (solvent)
      {
            NSEnumerator *e_res = [solvent objectEnumerator];
            id res;
            while ((res = [e_res nextObject]))
            {
                  [res setChain:nil];
            }
            [solvent removeAllObjects];
            RELEASE(solvent);
      }
      if (heterogens)
      {
            NSEnumerator *e_res = [heterogens objectEnumerator];
            id res;
            while ((res = [e_res nextObject]))
            {
                  [res setChain:nil];
            }
            [heterogens removeAllObjects];
            RELEASE(heterogens);
      }
      if (residues)
      {
            NSEnumerator *e_res = [residues objectEnumerator];
            id res;
            while ((res = [e_res nextObject]))
            {
                  [res setChain:nil];
            }
            [residues removeAllObjects];
            RELEASE(residues);
      }
      /* release fields */
      if (source)
      {
            RELEASE(source);
      }
      if (compound)
      {
            RELEASE(compound);
      }
      if (eccode)
      {
            RELEASE(eccode);
      }
      [super dealloc];
}


/*
 *   apply transformation matrix to this chain
 */
-(id)transformBy:(Matrix53*)m
{
      Residue *res;
      //printf("chain-transformBy %@\n",self);
      NSEnumerator *e_res = [solvent objectEnumerator];
      while ((res = [e_res nextObject]))
      {
            [res transformBy: m];
      }
      e_res = [heterogens objectEnumerator];
      while ((res = [e_res nextObject]))
      {
            [res transformBy: m];
      }
      e_res = [residues objectEnumerator];
      while ((res = [e_res nextObject]))
      {
            [res transformBy: m];
      }

      return self;
}


/*
 *   return character code of this chain
 */
-(char)code
{
      return code;
}


/*
 *   return numerical character code of this chain
 */
-(NSNumber*)codeNumber
{
      return [NSNumber numberWithChar:code];
}


/*
 *   return code of this chain as string
 */
-(NSString*)name
{
      char buffer;
      buffer=code;
      buffer='\0';
      return [NSString stringWithCString: buffer];
}


/* 
 *   see @method(Chain,-name)
 */
-(NSString*)description
{
      return [self name];
}


/*
 *   return source organism
 */
-(NSString*)source
{
      return source;
}


/*
 *   return compound name
 */
-(NSString*)compound
{
      return compound;
}


/*
 *   return E.C. code (if possible)
 */
-(NSString*)eccode
{
      return eccode;
}


/*
 *   give access to the @Structure we are part of
 */
-(Structure*)structure
{
      return strx;
}


/*
 *   return concatenation of the structure's code and this one's code
 */
-(NSString*)fullPDBCode
{
      if (strx)
      {
            return [NSString stringWithFormat:@"%@%c", [strx pdbcode], code];
      } else {
            return [NSString stringWithFormat:@"0000%c", code];
      }
}


/*
 *   return enumerator over all residues
 */
-(NSEnumerator*)allResidues
{
      return [residues objectEnumerator];
}


/*
 *   return enumerator over all hetero groups
 */
-(NSEnumerator*)allHeterogens
{
      //return [[heterogens allKeys] reverseObjectEnumerator];
      return [heterogens objectEnumerator];
}


/*
 *   return enumerator over all solvent
 */
-(NSEnumerator*)allSolvent
{
      //return [[solvent allKeys] reverseObjectEnumerator];
      return [solvent objectEnumerator];
}


/*
 *   count residues
 */
-(int)countResidues
{
      if (!residues)
      {
            return 0;
      }
      return [residues count];
}


/*
 *   count residues which are of the 20 standard amino acids
 */
-(int)countStandardAminoAcids
{
      NSEnumerator *resenum = [self allResidues];
      Residue *residue;
      int count=0;
      while ((residue = [resenum nextObject]))
      {
            if ([residue isStandardAminoAcid])
            {
                  count++;
            }
      }
      return count;
}


/*
 *   count hetero groups
 */
-(int)countHeterogens
{
      return [heterogens count];
}


/*
 *   count solvent
 */
-(int)countSolvent
{
      return [solvent count];
}


/*
 *   find residue for key
 */
-(Residue*)getResidue:(NSString*)nr
{
      return [residueKeys objectForKey: nr];
}


/*
 *   find hetero for key
 */
-(Residue*)getHeterogen:(NSString*)nr
{
      return [heterogens objectForKey: nr];
}


/*
 *   find solvent for key
 */
-(Residue*)getSolvent:(NSString*)nr
{
      return [solvent objectForKey: nr];
}


/*
 *   add residue 
 */
-(id)addResidue:(Residue*)res
{
      [residues addObject: res];
      [residueKeys setObject: res forKey: [res key]];
      [res setChain:self];
      return self;
}


/*
 *   add hetero group
 */
-(id)addHeterogen:(Residue*)het
{
      [heterogens setObject: het forKey: [het key]];
      [het setChain:self];
      return self;
}


/*
 *   add solvent
 */
-(id)addSolvent:(Residue*)sol
{
      [solvent setObject: sol forKey: [sol key]];
      [sol setChain:self];
      return self;
}


/*
 *   remove residue (either real residue, or hetero group, or solvent)
 */
-(void)removeResidue:(Residue*)p_res
{
      NSLog(@"Chain_remove residue: %@ (%@)\n",p_res, [p_res class]);
      NSString *p_num = [p_res key];
      if ([heterogens objectForKey:p_num])
      {
            [heterogens removeObjectForKey:p_num];
            return;
      }
      if ([solvent objectForKey:p_num])
      {
            [solvent removeObjectForKey:p_num];
            return;
      }
      NSArray *keys = [residueKeys allKeysForObject:p_res];
      if (keys && ([keys count] > 0))
      {
            [residueKeys removeObjectsForKeys:keys];
      }
      [residues removeObject:p_res];
}

  
/*
 *   return parsed SEQRES sequence
 */
-(NSString*)getSEQRES
{
      return seqres;
}


/*
 *   derive amino acid sequence from coordinates. 
 *   wherever there is a gap in the mainchain and in the number of the residues,
 *   insert as many neutral elements 'X' as necessary.
 */
-(NSString*)getSequence
{
      NSString *res = nil;
      char buffer;
      /* extract FASTA sequence */
      NSEnumerator *allresidues = [self allResidues];
      Residue *residue = nil;
      Residue *lastres = nil;
      int lastnumber=0;
      int diff;
      int seqpos = 0;
      float distanceCA;
      while (allresidues && (residue = [allresidues nextObject]))
      {
            diff = [[residue number] intValue] - lastnumber - 1;
            if (diff<0)
            {
                  diff = 0;
            }
            if (lastnumber>0 && diff>0)
            {
                  if (lastres)
                  {
                        distanceCA = [residue distanceCATo:lastres];
                  } else {
                        distanceCA = 0.0f;
                  }
                  //printf("GAP last=%d diff=%d distance=%1.1f\n",lastnumber,diff,distanceCA);
                  if (distanceCA > 4.5f)
                  {
                        // limit to 99 residues in insertion
                        if (diff>99)
                        {
                              diff = 99;
                        }
                        while (diff>0)
                        {
                              /* write missing residues */
                              buffer= 'X'; seqpos++; diff--;
                        }
                  }
            }
            buffer=(char)[[residue oneLetterCode] characterAtIndex:0];
            [residue setSeqNum:(seqpos+1)];
            lastnumber = [[residue number] intValue];
            lastres = residue;
            seqpos++;
      }
      if (seqpos>0)
      {
            buffer='\0';
            res = [NSString stringWithCString: buffer];
      }
      return res;
}


/*
 *   derive sequence of present standard amino acids
 */
-(NSString*)get3DSequence
{
      NSString *res = nil;
      char buffer;
      /* extract FASTA sequence */
      NSEnumerator *allresidues = [self allResidues];
      Residue *residue = nil;
      int seqpos = 0;
      while (allresidues && (residue = [allresidues nextObject]))
      {
            if ([residue isStandardAminoAcid])
            {
                  buffer=(char)[[residue oneLetterCode] characterAtIndex:0];
                  [residue setSeqNum:(seqpos+1)];
                  seqpos++;
            }
      }
      if (seqpos>0)
      {
            buffer='\0';
            res = [NSString stringWithCString: buffer];
      }
      return res;
}


/*
 *   select all residues in both chains which are at most /maxdist/ separated, return in selections
 */
-(int)selectResiduesCloseTo:(Chain*)other maxDistance:(float)maxdist sel1:(Selection**)p_sel1 sel2:(Selection**)p_sel2
{
      *p_sel1 = [Selection selectionWithChain: self];
      *p_sel2 = [Selection selectionWithChain: other];
      Residue *here, *there;
      Residue *minthere;
      float dist, mindist;
      NSEnumerator *residues1 = [self allResidues];
      while (residues1 && (here = [residues1 nextObject]))
      {
            mindist = 2*maxdist;
            minthere = nil;
            NSEnumerator *residues2 = [other allResidues];
            while (residues2 && (there = [residues2 nextObject]))
            {
                  dist = [here distanceCATo: there];
                  if (dist>=0.0 && dist <= mindist)
                  {
                        mindist = dist;
                        minthere = there;
                  }
            }
            if (mindist<=maxdist && minthere)
            {
                  [*p_sel1 addResidue: here];
                  [*p_sel2 addResidue: minthere];
            }
      }
      return [*p_sel1 count];
}


/*
 *   rather used for testing 
 */
-(void)writeDistanceMapTo:(Chain*)other toFile:(NSString*)fpath
{
      Residue *here, *there;
      float dist;
      FILE *mapfile = fopen([fpath cString],"w");
      if (!mapfile)
      {
            NSLog(@"writeDistanceMapTo:toFile: failed to open file: %@",fpath);
            return;
      }
      NSEnumerator *residues1 = [self allResidues];
      while (residues1 && (here = [residues1 nextObject]))
      {
            NSEnumerator *residues2 = [other allResidues];
            while (residues2 && (there = [residues2 nextObject]))
            {
                  dist = [here distanceCATo: there];
                  fprintf(mapfile,"%2f ",dist);
            }
            fprintf(mapfile,"\n");
      }
      fclose(mapfile);
}


/*
 *   computes the hash of all residues
 */
-(void)prepareResidueHash:(float)gridsize
{
      if (residuehash)
      {
            [residuehash removeAllObjects];
            RELEASE(residuehash);
      }
      hashingbits = (unsigned int)(log(MAX_RANGE / gridsize) / log(2.0) + 0.5);
      //printf("input: %1.1f  hashingbits: %d\n",gridsize,hashingbits);
      hash_value_offset = (MAX_RANGE / (double)(1UL << hashingbits))/2.0;
      //printf("resolution offset: %1.1f\n",hash_value_offset);

      if (hashingbits > 10 || hashingbits < 2)
      {
            NSLog(@"Chain_prepareResidueHash: gridsize is either too big or too small! Cannot compute residue hash.");
            return;
      }
      residuehash = RETAIN([NSMutableDictionary new]);
      Residue *t_res;
      Atom *t_atm;
      NSEnumerator *allres = [self allResidues];
      NSEnumerator *allatm;
      //NSNumber *hashvalue;
      //NSMutableArray *t_arr;
      while ((t_res = [allres nextObject]))
      {
            allatm = [t_res allAtoms];
            while ((t_atm = [allatm nextObject]))
            {
                  [self enterHashAtom:(Coordinates*)t_atm for:t_res];
            }
      }
      allres = [self allHeterogens];
      while ((t_res = [allres nextObject]))
      {
            allatm = [t_res allAtoms];
            while ((t_atm = [allatm nextObject]))
            {
                  [self enterHashAtom:(Coordinates*)t_atm for:t_res];
            }
      }
}


/*
 *   returns a list of residues which are close to the given coordinates<br>
 *   must compute the hash table first using: @method(Chain, -prepareResidueHash:)
 */
-(NSArray*)findResiduesCloseTo:(Coordinates*)p_coords
{
      NSNumber *hashvalue = [self mkCoordinatesHashX:[p_coords x] Y:[p_coords y] Z:[p_coords z]];
      NSArray *t_arr = [residuehash objectForKey: hashvalue];
      return t_arr;
}



-(NSNumber*)mkCoordinatesHashX:(double)p_x Y:(double)p_y Z:(double)p_z
{

// input:-320          0       +320
//          0         320       640
//          |----------|---------|
// output:  0          32        64 (6 bits) // (HASH_GRID_SIZE=10.0)
// output:  0          64       128 (7 bits) // (HASH_GRID_SIZE=5.0)
// output:  0         128       256 (8 bits) // (HASH_GRID_SIZE=2.5)
// output:  0         256       512 (9 bits) // (HASH_GRID_SIZE=1.25)
// output:  0         512      1024 (10 bits) // (HASH_GRID_SIZE=0.625)

      double x,y,z;
      unsigned long hashv;
      unsigned long mask = (1UL << hashingbits) - 1;
      double factor = (double)(1UL << hashingbits);
      x = (p_x+(MAX_RANGE/2.0)) / MAX_RANGE * factor;
      if (x<0.0) { x=0.0; printf("clipping -x\n"); }
      if (x>factor) { x=factor; printf("clipping +x\n"); }
      y = (p_y+(MAX_RANGE/2.0)) / MAX_RANGE * factor;
      if (y<0.0) { y=0.0; printf("clipping -y\n"); }
      if (y>factor) { y=factor; printf("clipping +y\n"); }
      z = (p_z+(MAX_RANGE/2.0)) / MAX_RANGE * factor;
      if (z<0.0) { z=0.0; printf("clipping -z\n"); }
      if (z>factor) { z=factor; printf("clipping +z\n"); }
      hashv = (unsigned long)(x+0.5) & mask;
      hashv |= ((unsigned long)(y+0.5) & mask)<<hashingbits;
      hashv |= (((unsigned long)(z+0.5) & mask)<<hashingbits)<<hashingbits;
      // hashv is now a 3 times HASHING_BITS number encoding coordinates
      //printf("mask:%lx factor:%1.1f hashv:%lx x:%1.1f y:%1.1f z:%1.1f -> x:%1.1f y:%1.1f z:%1.1f\n",mask,factor,hashv,p_x,p_y,p_z,x,y,z);
      return [NSNumber numberWithUnsignedLong:hashv];
}

@end


Generated by  Doxygen 1.6.0   Back to index