source: branches/alilink/READSEQ/ureadasn.c

/* ureadasn.c
  -- parse, mangle and otherwise rewrite ASN1 file/entries for readseq reading
  -- from NCBI toolkit (ncbi.nlm.nih.gov:/toolkit)
*/

#ifdef NCBI

#include <stdio.h>
#include <ctype.h>
#include <string.h>

/* NCBI toolkit :include: must be on lib path */
#include <ncbi.h>
#include <seqport.h>

#define UREADASN
#include "ureadseq.h"

#pragma segment ureadasn

/* this stuff is hacked up from tofasta.c of ncbitools */
#define   kBaseAny   0
#define   kBaseNucleic 1
#define   kBaseAmino   2

typedef struct tofasta {
    Boolean idonly;
    short   *seqnum;
    short   whichSeq;
    char    **seq, **seqid;
    long    *seqlen;
} FastaDat, PNTR FastaPtr;


void BioseqRawToRaw(BioseqPtr bsp, Boolean idonly,
              short whichSeq, short *seqnum,
              char **seq, char **seqid, long *seqlen)
{
  SeqPortPtr spp;
  SeqIdPtr bestid;
  Uint1 repr, code, residue;
  CharPtr tmp, title;
  long  outlen, outmax;
  char  localid[256], *sp;

  /* !!! this may be called several times for a single sequence
    because SeqEntryExplore looks for parts and joins them...
    assume seq, seqid, seqlen may contain data (or NULL)
  */
  if (bsp == NULL) return;
  repr = Bioseq_repr(bsp);
  if (!(repr == Seq_repr_raw || repr == Seq_repr_const)) return;

  (*seqnum)++;
  if (!(whichSeq == *seqnum || whichSeq == 0)) return;

  bestid = SeqIdFindBest(bsp->id, (Uint1) 0);
  title = BioseqGetTitle(bsp);
  if (idonly) {
    sprintf(localid, " %d)  ", *seqnum);
    tmp= localid + strlen(localid)-1;
    }
  else {
    strcpy(localid," ");
    tmp= localid;
    }
  tmp = SeqIdPrint(bestid, tmp, PRINTID_FASTA_SHORT);
  tmp = StringMove(tmp, " ");
  StringNCpy(tmp, title, 200);
/* fprintf(stderr,"BioseqRawToRaw: localid='%s'\n",localid); */

          /* < seqid is fixed storage */
  /* strcpy( *seqid, localid);  */
          /* < seqid is variable sized */
  outmax= strlen(localid) + 3;
  if (*seqid==NULL) {
    *seqid= (char*) malloc(outmax);
    if (*seqid==NULL) return;
    strcpy(*seqid, localid);
    }
  else {
    outmax += strlen(*seqid) + 2;
    *seqid= (char*) realloc( *seqid, outmax);
    if (*seqid==NULL) return;
    if (!idonly) strcat(*seqid, "; ");
    strcat(*seqid, localid);
    }

  if (idonly) {
    strcat(*seqid,"\n");
    return;
    }

  if (ISA_na(bsp->mol)) code = Seq_code_iupacna;
  else code = Seq_code_iupacaa;
  spp = SeqPortNew(bsp, 0, -1, 0, code);
  SeqPortSeek(spp, 0, SEEK_SET);

  sp= *seq;
  if (sp==NULL) {  /* this is always true now !? */
    outlen= 0;
    outmax= 500;
    sp= (char*) malloc(outmax);
    }
  else {
    outlen= strlen(sp);
    outmax= outlen + 500;
    sp= (char*) realloc( sp, outmax);
    }
  if (sp==NULL) return;

  while ((residue = SeqPortGetResidue(spp)) != SEQPORT_EOF) {
    if (outlen>=outmax) {
      outmax= outlen + 500;
      sp= (char*) realloc(sp, outmax);
      if (sp==NULL) return;
      }
    sp[outlen++] = residue;
    }
  sp= (char*) realloc(sp, outlen+1);
  if (sp!=NULL) sp[outlen]= '\0';
  *seq= sp;
  *seqlen= outlen;
  SeqPortFree(spp);
  return;
}


static void SeqEntryRawseq(SeqEntryPtr sep, Pointer data, Int4 index, Int2 indent)
{
  FastaPtr tfa;
  BioseqPtr bsp;

  if (!IS_Bioseq(sep)) return;
  bsp = (BioseqPtr)sep->data.ptrvalue;
  tfa = (FastaPtr) data;
  BioseqRawToRaw(bsp, tfa->idonly, tfa->whichSeq, tfa->seqnum,
                  tfa->seq, tfa->seqid, tfa->seqlen);
}

void SeqEntryToRaw(SeqEntryPtr sep, Boolean idonly, short whichSeq, short *seqnum,
                        char **seq, char **seqid, long *seqlen)
{
  FastaDat tfa;

  if (sep == NULL) return;
  tfa.idonly= idonly;
  tfa.seqnum= seqnum;
  tfa.whichSeq= whichSeq;
  tfa.seq   = seq;
  tfa.seqid = seqid;
  tfa.seqlen= seqlen;
  SeqEntryExplore(sep, (Pointer)&tfa, SeqEntryRawseq);
}




char *listASNSeqs(const char *filename, const long skiplines,
                  const short format,   /* note: this is kASNseqentry or kASNseqset */
                  short *nseq, short *error )
{
  AsnIoPtr aip = NULL;
  SeqEntryPtr the_set;
  AsnTypePtr atp, atp2;
  AsnModulePtr amp;
  Boolean inIsBinary= FALSE; /* damn, why can't asn routines test this? */
  char  *seq = NULL;
  char  *seqid = NULL, stemp[256];
  long  seqlen;
  int   i, count;

  *nseq= 0;
  *error= 0;

    /* asn dictionary setups */
/*fprintf(stderr,"listASNSeqs: SeqEntryLoad\n");*/
  if (! SeqEntryLoad()) goto errxit; /*  sequence alphabets (and sequence parse trees) */
  amp = AsnAllModPtr();   /* get pointer to all loaded ASN.1 modules */
  if (amp == NULL) goto errxit;
  atp = AsnFind("Bioseq-set");    /* get the initial type pointers */
  if (atp == NULL) goto errxit;
  atp2 = AsnFind("Bioseq-set.seq-set.E");
  if (atp2 == NULL) goto errxit;

/*fprintf(stderr,"listASNSeqs: AsnIoOpen\n");*/
      /* open the ASN.1 input file in the right mode */
      /* !!!! THIS FAILS when filename has MAC PATH (& other paths?) (:folder:filename) */
  if ((aip = AsnIoOpen(filename, inIsBinary?"rb":"r")) == NULL) goto errxit;
  for (i=0; i<skiplines; i++) fgets( stemp, 255, aip->fp);  /* this may mess up asn routines... */

  if (! ErrSetLog ("stderr"))  goto errxit;
  else ErrSetOpts(ERR_CONTINUE, ERR_LOG_ON);    /*??  log errors instead of die */

  if (format == kASNseqentry) {  /* read one Seq-entry */
/*fprintf(stderr,"listASNSeqs: SeqEntryAsnRead\n");*/
    the_set = SeqEntryAsnRead(aip, NULL);
    SeqEntryToRaw(the_set, true,  0, nseq, &seq, &seqid, &seqlen);
    if (seq) free(seq); seq= NULL;
    SeqEntryFree(the_set);
    }
  else   {                   /* read Seq-entry's from a Bioseq-set */
    count = 0;
/*fprintf(stderr,"listASNSeqs: AsnReadId\n");*/
    while ((atp = AsnReadId(aip, amp, atp)) != NULL) {
      if (atp == atp2)  {  /* top level Seq-entry */
        the_set = SeqEntryAsnRead(aip, atp);
        SeqEntryToRaw(the_set, true, 0, nseq, &seq, &seqid, &seqlen);
        SeqEntryFree(the_set);
        if (seq) free(seq); seq= NULL;
        }
      else
        AsnReadVal(aip, atp, NULL);
      count++;
      }
    }

  AsnIoClose(aip);
  *error= 0;
  return seqid;

errxit:
  AsnIoClose(aip);
  if (seqid) free(seqid);
  *error= eASNerr;
  return NULL;
}


char *readASNSeq(const short whichEntry, const char *filename,
                const long skiplines,
                const short format,     /* note: this is kASNseqentry or kASNseqset */
                long *seqlen, short *nseq,
                short *error, char **seqid )
{
  AsnIoPtr aip = NULL;
  SeqEntryPtr the_set;
  AsnTypePtr atp, atp2;
  AsnModulePtr amp;
  Boolean inIsBinary= FALSE; /* damn, why can't asn routines test this? */
  char  *seq, stemp[200];
  int   i, count;

  *seqlen= 0;
  *nseq= 0;
  *error= 0;
  seq= NULL;

/*fprintf(stderr,"readASNseq: SeqEntryLoad\n");*/
    /* asn dictionary setups */
  if (! SeqEntryLoad()) goto errxit; /*  sequence alphabets (and sequence parse trees) */
  amp = AsnAllModPtr();   /* get pointer to all loaded ASN.1 modules */
  if (amp == NULL) goto errxit;
  atp = AsnFind("Bioseq-set");    /* get the initial type pointers */
  if (atp == NULL) goto errxit;
  atp2 = AsnFind("Bioseq-set.seq-set.E");
  if (atp2 == NULL) goto errxit;

      /* open the ASN.1 input file in the right mode */
/*fprintf(stderr,"readASNseq: AsnIoOpen(%s)\n", filename);*/
  if ((aip = AsnIoOpen(filename, inIsBinary?"rb":"r")) == NULL) goto errxit;
  for (i=0; i<skiplines; i++) fgets( stemp, 255, aip->fp);  /* this may mess up asn routines... */

  if (! ErrSetLog ("stderr"))  goto errxit;
  else ErrSetOpts(ERR_CONTINUE, ERR_LOG_ON);    /*??  log errors instead of die */

  seq= NULL;
  if (format == kASNseqentry) {  /* read one Seq-entry */
/*fprintf(stderr,"readASNseq: SeqEntryAsnRead\n");*/
    the_set = SeqEntryAsnRead(aip, NULL);
    SeqEntryToRaw(the_set, false, whichEntry, nseq, &seq, seqid, seqlen);
    SeqEntryFree(the_set);
    goto goodexit;
    }

  else   {                   /* read Seq-entry's from a Bioseq-set */
    count = 0;
/*fprintf(stderr,"readASNseq: AsnReadId\n");*/
    while ((atp = AsnReadId(aip, amp, atp)) != NULL) {
      if (atp == atp2)  {  /* top level Seq-entry */
        the_set = SeqEntryAsnRead(aip, atp);
        SeqEntryToRaw(the_set, false, whichEntry, nseq, &seq, seqid, seqlen);
        SeqEntryFree(the_set);
        if (*nseq >= whichEntry) goto goodexit;
        }
      else
        AsnReadVal(aip, atp, NULL);
      count++;
      }
    }

goodexit:
  AsnIoClose(aip);
  *error= 0;
  return seq;

errxit:
  AsnIoClose(aip);
  *error= eASNerr;
  if (seq) free(seq);
  return NULL;
}


#endif /*NCBI*/