forester/archive/RIO/others/hmmer/squid/seqsplit_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
  12 /* seqsplit_main.c
  13  * SRE, Mon Sep 25 11:43:58 2000
  14  *
  15  * Split sequences into smaller chunks of defined size and overlap;
  16  * output a FASTA file.
  17  *
  18  * Limitations:
  19  *   still working in 32 bits -- no sequence can be more than 2 GB
  20  *   in size.
  21  * CVS $Id: seqsplit_main.c,v 1.1.1.1 2005/03/22 08:34:26 cmzmasek Exp $
  22  */
  23
  24 #include <stdio.h>
  25 #include <string.h>
  26 #include "squid.h"
  27 #include "msa.h"
  28
  29 static char banner[] = "seqsplit - split seqs into chunks of defined size and overlap";
  30
  31 static char usage[]  = "\
  32 Usage: seqsplit [-options] <seqfile>\n\
  33   Available options:\n\
  34   -h        : help; display usage and version\n\
  35   -o <file> : output the new FASTA file to <file>\n\
  36 ";
  37
  38 static char experts[] = "\
  39   --informat <s> : specify sequence file format <s>\n\
  40   --length <n>   : set max length of each unique seq frag to <n>\n\
  41   --overlap <n>  : set overlap length to <n> (total frag size = length+overlap)\n\
  42 ";
  43
  44 struct opt_s OPTIONS[] = {
  45   { "-h", TRUE, sqdARG_NONE },
  46   { "-o", TRUE, sqdARG_STRING },
  47   { "--informat", FALSE, sqdARG_STRING },
  48   { "--length",   FALSE, sqdARG_INT },
  49   { "--overlap",  FALSE, sqdARG_INT },
  50 };
  51 #define NOPTIONS (sizeof(OPTIONS) / sizeof(struct opt_s))
  52
  53
  54 int
  55 main(int argc, char **argv)
  56 {
  57   char     *seqfile;            /* name of sequence file     */
  58   char     *outfile;            /* name of output file       */
  59   SQFILE   *dbfp;               /* open sequence file        */
  60   FILE     *ofp;                /* open output file          */
  61   int       fmt;                /* format of seqfile         */
  62   char     *seq;                /* sequence                  */
  63   SQINFO    sqinfo;             /* extra info about sequence */
  64   char     *seqfrag;            /* space for a seq fragment  */
  65   int       fraglength;         /* length of unique seq per frag     */
  66   int       overlap;            /* length of overlap. frags are fraglength+overlap*/
  67   char      seqname[256];       /* renamed fragment, w/ coord info */
  68   int       num;                /* number of this fragment   */
  69   int       pos;                /* position in a sequence */
  70   int       len;                /* length of a fragment   */
  71   char     *desc;
  72
  73   int       nseqs;              /* total number of sequences */
  74   int       nsplit;             /* number of seqs that get split */
  75   int       nnewfrags;          /* total number of new fragments */
  76
  77   char  *optname;
  78   char  *optarg;
  79   int    optind;
  80
  81   /***********************************************
  82    * Parse command line
  83    ***********************************************/
  84
  85   fmt       = SQFILE_UNKNOWN;   /* default: autodetect      */
  86   fraglength = 100000;
  87   overlap   = 1000;
  88   outfile   = NULL;
  89
  90   while (Getopt(argc, argv, OPTIONS, NOPTIONS, usage,
  91                 &optind, &optname, &optarg))
  92     {
  93       if      (strcmp(optname, "-o")         == 0)  outfile    = optarg;
  94       else if (strcmp(optname, "--length")   == 0)  fraglength = atoi(optarg);
  95       else if (strcmp(optname, "--overlap")  == 0)  overlap    = atoi(optarg);
  96       else if (strcmp(optname, "--informat") == 0) {
  97         fmt = String2SeqfileFormat(optarg);
  98         if (fmt == SQFILE_UNKNOWN)
  99           Die("unrecognized sequence file format \"%s\"", optarg);
 100       }
 101       else if (strcmp(optname, "-h") == 0) {
 102         Banner(stdout, banner);
 103         puts(usage);
 104         puts(experts);
 105         exit(EXIT_SUCCESS);
 106       }
 107     }
 108
 109   if (argc - optind != 1) Die("%s\n", usage);
 110   seqfile = argv[argc-1];
 111
 112   seqfrag = MallocOrDie(sizeof(char) * (fraglength+overlap));
 113   seqfrag[fraglength+overlap] = '\0';
 114
 115   /***********************************************
 116    * Read the file.
 117    ***********************************************/
 118
 119   if (outfile == NULL)  ofp = stdout;
 120   else {
 121     if ((ofp = fopen(outfile, "w")) == NULL)
 122       Die("Failed to open output sequence file %s for writing", outfile);
 123   }
 124
 125   if ((dbfp = SeqfileOpen(seqfile, fmt, NULL)) == NULL)
 126     Die("Failed to open sequence file %s for reading", seqfile);
 127
 128   nseqs = nsplit = nnewfrags = 0;
 129   while (ReadSeq(dbfp, dbfp->format, &seq, &sqinfo))
 130     {
 131       nseqs++;
 132       if (sqinfo.flags & SQINFO_DESC) desc = sqinfo.desc;
 133       else desc = NULL;
 134
 135       if (sqinfo.len <= fraglength+overlap) {
 136         WriteSimpleFASTA(ofp, seq, sqinfo.name, desc);
 137         continue;
 138       }
 139
 140       num = 1;
 141       nsplit++;
 142       for (pos = 0; pos < sqinfo.len; pos += fraglength)
 143         {
 144           if (sqinfo.len - pos <= overlap) continue;
 145           strncpy(seqfrag, seq+pos, fraglength+overlap);
 146           len = strlen(seqfrag);
 147           sprintf(seqname, "%s/frag%d/%d-%d",
 148                   sqinfo.name, num, pos+1, pos+len);
 149           WriteSimpleFASTA(ofp, seqfrag, seqname, desc);
 150           nnewfrags++;
 151           num ++;
 152         }
 153       FreeSequence(seq, &sqinfo);
 154     }
 155   SeqfileClose(dbfp);
 156   if (outfile != NULL) fclose(ofp);
 157
 158   printf("Total # of seqs:         %d\n", nseqs);
 159   printf("Affected by splitting:   %d\n", nsplit);
 160   printf("New # of seqs:           %d\n", nseqs-nsplit + nnewfrags);
 161
 162   return 0;
 163 }