← Back to branch summary

~ubuntu-branches/ubuntu/jaunty/ncbi-tools6/jaunty

« back to all changes in this revision

Viewing changes to demo/nps2gps.c

  • Committer: Bazaar Package Importer
  • Author(s): Aaron M. Ucko
  • Date: 2008-07-14 19:43:15 UTC
  • mfrom: (2.1.12 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080714194315-ed44u9ek7txva2rz
Tags: 6.1.20080302-3
http://bugs.debian.org/490437
tools/readdb.c: enable madvise()-based code on all glibc (hence all
Debian) systems, not just Linux.  (Closes: #490437.)

Show diffs side-by-side

added added

removed removed

Lines of Context:
demo/nps2gps.c
29
29
*
30
30
* Version Creation Date:   5/12/05
31
31
*
32
 
* $Revision: 1.8 $
 
32
* $Revision: 1.13 $
33
33
*
34
34
* File Description:
35
35
*
50
50
#include <toasn3.h>
51
51
#include <pmfapi.h>
52
52
 
53
 
#define NPS2GPSAPP_VER "1.6"
 
53
#define NPS2GPSAPP_VER "2.1"
54
54
 
55
55
CharPtr NPS2GPSAPPLICATION = NPS2GPSAPP_VER;
56
56
 
59
59
  CharPtr  outfile;
60
60
  Boolean  failure;
61
61
  Boolean  lock;
 
62
  Boolean  byFeatID;
 
63
  Boolean  useProtID;
62
64
} N2GData, PNTR N2GPtr;
63
65
 
64
66
typedef struct npsseqs {
152
154
  }
153
155
  SetSeqFeatProduct (sfp, ns.prot);
154
156
 
155
 
  /* paste fields from temp copy of original cds */
 
157
  /* paste fields from temp copy of original cds - but not feature ID */
156
158
  crp = (CdRegionPtr) temp->data.value.ptrvalue;
157
159
  sfp->data.value.ptrvalue = (Pointer) crp;
158
160
 
178
180
    stop = GetOffsetInLoc (cbp->loc, mrna->location, SEQLOC_STOP);
179
181
    if (start < 0 || start >= ns.nuc->length ||
180
182
        stop < 0 || stop >= ns.nuc->length) continue;
181
 
        cbp->loc = SeqLocFree (cbp->loc);
182
 
        cbp->loc = AddIntervalToLocation (NULL, sip, start, stop, partial5, partial3);;
 
183
    cbp->loc = SeqLocFree (cbp->loc);
 
184
    cbp->loc = AddIntervalToLocation (NULL, sip, start, stop, partial5, partial3);;
183
185
  }
 
186
 
 
187
  /* remove feature ID cross-references of original cds */
 
188
  ClearFeatIDXrefs (sfp);
184
189
}
185
190
 
186
191
/* copy gene from contig onto nuc-prot, single interval on cdna bioseq */
338
343
  SeqDescrAddPointer (&(bsp->descr), Seq_descr_title, (Pointer) str);
339
344
}
340
345
 
341
 
static SeqIdPtr MakeIdFromLocusTag (SeqFeatPtr mrna)
 
346
static SeqIdPtr MakeIdFromLocusTag (
 
347
  SeqFeatPtr mrna
 
348
)
342
349
 
343
350
{
344
351
  Char        buf [64], suffix [8];
345
352
  SeqFeatPtr  gene;
346
353
  GeneRefPtr  grp;
 
354
  Int4Ptr     iptr;
347
355
 
348
356
  if (mrna == NULL) return NULL;
349
357
  suffix [0] = '\0';
355
363
    gene = SeqMgrGetOverlappingGene (mrna->location, NULL);
356
364
    if (gene != NULL) {
357
365
      grp = (GeneRefPtr) gene->data.value.ptrvalue;
358
 
      if (gene->id.choice == 1) {
359
 
        (gene->id.value.intvalue)++;
360
 
        sprintf (suffix, "_%ld", (long) gene->id.value.intvalue);
 
366
      iptr = (Int4Ptr) gene->idx.scratch;
 
367
      if (iptr != NULL) {
 
368
        (*iptr)++;
 
369
        sprintf (suffix, "_%ld", (long) (*iptr));
361
370
      }
362
371
    }
363
372
  }
373
382
  return NULL;
374
383
}
375
384
 
376
 
static void InstantiateMrnaIntoProt (SeqFeatPtr cds, SeqFeatPtr mrna, Int2Ptr ctrp, Boolean ambig)
 
385
static SeqIdPtr MakeIdFromProtein (
 
386
  BioseqPtr pbsp
 
387
)
 
388
 
 
389
{
 
390
  Char      buf [64], tmp [64];
 
391
  SeqIdPtr  sip;
 
392
 
 
393
  if (pbsp == NULL) return NULL;
 
394
  sip = SeqIdFindBestAccession (pbsp->id);
 
395
  if (sip == NULL) return NULL;
 
396
  SeqIdWrite (sip, buf, PRINTID_TEXTID_ACC_VER, sizeof (buf) - 1);
 
397
  if (StringHasNoText (buf)) return NULL;
 
398
  StringCpy (tmp, "gnl|MTRACK|");
 
399
  StringCat (tmp, buf);
 
400
  StringCat (tmp, "_mrna");
 
401
  return MakeSeqID (tmp);
 
402
}
 
403
 
 
404
static void InstantiateMrnaIntoProt (
 
405
  SeqFeatPtr cds,
 
406
  SeqFeatPtr mrna,
 
407
  Int2Ptr ctrp,
 
408
  N2GPtr ngp
 
409
)
377
410
 
378
411
{
379
412
  ByteStorePtr  bs;
385
418
  CharPtr       rnaseq;
386
419
  ValNodePtr    vnp;
387
420
 
388
 
  if (cds == NULL || mrna == NULL) return;
 
421
  if (cds == NULL || mrna == NULL || ngp == NULL) return;
389
422
  if (cds->product == NULL || mrna->product != NULL) return;
390
423
 
391
424
  pbsp = BioseqFindFromSeqLoc (cds->product);
392
 
  if (pbsp  == NULL) return;
 
425
  if (pbsp == NULL) return;
393
426
  psep = SeqMgrGetSeqEntryForData (pbsp);
394
427
  if (psep == NULL) return;
395
428
 
407
440
  mbsp->repr = Seq_repr_raw;
408
441
  mbsp->mol = Seq_mol_rna;
409
442
  mbsp->seq_data_type = Seq_code_iupacna;
410
 
  mbsp->seq_data = bs;
 
443
  mbsp->seq_data = (SeqDataPtr) bs;
411
444
  mbsp->length = BSLen (bs);
412
445
  BioseqPack (mbsp);
413
446
 
 
447
  /* now adds _# suffix to general Seq-id if ambiguous - but messed up Drosophila record, so only use feature ID */
414
448
  /*
415
 
  if (! ambig) {
416
 
    mbsp->id = MakeIdFromLocusTag (mrna);
417
 
  }
418
 
  */
419
 
  /* now adds _# suffix to general Seq-id if ambiguous */
420
449
  mbsp->id = MakeIdFromLocusTag (mrna);
 
450
  */
 
451
  if (ngp->useProtID) {
 
452
    mbsp->id = MakeIdFromProtein (pbsp);
 
453
  }
421
454
  if (mbsp->id == NULL) {
422
455
    mbsp->id = MakeNewProteinSeqIdEx (mrna->location, NULL, NULL, ctrp);
423
456
  }
452
485
  AddSeqEntryToSeqEntry (psep, msep, FALSE);
453
486
}
454
487
 
455
 
static void RemoveFeatIDs (SeqFeatPtr sfp, Pointer userdata)
 
488
static void RemoveFeatScratch (
 
489
  SeqFeatPtr sfp,
 
490
  Pointer userdata
 
491
)
456
492
 
457
493
{
458
494
  if (sfp == NULL) return;
459
 
  if (sfp->id.choice == 1) {
460
 
    sfp->id.choice = 0;
461
 
    sfp->id.value.intvalue = 0;
462
 
  }
 
495
  if (sfp->idx.scratch == NULL) return;
 
496
  sfp->idx.scratch = MemFree (sfp->idx.scratch);
463
497
}
464
498
 
465
499
typedef struct loopdata {
466
 
  Boolean     ambig;
467
500
  Int2        count;
468
501
  SeqFeatPtr  mrna;
469
502
} LoopData, PNTR LoopDataPtr;
478
511
 
479
512
  ldp = (LoopDataPtr) context->userdata;
480
513
 
481
 
  if (sfp->id.choice == 0) {
 
514
  if (sfp->idx.scratch == NULL) {
482
515
    (ldp->count)++;
483
516
    ldp->mrna = sfp;
484
 
  } else {
485
 
    ldp->ambig = TRUE;
486
517
  }
487
518
 
488
519
  return TRUE;
489
520
}
490
521
 
491
 
static Boolean CdsIsPseudo (SeqFeatPtr cds)
 
522
static Boolean CdsIsPseudo (
 
523
  SeqFeatPtr cds
 
524
)
492
525
 
493
526
{
494
527
  SeqMgrFeatContext  gcontext;
510
543
  return grp->pseudo;
511
544
}
512
545
 
513
 
static void LoopThroughCDSs (BioseqPtr bsp, Int2Ptr ctrp, N2GPtr ngp)
 
546
static void LoopThroughCDSs (
 
547
  BioseqPtr bsp,
 
548
  Int2Ptr ctrp,
 
549
  N2GPtr ngp
 
550
)
514
551
 
515
552
{
516
553
  Int2               count;
527
564
    goOn = FALSE;
528
565
    sfp = SeqMgrGetNextFeature (bsp, NULL, SEQFEAT_CDREGION, 0, &fcontext);
529
566
    while (sfp != NULL) {
530
 
      if (sfp->id.choice == 0) {
531
 
        ld.ambig = FALSE;
 
567
      if (sfp->idx.scratch == NULL) {
532
568
        ld.count = 0;
533
569
        ld.mrna = NULL;
534
570
        count = SeqMgrGetAllOverlappingFeatures (sfp->location, FEATDEF_mRNA, NULL, 0,
535
571
                                                 CHECK_INTERVALS, (Pointer) &ld, FindSingleMrnaProc);
536
572
        if (ld.count == 1 && ld.mrna != NULL) {
537
 
          InstantiateMrnaIntoProt (sfp, ld.mrna, ctrp, ld.ambig);
538
 
          sfp->id.choice = 1;
539
 
          ld.mrna->id.choice = 1;
 
573
          InstantiateMrnaIntoProt (sfp, ld.mrna, ctrp, ngp);
 
574
          sfp->idx.scratch = (Pointer) MemNew (sizeof (Int4));
 
575
          ld.mrna->idx.scratch = (Pointer) MemNew (sizeof (Int4));
540
576
          goOn = TRUE;
541
577
        } else {
542
578
          str = fcontext.label;
546
582
          Message (MSG_POSTERR, " CDS '%s' has %d overlapping mRNA", str, (int) ld.count);
547
583
        }
548
584
      }
 
585
 
549
586
      sfp = SeqMgrGetNextFeature (bsp, sfp, SEQFEAT_CDREGION, 0, &fcontext);
550
587
    }
551
588
  }
554
591
 
555
592
  sfp = SeqMgrGetNextFeature (bsp, NULL, SEQFEAT_CDREGION, 0, &fcontext);
556
593
  while (sfp != NULL) {
557
 
    if (sfp->id.choice == 0) {
 
594
    if (sfp->idx.scratch == NULL) {
558
595
      if (sfp->product == NULL && CdsIsPseudo (sfp)) {
559
596
        /* do not complain about pseudo CDS not having product or matching mRNA */
560
597
      } else {
572
609
  }
573
610
}
574
611
 
575
 
static void NPStoGPS (Uint2 entityID, CharPtr filename, N2GPtr ngp, SeqDescrPtr descr)
 
612
static SeqFeatPtr GetmRNAByFeatureID (
 
613
  SeqFeatPtr cds
 
614
)
 
615
 
 
616
{
 
617
  SeqFeatPtr         mRNA = NULL;
 
618
  SeqFeatXrefPtr     xref;
 
619
  ObjectIdPtr        oip;
 
620
  Char               buf [32];
 
621
 
 
622
  if (cds == NULL) return NULL;
 
623
  
 
624
  xref = cds->xref;
 
625
  while (xref != NULL && mRNA == NULL) {
 
626
    if (xref->id.choice == 3) {
 
627
      oip = (ObjectIdPtr) xref->id.value.ptrvalue;
 
628
      if (oip != NULL) {
 
629
        if (StringDoesHaveText (oip->str)) {
 
630
          mRNA = SeqMgrGetFeatureByFeatID (cds->idx.entityID, NULL, oip->str, NULL, NULL);
 
631
        } else {
 
632
          sprintf (buf, "%ld", (long) oip->id);
 
633
          mRNA = SeqMgrGetFeatureByFeatID (cds->idx.entityID, NULL, buf, NULL, NULL);
 
634
        }
 
635
      }
 
636
    }
 
637
    xref = xref->next;
 
638
  }
 
639
  return mRNA;
 
640
}
 
641
 
 
642
static void NPStoGPS (
 
643
  Uint2 entityID,
 
644
  CharPtr filename,
 
645
  N2GPtr ngp,
 
646
  SeqDescrPtr descr
 
647
)
576
648
 
577
649
{
578
650
  BioSourcePtr       biop;
579
651
  BioseqPtr          bsp;
580
652
  BioseqSetPtr       bssp;
581
653
  Int2               ctr = 1;
582
 
  SeqMgrFeatContext  gcontext, mcontext;
 
654
  SeqMgrFeatContext  fcontext, gcontext, mcontext;
583
655
  SeqFeatPtr         gene, mrna, sfp, lastsfp;
584
656
  GeneRefPtr         grp;
 
657
  Int4Ptr            iptr;
585
658
  SeqEntryPtr        old, top;
586
659
  CharPtr            organism;
587
660
  OrgRefPtr          orp;
619
692
  annot = bsp->annot;
620
693
  if (annot == NULL) {
621
694
    bsp->annot = sap;
 
695
  } else if (sap == NULL) {
622
696
  } else if (sap->next == NULL && annot->next == NULL &&
623
697
             sap->type == 1 && annot->type == 1 &&
624
698
             sap->data != NULL && annot->data != NULL) {
638
712
    lastsap->next = sap;
639
713
  }
640
714
 
641
 
  VisitFeaturesInSep (top, NULL, RemoveFeatIDs);
642
 
 
643
715
  old = SeqEntrySetScope (top);
644
716
 
645
717
  ctr = (Int2) VisitBioseqsInSep (top, NULL, NULL) + 1;
652
724
    if (grp == NULL) {
653
725
      gene = SeqMgrGetOverlappingGene (mrna->location, NULL);
654
726
      if (gene != NULL) {
655
 
        if (gene->id.choice == 0) {
656
 
          gene->id.choice = 1;
657
 
          gene->id.value.intvalue = 1;
658
 
        } else if (gene->id.choice == 1) {
659
 
          (gene->id.value.intvalue)++;
 
727
        iptr = (Int4Ptr) gene->idx.scratch;
 
728
        if (iptr == NULL) {
 
729
          iptr = (Int4Ptr) MemNew (sizeof (Int4));
 
730
          gene->idx.scratch = (Pointer) iptr;
 
731
        }
 
732
        if (iptr != NULL) {
 
733
          (*iptr)++;
660
734
        }
661
735
      }
662
736
    }
667
741
 
668
742
  gene = SeqMgrGetNextFeature (bsp, NULL, SEQFEAT_GENE, 0, &gcontext);
669
743
  while (gene != NULL) {
670
 
    if (gene->id.choice == 1) {
671
 
      if (gene->id.value.intvalue == 1) {
672
 
        /* if count was 1, clear id */
673
 
        gene->id.value.intvalue = 0;
674
 
        gene->id.choice = 0;
 
744
    iptr = (Int4Ptr) gene->idx.scratch;
 
745
    if (iptr != NULL) {
 
746
      if (*iptr == 1) {
 
747
        /* if count was 1, clear scratch */
 
748
        gene->idx.scratch = MemFree (gene->idx.scratch);
675
749
      } else {
676
 
        /* if count was > 1, just reset count, do not clear id */
677
 
        gene->id.value.intvalue = 0;
 
750
        /* if count was > 1, just reset count, do not clear scratch */
 
751
        *iptr = 0;
678
752
      }
679
753
    }
680
754
    gene = SeqMgrGetNextFeature (bsp, gene, SEQFEAT_GENE, 0, &gcontext);
681
755
  }
682
756
 
683
 
  /*
684
 
  sfp = SeqMgrGetNextFeature (bsp, NULL, SEQFEAT_CDREGION, 0, &fcontext);
685
 
  while (sfp != NULL) {
686
 
    mrna = SeqMgrGetOverlappingFeature (sfp->location, FEATDEF_mRNA, NULL, 0,
687
 
                                        NULL, CHECK_INTERVALS, &mcontext);
688
 
    if (mrna != NULL) {
689
 
      InstantiateMrnaIntoProt (sfp, mrna, &ctr);
 
757
  if (ngp->byFeatID) {
 
758
 
 
759
    /* trust feature ID cross-references */
 
760
 
 
761
    sfp = SeqMgrGetNextFeature (bsp, NULL, SEQFEAT_CDREGION, 0, &fcontext);
 
762
    while (sfp != NULL) {
 
763
      mrna = GetmRNAByFeatureID (sfp);
 
764
      if (mrna != NULL) {
 
765
        InstantiateMrnaIntoProt (sfp, mrna, &ctr, ngp);
 
766
      }
 
767
      sfp = SeqMgrGetNextFeature (bsp, sfp, SEQFEAT_CDREGION, 0, &fcontext);
690
768
    }
691
 
    sfp = SeqMgrGetNextFeature (bsp, sfp, SEQFEAT_CDREGION, 0, &fcontext);
 
769
 
 
770
  } else {
 
771
 
 
772
    /* make cDNA bioseq from mRNA feature, package with protein product */
 
773
 
 
774
    LoopThroughCDSs (bsp, &ctr, ngp);
692
775
  }
693
 
  */
694
 
 
695
 
  /* make cDNA bioseq from mRNA feature, package with protein product */
696
 
 
697
 
  LoopThroughCDSs (bsp, &ctr, ngp);
 
776
 
 
777
  VisitFeaturesInSep (top, NULL, RemoveFeatScratch);
698
778
 
699
779
  SeqMgrLinkSeqEntry (top, parenttype, parentptr);
700
780
 
702
782
 
703
783
  SeqMgrClearFeatureIndexes (bssp->idx.entityID, NULL);
704
784
 
705
 
  VisitFeaturesInSep (top, NULL, RemoveFeatIDs);
706
 
 
707
785
  /* need to reindex to get mRNA and CDS features from cDNA and protein */
708
786
  SeqMgrIndexFeatures (bssp->idx.entityID, NULL);
709
787
  VisitSetsInSet (bssp, NULL, LclMakeNucProtCDS);
906
984
#define x_argSuffix      5
907
985
#define R_argRemote      6
908
986
#define L_argLockFar     7
 
987
#define F_argUseFeatID   8
 
988
#define P_argUseProtID   9
909
989
 
910
990
 
911
991
Args myargs [] = {
925
1005
    TRUE, 'R', ARG_BOOLEAN, 0.0, 0, NULL},
926
1006
  {"Lock Components in Advance", "F", NULL, NULL,
927
1007
    TRUE, 'L', ARG_BOOLEAN, 0.0, 0, NULL},
 
1008
  {"Map by Feature ID", "F", NULL, NULL,
 
1009
    TRUE, 'F', ARG_BOOLEAN, 0.0, 0, NULL},
 
1010
  {"mRNA ID from Protein", "F", NULL, NULL,
 
1011
    TRUE, 'P', ARG_BOOLEAN, 0.0, 0, NULL},
928
1012
};
929
1013
 
930
1014
Int2 Main (void)
976
1060
  ngp = &ngd;
977
1061
  ngd.failure = FALSE;
978
1062
  ngd.lock = (Boolean) myargs [L_argLockFar].intvalue;
 
1063
  ngd.byFeatID = (Boolean) myargs [F_argUseFeatID].intvalue;
 
1064
  ngd.useProtID = (Boolean) myargs [P_argUseProtID].intvalue;
979
1065
 
980
1066
  directory = (CharPtr) myargs [p_argInputPath].strvalue;
981
1067
  results = (CharPtr) myargs [r_argOutputPath].strvalue;