← Back to branch summary

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

« back to all changes in this revision

Viewing changes to algo/blast/core/blast_query_info.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/core/blast_query_info.c
 
1
/* $Id: blast_query_info.c,v 1.3 2006/04/25 19:06:15 camacho Exp $
 
2
 * ===========================================================================
 
3
 *
 
4
 *                            PUBLIC DOMAIN NOTICE
 
5
 *               National Center for Biotechnology Information
 
6
 *
 
7
 *  This software/database is a "United States Government Work" under the
 
8
 *  terms of the United States Copyright Act.  It was written as part of
 
9
 *  the author's official duties as a United States Government employee and
 
10
 *  thus cannot be copyrighted.  This software/database is freely available
 
11
 *  to the public for use. The National Library of Medicine and the U.S.
 
12
 *  Government have not placed any restriction on its use or reproduction.
 
13
 *
 
14
 *  Although all reasonable efforts have been taken to ensure the accuracy
 
15
 *  and reliability of the software and data, the NLM and the U.S.
 
16
 *  Government do not and cannot warrant the performance or results that
 
17
 *  may be obtained by using this software or data. The NLM and the U.S.
 
18
 *  Government disclaim all warranties, express or implied, including
 
19
 *  warranties of performance, merchantability or fitness for any particular
 
20
 *  purpose.
 
21
 *
 
22
 *  Please cite the author in any work or product based on this material.
 
23
 *
 
24
 * ===========================================================================
 
25
 *
 
26
 * Author: Christiam Camacho
 
27
 *
 
28
 */
 
29
 
 
30
/** @file blast_query_info.c
 
31
 * Functions to manipulate the BlastQueryInfo structure
 
32
 */
 
33
 
 
34
 
 
35
#ifndef SKIP_DOXYGEN_PROCESSING
 
36
static char const rcsid[] = 
 
37
    "$Id: blast_query_info.c,v 1.3 2006/04/25 19:06:15 camacho Exp $";
 
38
#endif /* SKIP_DOXYGEN_PROCESSING */
 
39
 
 
40
#include <algo/blast/core/blast_query_info.h>
 
41
#include <algo/blast/core/blast_util.h>
 
42
#include <algo/blast/core/pattern.h>
 
43
 
 
44
Int4 
 
45
Blast_GetQueryIndexFromContext(Int4 context, EBlastProgramType program)
 
46
{
 
47
   Int4 index = -1;     /* -1 is used to indicate error */
 
48
   switch (program) {
 
49
   case eBlastTypeBlastn:
 
50
      index = context/NUM_STRANDS; 
 
51
      break;
 
52
   case eBlastTypeBlastp: 
 
53
   case eBlastTypeTblastn: 
 
54
   case eBlastTypeRpsBlast: 
 
55
   case eBlastTypePsiBlast:
 
56
   case eBlastTypeRpsTblastn:
 
57
   case eBlastTypePhiBlastn:
 
58
   case eBlastTypePhiBlastp:
 
59
      index = context; 
 
60
      break;
 
61
   case eBlastTypeBlastx: 
 
62
   case eBlastTypeTblastx:
 
63
      index = context/NUM_FRAMES; 
 
64
      break;
 
65
   default:
 
66
      break;
 
67
   }
 
68
   return index;
 
69
}
 
70
 
 
71
BlastQueryInfo* BlastQueryInfoNew(EBlastProgramType program, int num_queries)
 
72
{
 
73
    const unsigned int kNumContexts = BLAST_GetNumberOfContexts(program);
 
74
    BlastQueryInfo* retval = NULL;
 
75
    
 
76
    if (num_queries <= 0) {
 
77
        return retval;
 
78
    }
 
79
    ASSERT(kNumContexts != 0);
 
80
 
 
81
    retval = (BlastQueryInfo*) calloc(1, sizeof(BlastQueryInfo));
 
82
    if ( !retval ) {
 
83
        return BlastQueryInfoFree(retval);
 
84
    }
 
85
 
 
86
    retval->num_queries = num_queries;
 
87
 
 
88
    retval->first_context = 0;
 
89
    retval->last_context = retval->num_queries * kNumContexts - 1;
 
90
 
 
91
    retval->contexts = (BlastContextInfo*) calloc(retval->last_context + 1,
 
92
                                                  sizeof(BlastContextInfo));
 
93
 
 
94
    if ( !retval->contexts ) {
 
95
        return BlastQueryInfoFree(retval);
 
96
    } else {
 
97
        int i;
 
98
        for (i = 0; i < retval->last_context + 1; i++) {
 
99
            retval->contexts[i].query_index =
 
100
                Blast_GetQueryIndexFromContext(i, program);
 
101
            ASSERT(retval->contexts[i].query_index != -1);
 
102
 
 
103
            retval->contexts[i].frame = BLAST_ContextToFrame(program,  i);
 
104
            ASSERT(retval->contexts[i].frame != INT1_MAX);
 
105
 
 
106
            retval->contexts[i].is_valid = TRUE;
 
107
        }
 
108
    }
 
109
 
 
110
    return retval;
 
111
}
 
112
 
 
113
BlastQueryInfo* BlastQueryInfoFree(BlastQueryInfo* query_info)
 
114
{
 
115
    if (query_info) {
 
116
        sfree(query_info->contexts);
 
117
        query_info->pattern_info = 
 
118
            SPHIQueryInfoFree(query_info->pattern_info);
 
119
        sfree(query_info);
 
120
    }
 
121
    return NULL;
 
122
}
 
123
 
 
124
BlastQueryInfo* BlastQueryInfoDup(BlastQueryInfo* query_info)
 
125
{
 
126
   BlastQueryInfo* retval = BlastMemDup(query_info, sizeof(BlastQueryInfo));
 
127
   Int4 num_contexts = query_info->last_context + 1;
 
128
   
 
129
   retval->contexts =
 
130
       BlastMemDup(query_info->contexts, num_contexts * sizeof(BlastContextInfo));
 
131
   
 
132
   if (query_info->pattern_info) {
 
133
       retval->pattern_info = 
 
134
           SPHIQueryInfoCopy(query_info->pattern_info);
 
135
   }
 
136
 
 
137
   return retval;
 
138
}
 
139
 
 
140
/** Calculates length of the DNA query from the BlastQueryInfo structure that 
 
141
 * contains context information for translated frames for a set of queries.
 
142
 * @param query_info Query information containing data for all contexts [in]
 
143
 * @param query_index Which query to find DNA length for?
 
144
 * @return DNA length of the query, calculated as sum of 3 protein frame lengths, 
 
145
 *         plus 2, because 2 last nucleotide residues do not have a 
 
146
 *         corresponding codon.
 
147
 */
 
148
static Int4 
 
149
s_GetTranslatedQueryDNALength(const BlastQueryInfo* query_info, Int4 query_index)
 
150
{
 
151
    Int4 start_context = NUM_FRAMES*query_index;
 
152
    Int4 dna_length = 2;
 
153
    Int4 index;
 
154
 
 
155
    /* Make sure that query index is within appropriate range, and that this is
 
156
       really a translated search */
 
157
    ASSERT(query_index < query_info->num_queries);
 
158
    ASSERT(start_context < query_info->last_context);
 
159
 
 
160
    /* If only reverse strand is searched, then forward strand contexts don't 
 
161
       have lengths information */
 
162
    if (query_info->contexts[start_context].query_length == 0)
 
163
        start_context += 3;
 
164
 
 
165
    for (index = start_context; index < start_context + 3; ++index)
 
166
        dna_length += query_info->contexts[index].query_length;
 
167
 
 
168
    return dna_length;
 
169
}
 
170
 
 
171
Int4 BlastQueryInfoGetQueryLength(const BlastQueryInfo* qinfo,
 
172
                                  EBlastProgramType program,
 
173
                                  Int4 query_index)
 
174
{
 
175
    const Uint4 kNumContexts = BLAST_GetNumberOfContexts(program);
 
176
    ASSERT(query_index < qinfo->num_queries);
 
177
 
 
178
    if (Blast_QueryIsTranslated(program)) {
 
179
        return s_GetTranslatedQueryDNALength(qinfo, query_index);
 
180
    } else if (program == eBlastTypeBlastn) {
 
181
        Int4 retval = qinfo->contexts[query_index*kNumContexts].query_length;
 
182
        if (retval <= 0) {
 
183
            retval = qinfo->contexts[query_index*kNumContexts+1].query_length;
 
184
        }
 
185
        return retval;
 
186
    } else {
 
187
        return qinfo->contexts[query_index*kNumContexts].query_length;
 
188
    }
 
189
}
 
190
 
 
191
/* FIXME: should the EBlastProgramType be added as a member of the
 
192
 * BlastQueryInfo structure? Without it, there's many operations that can't be
 
193
 * done, so it doesn't make sense to have them separate... */
 
194
Int8
 
195
BlastQueryInfoGetEffSearchSpace(const BlastQueryInfo* qinfo,
 
196
                                EBlastProgramType program,
 
197
                                Int4 query_index)
 
198
{
 
199
    Int8 retval = 0;
 
200
    Int4 i = 0;
 
201
    const Uint4 kNumContexts = BLAST_GetNumberOfContexts(program);
 
202
    ASSERT(query_index < qinfo->num_queries);
 
203
 
 
204
    for (i = query_index*kNumContexts; i < (query_index+1)*kNumContexts; i++) {
 
205
        if ( (retval = qinfo->contexts[i].eff_searchsp) != 0) {
 
206
            break;
 
207
        }
 
208
    }
 
209
    return retval;
 
210
}
 
211
 
 
212
void
 
213
BlastQueryInfoSetEffSearchSpace(BlastQueryInfo* qinfo,
 
214
                                EBlastProgramType program,
 
215
                                Int4 query_index,
 
216
                                Int8 eff_searchsp)
 
217
{
 
218
    Int4 i = 0;
 
219
    const Uint4 kNumContexts = BLAST_GetNumberOfContexts(program);
 
220
    ASSERT(query_index < qinfo->num_queries);
 
221
 
 
222
    for (i = query_index*kNumContexts; i < (query_index+1)*kNumContexts; i++) {
 
223
        qinfo->contexts[i].eff_searchsp = eff_searchsp;
 
224
    }
 
225
}
 
226
 
 
227
Int4 BSearchContextInfo(Int4 n, const BlastQueryInfo * A)
 
228
{
 
229
    Int4 m=0, b=0, e=0, size=0;
 
230
    
 
231
    size = A->last_context+1;
 
232
    
 
233
    b = 0;
 
234
    e = size;
 
235
    while (b < e - 1) {
 
236
        m = (b + e) / 2;
 
237
        if (A->contexts[m].query_offset > n)
 
238
            e = m;
 
239
        else
 
240
            b = m;
 
241
    }
 
242
    return b;
 
243
}
 
244
 
 
245
Uint4
 
246
QueryInfo_GetSeqBufLen(const BlastQueryInfo* qinfo)
 
247
{
 
248
    BlastContextInfo * cinfo = & qinfo->contexts[qinfo->last_context];
 
249
    return cinfo->query_offset + cinfo->query_length + (cinfo->query_length ? 2 : 1);
 
250
}
 
251
 
 
252
Int4 *
 
253
ContextOffsetsToOffsetArray(BlastQueryInfo* info)
 
254
{
 
255
    /* The Many Values of 'Length'
 
256
     *
 
257
     * 1. info->last_context: the index of the last query offset.
 
258
     *
 
259
     * 2. count: the number of query offsets.
 
260
     *
 
261
     * 3. count + 1: the size of the output array (has an 'extra'
 
262
     *    member so as to communicate the last sequence length).
 
263
     *
 
264
     * 4. sz: the size of the context object array
 
265
     */
 
266
    
 
267
    Uint4   count  = (info->last_context + 1);
 
268
    Uint4   sz     = sizeof(Int4) * (count+1);
 
269
    Uint4   frame  = 0;
 
270
    Int4  * result = 0;
 
271
    
 
272
    ASSERT(info);
 
273
    ASSERT(info->contexts);
 
274
    
 
275
    result = malloc(sz);
 
276
    memset(result, 0, sz);
 
277
    
 
278
    for(frame = 0; frame < count; frame++) {
 
279
        result[frame] = info->contexts[frame].query_offset;
 
280
    }
 
281
    
 
282
    /* One more entry, provides length info for last element. */
 
283
    
 
284
    result[count] = info->contexts[count-1].query_offset;
 
285
    
 
286
    if (info->contexts[count-1].query_length) {
 
287
        result[count] += info->contexts[count-1].query_length + 1;
 
288
    }
 
289
    
 
290
    return result;
 
291
}
 
292
 
 
293
void
 
294
OffsetArrayToContextOffsets(BlastQueryInfo    * info,
 
295
                            Int4              * new_offsets,
 
296
                            EBlastProgramType   prog)
 
297
{
 
298
    Uint4 count = (info->last_context + 1);
 
299
    Uint4 i     = 0;
 
300
    
 
301
    ASSERT(info);
 
302
    ASSERT(new_offsets);
 
303
    
 
304
    if (! info->contexts) {
 
305
        info->contexts = calloc(count, sizeof(BlastContextInfo));
 
306
    }
 
307
    
 
308
    for(i = 0; i < count; i++) {
 
309
        Int4 distance = 0;
 
310
        
 
311
        info->contexts[i].query_offset = new_offsets[i];
 
312
        
 
313
        distance = new_offsets[i+1] - new_offsets[i];
 
314
        info->contexts[i].query_length = distance ? distance-1 : 0;
 
315
        
 
316
        /* Set the frame and query index */
 
317
        
 
318
        info->contexts[i].frame =
 
319
            BLAST_ContextToFrame(prog, i);
 
320
        
 
321
        info->contexts[i].query_index =
 
322
            Blast_GetQueryIndexFromContext(i, prog);
 
323
    }
 
324
}
 
325
 
 
326
Int2
 
327
Blast_GetOneQueryStructs(BlastQueryInfo** one_query_info_ptr, 
 
328
                         BLAST_SequenceBlk** one_query_ptr,
 
329
                         const BlastQueryInfo* query_info, 
 
330
                         BLAST_SequenceBlk* query, Int4 query_index)
 
331
{
 
332
    Int4 num_frames;
 
333
    Int4 index;
 
334
    Int4 first_context;
 
335
    Int4 query_offset;
 
336
    BlastQueryInfo* one_query_info = NULL;
 
337
    BLAST_SequenceBlk* one_query = NULL;
 
338
 
 
339
    if (!one_query_info_ptr || !one_query_ptr || !query_info || !query ||
 
340
        query_index >= query_info->num_queries)
 
341
        return -1;
 
342
 
 
343
    num_frames = (query_info->last_context / query_info->num_queries) + 1;
 
344
    first_context = query_index*num_frames;
 
345
    query_offset = query_info->contexts[first_context].query_offset;
 
346
 
 
347
    one_query_info = *one_query_info_ptr;
 
348
    /* If this hasn't been already done, allocate new query information 
 
349
       structure. */
 
350
    if (!one_query_info) {
 
351
        one_query_info = (BlastQueryInfo*) calloc(1, sizeof(BlastQueryInfo));
 
352
        *one_query_info_ptr = one_query_info;
 
353
        one_query_info->contexts = (BlastContextInfo*) calloc(num_frames, sizeof(BlastContextInfo));
 
354
    }
 
355
    one_query = *one_query_ptr;
 
356
    /* If this hasn't been already done, allocate new sequence block. */
 
357
    if (!one_query) {
 
358
        one_query = (BLAST_SequenceBlk*) calloc(1, sizeof(BLAST_SequenceBlk));
 
359
        *one_query_ptr = one_query;
 
360
    }
 
361
    if (!one_query_info || !one_query)
 
362
        return -1;
 
363
 
 
364
    one_query_info->num_queries = 1;
 
365
    one_query_info->last_context = num_frames - 1;
 
366
    
 
367
    memcpy(one_query_info->contexts,
 
368
           &query_info->contexts[first_context],
 
369
           num_frames * sizeof(BlastContextInfo));
 
370
    
 
371
    /* Make context offsets relative to this query. */
 
372
    for (index = 0; index < num_frames; ++index) {
 
373
        one_query_info->contexts[index].query_offset -= query_offset;
 
374
    }
 
375
    
 
376
    /* Fill the sequence block information for this one query. */
 
377
    memset(one_query, 0, sizeof(BLAST_SequenceBlk));
 
378
    one_query->sequence = &query->sequence[query_offset];
 
379
    one_query->length =
 
380
        one_query_info->contexts[num_frames-1].query_offset +
 
381
        one_query_info->contexts[num_frames-1].query_length;
 
382
    one_query->sequence_allocated = FALSE;
 
383
    one_query->oid = query_index;
 
384
 
 
385
    return 0;
 
386
}
 
387