← Back to branch summary

~ubuntu-branches/ubuntu/precise/ncbi-tools6/precise

« back to all changes in this revision

Viewing changes to tools/rpsutil.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:
tools/rpsutil.c
1
 
/* $Id: rpsutil.c,v 6.44 2001/12/21 20:41:13 bauer Exp $
 
1
static char const rcsid[] = "$Id: rpsutil.c,v 6.71 2004/09/21 13:59:31 dondosha Exp $";
 
2
 
 
3
/* $Id: rpsutil.c,v 6.71 2004/09/21 13:59:31 dondosha Exp $
2
4
* ===========================================================================
3
5
*
4
6
*                            PUBLIC DOMAIN NOTICE
29
31
*
30
32
* Initial Version Creation Date: 12/14/1999
31
33
*
32
 
* $Revision: 6.44 $
 
34
* $Revision: 6.71 $
33
35
*
34
36
* File Description:
35
37
*         Reversed PSI BLAST utilities file
36
38
*
37
39
* $Log: rpsutil.c,v $
 
40
* Revision 6.71  2004/09/21 13:59:31  dondosha
 
41
* Use query length in protein scale for search space calculation for RPS tblastn; do not multiply search space by extra 6
 
42
*
 
43
* Revision 6.70  2004/05/13 16:58:28  kans
 
44
* in AnnotateRegionsFromCDD, do not put cdd->ShortName into comment if same as cdd->Definition
 
45
*
 
46
* Revision 6.69  2004/03/18 15:09:22  papadopo
 
47
* use the score as a tiebreaker during final sort of seqaligns in RPS blast
 
48
*
 
49
* Revision 6.68  2003/10/07 15:24:10  kans
 
50
* added FreeCDDAligns
 
51
*
 
52
* Revision 6.67  2003/09/23 15:59:20  dondosha
 
53
* Initialize global variable search_is_done; allow to reset it
 
54
*
 
55
* Revision 6.66  2003/09/22 16:24:17  merezhuk
 
56
* made RPSBlastSearchEx public.
 
57
*
 
58
* Revision 6.65  2003/09/15 18:34:33  merezhuk
 
59
* made RPSBlastSearchLight public  and declared in rpsutil.h
 
60
*
 
61
* Revision 6.64  2003/09/12 20:07:05  kans
 
62
* AnnotateRegionsFromCDD puts CDDs into their own SeqAnnot with a create date AnnotDescr
 
63
*
 
64
* Revision 6.63  2003/05/30 17:25:37  coulouri
 
65
* add rcsid
 
66
*
 
67
* Revision 6.62  2003/05/13 16:02:54  coulouri
 
68
* make ErrPostEx(SEV_FATAL, ...) exit with nonzero status
 
69
*
 
70
* Revision 6.61  2003/01/24 22:26:03  camacho
 
71
* RPSInit is deprecated, use RPSInitEx instead
 
72
*
 
73
* Revision 6.60  2003/01/22 17:09:33  camacho
 
74
* Bug fix in RPSInitEx
 
75
*
 
76
* Revision 6.59  2003/01/21 14:18:45  merezhuk
 
77
* add RPSBlastSearchLight  function to bypass RPS search and use passed results structure.
 
78
*
 
79
* Revision 6.58  2003/01/07 19:45:41  merezhuk
 
80
* add RPSBlastSearchEx to support RPS searches with Net-result-struct type of results
 
81
*
 
82
* Revision 6.57  2002/11/27 22:23:09  camacho
 
83
* Fix to allow multiple sequences to be searched
 
84
*
 
85
* Revision 6.56  2002/10/30 16:35:38  camacho
 
86
* If rpsinfo->query_slp is uninitialized, assume whole query is searched
 
87
*
 
88
* Revision 6.55  2002/10/17 20:36:00  camacho
 
89
* Disallow -L option for tblastn
 
90
*
 
91
* Revision 6.54  2002/10/10 14:49:43  camacho
 
92
*
 
93
* 1. Removed irrelevant options: -E, -G, -S, -H
 
94
* 2. Added -L option to provide range restriction on query sequence
 
95
*
 
96
* Revision 6.53  2002/09/18 20:23:20  camacho
 
97
* Added BLASTCalculateSearchSpace
 
98
*
 
99
* Revision 6.52  2002/09/17 22:18:34  bauer
 
100
* added CDD database class in RPSBgetCddHits
 
101
*
 
102
* Revision 6.51  2002/09/13 19:11:43  camacho
 
103
* Added rps_qlen field
 
104
*
 
105
* Revision 6.50  2002/09/05 15:33:34  camacho
 
106
* Fix to use the effective database length
 
107
*
 
108
* Revision 6.49  2002/08/29 13:56:55  camacho
 
109
* Restore K parameter after impala statistical correction
 
110
*
 
111
* Revision 6.48  2002/08/26 16:55:52  madden
 
112
* Fix for scaling with translated searches
 
113
*
 
114
* Revision 6.47  2002/06/11 15:13:35  dondosha
 
115
* Syntax error fixed in previous commit
 
116
*
 
117
* Revision 6.46  2002/06/11 14:44:48  dondosha
 
118
* Return status from some functions instead of search block pointer
 
119
*
 
120
* Revision 6.45  2002/03/26 16:48:54  madden
 
121
* Correctly calculate effective lengths, use correct Karlin-Altschul parameters
 
122
*
38
123
* Revision 6.44  2001/12/21 20:41:13  bauer
39
124
* changed string parsing for Pfam in RPSBgetCddHits to deal with longer Deflines
40
125
*
188
273
typedef struct _RPSSap {
189
274
    SeqAlignPtr sap;
190
275
    Nlm_FloatHi e_value;
 
276
    Int4 score;
191
277
} RPSap, PNTR RPSapPtr;
192
278
 
193
279
typedef struct _RPSapSort {
205
291
    VoidPtr print_user_data;
206
292
} RPS_MTDataStruct, PNTR RPS_MTDataStructPtr;
207
293
 
208
 
static Boolean search_is_done;
 
294
static Boolean search_is_done=FALSE;
209
295
static Int4 glb_sequence_count=0;
210
296
 
211
 
void RPSFreeLookup(RPSLookupPtr lookup)
 
297
void ResetRPS(void)
 
298
{
 
299
   search_is_done = FALSE;
 
300
}
 
301
 
 
302
int RPSCalcEValue(BlastSearchBlkPtr search, RPSInfoPtr rpsinfo,
 
303
               Uint1Ptr subject_seq, SeqLocPtr slp, BioseqPtr subject_bsp);
 
304
 
 
305
 
 
306
 
 
307
static void RPSFreeLookup(RPSLookupPtr lookup)
212
308
{
213
309
    Nlm_MemMapFini(lookup->mmLookup);
214
310
    MemFree(lookup);
216
312
    return;
217
313
}
218
314
 
219
 
RPSLookupPtr RPSInitLookup(CharPtr LookupFile)
 
315
static RPSLookupPtr RPSInitLookup(CharPtr LookupFile)
220
316
{
221
317
    RPSLookupPtr lookup;
222
318
    
223
319
    lookup = MemNew(sizeof(RPSLookup));
224
320
 
225
321
    if((lookup->mmLookup = Nlm_MemMapInit(LookupFile)) == NULL) {
226
 
        ErrPostEx(SEV_FATAL, 0, 0, "RPSInit: mmap of lookup failed");
 
322
        ErrPostEx(SEV_FATAL, 1, 0, "RPSInit: mmap of lookup failed");
227
323
        return NULL;
228
324
    }
229
325
    
232
328
    if(lookup->header[0] != RPS_MAGIC_NUMBER) {
233
329
        /* Checking for wrong architecture */
234
330
        if(Nlm_SwapUint4(lookup->header[0]) == RPS_MAGIC_NUMBER) {
235
 
            ErrPostEx(SEV_FATAL, 0, 0, "RPSInit: database formated on "
 
331
            ErrPostEx(SEV_FATAL, 1, 0, "RPSInit: database formated on "
236
332
                      "opposite big/little endian architecture. Sorry.");
237
333
            return NULL;
238
334
        } else {
239
 
            ErrPostEx(SEV_FATAL, 0, 0, "RPSInit: invalid lookup file");
 
335
            ErrPostEx(SEV_FATAL, 1, 0, "RPSInit: invalid lookup file");
240
336
            return NULL;
241
337
        }
242
338
    }
270
366
 
271
367
/*get information from matrixAuxiliaryFile
272
368
*/
273
 
static Boolean getAuxInformation(CharPtr auxFile, Nlm_FloatHiPtr scalingFactor)
 
369
static Boolean getAuxInformation(CharPtr auxFile, Nlm_FloatHiPtr scalingFactor,
 
370
        RPSInfoPtr rpsinfo, Int4Ptr gap_open, Int4Ptr gap_extend, CharPtr matrix_name)
274
371
{
275
372
    Char underlyingMatrixName[PATH_MAX]; /*name of matrix to read*/
276
373
    Nlm_FloatHi privateScalingFactor;
277
374
    Int4 maxLength; /*maximum length of sequence*/
278
375
    Int4 dbLength; /*length of database*/
 
376
    Int4 index; /* index for "for" loop. */
279
377
    Nlm_FloatHi Kungapped, Hungapped; /*two values to read*/
280
 
    Int4 gap_open, gap_extend; /*gap costs to skip over in reading*/
281
378
    FILE *matrixAuxiliaryFile;
282
379
 
283
380
    
289
386
 
290
387
    matrixAuxiliaryFile = FileOpen(auxFile, "r"); 
291
388
    fscanf(matrixAuxiliaryFile,"%s",underlyingMatrixName);
292
 
    fscanf(matrixAuxiliaryFile,"%d\n", &gap_open);
293
 
    fscanf(matrixAuxiliaryFile,"%d\n", &gap_extend);
 
389
    if (underlyingMatrixName != NULL)
 
390
        StringCpy(matrix_name, underlyingMatrixName);
 
391
    else
 
392
        matrix_name[0] = NULLB;
 
393
 
 
394
    fscanf(matrixAuxiliaryFile,"%d\n", gap_open);
 
395
    fscanf(matrixAuxiliaryFile,"%d\n", gap_extend);
294
396
    fscanf(matrixAuxiliaryFile, "%le", &Kungapped);
295
397
    fscanf(matrixAuxiliaryFile, "%le", &Hungapped);
296
398
    fscanf(matrixAuxiliaryFile, "%d", &maxLength);
297
399
    fscanf(matrixAuxiliaryFile, "%d", &dbLength);
298
400
    fscanf(matrixAuxiliaryFile, "%lf", &privateScalingFactor);
 
401
 
 
402
    rpsinfo->karlinK = MemNew(rpsinfo->matrixCount*sizeof(FloatHi));
 
403
    for (index=0; index<rpsinfo->matrixCount; index++)
 
404
    {
 
405
        fscanf(matrixAuxiliaryFile, "%d", &maxLength);
 
406
        fscanf(matrixAuxiliaryFile, "%le", &(rpsinfo->karlinK[index]));
 
407
    }
299
408
    FileClose(matrixAuxiliaryFile);
300
409
 
301
410
    if (scalingFactor)
302
411
        *scalingFactor = privateScalingFactor;
303
 
   
304
412
 
305
413
    return TRUE; 
306
414
}
307
415
 
 
416
/*** This function is DEPRECATED. Please use RPSInitEx instead (01/24/03) ***/
308
417
RPSInfoPtr RPSInit(CharPtr dbname, Int4 query_is_prot)
309
418
{
310
419
        return RPSInitEx(dbname, query_is_prot, NULL);
314
423
{
315
424
    Nlm_FloatHi scalingFactor;
316
425
    RPSInfoPtr rpsinfo;
 
426
    Int4 gap_open, gap_extend;
317
427
    Int4Ptr header;
318
 
    Char rps_matrix[PATH_MAX], rps_lookup[PATH_MAX], rps_aux[PATH_MAX];
 
428
    Char rps_matrix[PATH_MAX], rps_lookup[PATH_MAX], rps_aux[PATH_MAX], matrix_name[PATH_MAX];
 
429
    
319
430
 
320
431
    rpsinfo = MemNew(sizeof(RPSInfo));
321
432
    
333
444
    
334
445
    
335
446
    if((rpsinfo->mmMatrix = Nlm_MemMapInit(rps_matrix)) == NULL) {
336
 
        ErrPostEx(SEV_FATAL, 0, 0, "RPSInit: mmap of matrix failed");
 
447
        ErrPostEx(SEV_FATAL, 1, 0, "RPSInit: mmap of matrix failed");
337
448
        return (NULL);
338
449
    }
339
450
 
341
452
    
342
453
    if(header[0] != RPS_MAGIC_NUMBER) {
343
454
        if(Nlm_SwapUint4(header[0]) == RPS_MAGIC_NUMBER) {
344
 
            ErrPostEx(SEV_FATAL, 0, 0, "RPSInit: database formated on "
 
455
            ErrPostEx(SEV_FATAL, 1, 0, "RPSInit: database formated on "
345
456
                      "opposite big/little endian architecture. Sorry.");
346
457
            return NULL;
347
458
        } else {
348
 
            ErrPostEx(SEV_FATAL, 0, 0, "RPSInit: invalid matrix memmap file");
 
459
            ErrPostEx(SEV_FATAL, 1, 0, "RPSInit: invalid matrix memmap file");
349
460
            return NULL;
350
461
        }
351
462
    }
363
474
    rpsinfo->lookup = RPSInitLookup(rps_lookup);
364
475
 
365
476
    /* Now look in aux file. */
366
 
    getAuxInformation(rps_aux, &scalingFactor);
 
477
    getAuxInformation(rps_aux, &scalingFactor, rpsinfo, &gap_open, &gap_extend, matrix_name);
367
478
 
368
479
    if (blast_options && scalingFactor != 0.0)
369
 
        blast_options->scalingFactor = scalingFactor;
 
480
    {
 
481
        blast_options->scalingFactor = scalingFactor;
 
482
        blast_options->gap_open = gap_open;
 
483
        blast_options->gap_extend = gap_extend;
 
484
        if (matrix_name[0] != NULLB)
 
485
        {
 
486
            blast_options->matrix = MemFree(blast_options->matrix);
 
487
            blast_options->matrix = StringSave(matrix_name);
 
488
        }
 
489
    }
 
490
    if (blast_options) {
 
491
        rpsinfo->start = blast_options->required_start;
 
492
        rpsinfo->stop  = blast_options->required_end;
 
493
        /* reset kludge */
 
494
        blast_options->required_start = 0;
 
495
        blast_options->required_end = 0; 
 
496
    }
370
497
    
371
498
    return rpsinfo;
372
499
}
379
506
 
380
507
    RPSFreeLookup(rpsinfo->lookup);
381
508
 
 
509
    MemFree(rpsinfo->karlinK);
 
510
 
382
511
    MemFree(rpsinfo);
383
512
    
384
513
    ReadDBBioseqFetchDisable();
405
534
    return b;
406
535
}
407
536
 
408
 
Int4 RPSFindSequence(ReadDBFILEPtr rdfp, Int4 start, Int4Ptr sequence_start,
409
 
                     Int4Ptr seq_length)
 
537
static Int4 RPSFindSequence(ReadDBFILEPtr rdfp, Int4 start, Int4Ptr sequence_start,
 
538
                            Int4Ptr seq_length)
410
539
{
411
540
    Int4 seq_num, num_seqs;
412
541
    
421
550
    return seq_num;
422
551
}
423
552
 
424
 
RPSequencePtr RPSGetSequenceEx(RPSInfoPtr rpsinfo, Int4 seqnum, Boolean all_seqs)
 
553
static RPSequencePtr RPSGetSequenceEx(RPSInfoPtr rpsinfo, Int4 seqnum, Boolean all_seqs)
425
554
{
426
555
    RPSequencePtr rpsp;
427
556
    Int4 offset, seqlength, i, row, num_seqs;
481
610
    return RPSGetSequenceEx(rpsinfo, seqnum, FALSE);
482
611
}
483
612
 
484
 
RPSequencePtr RPSGetBIGSequence(RPSInfoPtr rpsinfo, BioseqPtr PNTR bsp_out)
 
613
static RPSequencePtr RPSGetBIGSequence(RPSInfoPtr rpsinfo, BioseqPtr PNTR bsp_out)
485
614
{
486
615
    RPSequencePtr rpsp;
 
616
#if 0
487
617
    ValNodePtr vnp;
488
618
    BioseqPtr bsp;
 
619
#endif
489
620
 
490
621
    rpsp = RPSGetSequenceEx(rpsinfo, 0, TRUE);
491
622
    
536
667
 
537
668
void RPSUpdateDbSize(BLAST_OptionsBlkPtr options, RPSInfoPtr rpsinfo, Int4 query_length)
538
669
{
539
 
 
 
670
    Int8 dblen;
 
671
    FloatHi searchsp;
 
672
 
 
673
    if (options == NULL)
 
674
        return;
 
675
 
 
676
    dblen = rpsinfo->offsets[rpsinfo->matrixCount] - (rpsinfo->matrixCount); 
 
677
    if (options->db_length != 0) 
 
678
        dblen = options->db_length;
540
679
    options->dbseq_num = rpsinfo->matrixCount;
541
680
 
542
 
    if(options->db_length == 0) {
543
 
        options->db_length = rpsinfo->offsets[rpsinfo->matrixCount];
544
 
    }
545
 
 
546
 
    if(options->searchsp_eff == 0) {
547
 
        if(rpsinfo->query_is_prot)
548
 
            options->searchsp_eff = query_length * options->db_length;
549
 
        else
550
 
            options->searchsp_eff = query_length * 6 * options->db_length;
551
 
    }
552
 
    
 
681
    /* For nucleotide query, each translated frame's length is approximately
 
682
       1/3 of the nucleotide length. */
 
683
    if (!rpsinfo->query_is_prot)
 
684
       query_length /= 3;
 
685
 
 
686
    searchsp = BLASTCalculateSearchSpace(options, options->dbseq_num, dblen,
 
687
            query_length);
 
688
 
 
689
    if(options->db_length == 0)
 
690
        options->db_length = dblen;
 
691
 
 
692
    /* Save the search space, adjusting for tblastn */
 
693
    if (options->searchsp_eff == 0) {
 
694
        options->searchsp_eff = searchsp;
 
695
        if (!rpsinfo->query_is_prot)
 
696
            options->searchsp_eff *= 3;
 
697
    }
 
698
 
553
699
    return;
554
700
}
555
701
 
584
730
    
585
731
    rpsp->sap = sap;
586
732
        
587
 
    /* Extracting e_value from SeqAlignPtr */
 
733
    /* Extracting score and e_value from SeqAlignPtr */
 
734
    rpsp->e_value = 10.0;
 
735
    rpsp->score = 0;
588
736
    thisScorePtr = sap->score;
589
 
    while ((thisScorePtr != NULL) &&
590
 
           (StringICmp(thisScorePtr->id->str, "e_value") != 0) &&
591
 
           (StringICmp(thisScorePtr->id->str, "sum_e") != 0)) {
592
 
        thisScorePtr = thisScorePtr->next;
 
737
    while (thisScorePtr != NULL) {
 
738
       if ((StringICmp(thisScorePtr->id->str, "e_value") == 0) ||
 
739
           (StringICmp(thisScorePtr->id->str, "sum_e") == 0)) {
 
740
          rpsp->e_value = (Nlm_FloatHi) (thisScorePtr->value.realvalue);
 
741
       } else {
 
742
          if (StringICmp(thisScorePtr->id->str, "score") == 0) {
 
743
             rpsp->score = thisScorePtr->value.intvalue;
 
744
          }
 
745
       }
 
746
       thisScorePtr = thisScorePtr->next;
593
747
    }
594
748
 
595
 
    if(NULL == thisScorePtr)
596
 
        rpsp->e_value = 10.0;
597
 
    else
598
 
        rpsp->e_value = (Nlm_FloatHi) (thisScorePtr->value.realvalue);
599
 
 
600
749
    ssp->rpsap[ssp->count] = rpsp;
601
750
    ssp->count++;
602
751
 
615
764
    if (rpsp1->e_value < rpsp2->e_value)
616
765
        return (-1);
617
766
    
 
767
    if (rpsp1->score < rpsp2->score)
 
768
        return (1);
 
769
    
 
770
    if (rpsp1->score > rpsp2->score)
 
771
        return (-1);
 
772
    
618
773
    return (0);    
619
774
}
620
775
 
654
809
    return;
655
810
}
656
811
 
657
 
Boolean RPSReturnQuery(BlastSearchBlkPtr search, BioseqPtr query_bsp,
658
 
                       Uint1Ptr query_seq_start)
 
812
static Boolean RPSReturnQuery(BlastSearchBlkPtr search, BioseqPtr query_bsp,
 
813
                              Uint1Ptr query_seq_start)
659
814
{
660
815
    Int4 query_length;
661
816
  
682
837
    return TRUE;
683
838
}
684
839
 
685
 
Boolean RPSubstituteQueryLookup(BlastSearchBlkPtr search, 
686
 
                                RPSequencePtr rpseq, Boolean update_lookup)
 
840
static Boolean RPSubstituteQueryLookup(BlastSearchBlkPtr search, 
 
841
                                       RPSequencePtr rpseq, Boolean update_lookup)
687
842
{
688
843
    Int4 hitlist_max;
689
844
    LookupTablePtr lookup;
696
851
    
697
852
    if(search->allocated & BLAST_SEARCH_ALLOC_QUERY_SLP)
698
853
        search->query_slp = SeqLocFree(search->query_slp);
699
 
    else 
 
854
    else {
 
855
        if (search->query_slp != NULL)
 
856
            search->query_slp = SeqLocFree(search->query_slp);
700
857
        search->allocated += BLAST_SEARCH_ALLOC_QUERY_SLP;
 
858
    }
701
859
 
702
860
    ValNodeAddPointer(&search->query_slp, SEQLOC_WHOLE, 
703
861
                      SeqIdDup(SeqIdFindBest(rpseq->seqid, SEQID_GI)));
757
915
    return TRUE;
758
916
}
759
917
 
760
 
void RPSLookupCleanUp(LookupTablePtr lookup)
 
918
static void RPSLookupCleanUp(LookupTablePtr lookup)
761
919
{
762
920
    MemFree(lookup->pv_array);
763
921
    lookup->pv_array = NULL;
778
936
 
779
937
    return;
780
938
}
781
 
void RPSExchangeInt(Int4 *a, Int4 *b)
 
939
static void RPSExchangeInt(Int4 *a, Int4 *b)
782
940
{
783
941
    Int4 value;
784
942
    
789
947
    return;
790
948
}
791
949
 
792
 
Boolean RPSUpdateCoordinates(ReadDBFILEPtr rdfp, BLASTResultHitlistPtr result,
793
 
                             Boolean reverse)
 
950
static Boolean RPSUpdateCoordinates(ReadDBFILEPtr rdfp, BLASTResultHitlistPtr result,
 
951
                                    Boolean reverse)
794
952
{
795
953
    Int4 hspcnt, index;
796
954
    BLASTResultHspPtr hsp_array;
889
1047
    return TRUE;
890
1048
}
891
1049
 
892
 
BLASTResultHitlistPtr RPSExtractNewResult(BLASTResultHitlistPtr result)
 
1050
static BLASTResultHitlistPtr RPSExtractNewResult(BLASTResultHitlistPtr result)
893
1051
{
894
1052
    BLASTResultHitlistPtr new_result;
895
1053
    Int4 i, index, hspcnt, seq_no;
960
1118
        return((Nlm_FloatHi) (thisScorePtr->value.realvalue));
961
1119
}
962
1120
 
963
 
Boolean RPSUpdateResult(BlastSearchBlkPtr search, ReadDBFILEPtr rdfp)
 
1121
 
 
1122
static Boolean RPSUpdateResult(BlastSearchBlkPtr search, ReadDBFILEPtr rdfp)
964
1123
{
965
1124
    BLASTResultsStructPtr result_struct, result_struct_new;
966
1125
    Int4 index;
1081
1240
    return TRUE;
1082
1241
}
1083
1242
 
1084
 
SeqAlignPtr RPSAlignTraceBack(BlastSearchBlkPtr search, RPSInfoPtr rpsinfo,
1085
 
                              Uint1Ptr subject_seq, SeqLocPtr slp, BioseqPtr
1086
 
                              subject_bsp)
 
1243
/* Multiplier used by Impala, should use Method of Altshul, Bundschuh, etc.? */
 
1244
#define PRO_K_MULTIPLIER 1.2
 
1245
/*
 
1246
 * modify BlastResults structure.
 
1247
 *
 
1248
 */
 
1249
static int RPSSearchPartialPrepare( BlastSearchBlkPtr search, 
 
1250
                         RPSInfoPtr rpsinfo, 
 
1251
                         Uint1Ptr subject_seq,  
 
1252
                         SeqLocPtr slp,
 
1253
                         BioseqPtr subject_bsp){
 
1254
  Int4 subject_length = 0;
 
1255
 
 
1256
  if(search->result_struct->hitlist_count == 0 || 
 
1257
     search->result_struct->results[0]->hspcnt == 0) /* No brain no pain */
 
1258
    return 0;
 
1259
  
 
1260
  if(search->result_struct->hitlist_count != 1) {
 
1261
    ErrPostEx(SEV_ERROR, 0,0, "RPSSearchPartialPrepare: This function works only for hitlist_count == 1 (%d)", 
 
1262
              search->result_struct->hitlist_count);
 
1263
    return 0;
 
1264
  }
 
1265
 
 
1266
  HeapSort(search->result_struct->results[0]->hsp_array, 
 
1267
           search->result_struct->results[0]->hspcnt, 
 
1268
           sizeof(BLASTResultHsp), 
 
1269
           RPSResultHspScoreCmp);
 
1270
    
 
1271
  subject_length = SeqLocLen(slp);
 
1272
 
 
1273
 
 
1274
  search->subject_info = BLASTSubjectInfoNew(
 
1275
                                             SeqIdDup(SeqIdFindBest(subject_bsp->id, SEQID_GI)), 
 
1276
                                             StringSave(BioseqGetTitle(subject_bsp)), subject_length);    
 
1277
 
 
1278
  RPSCalcEValue( search, rpsinfo,subject_seq, slp, subject_bsp );
 
1279
 
 
1280
 
 
1281
  return 0;
 
1282
}
 
1283
 
 
1284
/*
 
1285
 *  modified version, everything moved out of loop
 
1286
 */
 
1287
int RPSCalcEValue(BlastSearchBlkPtr search, RPSInfoPtr rpsinfo,
 
1288
               Uint1Ptr subject_seq, SeqLocPtr slp, BioseqPtr subject_bsp)
 
1289
{
 
1290
  /*   BLASTResultHspPtr hsp_array, hsp; */
 
1291
  /* BLASTResultHsp *hsp_array; */
 
1292
  BLASTResultHsp *hsp;
 
1293
  Int4 hspcnt, index /* , index1 */ ;
 
1294
  BLAST_KarlinBlkPtr PNTR kbp;
 
1295
  Nlm_FloatHi searchsp_eff;
 
1296
  RPSequencePtr rpseq;  
 
1297
  Nlm_FloatHi e_value, lambda_ratio;
 
1298
  Int4 subject_length, subject_id;
 
1299
  Int4Ptr PNTR posMatrix;
 
1300
  Char tmp[128];  
 
1301
  
 
1302
  /*   hsp_array     = search->result_struct->results[0]->hsp_array; */
 
1303
  hspcnt  = search->result_struct->results[0]->hspcnt;
 
1304
 
 
1305
  
 
1306
  kbp = search->sbp->kbp_gap;
 
1307
  searchsp_eff = search->searchsp_eff;
 
1308
 
 
1309
 
 
1310
 
 
1311
  subject_length = SeqLocLen(slp);
 
1312
  subject_id = search->result_struct->results[0]->subject_id;
 
1313
 
 
1314
  /* #0.2 */
 
1315
  search->result_struct->results[0]->subject_info=
 
1316
    BLASTSubjectInfoNew(SeqIdDup(SeqIdFindBest(subject_bsp->id, SEQID_GI)), 
 
1317
                        StringSave(BioseqGetTitle(subject_bsp)), subject_length);    
 
1318
  /* #0.3 */
 
1319
  rpseq = RPSGetSequence(rpsinfo, subject_id);
 
1320
 
 
1321
  /* #0.4 */
 
1322
  /* (!!!!!) It looks like this function is called second time 
 
1323
  RPSubstituteQueryLookup(search, rpseq, FALSE);
 
1324
  */
 
1325
 
 
1326
  /* #0.5 */
 
1327
  if (rpsinfo->query_is_prot == TRUE)  {
 
1328
    if(!RPSimpalaStatCorrections(rpseq, &lambda_ratio, search->pbp->scalingFactor, 
 
1329
                                 subject_length, subject_seq)) {
 
1330
 
 
1331
      SeqIdWrite(rpseq->seqid, tmp, PRINTID_FASTA_LONG, sizeof(tmp));   
 
1332
      ErrPostEx(SEV_WARNING, 0, rpseq->number, 
 
1333
                "BAD matrix for the sequence %s %s",
 
1334
                tmp, rpseq->description);
 
1335
      RPSequenceFree(rpseq);      
 
1336
      return 1; 
 
1337
    }
 
1338
    /* #0.5.1 */
 
1339
    if (rpsinfo->karlinK) {
 
1340
      search->sbp->kbp_gap[0]->K = PRO_K_MULTIPLIER*rpsinfo->karlinK[ subject_id ];
 
1341
      search->sbp->kbp_gap[0]->logK = log(search->sbp->kbp_gap[0]->K);
 
1342
      /* scalingFactor should only be different than one (or zero) if
 
1343
         the aux file was read and rpsinfo->karlinK was filled in. */
 
1344
    }
 
1345
    /* #ifdef DEALLOC_MATRIX *
 
1346
    /* #0.5.2 */
 
1347
    if (rpseq->copyMatrix) {
 
1348
      posMatrix = search->sbp->posMatrix;
 
1349
      search->sbp->posMatrix = rpseq->copyMatrix;
 
1350
    }
 
1351
    /* #endif */
 
1352
  }
 
1353
 
 
1354
 
 
1355
  if (search->pbp->scalingFactor != 0) search->sbp->kbp_gap[0]->Lambda /= search->pbp->scalingFactor;
 
1356
  kbp = search->sbp->kbp_gap;
 
1357
 
 
1358
  for (index=0; index<hspcnt; index++) {
 
1359
    /* #1 */
 
1360
    hsp = &(search->result_struct->results[0]->hsp_array[ index] ); /* sp_array[index]; */
 
1361
 
 
1362
    /*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
 
1363
    /* CALCULATE E-VALUE HERE */    
 
1364
    e_value = BlastKarlinStoE_simple(hsp->score, kbp[0],   searchsp_eff);
 
1365
    hsp->e_value = e_value;
 
1366
    /*     fprintf(stderr,"score: %d e-value: %g\n",hsp->score,hsp->e_value); */
 
1367
    /*\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
 
1368
    
 
1369
  } /* end of for cycle */
 
1370
  
 
1371
  /* RESTORE DATA */
 
1372
  if (rpsinfo->query_is_prot == TRUE)  {
 
1373
    /* #ifdef DEALLOC_MATRIX */
 
1374
    /* Deallocate copyMatrix to prevent a leak. */
 
1375
    if (rpseq->copyMatrix) {
 
1376
/*       for (index1=0; index1<(1+rpseq->seqlen); index1++) */
 
1377
/*      MemFree(rpseq->copyMatrix[index1]); */
 
1378
 
 
1379
      rpseq->copyMatrix = MemFree(rpseq->copyMatrix);
 
1380
      search->sbp->posMatrix = posMatrix;
 
1381
    }
 
1382
    /* #endif */
 
1383
    if (rpsinfo->karlinK) {
 
1384
      /* Restore K */
 
1385
      search->sbp->kbp_gap[0]->K /= PRO_K_MULTIPLIER;
 
1386
      search->sbp->kbp_gap[0]->logK = log(search->sbp->kbp_gap[0]->K);
 
1387
    }    
 
1388
  }
 
1389
 
 
1390
  if (search->pbp->scalingFactor != 0) search->sbp->kbp_gap[0]->Lambda *= search->pbp->scalingFactor;  
 
1391
 
 
1392
 
 
1393
  RPSequenceFree(rpseq);      
 
1394
 
 
1395
  return 0;
 
1396
}
 
1397
 
 
1398
 
 
1399
 
 
1400
static SeqAlignPtr RPSAlignTraceBack(BlastSearchBlkPtr search, RPSInfoPtr rpsinfo,
 
1401
                                     Uint1Ptr subject_seq, SeqLocPtr slp, BioseqPtr
 
1402
                                     subject_bsp)
1087
1403
{
1088
1404
    SeqAlignPtr seqalign;
1089
1405
    BLASTResultsStructPtr result_struct, result_struct_new;
1156
1472
    }
1157
1473
 
1158
1474
    /* ---------------------- */
 
1475
 
1159
1476
    for(index = 0; index < result_struct_new->hitlist_max; index++) {
1160
1477
        
1161
1478
        new_result = RPSExtractNewResult(result_struct->results[0]);
1193
1510
                continue;
1194
1511
            }
1195
1512
 
 
1513
            if (rpsinfo->karlinK)
 
1514
            {
 
1515
                search->sbp->kbp_gap[0]->K = PRO_K_MULTIPLIER*rpsinfo->karlinK[new_result->subject_id];
 
1516
                search->sbp->kbp_gap[0]->logK = log(search->sbp->kbp_gap[0]->K);
 
1517
                /* scalingFactor should only be different than one (or zero) if
 
1518
                        the aux file was read and rpsinfo->karlinK was filled in. */
 
1519
            }
1196
1520
 
1197
1521
            if (rpseq->copyMatrix) {
1198
1522
                posMatrix = search->sbp->posMatrix;
1200
1524
            }
1201
1525
        }
1202
1526
 
 
1527
        if (search->pbp->scalingFactor != 0) 
 
1528
                search->sbp->kbp_gap[0]->Lambda /= search->pbp->scalingFactor;
 
1529
 
1203
1530
        if(rpsinfo->query_is_prot) {
1204
1531
 
1205
1532
            if(search->pbp->gapped_calculation == TRUE) {
1225
1552
                for (index1=0; index1<(1+rpseq->seqlen); index1++)
1226
1553
                        MemFree(rpseq->copyMatrix[index1]);
1227
1554
                rpseq->copyMatrix = MemFree(rpseq->copyMatrix);
1228
 
                search->sbp->posMatrix = posMatrix;
 
1555
        search->sbp->posMatrix = posMatrix;
 
1556
        if (rpsinfo->karlinK) {
 
1557
            /* Restore K */
 
1558
            search->sbp->kbp_gap[0]->K /= PRO_K_MULTIPLIER;
 
1559
            search->sbp->kbp_gap[0]->logK = log(search->sbp->kbp_gap[0]->K);
 
1560
        }
1229
1561
        }
1230
1562
        
1231
1563
        if(seqalign != NULL && (e_value = getEvalueFromSeqAlign(seqalign)) < 
1313
1645
    return;
1314
1646
}
1315
1647
#endif
 
1648
 
 
1649
/* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
 
1650
#ifdef USE_REMOVED
 
1651
SeqAlignPtr 
 
1652
RPS_FIX(BlastSearchBlkPtr search, RPSInfoPtr rpsinfo,
 
1653
        Uint1Ptr subject_seq, SeqLocPtr slp, BioseqPtr
 
1654
        subject_bsp)
 
1655
{ /* ^^^A */
 
1656
    SeqAlignPtr seqalign;
 
1657
    /* BLASTResultsStructPtr result_struct; */
 
1658
    BLASTResultsStructPtr result_struct_new;
 
1659
    Int4 index, index1, subject_length, rev_subject_length;
 
1660
    Int4Ptr PNTR posMatrix;
 
1661
    BLASTResultHitlistPtr new_result;
 
1662
    RPSapSortPtr rpssp;
 
1663
    RPSequencePtr rpseq;
 
1664
    Nlm_FloatHi e_value, lambda_ratio;
 
1665
    Uint1Ptr subject = NULL, rev_subject = NULL;
 
1666
    SeqPortPtr spp;
 
1667
 
 
1668
 
 
1669
    fprintf(stderr,"[^^^] RPS_FIX -- START\n");
 
1670
 
 
1671
    /*     result_struct = search->result_struct; */
 
1672
 
 
1673
    if(search->result_struct->hitlist_count == 0 || search-> result_struct->results[0]->hspcnt == 0) /* No brain no pain */
 
1674
      return NULL;
 
1675
 
 
1676
    
 
1677
    result_struct_new = BLASTResultsStructNew( search->result_struct->results[0]->hspcnt, search->pbp->max_pieces, search->pbp->hsp_range_max);
 
1678
    /* ?moved up ? */
 
1679
    if(search->result_struct->hitlist_count != 1) {
 
1680
        ErrPostEx(SEV_ERROR, 0,0, "RPSUpdateResult: This function works only for hitlist_count == 1 (%d)", search->result_struct->hitlist_count);
 
1681
        return NULL;
 
1682
    }
 
1683
 
 
1684
    HeapSort(search->result_struct->results[0]->hsp_array, 
 
1685
             search->result_struct->results[0]->hspcnt, 
 
1686
             sizeof(BLASTResultHsp), RPSResultHspScoreCmp);
 
1687
    
 
1688
    subject_length = SeqLocLen(slp);
 
1689
 
 
1690
    search->result_struct = result_struct_new;
 
1691
 
 
1692
    search->subject_info = BLASTSubjectInfoNew(SeqIdDup(SeqIdFindBest(subject_bsp->id, SEQID_GI)), StringSave(BioseqGetTitle(subject_bsp)), subject_length);    
 
1693
 
 
1694
 
 
1695
 
 
1696
    rpssp =  RPSapSortInit();
 
1697
    
 
1698
    /* First we have to update coordinates */
 
1699
    RPSUpdateCoordinates(rpsinfo->rdfp, search->result_struct->results[0], FALSE);
 
1700
    
 
1701
    /* Now we will extract correct "results" one by one 
 
1702
       and make alignment then */
 
1703
 
 
1704
    result_struct_new->hitlist_count = 1; /* Always will be 1 */
 
1705
    
 
1706
    /* We have to create 4na buffer for our query, that is subject */
 
1707
    if(!rpsinfo->query_is_prot) {
 
1708
        spp = SeqPortNew(subject_bsp, 0, -1, Seq_strand_plus, 
 
1709
                         Seq_code_ncbi4na);
 
1710
        /* make one longer to "protect" ALIGN. */
 
1711
        subject = MemNew((1+subject_bsp->length)*sizeof(Uint1));
 
1712
        for (index1 = 0; index1 < subject_bsp->length; index1++) {
 
1713
            subject[index1] = SeqPortGetResidue(spp);
 
1714
        }
 
1715
        
 
1716
        /* Gap character in last space. */
 
1717
        subject[subject_bsp->length] = NULLB;
 
1718
        subject_length = subject_bsp->length;
 
1719
        spp = SeqPortFree(spp);
 
1720
        
 
1721
        spp = SeqPortNew(subject_bsp, 0, -1, Seq_strand_minus, Seq_code_ncbi4na);
 
1722
        /* make one longer to "protect" ALIGN. */
 
1723
        rev_subject = MemNew((1+subject_bsp->length)*sizeof(Uint1));
 
1724
        for (index1=0; index1<subject_bsp->length; index1++) {
 
1725
            rev_subject[index1] = SeqPortGetResidue(spp);
 
1726
        }
 
1727
        /* Gap character in last space. */
 
1728
        rev_subject[subject_bsp->length] = NULLB;
 
1729
        rev_subject_length = subject_bsp->length;
 
1730
        spp = SeqPortFree(spp);
 
1731
    }
 
1732
 
 
1733
    /* ---------------------- */
 
1734
 
 
1735
    for(index = 0; index < result_struct_new->hitlist_max; index++) {
 
1736
        
 
1737
        new_result = RPSExtractNewResult(search->result_struct->results[0]);
 
1738
 
 
1739
        if(new_result == NULL)    break;
 
1740
 
 
1741
        new_result->subject_info = BLASTSubjectInfoNew(SeqIdDup(SeqIdFindBest(subject_bsp->id, SEQID_GI)), 
 
1742
                                                       StringSave(BioseqGetTitle(subject_bsp)), subject_length);        
 
1743
        result_struct_new->results[0] = new_result;
 
1744
        
 
1745
        /* subject_id == query_id in rdfp database */
 
1746
        rpseq = RPSGetSequence(rpsinfo, new_result->subject_id);
 
1747
        
 
1748
        RPSubstituteQueryLookup(search, rpseq, FALSE);
 
1749
        /* Correct statistics for proteins, should also be doen for nucleotides. */
 
1750
        if (rpsinfo->query_is_prot == TRUE) {
 
1751
            if(!RPSimpalaStatCorrections(rpseq, &lambda_ratio, search->pbp->scalingFactor, 
 
1752
                                         subject_length, subject_seq)) {
 
1753
                Char tmp[64];
 
1754
 
 
1755
                /* If this function failed - and returned FALSE - matrix
 
1756
                   used for this sequence is incorrect and results may not
 
1757
                   be trusted. We will post WARNING message and skip this
 
1758
                   profile completatly */
 
1759
 
 
1760
                SeqIdWrite(rpseq->seqid, tmp, PRINTID_FASTA_LONG, sizeof(tmp));
 
1761
 
 
1762
                ErrPostEx(SEV_WARNING, 0, rpseq->number, 
 
1763
                          "BAD matrix for the sequence %s %s",
 
1764
                          tmp, rpseq->description);
 
1765
                
 
1766
                seqalign = NULL; 
 
1767
                RPSequenceFree(rpseq);
 
1768
                BLASTResultHitlistFree(new_result);
 
1769
                result_struct_new->results[0] = NULL;
 
1770
                continue;
 
1771
            }
 
1772
 
 
1773
            if (rpsinfo->karlinK)
 
1774
            {
 
1775
                search->sbp->kbp_gap[0]->K = PRO_K_MULTIPLIER*rpsinfo->karlinK[new_result->subject_id];
 
1776
                search->sbp->kbp_gap[0]->logK = log(search->sbp->kbp_gap[0]->K);
 
1777
                /* scalingFactor should only be different than one (or zero) if
 
1778
                        the aux file was read and rpsinfo->karlinK was filled in. */
 
1779
            }
 
1780
 
 
1781
            if (rpseq->copyMatrix) {
 
1782
                posMatrix = search->sbp->posMatrix;
 
1783
                search->sbp->posMatrix = rpseq->copyMatrix;
 
1784
            }
 
1785
        }
 
1786
 
 
1787
        if (search->pbp->scalingFactor != 0) 
 
1788
                search->sbp->kbp_gap[0]->Lambda /= search->pbp->scalingFactor;
 
1789
 
 
1790
        if(rpsinfo->query_is_prot) {
 
1791
 
 
1792
            if(search->pbp->gapped_calculation == TRUE) {
 
1793
                seqalign = BlastGetGapAlgnTbck(search, 0, TRUE, FALSE, 
 
1794
                                               subject_seq, 
 
1795
                                               subject_length, NULL, 0);
 
1796
                
 
1797
            } else {
 
1798
                seqalign = GetSeqAlignForResultHitList(search, TRUE, FALSE, 
 
1799
                                                       search->pbp->discontinuous, 
 
1800
                                                       TRUE, FALSE);
 
1801
            }
 
1802
        } else { /* tblastn */
 
1803
            
 
1804
            seqalign = BlastGetGapAlgnTbck (search, 0, TRUE, FALSE, 
 
1805
                                            subject, subject_length, 
 
1806
                                            rev_subject, rev_subject_length);
 
1807
        }
 
1808
 
 
1809
        /* Deallocate copyMatrix to prevent a leak. */
 
1810
        if (rpseq->copyMatrix)
 
1811
        {
 
1812
                for (index1=0; index1<(1+rpseq->seqlen); index1++)
 
1813
                        MemFree(rpseq->copyMatrix[index1]);
 
1814
                rpseq->copyMatrix = MemFree(rpseq->copyMatrix);
 
1815
        search->sbp->posMatrix = posMatrix;
 
1816
        if (rpsinfo->karlinK) {
 
1817
            /* Restore K */
 
1818
            search->sbp->kbp_gap[0]->K /= PRO_K_MULTIPLIER;
 
1819
            search->sbp->kbp_gap[0]->logK = log(search->sbp->kbp_gap[0]->K);
 
1820
        }
 
1821
        }
 
1822
        
 
1823
        if(seqalign != NULL && (e_value = getEvalueFromSeqAlign(seqalign)) < 
 
1824
           search->pbp->cutoff_e) {
 
1825
            AdjustOffSetsInSeqAlign(seqalign, slp, search->query_slp);
 
1826
            RPSAddSap(rpssp, seqalign);
 
1827
        } else {
 
1828
            SeqAlignSetFree(seqalign);
 
1829
        }
 
1830
        
 
1831
        seqalign = NULL; 
 
1832
        
 
1833
        RPSequenceFree(rpseq);
 
1834
        BLASTResultHitlistFree(new_result);
 
1835
       
 
1836
        result_struct_new->results[0] = NULL;
 
1837
    }
 
1838
 
 
1839
    result_struct_new->hitlist_count = 0;
 
1840
    
 
1841
    /* ??? -- BLASTResultsStructDelete(result_struct); */
 
1842
 
 
1843
    seqalign = RPSReadSapSort(rpssp); 
 
1844
    
 
1845
    MemFree(rev_subject);
 
1846
    MemFree(subject);
 
1847
 
 
1848
    RPSapSortFree(rpssp);
 
1849
 
 
1850
  fprintf(stderr,"[^^^] RPS_FIX -- END\n");
 
1851
    return seqalign;
 
1852
 
 
1853
 
 
1854
 
 
1855
 
 
1856
}
 
1857
#endif
 
1858
 
 
1859
/* \\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\\ */
 
1860
/*
 
1861
 * support RPS searches with Net-result-struct type of results
 
1862
 */
 
1863
SeqAlignPtr RPSBlastSearchEx (BlastSearchBlkPtr search,
 
1864
                                    BioseqPtr query_bsp, RPSInfoPtr rpsinfo, Boolean partial)
 
1865
{
 
1866
    Int2 status;
 
1867
    Int4 index;
 
1868
    SeqAlignPtr seqalign=NULL, head = NULL, seqalign_var = NULL;
 
1869
    SeqPortPtr spp;
 
1870
    Uint1Ptr subject_seq, subject_seq_start;
 
1871
    Uint1 residue;
 
1872
    Int4 subject_length;
 
1873
    SeqLocPtr slp = NULL;
 
1874
    BioseqPtr subject_bsp;
 
1875
    RPSequencePtr rpseq;
 
1876
    BioseqPtr bsp = NULL;
 
1877
    SPCompressPtr spc = NULL;
 
1878
 
 
1879
    /* call normal RPS search if this is not SPLITDB */
 
1880
    if( !partial ) return  RPSBlastSearch(search,query_bsp,rpsinfo);
 
1881
 
 
1882
 
 
1883
    
 
1884
    if (search == NULL || search->query_invalid)
 
1885
        return NULL;
 
1886
    
 
1887
 
 
1888
    subject_bsp = query_bsp;    /* Reversed ... */
 
1889
    
 
1890
    if (rpsinfo->query_slp == NULL || 
 
1891
        SeqLocLen(rpsinfo->query_slp) == subject_bsp->length ||
 
1892
        ISA_na(subject_bsp->mol)) {
 
1893
        ValNodeAddPointer(&slp, SEQLOC_WHOLE, SeqIdDup(subject_bsp->id));
 
1894
    } else {
 
1895
        slp = SeqLocIntNew(SeqLocStart(rpsinfo->query_slp),
 
1896
                SeqLocStop(rpsinfo->query_slp),
 
1897
                SeqLocStrand(rpsinfo->query_slp),
 
1898
                SeqIdDup(subject_bsp->id));
 
1899
    }
 
1900
    subject_length = SeqLocLen(slp);
 
1901
    
 
1902
    search->subject_info = BLASTSubjectInfoDestruct(search->subject_info);
 
1903
 
 
1904
    /* Extracting sequence from the Bioseq */
 
1905
    
 
1906
    if(rpsinfo->query_is_prot) {         /* For blastp search */
 
1907
        subject_seq_start = subject_seq = NULL;
 
1908
        
 
1909
        subject_seq_start = (Uint1Ptr) MemNew(((subject_bsp->length)+2)*sizeof(Uint1));
 
1910
        /* The first residue is the sentinel. */
 
1911
        subject_seq_start[0] = NULLB;
 
1912
        subject_seq = subject_seq_start+1;
 
1913
        index = 0;
 
1914
        spp = SeqPortNewByLoc(slp, Seq_code_ncbistdaa);
 
1915
        while ((residue=SeqPortGetResidue(spp)) != SEQPORT_EOF) {
 
1916
            if (IS_residue(residue)) {
 
1917
                subject_seq[index++] = residue;
 
1918
            }
 
1919
        }
 
1920
        subject_seq[index] = NULLB;
 
1921
        spp = SeqPortFree(spp);
 
1922
    } else { /* tblastn */
 
1923
        spp = SeqPortNewByLoc(slp, Seq_code_ncbi4na);
 
1924
        subject_bsp = BioseqFindCore(SeqLocId(slp));
 
1925
        if (subject_bsp != NULL && subject_bsp->repr == Seq_repr_delta)
 
1926
            SeqPortSet_do_virtual(spp, TRUE);
 
1927
        spc = SPCompressDNA(spp);
 
1928
        if (spc == NULL)
 
1929
            return NULL;
 
1930
        subject_seq_start = subject_seq = spc->buffer;
 
1931
        spp = SeqPortFree(spp);
 
1932
 
 
1933
        /* Adjusting translation buffer */
 
1934
        spc->lbytes = MemFree(spc->lbytes);
 
1935
        MemFree(spc);
 
1936
        MemFree(search->translation_buffer);
 
1937
        search->translation_buffer = MemNew((3+(subject_length/3))*sizeof(Uint1));
 
1938
        search->translation_buffer_size = 1 + (subject_length/3);
 
1939
    }
 
1940
    
 
1941
    if(search->context[0].query->sequence_start != NULL) {
 
1942
        MemFree(search->context[0].query->sequence_start);
 
1943
        search->context[0].query->sequence_start = NULL;
 
1944
    }
 
1945
 
 
1946
    /* Cleaning up word finder */
 
1947
 
 
1948
/*
 
1949
    RPSLookupCleanUp(search->wfp->lookup);
 
1950
*/
 
1951
    
 
1952
 
 
1953
    search->subject_info = BLASTSubjectInfoNew(SeqIdDup(SeqIdFindBest(subject_bsp->id, SEQID_GI)), StringSave(BioseqGetTitle(subject_bsp)), subject_length);        
 
1954
    rpseq = RPSGetBIGSequence(rpsinfo, &bsp);
 
1955
 
 
1956
    /*    RRRPrintStatistics(rpseq); */
 
1957
  
 
1958
    RPSubstituteQueryLookup(search, rpseq, TRUE);
 
1959
 
 
1960
    /* Filter the subject sequence, that is actually query */    
 
1961
    if(rpsinfo->query_is_prot && search->mask != NULL) {
 
1962
        BlastMaskTheResidues(subject_seq, subject_length, 21, 
 
1963
                             search->mask->data.ptrvalue, FALSE, 0);
 
1964
    }
 
1965
 
 
1966
    /* Performing real search here */    
 
1967
    status = BLASTPerformSearch(search, subject_length, subject_seq);
 
1968
 
 
1969
    if (status) {
 
1970
       /* This is a fatal error, must exit now */
 
1971
       search = BlastSearchBlkDestruct(search);
 
1972
       return NULL;
 
1973
    }
 
1974
 
 
1975
    if(!rpsinfo->query_is_prot) {
 
1976
        /* status = BlastLinkHsps(search); */
 
1977
    } else if (search->pbp->gapped_calculation == FALSE) {
 
1978
        if (search->pbp->do_sum_stats == TRUE)
 
1979
            status = BlastLinkHsps(search);
 
1980
        else
 
1981
            status = BlastGetNonSumStatsEvalue(search);
 
1982
    } 
 
1983
    
 
1984
    status = BlastReapHitlistByEvalue(search);
 
1985
    
 
1986
    BlastSaveCurrentHitlist(search);
 
1987
 
 
1988
    /* call this function to re-arrange results to be more
 
1989
       compatible with BLASTResultHitlistToNetResultHitlist
 
1990
     */
 
1991
    search->rdfp = rpsinfo->rdfp;
 
1992
    RPSSearchPartialPrepare( search, rpsinfo, subject_seq,    slp, subject_bsp);   
 
1993
 
 
1994
 
 
1995
    return seqalign;
 
1996
 
 
1997
}
 
1998
/* 
 
1999
 * bypass RPS search and use passed results structure. 
 
2000
 * used in SPLITD
 
2001
 */
 
2002
SeqAlignPtr RPSBlastSearchLight (BlastSearchBlkPtr search,
 
2003
                            BioseqPtr query_bsp, RPSInfoPtr rpsinfo, 
 
2004
                                 BLASTResultsStructPtr   ready_result_struct){
 
2005
    Int2 status;
 
2006
    Int4 index;
 
2007
    SeqAlignPtr seqalign=NULL, head = NULL, seqalign_var = NULL;
 
2008
    SeqPortPtr spp;
 
2009
    Uint1Ptr subject_seq, subject_seq_start;
 
2010
    Uint1 residue;
 
2011
    Int4 subject_length;
 
2012
    SeqLocPtr slp = NULL;
 
2013
    BioseqPtr subject_bsp;
 
2014
    RPSequencePtr rpseq;
 
2015
    BioseqPtr bsp = NULL;
 
2016
    SPCompressPtr spc = NULL;
 
2017
    ValNodePtr vnp;
 
2018
    SeqLocPtr filter_slp;
 
2019
    
 
2020
    if (search == NULL || search->query_invalid)
 
2021
        return NULL;
 
2022
    
 
2023
 
 
2024
    subject_bsp = query_bsp;    /* Reversed ... */
 
2025
    
 
2026
    if (rpsinfo->query_slp == NULL || 
 
2027
        SeqLocLen(rpsinfo->query_slp) == subject_bsp->length ||
 
2028
        ISA_na(subject_bsp->mol)) {
 
2029
        ValNodeAddPointer(&slp, SEQLOC_WHOLE, SeqIdDup(subject_bsp->id));
 
2030
    } else {
 
2031
        slp = SeqLocIntNew(SeqLocStart(rpsinfo->query_slp),
 
2032
                SeqLocStop(rpsinfo->query_slp),
 
2033
                SeqLocStrand(rpsinfo->query_slp),
 
2034
                SeqIdDup(subject_bsp->id));
 
2035
    }
 
2036
    subject_length = SeqLocLen(slp);
 
2037
    
 
2038
    search->subject_info = BLASTSubjectInfoDestruct(search->subject_info);
 
2039
 
 
2040
    /* Extracting sequence from the Bioseq */
 
2041
    
 
2042
    if(rpsinfo->query_is_prot) {         /* For blastp search */
 
2043
        subject_seq_start = subject_seq = NULL;
 
2044
        
 
2045
        subject_seq_start = (Uint1Ptr) MemNew(((subject_bsp->length)+2)*sizeof(Uint1));
 
2046
        /* The first residue is the sentinel. */
 
2047
        subject_seq_start[0] = NULLB;
 
2048
        subject_seq = subject_seq_start+1;
 
2049
        index = 0;
 
2050
        spp = SeqPortNewByLoc(slp, Seq_code_ncbistdaa);
 
2051
        while ((residue=SeqPortGetResidue(spp)) != SEQPORT_EOF) {
 
2052
            if (IS_residue(residue)) {
 
2053
                subject_seq[index++] = residue;
 
2054
            }
 
2055
        }
 
2056
        subject_seq[index] = NULLB;
 
2057
        spp = SeqPortFree(spp);
 
2058
    } else { /* tblastn */
 
2059
        spp = SeqPortNewByLoc(slp, Seq_code_ncbi4na);
 
2060
        subject_bsp = BioseqFindCore(SeqLocId(slp));
 
2061
        if (subject_bsp != NULL && subject_bsp->repr == Seq_repr_delta)
 
2062
            SeqPortSet_do_virtual(spp, TRUE);
 
2063
        spc = SPCompressDNA(spp);
 
2064
        if (spc == NULL)
 
2065
            return NULL;
 
2066
        subject_seq_start = subject_seq = spc->buffer;
 
2067
        spp = SeqPortFree(spp);
 
2068
 
 
2069
        /* Adjusting translation buffer */
 
2070
        spc->lbytes = MemFree(spc->lbytes);
 
2071
        MemFree(spc);
 
2072
        MemFree(search->translation_buffer);
 
2073
        search->translation_buffer = MemNew((3+(subject_length/3))*sizeof(Uint1));
 
2074
        search->translation_buffer_size = 1 + (subject_length/3);
 
2075
    }
 
2076
    
 
2077
    if(search->context[0].query->sequence_start != NULL) {
 
2078
        MemFree(search->context[0].query->sequence_start);
 
2079
        search->context[0].query->sequence_start = NULL;
 
2080
    }
 
2081
 
 
2082
    /* Cleaning up word finder */
 
2083
 
 
2084
/*
 
2085
    RPSLookupCleanUp(search->wfp->lookup);
 
2086
*/
 
2087
    
 
2088
    search->subject_info = BLASTSubjectInfoNew(SeqIdDup(SeqIdFindBest(subject_bsp->id, SEQID_GI)), StringSave(BioseqGetTitle(subject_bsp)), subject_length);        
 
2089
    rpseq = RPSGetBIGSequence(rpsinfo, &bsp);
 
2090
 
 
2091
    /*    RRRPrintStatistics(rpseq); */
 
2092
  
 
2093
    RPSubstituteQueryLookup(search, rpseq, TRUE);
 
2094
 
 
2095
    /* Filter the subject sequence, that is actually query */    
 
2096
    if(rpsinfo->query_is_prot && search->mask != NULL) {
 
2097
        BlastMaskTheResidues(subject_seq, subject_length, 21, 
 
2098
                             search->mask->data.ptrvalue, FALSE, 0);
 
2099
    }
 
2100
 
 
2101
    /* NOT Performing real search here, but assigning ready results  */    
 
2102
    search->result_struct = ready_result_struct ;
 
2103
#ifdef USE_REMOVED
 
2104
    status = BLASTPerformSearch(search, subject_length, subject_seq) ; 
 
2105
    
 
2106
    if ( status ) {
 
2107
       /* This is a fatal error, must exit now */
 
2108
       search = BlastSearchBlkDestruct(search);
 
2109
       return NULL;
 
2110
    }
 
2111
#endif
 
2112
 
 
2113
    if(!rpsinfo->query_is_prot) {
 
2114
        /* status = BlastLinkHsps(search); */
 
2115
    } else if (search->pbp->gapped_calculation == FALSE) {
 
2116
        if (search->pbp->do_sum_stats == TRUE)
 
2117
            status = BlastLinkHsps(search);
 
2118
        else
 
2119
            status = BlastGetNonSumStatsEvalue(search);
 
2120
    }
 
2121
    
 
2122
    status = BlastReapHitlistByEvalue(search);
 
2123
    
 
2124
    BlastSaveCurrentHitlist(search);
 
2125
 
 
2126
    seqalign = RPSAlignTraceBack(search, rpsinfo, subject_seq, 
 
2127
                                 slp, subject_bsp); 
 
2128
    
 
2129
    /* Now converting protein filter SeqLoc-s into DNA seqlocs */
 
2130
    if(!rpsinfo->query_is_prot) {  
 
2131
        for(vnp = search->mask; vnp != NULL; vnp = vnp->next) {
 
2132
            filter_slp = (SeqLocPtr) vnp->data.ptrvalue;
 
2133
            BlastConvertProteinSeqLoc(filter_slp, DefineToFrame(vnp->choice), 
 
2134
                                      subject_bsp->length);
 
2135
        }
 
2136
    }
 
2137
 
 
2138
    search->context[0].query->sequence_start = NULL;
 
2139
    search->context[0].query->sequence = NULL;
 
2140
 
 
2141
    slp = SeqLocFree(slp);
 
2142
 
 
2143
    MemFree(subject_seq_start);
 
2144
    RPSequenceFree(rpseq);
 
2145
    BioseqFree(bsp);
 
2146
 
 
2147
    search->sbp->posMatrix = NULL;
 
2148
    search->wfp->lookup->mod_lt = NULL;
 
2149
    search->wfp->lookup->mod_lookup_table_memory = NULL;
 
2150
 
 
2151
    return seqalign;  
 
2152
}
 
2153
 
 
2154
 
1316
2155
SeqAlignPtr RPSBlastSearch (BlastSearchBlkPtr search,
1317
2156
                            BioseqPtr query_bsp, RPSInfoPtr rpsinfo)
1318
2157
{
1323
2162
    Uint1Ptr subject_seq, subject_seq_start;
1324
2163
    Uint1 residue;
1325
2164
    Int4 subject_length;
1326
 
    SeqLocPtr slp;
 
2165
    SeqLocPtr slp = NULL;
1327
2166
    BioseqPtr subject_bsp;
1328
2167
    RPSequencePtr rpseq;
1329
2168
    BioseqPtr bsp = NULL;
1337
2176
 
1338
2177
    subject_bsp = query_bsp;    /* Reversed ... */
1339
2178
    
1340
 
    slp = NULL;
1341
 
    ValNodeAddPointer(&slp, SEQLOC_WHOLE, 
1342
 
                      SeqIdDup(SeqIdFindBest(subject_bsp->id, SEQID_GI)));
1343
 
    
 
2179
    if (rpsinfo->query_slp == NULL || 
 
2180
        SeqLocLen(rpsinfo->query_slp) == subject_bsp->length ||
 
2181
        ISA_na(subject_bsp->mol)) {
 
2182
        ValNodeAddPointer(&slp, SEQLOC_WHOLE, SeqIdDup(subject_bsp->id));
 
2183
    } else {
 
2184
        slp = SeqLocIntNew(SeqLocStart(rpsinfo->query_slp),
 
2185
                SeqLocStop(rpsinfo->query_slp),
 
2186
                SeqLocStrand(rpsinfo->query_slp),
 
2187
                SeqIdDup(subject_bsp->id));
 
2188
    }
1344
2189
    subject_length = SeqLocLen(slp);
1345
2190
    
1346
2191
    search->subject_info = BLASTSubjectInfoDestruct(search->subject_info);
1358
2203
        spp = SeqPortNewByLoc(slp, Seq_code_ncbistdaa);
1359
2204
        while ((residue=SeqPortGetResidue(spp)) != SEQPORT_EOF) {
1360
2205
            if (IS_residue(residue)) {
1361
 
                subject_seq[index] = residue;
1362
 
                index++;
 
2206
                subject_seq[index++] = residue;
1363
2207
            }
1364
2208
        }
1365
2209
        subject_seq[index] = NULLB;
1408
2252
    }
1409
2253
 
1410
2254
    /* Performing real search here */    
1411
 
    search = BLASTPerformSearch(search, subject_length, subject_seq);
 
2255
    status = BLASTPerformSearch(search, subject_length, subject_seq);
 
2256
 
 
2257
    if (status) {
 
2258
       /* This is a fatal error, must exit now */
 
2259
       search = BlastSearchBlkDestruct(search);
 
2260
       return NULL;
 
2261
    }
1412
2262
 
1413
2263
    if(!rpsinfo->query_is_prot) {
1414
2264
        /* status = BlastLinkHsps(search); */
1422
2272
    status = BlastReapHitlistByEvalue(search);
1423
2273
    
1424
2274
    BlastSaveCurrentHitlist(search);
1425
 
    
 
2275
 
1426
2276
    seqalign = RPSAlignTraceBack(search, rpsinfo, subject_seq, 
1427
2277
                                 slp, subject_bsp); 
1428
2278
    
1438
2288
    search->context[0].query->sequence_start = NULL;
1439
2289
    search->context[0].query->sequence = NULL;
1440
2290
 
 
2291
    slp = SeqLocFree(slp);
 
2292
 
1441
2293
    MemFree(subject_seq_start);
1442
2294
    RPSequenceFree(rpseq);
1443
2295
    BioseqFree(bsp);
1446
2298
    search->wfp->lookup->mod_lt = NULL;
1447
2299
    search->wfp->lookup->mod_lookup_table_memory = NULL;
1448
2300
 
1449
 
    SeqLocFree(slp);
1450
 
    
1451
2301
    return seqalign;
1452
2302
}
1453
2303
 
1549
2399
                
1550
2400
                dbtag = bsp->id->data.ptrvalue;
1551
2401
                dbname = StringSave(dbtag->db);
1552
 
                cdhThis->CDDid = StringSave(dbtag->tag->str);
 
2402
                if (dbtag->tag->str) 
 
2403
                  cdhThis->CDDid = StringSave(dbtag->tag->str);
1553
2404
        }
1554
2405
        if (StrCmp(dbname,"Pfam") == 0) {
1555
2406
          cdhThis->ShortName = SaveUntilChar(title, &tmp, ',');
1568
2419
          cdhThis->ShortName = SaveUntilChar(title, &tmp, ',');
1569
2420
          if (tmp)
1570
2421
                cdhThis->Definition = SaveUntilChar(tmp+1, &tmp, ',');
 
2422
        } else if (StrCmp(dbname,"CDD") == 0) {
 
2423
          cdhThis->CDDid = SaveUntilChar(title, &tmp, ',');
 
2424
          if (tmp) cdhThis->ShortName = SaveUntilChar(tmp+1, &tmp, ',');
 
2425
          if (tmp) cdhThis->Definition = SaveUntilChar(tmp+1, &tmp, ';');
1571
2426
        } else {
1572
2427
          cdhThis->ShortName = StringSave(cdhThis->CDDid);
1573
2428
          cdhThis->Definition = StringSave(title);
1626
2481
    SeqEntryPtr sep;
1627
2482
    Char buffer[64];
1628
2483
    static TNlmMutex print_mutex;
1629
 
    SeqLocPtr query_lcase_mask = NULL;
 
2484
    SeqLocPtr query_lcase_mask = NULL, query_slp = NULL;
1630
2485
    RPSInfoPtr rpsinfo_local;
1631
 
    Int4  old_searchsp_eff;
 
2486
    Int4  old_searchsp_eff, qlen;
1632
2487
    static Int4 count_id = 0;
1633
2488
 
1634
2489
    if((mtdata = (RPS_MTDataStructPtr) data) == NULL) {
1660
2515
                        rpsbop->query_is_protein ? FindProt : FindNuc); 
1661
2516
        
1662
2517
        if (query_bsp == NULL) {
1663
 
            ErrPostEx(SEV_FATAL, 0, 0, "Unable to obtain bioseq\n");
 
2518
            ErrPostEx(SEV_FATAL, 1, 0, "Unable to obtain bioseq\n");
1664
2519
            return NULL;
1665
2520
        }  
1666
2521
 
1701
2556
        
1702
2557
        NlmMutexLock(print_mutex);
1703
2558
 
 
2559
        /* Create the query_slp */
 
2560
        if (rpsbop->query_is_protein && rpsinfo_local->query_slp == NULL) {
 
2561
            if (rpsinfo_local->start == 0 && rpsinfo_local->stop == 0) {
 
2562
                ValNodeAddPointer(&query_slp, SEQLOC_WHOLE,
 
2563
                        SeqIdDup(fake_bsp->id));
 
2564
            } else {
 
2565
 
 
2566
                rpsinfo_local->start = MAX(rpsinfo_local->start, 0);
 
2567
                if (rpsinfo_local->stop <= 0)
 
2568
                    rpsinfo_local->stop = fake_bsp->length - 1;
 
2569
                rpsinfo_local->stop = MIN(rpsinfo_local->stop, 
 
2570
                                          fake_bsp->length-1);
 
2571
 
 
2572
                if (rpsinfo_local->start >= fake_bsp->length) {
 
2573
                    ErrPostEx(SEV_ERROR, 0, 0, 
 
2574
                        "Incorrect range restriction (%ld-%ld) for query "
 
2575
                        "sequence (length %ld). Ignoring range restriction.", 
 
2576
                        rpsinfo_local->start, rpsinfo_local->stop, 
 
2577
                        fake_bsp->length);
 
2578
                    ValNodeAddPointer(&query_slp, SEQLOC_WHOLE,
 
2579
                            SeqIdDup(fake_bsp->id));
 
2580
                } else {
 
2581
                    query_slp = SeqLocIntNew(rpsinfo_local->start,
 
2582
                            rpsinfo_local->stop, Seq_strand_unknown,
 
2583
                            SeqIdDup(fake_bsp->id));
 
2584
                }
 
2585
            }
 
2586
            rpsinfo_local->query_slp = query_slp;
 
2587
        }
 
2588
 
 
2589
 
 
2590
        if (rpsbop->query_is_protein)
 
2591
            qlen = SeqLocLen(query_slp);
 
2592
        else
 
2593
            qlen = fake_bsp->length;
1704
2594
        old_searchsp_eff = rpsbop->options->searchsp_eff;
1705
 
        RPSUpdateDbSize(rpsbop->options, rpsinfo_local, fake_bsp->length);
 
2595
        RPSUpdateDbSize(rpsbop->options, rpsinfo_local, qlen);
1706
2596
        
1707
 
        search = BLASTSetUpSearch (bsp, rpsbop->query_is_protein ? 
1708
 
                                   "blastp" : "tblastn", 
1709
 
                                   bsp->length, 0, 
1710
 
                                   NULL, rpsbop->options, NULL);
 
2597
 
 
2598
        if (rpsbop->query_is_protein) {
 
2599
            search = BLASTSetUpSearchByLoc(query_slp, "blastp",
 
2600
                    SeqLocLen(query_slp), 0, NULL, rpsbop->options, NULL);
 
2601
        } else {
 
2602
            search = BLASTSetUpSearch(bsp, "tblastn", bsp->length, 0, NULL,
 
2603
                    rpsbop->options, NULL);
 
2604
        }
1711
2605
        rpsbop->options->searchsp_eff = old_searchsp_eff;
 
2606
        search->rps_qlen = qlen; /* save for reporting purposes */
 
2607
 
1712
2608
        NlmMutexUnlock(print_mutex); 
1713
2609
 
1714
2610
        /* External lower-case mask */
1717
2613
        
1718
2614
        NlmMutexLock(print_mutex);
1719
2615
        
1720
 
        other_returns = BlastOtherReturnsPrepare(search);
 
2616
        other_returns = BlastOtherReturnsPrepare(search); /* ,0); */
1721
2617
 
1722
2618
        mtdata->print_callback(fake_bsp, rpsbop, seqalign, 
1723
2619
                               other_returns, NULL, 
1729
2625
 
1730
2626
        if(!rpsbop->query_is_protein) 
1731
2627
            BioseqFree(bsp);
1732
 
        
 
2628
 
1733
2629
        SeqAlignSetFree(seqalign);
1734
2630
        search = BlastSearchBlkDestruct(search);
 
2631
        rpsinfo_local->query_slp = query_slp = NULL;
1735
2632
        
1736
2633
        if(!rpsbop->believe_query)
1737
2634
            fake_bsp = BlastDeleteFakeBioseq(fake_bsp);
1805
2702
            for (i=0; i<rpsbop->num_threads; i++) {
1806
2703
                NlmThreadJoin(thread_array[i], &status);
1807
2704
            }
 
2705
            thread_array = MemFree(thread_array);
1808
2706
 
1809
2707
        } else {
1810
2708
            RPSEngineThread((VoidPtr) mtdata);
1822
2720
}
1823
2721
 
1824
2722
/* These functions may be never be used ... */
1825
 
Int4Ptr PNTR RPSReadPSMatrix(CharPtr filename, Int4Ptr mat_len)
 
2723
static Int4Ptr PNTR RPSReadPSMatrix(CharPtr filename, Int4Ptr mat_len)
1826
2724
{
1827
2725
    Int4Ptr PNTR psmatrix;
1828
2726
    FILE *fd;
1864
2762
 
1865
2763
    return psmatrix;
1866
2764
}
1867
 
BioseqPtr RPSGetBioseqFromMatrix(Int4Ptr PNTR psmatrix, Int4 length)
 
2765
static BioseqPtr RPSGetBioseqFromMatrix(Int4Ptr PNTR psmatrix, Int4 length)
1868
2766
{
1869
2767
    BioseqPtr bsp;
1870
2768
    CharPtr buffer;
1985
2883
  }
1986
2884
}
1987
2885
 
 
2886
static SeqFeatPtr CreateNewFeatureOnSeqAnnot (SeqAnnotPtr sap, Uint1 choice, SeqLocPtr slp)
 
2887
 
 
2888
{
 
2889
  SeqFeatPtr  prev;
 
2890
  SeqFeatPtr  sfp;
 
2891
 
 
2892
  if (sap == NULL || slp == NULL) return NULL;
 
2893
  sfp = SeqFeatNew ();
 
2894
  if (sfp == NULL) return NULL;
 
2895
  if (sap->data != NULL) {
 
2896
    prev = sap->data;
 
2897
    while (prev->next != NULL) {
 
2898
      prev = prev->next;
 
2899
    }
 
2900
    prev->next = sfp;
 
2901
  } else {
 
2902
    sap->data = (Pointer) sfp;
 
2903
  }
 
2904
  sfp->data.choice = choice;
 
2905
  sfp->location = SeqLocFree (sfp->location);
 
2906
  sfp->location = AsnIoMemCopy (slp, (AsnReadFunc) SeqLocAsnRead,
 
2907
                                (AsnWriteFunc) SeqLocAsnWrite);
 
2908
  return sfp;
 
2909
}
 
2910
 
1988
2911
NLM_EXTERN void AnnotateRegionsFromCDD (
1989
2912
  BioseqPtr bsp,
1990
2913
  SeqAlignPtr salp,
1993
2916
 
1994
2917
{
1995
2918
  DbtagPtr       db;
 
2919
  AnnotDescrPtr  desc;
 
2920
  DatePtr        dp;
1996
2921
  CddHitPtr      cdd_head, cdd;
 
2922
  SeqAnnotPtr    lastsap;
1997
2923
  size_t         len;
1998
2924
  ObjectIdPtr    oip;
 
2925
  SeqAnnotPtr    sap = NULL;
1999
2926
  SeqFeatPtr     sfp;
2000
2927
  SeqInt         sint;
2001
2928
  SeqIdPtr       sip;
2029
2956
      vn.extended = 0;
2030
2957
      vn.data.ptrvalue = &sint;
2031
2958
      vn.next = NULL;
2032
 
      sfp = CreateNewFeatureOnBioseq (bsp, SEQFEAT_REGION, &vn);
 
2959
      if (sap == NULL) {
 
2960
        sap = SeqAnnotNew ();
 
2961
        if (sap != NULL) {
 
2962
          sap->type = 1;
 
2963
          dp = DateCurr ();
 
2964
          if (dp != NULL) {
 
2965
            desc = ValNodeNew (NULL);
 
2966
            if (desc != NULL) {
 
2967
              desc->choice = Annot_descr_create_date;
 
2968
              desc->data.ptrvalue = DateCurr ();
 
2969
              sap->desc = desc;
 
2970
            }
 
2971
          }
 
2972
          if (bsp->annot != NULL) {
 
2973
            lastsap = bsp->annot;
 
2974
            while (lastsap->next != NULL) {
 
2975
              lastsap = lastsap->next;
 
2976
            }
 
2977
            lastsap->next = sap;
 
2978
          } else {
 
2979
            bsp->annot = sap;
 
2980
          }
 
2981
        }
 
2982
      }
 
2983
      sfp = CreateNewFeatureOnSeqAnnot (sap, SEQFEAT_REGION, &vn);
2033
2984
      if (sfp != NULL) {
2034
2985
 
2035
2986
        if (! StringHasNoText (cdd->Definition)) {
2048
2999
          AddFieldToCddUserObject (uop, "evalue", NULL, 0, cdd->evalue);
2049
3000
          AddFieldToCddUserObject (uop, "bit_score", NULL, 0, cdd->bit_score);
2050
3001
        }
2051
 
        if (cdd->ShortName != NULL) {
 
3002
        if (cdd->ShortName != NULL && StringICmp (cdd->ShortName, cdd->Definition) != 0) {
2052
3003
          len = StringLen (cdd->ShortName) + 10;
2053
3004
          str = MemNew (len);
2054
3005
          if (str != NULL) {
2227
3178
  DeleteMarkedObjects (0, OBJ_SEQENTRY, topsep);
2228
3179
}
2229
3180
 
2230
 
static void FreeCDDProc (
 
3181
static void FreeCDDFeatProc (
2231
3182
  SeqFeatPtr sfp,
2232
3183
  Pointer userdata
2233
3184
)
2248
3199
)
2249
3200
 
2250
3201
{
2251
 
  VisitFeaturesInSep (topsep, NULL, FreeCDDProc);
2252
 
  DeleteMarkedObjects (0, OBJ_SEQENTRY, topsep);
2253
 
}
 
3202
  VisitFeaturesInSep (topsep, NULL, FreeCDDFeatProc);
 
3203
  DeleteMarkedObjects (0, OBJ_SEQENTRY, topsep);
 
3204
}
 
3205
 
 
3206
static void IsCDDSeqId (SeqIdPtr sip, Pointer userdata)
 
3207
 
 
3208
{
 
3209
  BoolPtr   cddIDp;
 
3210
  DbtagPtr  dbt;
 
3211
 
 
3212
  cddIDp = (BoolPtr) userdata;
 
3213
  if (sip == NULL || cddIDp == NULL) return;
 
3214
  if (sip->choice != SEQID_GENERAL) return;
 
3215
  dbt = (DbtagPtr) sip->data.ptrvalue;
 
3216
  if (dbt == NULL) return;
 
3217
  if (StringCmp (dbt->db, "CDD") != 0) return;
 
3218
  *cddIDp = TRUE;
 
3219
}
 
3220
 
 
3221
static void FreeCDDAlignProc (
 
3222
  SeqAlignPtr sap,
 
3223
  Pointer userdata
 
3224
)
 
3225
 
 
3226
{
 
3227
  Boolean  cddID = FALSE;
 
3228
 
 
3229
  VisitSeqIdsInSeqAlign (sap, (Pointer) &cddID, IsCDDSeqId);
 
3230
  if (cddID) {
 
3231
    sap->idx.deleteme = TRUE;
 
3232
  }
 
3233
}
 
3234
 
 
3235
NLM_EXTERN void FreeCDDAligns (
 
3236
  SeqEntryPtr topsep
 
3237
)
 
3238
 
 
3239
{
 
3240
  VisitAlignmentsInSep (topsep, NULL, FreeCDDAlignProc);
 
3241
  DeleteMarkedObjects (0, OBJ_SEQENTRY, topsep);
 
3242
}
 
3243
 
 
3244
 
2254
3245