forester/archive/RIO/others/hmmer/squid/afetch_main.c

   1 /*****************************************************************
   2  * HMMER - Biological sequence analysis with profile HMMs
   3  * Copyright (C) 1992-1999 Washington University School of Medicine
   4  * All Rights Reserved
   5  *
   6  *     This source code is distributed under the terms of the
   7  *     GNU General Public License. See the files COPYING and LICENSE
   8  *     for details.
   9  *****************************************************************/
  10
  11 /* afetch_main.c
  12  * SRE, Tue Nov  9 18:47:02 1999 [Saint Louis]
  13  *
  14  * afetch -- a program to extract alignments from the Pfam database
  15  *
  16  * CVS $Id: afetch_main.c,v 1.1.1.1 2005/03/22 08:34:30 cmzmasek Exp $
  17  */
  18
  19 #include <stdio.h>
  20 #include <string.h>
  21 #include "squid.h"
  22 #include "msa.h"
  23 #include "ssi.h"
  24
  25 static char banner[] = "afetch - retrieve an alignment from Pfam";
  26
  27 static char usage[] = "\
  28 Usage: afetch [-options] <alignment database> <name or accession>\n\
  29    or: afetch --index <alignment database>\n\
  30 \n\
  31    Get an alignment from a database.\n\
  32    Available options:\n\
  33    -h      : help; print version and usage info\n\
  34 ";
  35
  36 static char experts[] = "\
  37    --index : construct indices for the database\n\
  38 ";
  39
  40 struct opt_s OPTIONS[] = {
  41   { "-h",      TRUE,  sqdARG_NONE   },
  42   { "--index", FALSE, sqdARG_NONE   }
  43 };
  44 #define NOPTIONS (sizeof(OPTIONS) / sizeof(struct opt_s))
  45
  46 int
  47 main(int argc, char **argv)
  48 {
  49   char    *afile;               /* name of alignment file to read       */
  50   MSAFILE *afp;                 /* pointer to open index file           */
  51   char    *key;                 /* name/accession of alignment to fetch */
  52   MSA     *msa;                 /* the fetched alignment                */
  53   int      format;              /* format of afile */
  54   int      do_index;            /* TRUE to index instead of retrieve    */
  55
  56   char *optname;
  57   char *optarg;
  58   int   optind;
  59
  60   /***********************************************
  61    * Parse the command line
  62    ***********************************************/
  63
  64                                 /* initializations and defaults */
  65   format   = MSAFILE_STOCKHOLM; /* period. It's the only multi-MSA file format. */
  66   do_index = FALSE;
  67   key      = NULL;
  68
  69   while (Getopt(argc, argv, OPTIONS, NOPTIONS, usage,
  70                 &optind, &optname, &optarg))
  71     {
  72       if      (strcmp(optname, "--index") == 0) { do_index = TRUE; }
  73       else if (strcmp(optname, "-h")      == 0) {
  74         Banner(stdout, banner);
  75         puts(usage);
  76         puts(experts);
  77         exit(EXIT_SUCCESS);
  78       }
  79     }
  80
  81   if ((do_index && argc - optind != 1) || (! do_index && argc - optind != 2))
  82     Die("Incorrect number of command line arguments.\n%s\n", usage);
  83
  84   afile = argv[optind++];
  85   if (! do_index) key = argv[optind++];
  86
  87   if ((afp = MSAFileOpen(afile, format, NULL)) == NULL)
  88     Die("Alignment file %s could not be opened for reading", afile);
  89
  90   /***********************************************
  91    * Section 1. Alignment database indexing
  92    ***********************************************/
  93
  94   if (do_index) {
  95     int        mode;
  96     char      *ssifile;
  97     SSIINDEX  *si;
  98     int        fh;
  99     int        status;
 100     SSIOFFSET  offset;
 101     int        n = 0;
 102
 103     /* Not that we're expecting an alignment file so
 104      * large that it would require a 64-bit index, but...
 105      */
 106     if ((mode = SSIRecommendMode(afile)) == -1)
 107       Die("File %s doesn't exist, or is too large for your OS", afile);
 108
 109     ssifile = sre_strdup(afile, -1);
 110     sre_strcat(&ssifile, -1, ".ssi", -1);
 111
 112     if ((si = SSICreateIndex(mode)) == NULL)
 113       Die("Couldn't allocate/initialize the new SSI index");
 114     if (SSIAddFileToIndex(si, afile, afp->format, &fh) != 0)
 115       Die("SSIAddFileToIndex() failed");
 116
 117     status = SSIGetFilePosition(afp->f, mode, &offset);
 118     if (status != 0) Die("SSIGetFilePosition() failed");
 119
 120     while ((msa = MSAFileRead(afp)) != NULL)
 121       {
 122         if (msa->name == NULL)
 123           Die("SSI index requires that every MSA has a name");
 124
 125         status = SSIAddPrimaryKeyToIndex(si, msa->name, fh, &offset, NULL, 0);
 126         if (status != 0) Die("SSIAddPrimaryKeyToIndex() failed");
 127
 128         if (msa->acc != NULL) {
 129           status = SSIAddSecondaryKeyToIndex(si, msa->acc, msa->name);
 130           if (status != 0) Die("SSIAddSecondaryKeyToIndex() failed");
 131         }
 132
 133         status = SSIGetFilePosition(afp->f, mode, &offset);
 134         if (status != 0) Die("SSIGetFilePosition() failed");
 135
 136         n++;
 137         MSAFree(msa);
 138       }
 139
 140     status = SSIWriteIndex(ssifile, si);
 141     if (status != 0) Die("SSIWriteIndex() failed");
 142
 143     printf ("%d alignments indexed in SSI index %s\n", n, ssifile);
 144     free(ssifile);
 145     MSAFileClose(afp);
 146     SSIFreeIndex(si);
 147     SqdClean();
 148     exit (0);   /* exit indexing program here */
 149   }
 150
 151   /***********************************************
 152    * Section 2. Alignment retrieval
 153    ***********************************************/
 154
 155   /* Indexed retrieval:
 156    */
 157   if (afp->ssi != NULL) {
 158     if (! MSAFilePositionByKey(afp, key))
 159       Die("No such alignment %s found in file %s", key, afile);
 160     msa = MSAFileRead(afp);
 161   }
 162   /* Brute force retrieval:
 163    */
 164   else {
 165     while ((msa = MSAFileRead(afp)) != NULL)
 166       {
 167         if (strcmp(msa->name, key) == 0) break;
 168         if (strcmp(msa->acc,  key) == 0) break;
 169         MSAFree(msa);
 170       }
 171   }
 172
 173   if (msa == NULL) Die("Failed to retrieve %s from file %s", key, afile);
 174
 175   /* Output the alignment we retrieved
 176    */
 177   WriteStockholm(stdout, msa);
 178
 179   MSAFileClose(afp);
 180   MSAFree(msa);
 181   exit (0);
 182 }