← Back to branch summary

~ubuntu-branches/ubuntu/edgy/ncbi-tools6/edgy

« back to all changes in this revision

Viewing changes to api/seqport.c

  • Committer: Bazaar Package Importer
  • Author(s): Barry deFreese
  • Date: 2006-07-19 23:28:07 UTC
  • mfrom: (1.1.5 upstream)
  • Revision ID: james.westby@ubuntu.com-20060719232807-et3cdmcjgmnyleyx
Tags: 6.1.20060507-3ubuntu1
Re-merge with Debian

Show diffs side-by-side

added added

removed removed

Lines of Context:
api/seqport.c
29
29
*   
30
30
* Version Creation Date: 7/13/91
31
31
*
32
 
* $Revision: 6.144 $
 
32
* $Revision: 6.150 $
33
33
*
34
34
* File Description:  Ports onto Bioseqs
35
35
*
39
39
* -------  ----------  -----------------------------------------------------
40
40
*
41
41
* $Log: seqport.c,v $
 
42
* Revision 6.150  2006/03/22 15:31:32  kans
 
43
* SeqPortStreamSeqLoc gives unique message when bailing on gi 0 as opposed to failure after trying to load
 
44
*
 
45
* Revision 6.149  2006/03/07 21:34:28  kans
 
46
* checks for gi 0 now also check for negative value
 
47
*
 
48
* Revision 6.148  2006/03/07 20:02:01  kans
 
49
* SeqPortStreamSeqLoc immediately treats gi 0 as an error
 
50
*
 
51
* Revision 6.147  2006/01/23 13:01:41  bollin
 
52
* when converting sequences from raw to delta, adjust any alignments that the
 
53
* sequence may be part of.
 
54
*
 
55
* Revision 6.146  2005/12/16 20:19:56  bollin
 
56
* only allow reverse for delta sequences when the delta sequence has no far locations
 
57
*
 
58
* Revision 6.145  2005/12/15 19:45:24  bollin
 
59
* added functions to reverse and complement delta sequences
 
60
*
42
61
* Revision 6.144  2005/08/24 15:14:31  kans
43
62
* modified MolWtForLoc to use StreamCache, added MolWtForBsp and MolWtForStr
44
63
*
597
616
#include <subutil.h>
598
617
#include <tofasta.h>   /* for FastaSeqLineEx function */
599
618
#include <salutil.h> 
 
619
#include <alignmgr2.h> /* for correcting alignments when converting to delta */
600
620
 
601
621
 
602
622
NLM_EXTERN Boolean LIBCALL SeqPortAdjustLength (SeqPortPtr spp);
2967
2987
  sip = SeqLocId (slp);
2968
2988
  if (sip == NULL) return 0;
2969
2989
 
 
2990
  if (sip->choice == SEQID_GI && sip->data.intvalue <= 0) {
 
2991
 
 
2992
    /* gi 0 or negative is always a data error, just report and bail */
 
2993
 
 
2994
    SeqIdWrite (sip, buf, PRINTID_FASTA_SHORT, sizeof (buf) - 1);
 
2995
    if (parentID != NULL) {
 
2996
      SeqIdWrite (parentID, pid, PRINTID_FASTA_LONG, sizeof (pid) - 1);
 
2997
      ErrPostEx (SEV_ERROR, 0, 0, "SeqPortStream ignoring Bioseq %s component of %s", buf, pid);
 
2998
    } else {
 
2999
      ErrPostEx (SEV_ERROR, 0, 0, "SeqPortStream ignoring Bioseq %s", buf);
 
3000
    }
 
3001
    sdp->failed = TRUE;
 
3002
    return 0;
 
3003
  }
 
3004
 
2970
3005
  bsp = BioseqLockById (sip);
2971
3006
 
2972
3007
#ifdef OS_UNIX
4989
5024
        return retval;
4990
5025
}
4991
5026
 
4992
 
/*-------------- BioseqComplement () ---------------------------*/
4993
 
/***********************************************************************
4994
 
*   BioseqComplement:   Takes the nucleic acid sequence from Bioseq
4995
 
*       Entry and gives the complement sequence in place
4996
 
*       Does not change features.
4997
 
************************************************************************/
4998
 
NLM_EXTERN Boolean LIBCALL BioseqComplement (BioseqPtr bsp)
 
5027
static Boolean ComplementSeqData (Uint1 seqtype, Int4 seqlen, ByteStorePtr bysp)
4999
5028
{
5000
5029
        SeqCodeTablePtr sctp;
5001
 
        ByteStorePtr    bysp;
5002
 
        long            readbyte, bslen;
5003
 
        Int4            seqlen;
5004
 
        Uint1           seqtype, byte = 0, byte_to, newbyte = 0, residue;
 
5030
        long                    readbyte, bslen;
 
5031
        Uint1           byte = 0, byte_to, newbyte = 0, residue;
5005
5032
        Uint1           comp, bitctr, mask, lshift, rshift, bc;
5006
5033
        
5007
 
        if (bsp == NULL)
5008
 
        {
5009
 
                ErrPostEx(SEV_ERROR,0,0, "Error: not a BioseqPtr\n");
5010
 
                return FALSE;  
5011
 
        }
5012
 
                        
5013
 
        if (bsp->repr != Seq_repr_raw)
5014
 
        {
5015
 
                ErrPostEx(SEV_ERROR,0,0, "Error: not a raw sequence\n");
5016
 
                return FALSE;  
5017
 
        }
5018
 
                        
5019
 
        if (bsp->seq_data == NULL)
 
5034
        if (bysp == NULL)
5020
5035
        {
5021
5036
                ErrPostEx(SEV_ERROR,0,0, "Error:  no sequence data\n");
5022
 
                return FALSE;
 
5037
                return FALSE;     
5023
5038
        }
5024
5039
 
5025
 
        seqtype = bsp->seq_data_type;
5026
 
        if ( ISA_aa(bsp->mol)) {
5027
 
                ErrPostEx(SEV_ERROR,0,0, "Error:  cannot complement aa\n");
5028
 
                return FALSE;
5029
 
        }
5030
 
 
5031
5040
        if ((sctp = SeqCodeTableFind (seqtype)) == NULL)
5032
5041
        {
5033
5042
                ErrPostEx(SEV_ERROR,0,0, "Can't open table\n");
5056
5065
                        lshift = 0;
5057
5066
                        mask = 255;
5058
5067
                        break;
5059
 
                case Seq_code_iupacaa:
5060
 
                case Seq_code_ncbi8aa:
5061
 
                case Seq_code_ncbieaa:
5062
 
                case Seq_code_ncbipaa:
5063
 
                case Seq_code_iupacaa3:
5064
 
                case Seq_code_ncbistdaa:                        /* ignore amino acid */
5065
 
                    ErrPostEx(SEV_ERROR,0,0, "Error:  cannot complement aa ; No ->mol flag on Bioseq\n");
5066
 
                    return FALSE;
5067
 
                case Seq_code_ncbipna:
5068
 
                    ErrPostEx(SEV_WARNING,0,0, "Error: Don't yet know how to complement profile\n");
 
5068
    case Seq_code_iupacaa:
 
5069
    case Seq_code_ncbi8aa:
 
5070
    case Seq_code_ncbieaa:
 
5071
    case Seq_code_ncbipaa:
 
5072
    case Seq_code_iupacaa3:
 
5073
    case Seq_code_ncbistdaa:                    /* ignore amino acid */
 
5074
      ErrPostEx(SEV_ERROR,0,0, "Error:  cannot complement aa ; No ->mol flag on Bioseq\n");
 
5075
      return FALSE;
 
5076
    case Seq_code_ncbipna:
 
5077
      ErrPostEx(SEV_WARNING,0,0, "Error: Don't yet know how to complement profile\n");
 
5078
                        return FALSE;
5069
5079
                default:
5070
5080
                        return FALSE;
5071
5081
        }
5072
5082
 
5073
 
        seqlen = bsp->length;
5074
 
        bysp = bsp->seq_data;
5075
5083
        bslen = BSLen(bysp);
5076
5084
        bitctr = 0;
5077
5085
        readbyte = 0;
5110
5118
                }
5111
5119
        }
5112
5120
        return TRUE;
 
5121
  
 
5122
}
 
5123
 
 
5124
 
 
5125
static Boolean DeltaBioseqComplement (BioseqPtr bsp)
 
5126
{
 
5127
  DeltaSeqPtr dsp;
 
5128
  SeqLitPtr   slip;
 
5129
  Boolean     rval = FALSE;
 
5130
  
 
5131
  if (bsp == NULL || bsp->repr != Seq_repr_delta)
 
5132
  {
 
5133
    return rval;
 
5134
  }
 
5135
  
 
5136
  dsp = (DeltaSeqPtr) bsp->seq_ext;
 
5137
  while (dsp != NULL)
 
5138
  {
 
5139
    if (dsp->choice != 2) 
 
5140
    {
 
5141
      ErrPostEx(SEV_ERROR,0,0, "Error: Can't complement delta sequences with far locs\n");
 
5142
      return FALSE;  
 
5143
    }
 
5144
    dsp = dsp->next;
 
5145
  }
 
5146
  rval = TRUE;
 
5147
  dsp = (DeltaSeqPtr) bsp->seq_ext;
 
5148
  while (dsp != NULL)
 
5149
  {
 
5150
    slip = (SeqLitPtr) dsp->data.ptrvalue;
 
5151
    /* complement data */
 
5152
    if (slip->seq_data != NULL)
 
5153
    {
 
5154
      rval &= ComplementSeqData (slip->seq_data_type, slip->length, slip->seq_data);
 
5155
    }
 
5156
    dsp = dsp->next;
 
5157
  }
 
5158
  return rval;
 
5159
}
 
5160
 
 
5161
 
 
5162
/*-------------- BioseqComplement () ---------------------------*/
 
5163
/***********************************************************************
 
5164
*   BioseqComplement:   Takes the nucleic acid sequence from Bioseq
 
5165
*       Entry and gives the complement sequence in place
 
5166
*       Does not change features.
 
5167
************************************************************************/
 
5168
NLM_EXTERN Boolean LIBCALL BioseqComplement (BioseqPtr bsp)
 
5169
{
 
5170
        Boolean         rval = FALSE;
 
5171
        
 
5172
  if (bsp == NULL)
 
5173
  {
 
5174
    ErrPostEx(SEV_ERROR,0,0, "Error: not a BioseqPtr\n");
 
5175
    rval = FALSE;  
 
5176
  }
 
5177
  else if (ISA_aa(bsp->mol)) 
 
5178
  {
 
5179
    ErrPostEx(SEV_ERROR,0,0, "Error:  cannot complement aa\n");
 
5180
                rval = FALSE;  
 
5181
  }
 
5182
  else if (bsp->repr == Seq_repr_delta)
 
5183
  {
 
5184
    rval = DeltaBioseqComplement (bsp);
 
5185
  }
 
5186
  else if (bsp->repr == Seq_repr_raw)
 
5187
  {
 
5188
    rval = ComplementSeqData (bsp->seq_data_type, bsp->length, bsp->seq_data);
 
5189
  }
 
5190
  else
 
5191
  {
 
5192
    ErrPostEx(SEV_ERROR,0,0, "Error: not a raw or delta sequence\n");
 
5193
                rval = FALSE;  
 
5194
  }
 
5195
  return rval;                        
5113
5196
 
5114
5197
} /* BioseqComplement */
5115
5198
 
5116
 
           
5117
 
/*-------------- BioseqReverse () ---------------------------*/
5118
 
/***********************************************************************
5119
 
*   BioseqReverse:   Takes nucleic acid sequence from Bioseq Entry and 
5120
 
*       reverses the whole sequence in place
5121
 
*       Does not change features.
5122
 
************************************************************************/
5123
 
NLM_EXTERN Boolean LIBCALL BioseqReverse (BioseqPtr bsp)
 
5199
 
 
5200
static Boolean LIBCALL ReverseSeqData (Uint1 seqtype, Int4 seqlen, ByteStorePtr bysp1)
5124
5201
{
5125
 
        ByteStorePtr    bysp1 = '\0';
5126
5202
        ByteStorePtr    bysp2 = '\0';
5127
5203
        long            readbyte, bslen = 0;
5128
 
        Int4            seqlen, count = 0;
5129
 
        Uint1           seqtype, byte = 0, byte2, byte_to = 0, byte_to2, newbyte = 0;
 
5204
        Int4            count = 0;
 
5205
        Uint1   byte = 0, byte2, byte_to = 0, byte_to2, newbyte = 0;
5130
5206
        Uint1           newbyte2, finalbyte, residue, residue2, bitctr, bc2 = 0;
5131
5207
        Uint1           bitctr2, mask, mask2, lshift, rshift, bc = 0, jagged;
5132
5208
        
5133
 
        if (bsp == NULL)
5134
 
        {
5135
 
                ErrPostEx(SEV_ERROR,0,0, "Error: not a BioseqPtr\n");
5136
 
                return FALSE;  
5137
 
        }
5138
 
                        
5139
 
        if (bsp->repr != Seq_repr_raw)
5140
 
        {
5141
 
                ErrPostEx(SEV_ERROR,0,0, "Error: not a raw sequence\n");
5142
 
                return FALSE;  
5143
 
        }
5144
 
                        
5145
 
        if (bsp->seq_data == NULL)
5146
 
        {
5147
 
                ErrPostEx(SEV_ERROR,0,0, "Error:  No sequence data\n");
5148
 
                return FALSE;
5149
 
        }
 
5209
  if (bysp1 == NULL)
 
5210
  {
 
5211
    ErrPostEx(SEV_ERROR,0,0, "Error:  No sequence data\n");
 
5212
    return FALSE;
 
5213
  }
5150
5214
 
5151
 
        seqlen = bsp->length;
5152
 
        seqtype = bsp->seq_data_type;
5153
5215
        switch (seqtype){
5154
5216
                case Seq_code_ncbi2na:          /*bitshifts needed*/
5155
5217
                        mask = 192;
5219
5281
                default:                /*ignores amino acid sequence*/
5220
5282
                        return FALSE;
5221
5283
        }
5222
 
        bysp1 = bsp->seq_data;
5223
5284
        bysp2 = BSDup(bysp1);
5224
5285
        bslen = BSLen (bysp1);
5225
5286
        bitctr = bitctr2 = 0;
5314
5375
        }
5315
5376
        BSFree(bysp2);
5316
5377
        return TRUE;
 
5378
} /* ReverseSeqData */
 
5379
 
 
5380
 
 
5381
static Boolean DeltaBioseqReverse (BioseqPtr bsp)
 
5382
{
 
5383
  DeltaSeqPtr dsp, next_dsp, newchain = NULL;
 
5384
  SeqLitPtr   slip;
 
5385
  Boolean     rval = FALSE;
 
5386
  Boolean     split = FALSE;
 
5387
  
 
5388
  if (bsp == NULL || bsp->repr != Seq_repr_delta)
 
5389
  {
 
5390
    return rval;
 
5391
  }
 
5392
  
 
5393
  dsp = (DeltaSeqPtr) bsp->seq_ext;
 
5394
  while (dsp != NULL)
 
5395
  {
 
5396
    if (dsp->choice != 2) 
 
5397
    {
 
5398
      ErrPostEx(SEV_ERROR,0,0, "Error: Can't reverse delta sequences with far locs\n");
 
5399
      return FALSE;  
 
5400
    }
 
5401
    dsp = dsp->next;
 
5402
  }
 
5403
  
 
5404
  dsp = (DeltaSeqPtr) bsp->seq_ext;  
 
5405
  rval = TRUE;
 
5406
  while (dsp != NULL)
 
5407
  {
 
5408
    slip = (SeqLitPtr) dsp->data.ptrvalue;
 
5409
    /* reverse data */
 
5410
    if (slip->seq_data != NULL)
 
5411
    {
 
5412
      rval &= ReverseSeqData (slip->seq_data_type, slip->length, slip->seq_data);
 
5413
    }
 
5414
    
 
5415
    /* reverse the chain */
 
5416
    next_dsp = dsp->next;
 
5417
    dsp->next = newchain;
 
5418
    newchain = dsp;
 
5419
    
 
5420
    dsp = next_dsp;
 
5421
  }
 
5422
  bsp->seq_ext = newchain;
 
5423
  return rval;
 
5424
}
 
5425
           
 
5426
/*-------------- BioseqReverse () ---------------------------*/
 
5427
/***********************************************************************
 
5428
*   BioseqReverse:   Takes nucleic acid sequence from Bioseq Entry and 
 
5429
*       reverses the whole sequence in place
 
5430
*       Does not change features.
 
5431
************************************************************************/
 
5432
NLM_EXTERN Boolean LIBCALL BioseqReverse (BioseqPtr bsp)
 
5433
{
 
5434
        Boolean       rval;
 
5435
        
 
5436
  if (bsp == NULL)
 
5437
  {
 
5438
    ErrPostEx(SEV_ERROR,0,0, "Error: not a BioseqPtr\n");
 
5439
    rval = FALSE;  
 
5440
  }
 
5441
  else if (bsp->repr == Seq_repr_delta)
 
5442
  {
 
5443
    rval = DeltaBioseqReverse (bsp);
 
5444
  }                       
 
5445
  else if (bsp->repr == Seq_repr_raw)
 
5446
  {
 
5447
    rval = ReverseSeqData (bsp->seq_data_type, bsp->length, bsp->seq_data);
 
5448
  }
 
5449
  else
 
5450
  {
 
5451
    ErrPostEx(SEV_ERROR,0,0, "Error: not a raw or delta sequence\n");
 
5452
    rval = FALSE;  
 
5453
  }
 
5454
                        
 
5455
        return rval;
5317
5456
} /* BioseqReverse */
5318
5457
 
5319
5458
#define SPC_BUFF_CHUNK 1024
7817
7956
}
7818
7957
 
7819
7958
 
7820
 
static void FixGapLength (SeqIdPtr sip, Uint2 moltype, Int4 offset, Int4 diff)
 
7959
static void FixGapLength (BioseqPtr bsp, Int4 offset, Int4 diff)
7821
7960
{
7822
 
  CharPtr    extra_ns;
7823
 
  SeqLocPtr  slp;
7824
 
  
7825
 
  if (sip == NULL || diff == 0) return;
 
7961
  CharPtr     extra_ns;
 
7962
  SeqLocPtr   slp;
 
7963
  ValNodePtr  align_annot_list, vnp;
 
7964
  SeqAnnotPtr sanp;
 
7965
  
 
7966
  if (bsp == NULL || bsp->id == NULL || diff == 0) return;
 
7967
  
 
7968
  align_annot_list = FindAlignSeqAnnotsForBioseq (bsp);
 
7969
 
7826
7970
  if (diff > 0)
7827
7971
  {
7828
7972
    extra_ns = (CharPtr)MemNew ((diff + 1) * sizeof (Char));
7830
7974
    {
7831
7975
      MemSet (extra_ns, 'N', diff);
7832
7976
      extra_ns [diff] = 0;
7833
 
      insertchar (extra_ns, offset, sip, moltype, FALSE);
7834
 
    }
 
7977
      insertchar (extra_ns, offset, bsp->id, bsp->mol, FALSE);
 
7978
    }
 
7979
        slp = SeqLocIntNew (offset, offset + diff - 1, Seq_strand_plus, bsp->id);
 
7980
    for (vnp = align_annot_list; vnp != NULL; vnp = vnp->next)
 
7981
    {
 
7982
      sanp = vnp->data.ptrvalue;
 
7983
      if (sanp != NULL && sanp->type == 2)
 
7984
      {
 
7985
        sanp->data = SeqAlignInsertByLoc (slp, sanp->data);
 
7986
      }
 
7987
    }
 
7988
    SeqLocFree (slp);
7835
7989
  }
7836
7990
  else
7837
7991
  {
7838
 
        slp = SeqLocIntNew (offset, offset - diff - 1, Seq_strand_plus, sip);
 
7992
        slp = SeqLocIntNew (offset, offset - diff - 1, Seq_strand_plus, bsp->id);
7839
7993
    SeqDeleteByLoc (slp, TRUE, FALSE);
 
7994
 
 
7995
    for (vnp = align_annot_list; vnp != NULL; vnp = vnp->next)
 
7996
    {
 
7997
      sanp = vnp->data.ptrvalue;
 
7998
      if (sanp != NULL && sanp->type == 2)
 
7999
      {
 
8000
        sanp->data = SeqAlignDeleteByLoc (slp, sanp->data);
 
8001
      }
 
8002
    }
 
8003
 
7840
8004
    SeqLocFree (slp);
7841
8005
  }
7842
8006
}
7998
8162
          slp->fuzz = ifp;
7999
8163
          if (slp->length != 100)
8000
8164
          {
8001
 
            FixGapLength (bsp->id, bsp->mol, len, 100 - slp->length);
 
8165
            FixGapLength (bsp, len, 100 - slp->length);
8002
8166
                slp->length = 100;
8003
8167
          }
8004
8168
        }