← Back to branch summary

~ubuntu-branches/debian/experimental/ncbi-tools6/experimental

« back to all changes in this revision

Viewing changes to algo/blast/core/blast_psi.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2005-03-27 12:00:15 UTC
  • mfrom: (2.1.2 hoary)
  • Revision ID: james.westby@ubuntu.com-20050327120015-embhesp32nj73p9r
Tags: 6.1.20041020-3
http://bugs.debian.org/295110
* Fix FTBFS under GCC 4.0 caused by inconsistent use of "static" on
  functions.  (Closes: #295110.)
* Add a watch file, now that we can.  (Upstream's layout needs version=3.)

Show diffs side-by-side

added added

removed removed

Lines of Context:
algo/blast/core/blast_psi.c
 
1
static char const rcsid[] =
 
2
    "$Id: blast_psi.c,v 1.14 2004/10/18 14:33:14 camacho Exp $";
 
3
/* ===========================================================================
 
4
 *
 
5
 *                            PUBLIC DOMAIN NOTICE
 
6
 *               National Center for Biotechnology Information
 
7
 *
 
8
 *  This software/database is a "United States Government Work" under the
 
9
 *  terms of the United States Copyright Act.  It was written as part of
 
10
 *  the author's offical duties as a United States Government employee and
 
11
 *  thus cannot be copyrighted.  This software/database is freely available
 
12
 *  to the public for use. The National Library of Medicine and the U.S.
 
13
 *  Government have not placed any restriction on its use or reproduction.
 
14
 *
 
15
 *  Although all reasonable efforts have been taken to ensure the accuracy
 
16
 *  and reliability of the software and data, the NLM and the U.S.
 
17
 *  Government do not and cannot warrant the performance or results that
 
18
 *  may be obtained by using this software or data. The NLM and the U.S.
 
19
 *  Government disclaim all warranties, express or implied, including
 
20
 *  warranties of performance, merchantability or fitness for any particular
 
21
 *  purpose.
 
22
 *
 
23
 *  Please cite the author in any work or product based on this material.
 
24
 *
 
25
 * ===========================================================================
 
26
 *
 
27
 * Author:  Christiam Camacho
 
28
 *
 
29
 */
 
30
 
 
31
/** @file blast_psi.c
 
32
 * Implementation of the high level functions of PSI-BLAST's PSSM engine.
 
33
 */
 
34
    
 
35
#include <algo/blast/core/blast_stat.h>
 
36
#include <algo/blast/core/blast_encoding.h>
 
37
#include "blast_psi_priv.h"
 
38
 
 
39
/****************************************************************************/
 
40
/* Function prototypes */
 
41
 
 
42
/** Convenience function to deallocate data structures allocated in
 
43
 * PSICreatePssmWithDiagnostics.
 
44
 */
 
45
static void
 
46
_PSICreatePssmCleanUp(PSIMatrix** pssm,
 
47
                      _PSIMsa* msa,
 
48
                      _PSIAlignedBlock* aligned_block,
 
49
                      _PSISequenceWeights* seq_weights,
 
50
                      _PSIInternalPssmData* internal_pssm);
 
51
 
 
52
/** Copies pssm data from internal_pssm and sbp into pssm. None of its
 
53
 * parameters can be NULL. */
 
54
static void
 
55
_PSISavePssm(const _PSIInternalPssmData* internal_pssm,
 
56
             const BlastScoreBlk* sbp,
 
57
             PSIMatrix* pssm);
 
58
 
 
59
/****************************************************************************/
 
60
 
 
61
int
 
62
PSICreatePssm(const PSIMsa* msap,
 
63
              const PSIBlastOptions* options,
 
64
              BlastScoreBlk* sbp,
 
65
              PSIMatrix** pssm)
 
66
{
 
67
    return PSICreatePssmWithDiagnostics(msap, options, sbp, NULL,
 
68
                                        pssm, NULL);
 
69
}
 
70
 
 
71
int
 
72
PSICreatePssmWithDiagnostics(const PSIMsa* msap,                    /* [in] */
 
73
                             const PSIBlastOptions* options,        /* [in] */
 
74
                             BlastScoreBlk* sbp,                    /* [in] */
 
75
                             const PSIDiagnosticsRequest* request,  /* [in] */
 
76
                             PSIMatrix** pssm,                      /* [out] */
 
77
                             PSIDiagnosticsResponse** diagnostics)  /* [out] */
 
78
{
 
79
    _PSIMsa* msa = NULL;
 
80
    _PSIAlignedBlock* aligned_block = NULL;
 
81
    _PSISequenceWeights* seq_weights = NULL; 
 
82
    _PSIInternalPssmData* internal_pssm = NULL;
 
83
    int status = 0;
 
84
 
 
85
    if ( !msap || !options || !sbp || !pssm ) {
 
86
        return PSIERR_BADPARAM;
 
87
    }
 
88
 
 
89
    /*** Allocate data structures ***/
 
90
    msa = _PSIMsaNew(msap, (Uint4) sbp->alphabet_size);
 
91
    aligned_block = _PSIAlignedBlockNew(msa->dimensions->query_length);
 
92
    seq_weights = _PSISequenceWeightsNew(msa->dimensions, sbp);
 
93
    internal_pssm = _PSIInternalPssmDataNew(msa->dimensions->query_length,
 
94
                                            (Uint4) sbp->alphabet_size);
 
95
    *pssm = PSIMatrixNew(msa->dimensions->query_length, 
 
96
                         (Uint4) sbp->alphabet_size);
 
97
    if ( !msa || ! aligned_block || !seq_weights || !internal_pssm || !*pssm ) {
 
98
        _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights,
 
99
                              internal_pssm);
 
100
        return PSIERR_OUTOFMEM;
 
101
    }
 
102
 
 
103
    /*** Enable structure group customization if needed ***/
 
104
    if (options->ignore_consensus) {
 
105
        _PSIPurgeAlignedRegion(msa, kQueryIndex, 0,
 
106
                               msa->dimensions->query_length);
 
107
    }
 
108
 
 
109
    /*** Run the engine's stages ***/
 
110
    status = _PSIPurgeBiasedSegments(msa);
 
111
    if (status != PSI_SUCCESS) {
 
112
        _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights, 
 
113
                              internal_pssm);
 
114
        return status;
 
115
    }
 
116
 
 
117
    status = _PSIValidateMSA(msa, options->ignore_consensus);
 
118
    if (status != PSI_SUCCESS) {
 
119
        _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights, 
 
120
                              internal_pssm);
 
121
        return status;
 
122
    }
 
123
 
 
124
    status = _PSIComputeAlignmentBlocks(msa, aligned_block);
 
125
    if (status != PSI_SUCCESS) {
 
126
        _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights, 
 
127
                              internal_pssm);
 
128
        return status;
 
129
    }
 
130
 
 
131
    status = _PSIComputeSequenceWeights(msa, aligned_block, 
 
132
                                        options->ignore_consensus,
 
133
                                        seq_weights);
 
134
    if (status != PSI_SUCCESS) {
 
135
        _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights, 
 
136
                              internal_pssm);
 
137
        return status;
 
138
    }
 
139
 
 
140
    status = _PSIComputeResidueFrequencies(msa, seq_weights, sbp, 
 
141
                                           aligned_block, 
 
142
                                           options->pseudo_count, 
 
143
                                           internal_pssm);
 
144
    if (status != PSI_SUCCESS) {
 
145
        _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights, 
 
146
                              internal_pssm);
 
147
        return status;
 
148
    }
 
149
 
 
150
    status = _PSIConvertResidueFreqsToPSSM(internal_pssm, msa->query, sbp, 
 
151
                                           seq_weights->std_prob);
 
152
    if (status != PSI_SUCCESS) {
 
153
        _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights, 
 
154
                              internal_pssm);
 
155
        return status;
 
156
    }
 
157
 
 
158
    /* FIXME: instead of NULL pass options->scaling_factor */
 
159
    status = _PSIScaleMatrix(msa->query, msa->dimensions->query_length, 
 
160
                             seq_weights->std_prob, NULL, internal_pssm, sbp);
 
161
    if (status != PSI_SUCCESS) {
 
162
        _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights, 
 
163
                              internal_pssm);
 
164
        return status;
 
165
    }
 
166
 
 
167
    /*** Save the pssm outgoing parameter ***/
 
168
    _PSISavePssm(internal_pssm, sbp, *pssm);
 
169
 
 
170
    /*** Save diagnostics if required ***/
 
171
    if (request && diagnostics) {
 
172
        *diagnostics = PSIDiagnosticsResponseNew(msa->dimensions->query_length,
 
173
                                                 (Uint4) sbp->alphabet_size,
 
174
                                                 request);
 
175
        if ( !*diagnostics ) {
 
176
            /* FIXME: This could be changed to return a warning and not
 
177
             * deallocate PSSM data */
 
178
            _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights, 
 
179
                                  internal_pssm);
 
180
            return PSIERR_OUTOFMEM;
 
181
        }
 
182
        status = _PSISaveDiagnostics(msa, aligned_block, seq_weights, 
 
183
                                     internal_pssm, *diagnostics);
 
184
        if (status != PSI_SUCCESS) {
 
185
            *diagnostics = PSIDiagnosticsResponseFree(*diagnostics);
 
186
            _PSICreatePssmCleanUp(pssm, msa, aligned_block, seq_weights,
 
187
                                  internal_pssm);
 
188
            return status;
 
189
        }
 
190
    }
 
191
    _PSICreatePssmCleanUp(NULL, msa, aligned_block, seq_weights, internal_pssm);
 
192
 
 
193
    return PSI_SUCCESS;
 
194
}
 
195
 
 
196
/****************************************************************************/
 
197
 
 
198
static void
 
199
_PSICreatePssmCleanUp(PSIMatrix** pssm,
 
200
                      _PSIMsa* msa,
 
201
                      _PSIAlignedBlock* aligned_block,
 
202
                      _PSISequenceWeights* seq_weights,
 
203
                      _PSIInternalPssmData* internal_pssm)
 
204
{
 
205
    if (pssm) {
 
206
        *pssm = PSIMatrixFree(*pssm);
 
207
    }
 
208
    _PSIMsaFree(msa);
 
209
    _PSIAlignedBlockFree(aligned_block);
 
210
    _PSISequenceWeightsFree(seq_weights);
 
211
    _PSIInternalPssmDataFree(internal_pssm);
 
212
}
 
213
 
 
214
static void
 
215
_PSISavePssm(const _PSIInternalPssmData* internal_pssm,
 
216
             const BlastScoreBlk* sbp,
 
217
             PSIMatrix* pssm)
 
218
{
 
219
    ASSERT(internal_pssm);
 
220
    ASSERT(sbp);
 
221
    ASSERT(pssm);
 
222
 
 
223
    _PSICopyIntMatrix(pssm->pssm, internal_pssm->pssm,
 
224
                      pssm->ncols, pssm->nrows);
 
225
 
 
226
    pssm->lambda = sbp->kbp_gap_psi[0]->Lambda;
 
227
    pssm->kappa = sbp->kbp_gap_psi[0]->K;
 
228
    pssm->h = sbp->kbp_gap_psi[0]->H;
 
229
}
 
230
 
 
231
/****************************************************************************/
 
232
 
 
233
PSIMsa*
 
234
PSIMsaNew(const PSIMsaDimensions* dimensions)
 
235
{
 
236
    PSIMsa* retval = NULL;
 
237
 
 
238
    if ( !dimensions ) {
 
239
        return NULL;
 
240
    }
 
241
 
 
242
    retval = (PSIMsa*) malloc(sizeof(PSIMsa));
 
243
    if ( !retval ) {
 
244
        return PSIMsaFree(retval);
 
245
    }
 
246
 
 
247
    retval->dimensions = (PSIMsaDimensions*) malloc(sizeof(PSIMsaDimensions));
 
248
    if ( !retval->dimensions ) {
 
249
        return PSIMsaFree(retval);
 
250
    }
 
251
    memcpy((void*) retval->dimensions,
 
252
           (void*) dimensions, 
 
253
           sizeof(PSIMsaDimensions));
 
254
 
 
255
    retval->data = (PSIMsaCell**) _PSIAllocateMatrix(dimensions->num_seqs + 1,
 
256
                                                     dimensions->query_length,
 
257
                                                     sizeof(PSIMsaCell));
 
258
    if ( !retval->data ) {
 
259
        return PSIMsaFree(retval);
 
260
    }
 
261
    {
 
262
        Uint4 s = 0;    /* index on sequences */
 
263
        Uint4 p = 0;    /* index on positions */
 
264
 
 
265
        for (s = 0; s < dimensions->num_seqs + 1; s++) {
 
266
            for (p = 0; p < dimensions->query_length; p++) {
 
267
                /* FIXME: change to 0 when old code is dropped */
 
268
                retval->data[s][p].letter = (Uint1) -1;
 
269
                retval->data[s][p].is_aligned = FALSE;
 
270
            }
 
271
        }
 
272
    }
 
273
 
 
274
    return retval;
 
275
}
 
276
 
 
277
PSIMsa*
 
278
PSIMsaFree(PSIMsa* msa)
 
279
{
 
280
    if ( !msa ) {
 
281
        return NULL;
 
282
    }
 
283
 
 
284
    if ( msa->data && msa->dimensions ) {
 
285
        _PSIDeallocateMatrix((void**) msa->data,
 
286
                             msa->dimensions->num_seqs + 1);
 
287
        msa->data = NULL;
 
288
    }
 
289
 
 
290
    if ( msa->dimensions ) {
 
291
        sfree(msa->dimensions);
 
292
    }
 
293
 
 
294
    sfree(msa);
 
295
 
 
296
    return NULL;
 
297
}
 
298
 
 
299
PSIMatrix*
 
300
PSIMatrixNew(Uint4 query_length, Uint4 alphabet_size)
 
301
{
 
302
    PSIMatrix* retval = NULL;
 
303
 
 
304
    retval = (PSIMatrix*) malloc(sizeof(PSIMatrix));
 
305
    if ( !retval ) {
 
306
        return NULL;
 
307
    }
 
308
    retval->ncols = query_length;
 
309
    retval->nrows = alphabet_size;
 
310
 
 
311
    retval->pssm = (int**) _PSIAllocateMatrix(query_length, alphabet_size,
 
312
                                              sizeof(int));
 
313
    if ( !(retval->pssm) ) {
 
314
        return PSIMatrixFree(retval);
 
315
    }
 
316
 
 
317
    retval->lambda = 0.0;
 
318
    retval->kappa = 0.0;
 
319
    retval->h = 0.0;
 
320
 
 
321
    return retval;
 
322
}
 
323
 
 
324
PSIMatrix*
 
325
PSIMatrixFree(PSIMatrix* matrix)
 
326
{
 
327
    if ( !matrix ) {
 
328
        return NULL;
 
329
    }
 
330
 
 
331
    if (matrix->pssm) {
 
332
        _PSIDeallocateMatrix((void**) matrix->pssm, matrix->ncols);
 
333
    }
 
334
 
 
335
    sfree(matrix);
 
336
 
 
337
    return NULL;
 
338
}
 
339
 
 
340
PSIDiagnosticsResponse*
 
341
PSIDiagnosticsResponseNew(Uint4 query_length, Uint4 alphabet_size,
 
342
                          const PSIDiagnosticsRequest* wants)
 
343
{
 
344
    PSIDiagnosticsResponse* retval = NULL;
 
345
 
 
346
    if ( !wants ) {
 
347
        return NULL;
 
348
    }
 
349
 
 
350
    /* MUST use calloc to allocate structure because code that uses this
 
351
     * structure assumes that non-NULL members will require to be populated */
 
352
    retval = (PSIDiagnosticsResponse*) calloc(1, 
 
353
                                              sizeof(PSIDiagnosticsResponse));
 
354
    if ( !retval ) {
 
355
        return NULL;
 
356
    }
 
357
 
 
358
    retval->query_length = query_length;
 
359
    retval->alphabet_size = alphabet_size;
 
360
 
 
361
    if (wants->information_content) {
 
362
        retval->information_content = (double*) 
 
363
            calloc(query_length, sizeof(double));
 
364
        if ( !(retval->information_content) ) {
 
365
            return PSIDiagnosticsResponseFree(retval);
 
366
        }
 
367
    }
 
368
 
 
369
    if (wants->residue_frequencies) {
 
370
        retval->residue_freqs = (Uint4**) _PSIAllocateMatrix(query_length, 
 
371
                                                             alphabet_size, 
 
372
                                                             sizeof(Uint4));
 
373
        if ( !(retval->residue_freqs) ) {
 
374
            return PSIDiagnosticsResponseFree(retval);
 
375
        }
 
376
    }
 
377
 
 
378
    if (wants->weighted_residue_frequencies) {
 
379
        retval->weighted_residue_freqs = (double**) 
 
380
            _PSIAllocateMatrix(query_length, 
 
381
                               alphabet_size, 
 
382
                               sizeof(double));
 
383
        if ( !(retval->weighted_residue_freqs) ) {
 
384
            return PSIDiagnosticsResponseFree(retval);
 
385
        }
 
386
    }
 
387
 
 
388
    if (wants->frequency_ratios) {
 
389
        retval->frequency_ratios = (double**)
 
390
            _PSIAllocateMatrix(query_length, 
 
391
                               alphabet_size, 
 
392
                               sizeof(double));
 
393
        if ( !retval->frequency_ratios ) {
 
394
            return PSIDiagnosticsResponseFree(retval);
 
395
        }
 
396
    }
 
397
 
 
398
    if (wants->gapless_column_weights) {
 
399
        retval->gapless_column_weights = (double*) 
 
400
            calloc(query_length, sizeof(double));
 
401
        if ( !(retval->gapless_column_weights) ) {
 
402
            return PSIDiagnosticsResponseFree(retval);
 
403
        }
 
404
    }
 
405
 
 
406
    return retval;
 
407
}
 
408
 
 
409
PSIDiagnosticsResponse*
 
410
PSIDiagnosticsResponseFree(PSIDiagnosticsResponse* diags)
 
411
{
 
412
    if ( !diags )
 
413
        return NULL;
 
414
 
 
415
    if (diags->information_content) {
 
416
        sfree(diags->information_content);
 
417
    }
 
418
 
 
419
    if (diags->residue_freqs) {
 
420
        _PSIDeallocateMatrix((void**) diags->residue_freqs,
 
421
                             diags->query_length);
 
422
    }
 
423
 
 
424
    if (diags->weighted_residue_freqs) {
 
425
        _PSIDeallocateMatrix((void**) diags->weighted_residue_freqs,
 
426
                             diags->query_length);
 
427
    }
 
428
 
 
429
    if (diags->frequency_ratios) {
 
430
        _PSIDeallocateMatrix((void**) diags->frequency_ratios,
 
431
                             diags->query_length);
 
432
    }
 
433
 
 
434
    if (diags->gapless_column_weights) {
 
435
        sfree(diags->gapless_column_weights);
 
436
    }
 
437
 
 
438
    sfree(diags);
 
439
 
 
440
    return NULL;
 
441
}
 
442