forester/archive/RIO/others/hmmer/src/hmmalign.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 /* hmmalign.c
  12  * SRE, Thu Dec 18 16:05:29 1997 [St. Louis]
  13  *
  14  * main() for aligning a set of sequences to an HMM.
  15  * RCS $Id: hmmalign.c,v 1.1.1.1 2005/03/22 08:34:00 cmzmasek Exp $
  16  */
  17
  18 #include <stdio.h>
  19 #include <stdlib.h>
  20
  21 #include "structs.h"            /* data structures, macros, #define's   */
  22 #include "config.h"             /* compile-time configuration constants */
  23 #include "funcs.h"              /* function declarations                */
  24 #include "globals.h"            /* alphabet global variables            */
  25 #include "squid.h"              /* general sequence analysis library    */
  26 #include "msa.h"                /* squid's multiple alignment i/o       */
  27
  28 static char banner[] = "hmmalign - align sequences to an HMM profile";
  29
  30 static char usage[]  = "\
  31 Usage: hmmalign [-options] <hmm file> <sequence file>\n\
  32 Available options are:\n\
  33    -h     : help; print brief help on version and usage\n\
  34    -m     : only print symbols aligned to match states\n\
  35    -o <f> : save alignment in file <f> in SELEX format\n\
  36    -q     : quiet - suppress verbose banner\n\
  37 ";
  38
  39 static char experts[] = "\
  40    --informat <s>: sequence file is in format <s>, not FASTA\n\
  41    --mapali <f>  : include alignment in file <f> using map in HMM\n\
  42    --withali <f> : include alignment to (fixed) alignment in file <f>\n\
  43 \n";
  44
  45 static struct opt_s OPTIONS[] = {
  46   { "-h", TRUE, sqdARG_NONE   },
  47   { "-m", TRUE, sqdARG_NONE   } ,
  48   { "-o", TRUE, sqdARG_STRING },
  49   { "-q", TRUE, sqdARG_NONE   },
  50   { "--informat",FALSE, sqdARG_STRING },
  51   { "--mapali",  FALSE, sqdARG_STRING },
  52   { "--withali", FALSE, sqdARG_STRING },
  53 };
  54 #define NOPTIONS (sizeof(OPTIONS) / sizeof(struct opt_s))
  55
  56 static void include_alignment(char *seqfile, struct plan7_s *hmm, int do_mapped,
  57                               char ***rseq, char ***dsq, SQINFO **sqinfo,
  58                               struct p7trace_s ***tr, int *nseq);
  59
  60 int
  61 main(int argc, char **argv)
  62 {
  63   char            *hmmfile;     /* file to read HMMs from                  */
  64   HMMFILE         *hmmfp;       /* opened hmmfile for reading              */
  65   struct plan7_s  *hmm;         /* HMM to align to                         */
  66   char            *seqfile;     /* file to read target sequence from       */
  67   int              format;      /* format of seqfile                       */
  68   char           **rseq;        /* raw, unaligned sequences                */
  69   SQINFO          *sqinfo;      /* info associated with sequences          */
  70   char           **dsq;         /* digitized raw sequences                 */
  71   int              nseq;        /* number of sequences                     */
  72   float           *wgt;         /* weights to assign to alignment          */
  73   MSA             *msa;         /* alignment that's created                */
  74   int              i;
  75   struct p7trace_s **tr;        /* traces for aligned sequences            */
  76
  77   char *optname;                /* name of option found by Getopt()         */
  78   char *optarg;                 /* argument found by Getopt()               */
  79   int   optind;                 /* index in argv[]                          */
  80   int   be_quiet;               /* TRUE to suppress verbose banner          */
  81   int   matchonly;              /* TRUE to show only match state syms       */
  82   char *outfile;                /* optional alignment output file           */
  83   FILE *ofp;                    /* handle on alignment output file          */
  84   char *withali;                /* name of additional alignment file to align */
  85   char *mapali;                 /* name of additional alignment file to map   */
  86
  87   /***********************************************
  88    * Parse command line
  89    ***********************************************/
  90
  91   format    = SQFILE_UNKNOWN;   /* default: autodetect format */
  92   matchonly = FALSE;
  93   outfile   = NULL;
  94   be_quiet  = FALSE;
  95   withali   = NULL;
  96   mapali    = NULL;
  97
  98   while (Getopt(argc, argv, OPTIONS, NOPTIONS, usage,
  99                 &optind, &optname, &optarg))  {
 100     if      (strcmp(optname, "-m")        == 0) matchonly= TRUE;
 101     else if (strcmp(optname, "-o")        == 0) outfile  = optarg;
 102     else if (strcmp(optname, "-q")        == 0) be_quiet = TRUE;
 103     else if (strcmp(optname, "--mapali")  == 0) mapali   = optarg;
 104     else if (strcmp(optname, "--withali") == 0) withali  = optarg;
 105     else if (strcmp(optname, "--informat") == 0) {
 106         format = String2SeqfileFormat(optarg);
 107         if (format == SQFILE_UNKNOWN)
 108           Die("unrecognized sequence file format \"%s\"", optarg);
 109       }
 110     else if (strcmp(optname, "-h") == 0)
 111       {
 112         Banner(stdout, banner);
 113         puts(usage);
 114         puts(experts);
 115         exit(EXIT_SUCCESS);
 116       }
 117   }
 118   if (argc - optind != 2)
 119     Die("Incorrect number of arguments.\n%s\n", usage);
 120
 121   hmmfile = argv[optind++];
 122   seqfile = argv[optind++];
 123
 124  /***********************************************
 125   * Open HMM file (might be in HMMERDB or current directory).
 126   * Read a single HMM from it.
 127   *
 128   * Currently hmmalign disallows the J state and
 129   * only allows one domain per sequence. To preserve
 130   * the S/W entry information, the J state is explicitly
 131   * disallowed, rather than calling a Plan7*Config() function.
 132   * this is a workaround in 2.1 for the 2.0.x "yo!" bug.
 133   ***********************************************/
 134
 135   if ((hmmfp = HMMFileOpen(hmmfile, "HMMERDB")) == NULL)
 136     Die("Failed to open HMM file %s\n%s", hmmfile, usage);
 137   if (!HMMFileRead(hmmfp, &hmm))
 138     Die("Failed to read any HMMs from %s\n", hmmfile);
 139   HMMFileClose(hmmfp);
 140   if (hmm == NULL)
 141     Die("HMM file %s corrupt or in incorrect format? Parse failed", hmmfile);
 142   hmm->xt[XTE][MOVE] = 1.;            /* only 1 domain/sequence ("global" alignment) */
 143   hmm->xt[XTE][LOOP] = 0.;
 144   P7Logoddsify(hmm, TRUE);
 145                                 /* do we have the map we might need? */
 146   if (mapali != NULL && ! (hmm->flags & PLAN7_MAP))
 147     Die("HMMER: HMM file %s has no map; you can't use --mapali.", hmmfile);
 148
 149   /***********************************************
 150    * Open sequence file in current directory.
 151    * Read all seqs from it.
 152    ***********************************************/
 153
 154   if (! ReadMultipleRseqs(seqfile, format, &rseq, &sqinfo, &nseq))
 155     Die("Failed to read any sequences from file %s", seqfile);
 156
 157   /***********************************************
 158    * Show the banner
 159    ***********************************************/
 160
 161   if (! be_quiet)
 162     {
 163       Banner(stdout, banner);
 164       printf(   "HMM file:             %s\n", hmmfile);
 165       printf(   "Sequence file:        %s\n", seqfile);
 166       printf("- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -\n\n");
 167     }
 168
 169   /***********************************************
 170    * Do the work
 171    ***********************************************/
 172
 173   /* Allocations and initializations.
 174    */
 175   dsq = MallocOrDie(sizeof(char *) * nseq);
 176   tr  = MallocOrDie(sizeof(struct p7trace_s *) * nseq);
 177
 178   /* Align each sequence to the model, collect traces
 179    */
 180   for (i = 0; i < nseq; i++)
 181     {
 182       dsq[i] = DigitizeSequence(rseq[i], sqinfo[i].len);
 183
 184       if (P7ViterbiSize(sqinfo[i].len, hmm->M) <= RAMLIMIT)
 185         (void) P7Viterbi(dsq[i], sqinfo[i].len, hmm, &(tr[i]));
 186       else
 187         (void) P7SmallViterbi(dsq[i], sqinfo[i].len, hmm, &(tr[i]));
 188     }
 189
 190   /* Include an aligned alignment, if desired.
 191    */
 192   if (mapali != NULL)
 193     include_alignment(mapali, hmm, TRUE, &rseq, &dsq, &sqinfo, &tr, &nseq);
 194   if (withali != NULL)
 195     include_alignment(withali, hmm, FALSE, &rseq, &dsq, &sqinfo, &tr, &nseq);
 196
 197   /* Turn traces into a multiple alignment
 198    */
 199   wgt = MallocOrDie(sizeof(float) * nseq);
 200   FSet(wgt, nseq, 1.0);
 201   msa = P7Traces2Alignment(dsq, sqinfo, wgt, nseq, hmm->M, tr, matchonly);
 202
 203   /***********************************************
 204    * Output the alignment
 205    ***********************************************/
 206
 207   if (outfile != NULL && (ofp = fopen(outfile, "w")) != NULL)
 208     {
 209       WriteStockholm(ofp, msa);
 210       printf("Alignment saved in file %s\n", outfile);
 211       fclose(ofp);
 212     }
 213   else
 214     WriteStockholm(stdout, msa);
 215
 216   /***********************************************
 217    * Cleanup and exit
 218    ***********************************************/
 219
 220   for (i = 0; i < nseq; i++)
 221     {
 222       P7FreeTrace(tr[i]);
 223       FreeSequence(rseq[i], &(sqinfo[i]));
 224       free(dsq[i]);
 225     }
 226   MSAFree(msa);
 227   FreePlan7(hmm);
 228   free(sqinfo);
 229   free(rseq);
 230   free(dsq);
 231   free(wgt);
 232   free(tr);
 233
 234   SqdClean();
 235   return 0;
 236 }
 237
 238
 239 /* Function: include_alignment()
 240  * Date:     SRE, Sun Jul  5 15:25:13 1998 [St. Louis]
 241  *
 242  * Purpose:  Given the name of a multiple alignment file,
 243  *           align that alignment to the HMM, and add traces
 244  *           to an existing array of traces. If do_mapped
 245  *           is TRUE, we use the HMM's map file. If not,
 246  *           we use P7ViterbiAlignAlignment().
 247  *
 248  * Args:     seqfile  - name of alignment file
 249  *           hmm      - model to align to
 250  *           do_mapped- TRUE if we're to use the HMM's alignment map
 251  *           rsq      - RETURN: array of rseqs to add to
 252  *           dsq      - RETURN: array of dsq to add to
 253  *           sqinfo   - RETURN: array of SQINFO to add to
 254  *           tr       - RETURN: array of traces to add to
 255  *           nseq     - RETURN: number of seqs
 256  *
 257  * Returns:  new, realloc'ed arrays for rsq, dsq, sqinfo, tr; nseq is
 258  *           increased to nseq+ainfo.nseq.
 259  */
 260 static void
 261 include_alignment(char *seqfile, struct plan7_s *hmm, int do_mapped,
 262                   char ***rsq, char ***dsq, SQINFO **sqinfo,
 263                   struct p7trace_s ***tr, int *nseq)
 264 {
 265   int format;                   /* format of alignment file */
 266   MSA   *msa;                   /* alignment to align to    */
 267   MSAFILE *afp;
 268   SQINFO  *newinfo;             /* sqinfo array from msa */
 269   char **newdsq;
 270   char **newrseq;
 271   int   idx;                    /* counter over aseqs       */
 272   struct p7trace_s *master;     /* master trace             */
 273   struct p7trace_s **addtr;     /* individual traces for aseq */
 274
 275   format = MSAFILE_UNKNOWN;     /* invoke Babelfish */
 276   if ((afp = MSAFileOpen(seqfile, format, NULL)) == NULL)
 277     Die("Alignment file %s could not be opened for reading", seqfile);
 278   if ((msa = MSAFileRead(afp)) == NULL)
 279     Die("Failed to read an alignment from %s\n", seqfile);
 280   MSAFileClose(afp);
 281   for (idx = 0; idx < msa->nseq; idx++)
 282     s2upper(msa->aseq[idx]);
 283   newinfo = MSAToSqinfo(msa);
 284
 285                                 /* Verify checksums before mapping */
 286   if (do_mapped && GCGMultchecksum(msa->aseq, msa->nseq) != hmm->checksum)
 287     Die("The checksums for alignment file %s and the HMM alignment map don't match.",
 288         seqfile);
 289                                 /* Get a master trace */
 290   if (do_mapped) master = MasterTraceFromMap(hmm->map, hmm->M, msa->alen);
 291   else           master = P7ViterbiAlignAlignment(msa, hmm);
 292
 293                                 /* convert to individual traces */
 294   ImposeMasterTrace(msa->aseq, msa->nseq, master, &addtr);
 295                                 /* add those traces to existing ones */
 296   *tr = MergeTraceArrays(*tr, *nseq, addtr, msa->nseq);
 297
 298                                 /* additional bookkeeping: add to dsq, sqinfo */
 299   *rsq = ReallocOrDie((*rsq), sizeof(char *) * (*nseq + msa->nseq));
 300   DealignAseqs(msa->aseq, msa->nseq, &newrseq);
 301   for (idx = *nseq; idx < *nseq + msa->nseq; idx++)
 302     (*rsq)[idx] = newrseq[idx - (*nseq)];
 303   free(newrseq);
 304
 305   *dsq = ReallocOrDie((*dsq), sizeof(char *) * (*nseq + msa->nseq));
 306   DigitizeAlignment(msa, &newdsq);
 307   for (idx = *nseq; idx < *nseq + msa->nseq; idx++)
 308     (*dsq)[idx] = newdsq[idx - (*nseq)];
 309   free(newdsq);
 310                                 /* unnecessarily complex, but I can't be bothered... */
 311   *sqinfo = ReallocOrDie((*sqinfo), sizeof(SQINFO) * (*nseq + msa->nseq));
 312   for (idx = *nseq; idx < *nseq + msa->nseq; idx++)
 313     SeqinfoCopy(&((*sqinfo)[idx]), &(newinfo[idx - (*nseq)]));
 314
 315   *nseq = *nseq + msa->nseq;
 316
 317                                 /* Cleanup */
 318   P7FreeTrace(master);
 319   MSAFree(msa);
 320                                 /* Return */
 321   return;
 322 }
 323
 324
 325