extractseq.cpp

/*    extractseq: a program to extract sequence regions from
 *    chromosome files corresponding to genomic intervals specified in
 *    a BED file.
 *
 *    Copyright (C) 2009 University of Southern California and
 *                       Andrew D. Smith
 *
 *    Authors: Andrew D. Smith
 *
 *    This program is free software: you can redistribute it and/or modify
 *    it under the terms of the GNU General Public License as published by
 *    the Free Software Foundation, either version 3 of the License, or
 *    (at your option) any later version.
 *
 *    This program is distributed in the hope that it will be useful,
 *    but WITHOUT ANY WARRANTY; without even the implied warranty of
 *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *    GNU General Public License for more details.
 *
 *    You should have received a copy of the GNU General Public License
 *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
 */

#include <string>
#include <vector>
#include <iostream>
#include <fstream>

#include "OptionParser.hpp"
#include "smithlab_utils.hpp"
#include "smithlab_os.hpp"
#include "GenomicRegion.hpp"
#include "chromosome_utils.hpp"

using std::string;
using std::vector;
using std::cout;
using std::cerr;
using std::endl;
using std::ifstream;

int
main(int argc, const char **argv) {

  try {

    bool VERBOSE = false;

    string chrom_dir;
    string outfile;

    /****************** COMMAND LINE OPTIONS ********************/
    OptionParser opt_parse("extractseq", "program to extract sequence regions from "
         "chromosome files corresponding to genomic intervals "
         "specified in a BED file",
         "<bed-format-regions>");
    opt_parse.add_opt("output", 'o', "Name of output file (default: stdout)",
          false, outfile);
    opt_parse.add_opt("chrom", 'c', "directory with chrom files (FASTA format)",
          true , chrom_dir);
    opt_parse.add_opt("verbose", 'v', "print more run info", false, VERBOSE);
    vector<string> leftover_args;
    opt_parse.parse(argc, argv, leftover_args);
    if (argc == 1 || opt_parse.help_requested()) {
      cerr << opt_parse.help_message() << endl;
      return EXIT_SUCCESS;
    }
    if (opt_parse.about_requested()) {
      cerr << opt_parse.about_message() << endl;
      return EXIT_SUCCESS;
    }
    if (opt_parse.option_missing()) {
      cerr << opt_parse.option_missing_message() << endl;
      return EXIT_SUCCESS;
    }
    if (leftover_args.empty()) {
      cerr << opt_parse.help_message() << endl;
      return EXIT_SUCCESS;
    }
    const string regions_file = leftover_args.front();
    /****************** END COMMAND LINE OPTIONS *****************/

    if (VERBOSE)
      cerr << "loading regions" << endl;
    vector<GenomicRegion> regions;
    ReadBEDFile(regions_file, regions);

    if (VERBOSE)
      cerr << "extracting reference sequences" << endl;
    vector<string> region_sequences;
    extract_regions_fasta(chrom_dir, regions, region_sequences);

    std::ostream *out = (outfile.empty()) ? &cout :
      new std::ofstream(outfile.c_str());
    for (size_t i = 0; i < region_sequences.size(); ++i)
      *out << ">" << assemble_region_name(regions[i]) << endl
     << region_sequences[i] << endl;
    if (out != &cout) delete out;
  }
  catch (const SMITHLABException &e) {
    cerr << e.what() << endl;
    return EXIT_FAILURE;
  }
  catch (std::bad_alloc &ba) {
    cerr << "ERROR: could not allocate memory" << endl;
    return EXIT_FAILURE;
  }
  return EXIT_SUCCESS;
}