← Back to branch summary

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

« back to all changes in this revision

Viewing changes to algo/blast/api/blast_seq.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:
algo/blast/api/blast_seq.c
1
 
static char const rcsid[] = "$Id: blast_seq.c,v 1.73 2005/11/09 14:49:49 camacho Exp $";
 
1
static char const rcsid[] = "$Id: blast_seq.c,v 1.83 2006/05/04 15:53:12 camacho Exp $";
2
2
/*
3
3
* ===========================================================================
4
4
*
32
32
#include <seqport.h>
33
33
#include <sequtil.h>
34
34
#include <objloc.h>
 
35
#include <readdb.h>
35
36
#include <algo/blast/api/blast_seq.h>
36
37
#include <algo/blast/core/blast_filter.h>
37
38
#include <algo/blast/core/blast_util.h>
38
39
#include <algo/blast/core/blast_encoding.h>
39
40
#include <algo/blast/core/blast_setup.h> /* For BlastSeqLoc_RestrictToInterval */
 
41
#include <algo/blast/core/blast_inline.h>
40
42
 
41
43
/** @addtogroup CToolkitAlgoBlast
42
44
 *
43
45
 * @{
44
46
 */
45
47
 
 
48
 
 
49
/** Structure used for hash-based comparison of sequence IDs */
 
50
typedef struct SeqIdHash {
 
51
    SeqId *id;           /**< The ID of this entry */
 
52
    Int4 query_index;    /**< index of query with this ID */
 
53
    Int4 next_id;        /**< Offset of the next hash entry in a chain */
 
54
} SeqIdHash;
 
55
 
 
56
Boolean
 
57
BlastSeqlocsHaveDuplicateIDs(SeqLoc* query_seqlocs)
 
58
{
 
59
   Boolean retval = FALSE;
 
60
   const Int4 kNumSeqs = ValNodeLen(query_seqlocs);
 
61
   const Int4 kLog2HashSize = 11;
 
62
   SeqIdHash *id_entries;
 
63
   Uint4 *hashtable;
 
64
   Int4 curr_id_num;
 
65
   SeqLocPtr slp;
 
66
 
 
67
   if (kNumSeqs == 1)
 
68
      return FALSE;
 
69
 
 
70
   /* allocate hashtable */
 
71
   hashtable = (Uint4 *)calloc(1 << kLog2HashSize, sizeof(Uint4));
 
72
   id_entries = (SeqIdHash *)malloc((kNumSeqs + 1) * sizeof(SeqIdHash));
 
73
 
 
74
   for (slp = query_seqlocs, curr_id_num = 1; slp; slp = slp->next) {
 
75
 
 
76
       Uint4 hashval;
 
77
       SeqIdPtr id = SeqLocId(slp);
 
78
       Char buffer[64];
 
79
 
 
80
       /* hash the ID of the next query sequence */
 
81
       SeqIdLabel(id, buffer, sizeof(buffer), OM_LABEL_CONTENT);
 
82
       hashval = readdb_sequence_hash(buffer, strlen(buffer));
 
83
       hashval = hashval >> (32 - kLog2HashSize);
 
84
       if (hashtable[hashval] != 0) {
 
85
          Int4 offset = hashtable[hashval];
 
86
          /* check the chain of ID's in the resulting hash
 
87
             entry for a duplicate */
 
88
          while (offset != 0) {
 
89
             SeqIdHash *curr_entry = id_entries + offset;
 
90
             if (SeqIdMatch(id, curr_entry->id)) {
 
91
                 retval = TRUE;
 
92
                 goto clean_up;
 
93
             }
 
94
             offset = curr_entry->next_id;
 
95
          }
 
96
       }
 
97
 
 
98
       /* no duplicate found; add the ID to the hashtable */
 
99
       id_entries[curr_id_num].id = id;
 
100
       id_entries[curr_id_num].next_id = hashtable[hashval];
 
101
       hashtable[hashval] = curr_id_num++;
 
102
   }
 
103
 
 
104
clean_up:
 
105
   sfree(hashtable);
 
106
   sfree(id_entries);
 
107
   return retval;
 
108
}
 
109
 
46
110
/** Converts a SeqLocPtr to a BlastSeqLoc, used for formatting.
47
111
 * @param slp SeqLocPtr to be converted [in]
48
112
 * @param head_loc BlastSeqLoc returned from last call [in]
78
142
    BlastMaskLoc* retval = NULL;
79
143
    Int4 query_index = 0;
80
144
    const unsigned int kNumContexts = BLAST_GetNumberOfContexts(program_number);
81
 
    SeqLocPtr current_query_loc = NULL;
 
145
    const Int4 kLog2HashSize = 11;
 
146
    SeqIdHash *id_entries;
 
147
    Uint4 *hashtable;
 
148
    Int4 curr_id_num = 1;
 
149
    SeqLocPtr query_slp, mask_slp;
82
150
 
83
151
    if (!mask_seqlocs)
84
152
        return NULL;
85
153
 
86
154
    retval = BlastMaskLocNew(kNumSeqs*kNumContexts);
87
155
 
88
 
    for (current_query_loc = query_seqlocs, query_index = 0; 
89
 
         current_query_loc; 
90
 
         current_query_loc = current_query_loc->next, query_index++) {
 
156
    /* create hashtable for query IDs */
 
157
    hashtable = (Uint4 *)calloc(1 << kLog2HashSize, sizeof(Uint4));
 
158
    id_entries = (SeqIdHash *)malloc((kNumSeqs + 1) * sizeof(SeqIdHash));
 
159
 
 
160
    /* add the ID of each query sequence to the hashtable */
 
161
    for (query_slp = query_seqlocs; query_slp; query_slp = query_slp->next) {
 
162
        Uint4 hashval;
 
163
        SeqIdPtr seq_id = SeqLocId(query_slp);
 
164
        Char buffer[64];
 
165
 
 
166
        SeqIdLabel(seq_id, buffer, sizeof(buffer), OM_LABEL_CONTENT);
 
167
        hashval = readdb_sequence_hash(buffer, strlen(buffer));
 
168
        hashval = hashval >> (32 - kLog2HashSize);
 
169
 
 
170
        id_entries[curr_id_num].id = seq_id;
 
171
        id_entries[curr_id_num].query_index = query_index++;
 
172
        id_entries[curr_id_num].next_id = hashtable[hashval];
 
173
        hashtable[hashval] = curr_id_num++;
 
174
    }
 
175
 
 
176
    /* for each mask location, find the query sequence containing
 
177
       that mask and add to the list of filter locations for
 
178
       that query. Note that this assumes IDs for all query
 
179
       sequences are unique */
 
180
 
 
181
    for (mask_slp = mask_seqlocs; mask_slp; mask_slp = mask_slp->next) {
 
182
       SeqLocPtr current_mask = (SeqLocPtr) mask_slp->data.ptrvalue;
 
183
       Uint4 hashval;
 
184
       SeqIdPtr mask_id;
 
185
       Char buffer[64];
 
186
 
 
187
       if (current_mask == NULL)
 
188
           continue;
 
189
 
 
190
       mask_id = SeqLocId(current_mask);
 
191
       SeqIdLabel(mask_id, buffer, sizeof(buffer), OM_LABEL_CONTENT);
 
192
       hashval = readdb_sequence_hash(buffer, strlen(buffer));
 
193
       hashval = hashval >> (32 - kLog2HashSize);
 
194
 
 
195
       /* examine only the query IDs that hash to the same value */
 
196
       if (hashtable[hashval] != 0) {
 
197
          Int4 offset = hashtable[hashval];
 
198
          while (offset != 0) {
 
199
 
 
200
             SeqIdHash *q_entry = id_entries + offset;
 
201
 
 
202
             if (SeqIdMatch(mask_id, q_entry->id)) {
 
203
                Int4 context_idx = kNumContexts * q_entry->query_index;
 
204
                retval->seqloc_array[context_idx] = 
 
205
                              s_BlastSeqLocFromSeqLoc(current_mask,
 
206
                                         retval->seqloc_array[context_idx]);
 
207
                break;
 
208
             }
 
209
             offset = q_entry->next_id;
 
210
          }
 
211
       }
 
212
    }
 
213
 
 
214
    sfree(hashtable);
 
215
    sfree(id_entries);
 
216
 
 
217
    /* iterate through the query sequences and compute
 
218
       the complement of the filtering locations for each */
 
219
 
 
220
    for (query_slp = query_seqlocs, query_index = 0; 
 
221
         query_slp; 
 
222
         query_slp = query_slp->next, query_index++) {
91
223
 
92
224
        const int kCtxIndex = kNumContexts * query_index; /* context index */
93
 
        SeqLocPtr mask_slp = NULL;
94
 
 
95
 
        for (mask_slp = mask_seqlocs; mask_slp; mask_slp = mask_slp->next)
96
 
        {
97
 
           SeqLocPtr current_mask = (SeqLocPtr) mask_slp->data.ptrvalue;
98
 
           /* If mask is empty, advance to the next link in the mask chain.
99
 
              If mask Seq-id does not match sequence Seq-id, stay with this mask
100
 
              for the next link in the sequence Seq-loc chain. */
101
 
           if (current_mask && 
102
 
               SeqIdMatch(SeqLocId(current_mask), SeqLocId(current_query_loc))) 
103
 
           {
104
 
               retval->seqloc_array[kCtxIndex] = 
105
 
                     s_BlastSeqLocFromSeqLoc(current_mask,
106
 
                                             retval->seqloc_array[kCtxIndex]);
107
 
           }
108
 
        }
109
225
         
110
226
        if (retval->seqloc_array[kCtxIndex])
111
227
        {
113
229
                !Blast_QueryIsTranslated(program_number) &&
114
230
                !Blast_ProgramIsPhiBlast(program_number);
115
231
            BlastSeqLoc_RestrictToInterval(&retval->seqloc_array[kCtxIndex], 
116
 
                                           SeqLocStart(current_query_loc), 
117
 
                                           SeqLocStop(current_query_loc));
 
232
                                           SeqLocStart(query_slp), 
 
233
                                           SeqLocStop(query_slp));
118
234
            if (kIsNa) {
119
235
                /* N.B.: Unlike in the C++ APIs, this logic is only applied to
120
236
                 * non-translated nucleotide queries. See comment for
121
237
                 * BlastMaskLocDNAToProtein */
122
 
                Uint1 strand = SeqLocStrand(current_query_loc);
 
238
                Uint1 strand = SeqLocStrand(query_slp);
123
239
                if (strand == Seq_strand_minus) {
124
240
                    retval->seqloc_array[kCtxIndex+1] = 
125
241
                        retval->seqloc_array[kCtxIndex];
155
271
   Int4 index;
156
272
   const Boolean k_translate = Blast_QueryIsTranslated(program_number);
157
273
   const Uint1 k_num_frames = BLAST_GetNumberOfContexts(program_number);
 
274
   const Boolean kIsNucl = (program_number == eBlastTypeBlastn);
158
275
   SeqLoc* slp;
 
276
   Boolean all_minus = TRUE;
159
277
 
160
278
   if (mask_loc == NULL || mask_loc->seqloc_array == NULL)
161
279
      return NULL;
162
280
 
 
281
   for (slp = query_loc; slp; slp = slp->next)
 
282
   {
 
283
        Uint1 strand = SeqLocStrand(slp);
 
284
        if (strand != Seq_strand_minus)
 
285
        {
 
286
           all_minus = FALSE;
 
287
           break;
 
288
        }
 
289
   }
 
290
 
163
291
   for (index=0, slp = query_loc; slp; ++index, slp = slp->next)
164
292
   {
165
293
      const int kCtxIndex = k_num_frames * index; /* context index */
170
298
      {
171
299
         BlastSeqLoc* loc = NULL;
172
300
         SeqLocPtr mask_slp_head = NULL, mask_slp_tail = NULL;
173
 
         for (loc = mask_loc->seqloc_array[tmp_index]; loc; loc = loc->next)
 
301
         if (all_minus || BlastIsReverseStrand(kIsNucl , tmp_index) == FALSE)
174
302
         {
175
 
            SeqIntPtr si = SeqIntNew();
176
 
            si->from = loc->ssr->left + slp_from;
177
 
            si->to = loc->ssr->right + slp_from;
178
 
            si->id = SeqIdDup(seqid);
179
 
            /* Append the pointer, but also keep track of the tail of the list
180
 
             * so that appending to the list is a constant operation */
181
 
            mask_slp_tail = ValNodeAddPointer
182
 
                ( (mask_slp_tail ? &mask_slp_tail : &mask_slp_head), 
183
 
                  SEQLOC_INT, si);
 
303
            for (loc = mask_loc->seqloc_array[tmp_index]; loc; loc = loc->next)
 
304
            {
 
305
               SeqIntPtr si = SeqIntNew();
 
306
               si->from = loc->ssr->left + slp_from;
 
307
               si->to = loc->ssr->right + slp_from;
 
308
               si->id = SeqIdDup(seqid);
 
309
               /* Append the pointer, but also keep track of the tail of the list
 
310
                * so that appending to the list is a constant operation */
 
311
               mask_slp_tail = ValNodeAddPointer
 
312
                   ( (mask_slp_tail ? &mask_slp_tail : &mask_slp_head), 
 
313
                     SEQLOC_INT, si);
 
314
            }
184
315
         }
185
316
 
186
317
         if (mask_slp_head) {
212
343
 * @param qinfo  Query info structure containing contexts. [in/out]
213
344
 * @param index  Index of the context to fill. [in]
214
345
 * @param length Length of this context. [in]
215
 
 * @param prog   Program type of this search. [in]
216
346
 */
217
347
static void
218
348
s_QueryInfoSetContextInfo(BlastQueryInfo*   qinfo,
219
349
                          Uint4             index,
220
 
                          Uint4             length,
221
 
                          EBlastProgramType prog)
 
350
                          Uint4             length)
222
351
{
223
 
    qinfo->contexts[index].frame =
224
 
        BLAST_ContextToFrame(prog, index);
225
 
    
226
 
    qinfo->contexts[index].query_index =
227
 
        Blast_GetQueryIndexFromContext(index, prog);
228
 
    
229
352
    if (index) {
230
353
        Uint4 prev_loc = qinfo->contexts[index-1].query_offset;
231
354
        Uint4 prev_len = qinfo->contexts[index-1].query_length;
234
357
        
235
358
        qinfo->contexts[index].query_offset = prev_loc + shift;
236
359
        qinfo->contexts[index].query_length = length;
 
360
        if (length == 0)
 
361
            qinfo->contexts[index].is_valid = FALSE;
 
362
 
237
363
    } else {
238
364
        /* First context */
239
365
        qinfo->contexts[0].query_offset = 0;
240
366
        qinfo->contexts[0].query_length = length;
 
367
        if (length == 0)
 
368
            qinfo->contexts[0].is_valid = FALSE;
241
369
    }
242
370
}
243
371
 
260
388
   Uint1 strand;
261
389
   BlastQueryInfo* query_info;
262
390
   Int4 index;
263
 
   Int4 total_contexts;
264
391
   Uint4 max_length = 0;
265
392
 
266
393
   if (translate)
270
397
   else
271
398
      num_frames = 1;
272
399
 
273
 
   if ((query_info = (BlastQueryInfo*) calloc(1, sizeof(BlastQueryInfo)))
274
 
       == NULL)
 
400
   if ((query_info = BlastQueryInfoNew(program, ValNodeLen(slp))) == NULL) 
275
401
      return -1;
276
402
 
277
 
   query_info->first_context = 0;
278
 
   query_info->num_queries = ValNodeLen(slp);
279
 
   query_info->last_context = query_info->num_queries*num_frames - 1;
280
 
   total_contexts = query_info->last_context + 1;
281
 
 
282
403
   if ((strand = SeqLocStrand(slp)) == Seq_strand_minus) {
283
404
      if (translate)
284
405
         query_info->first_context = 3;
285
406
      else
286
407
         query_info->first_context = 1;
287
408
   }
288
 
 
289
 
   query_info->contexts = calloc(total_contexts, sizeof(BlastContextInfo));
290
409
   
291
410
   /* Fill the context offsets */
292
411
   for (index = 0; slp; slp = slp->next, index += num_frames) {
307
426
 
308
427
         /* Set the unused initial contexts if any */
309
428
         for (frame = 0; frame < first_frame; ++frame) {
310
 
             s_QueryInfoSetContextInfo(query_info, index+frame, 0, program);
 
429
             s_QueryInfoSetContextInfo(query_info, index+frame, 0);
311
430
         }
312
431
         
313
432
         for (frame = first_frame; frame <= last_frame; ++frame) {
316
435
 
317
436
            s_QueryInfoSetContextInfo(query_info,
318
437
                                     index+frame,
319
 
                                     protein_length,
320
 
                                     program);
 
438
                                     protein_length);
321
439
         }
322
440
 
323
441
         /* Set the unused trailing contexts if any */
324
442
         for (frame = last_frame + 1; frame < num_frames; ++frame) {
325
 
             s_QueryInfoSetContextInfo(query_info, index+frame, 0, program);
 
443
             s_QueryInfoSetContextInfo(query_info, index+frame, 0);
326
444
         }
327
445
      } else {
328
446
         max_length = MAX(max_length, length);
329
447
         
330
448
         if (is_na) {
331
449
            if (strand == Seq_strand_plus) {
332
 
                s_QueryInfoSetContextInfo(query_info, index,   length, program);
333
 
                s_QueryInfoSetContextInfo(query_info, index+1, 0,      program);
 
450
                s_QueryInfoSetContextInfo(query_info, index,   length);
 
451
                s_QueryInfoSetContextInfo(query_info, index+1, 0);
334
452
            } else if (strand == Seq_strand_minus) {
335
 
                s_QueryInfoSetContextInfo(query_info, index,   0,      program);
336
 
                s_QueryInfoSetContextInfo(query_info, index+1, length, program);
 
453
                s_QueryInfoSetContextInfo(query_info, index,   0);
 
454
                s_QueryInfoSetContextInfo(query_info, index+1, length);
337
455
            } else {
338
 
                s_QueryInfoSetContextInfo(query_info, index,   length, program);
339
 
                s_QueryInfoSetContextInfo(query_info, index+1, length, program);
 
456
                s_QueryInfoSetContextInfo(query_info, index,   length);
 
457
                s_QueryInfoSetContextInfo(query_info, index+1, length);
340
458
            }
341
459
         } else {
342
 
             s_QueryInfoSetContextInfo(query_info, index, length, program);
 
460
             s_QueryInfoSetContextInfo(query_info, index, length);
343
461
         }
344
462
      }
345
463
   }
369
487
 
370
488
   switch (encoding) {
371
489
   case eBlastEncodingProtein: 
372
 
      for (index = 0; index < size; index++) 
373
 
          buffer_var[index] = AMINOACID_TO_NCBISTDAA[buffer_var[index]];
 
490
      for (index = 0; index < size; index++) {
 
491
         Uint1 letter = buffer_var[index];
 
492
         if (letter == 'U' || letter == 'O' || letter == 'J')
 
493
            buffer_var[index] = AMINOACID_TO_NCBISTDAA['X'];
 
494
         else
 
495
            buffer_var[index] = AMINOACID_TO_NCBISTDAA[letter];
 
496
      }
374
497
      break;
375
498
   case eBlastEncodingNcbi4na:
376
499
      for (index = 0; index < size; index++) 
613
736
       query_info == NULL || query_blk == NULL)
614
737
      return -1;
615
738
 
616
 
 
617
739
   if ((status = s_QueryInfoSetUp(query_slp, program_number, query_info)))
618
740
      return status;
619
741
 
639
761
   /* Do not count the first and last sentinel bytes in the 
640
762
      query length */
641
763
   if ((status=BlastSetUp_SeqBlkNew(buffer, buffer_length-2, 
642
 
                                    0, query_blk, TRUE)))
 
764
                                    query_blk, TRUE)))
643
765
      return status;
644
766
 
645
767
   if (masking_locs) {
684
806
   /* Initialize the sequence block, saving the sequence buffer in 
685
807
      'sequence_start'. */
686
808
   if ((status=BlastSetUp_SeqBlkNew(subject_buffer, buffer_length,
687
 
                                    0, subject, TRUE)))
 
809
                                    subject, TRUE)))
688
810
      return status;
689
811
 
690
812
   /* If subject sequence is nucleotide, create compressed sequence buffer