src/main.cpp

#include <boost/unordered_map.hpp>
#include <boost/unordered_set.hpp>
#include <boost/filesystem.hpp>
#include <boost/program_options.hpp>
#include <boost/thread.hpp>
#include <iostream>
#include <iomanip>
#include <fstream>
#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include "main.h"
#include "bundles.h"
#include "targets.h"
#include "lengthdistribution.h"
#include "fragments.h"
#include "biascorrection.h"
#include "mismatchmodel.h"
#include "mapparser.h"
#include "threadsafety.h"
#include "robertsfilter.h"
#include "directiondetector.h"
#include "library.h"

#ifdef PROTO
  #include PROTO_ALIGNMENT_INCL
  #include PROTO_TARGET_INCL
  #include <boost/archive/iterators/base64_from_binary.hpp>
  #include <boost/archive/iterators/transform_width.hpp>
  #include <boost/archive/iterators/ostream_iterator.hpp>
#endif

#ifndef WIN32
  #include "update_check.h"
#endif

using namespace std;
namespace po = boost::program_options;
namespace fs = boost::filesystem;

Logger logger;

// the forgetting factor parameter controls the growth of the fragment mass
double ff_param = 0.85;

// the burn-in parameter determines how many reads are required before the
// error and bias models are applied to probabilistic assignment
size_t burn_in = 100000;
size_t burn_out = 5000000;
bool burned_out = false;

size_t max_read_len = 250;

size_t max_indel_size = 10;

size_t stop_at = 0;

// file location parameters
string output_dir = ".";
string fasta_file_name = "";
string in_map_file_names = "";
string param_file_name = "";
string haplotype_file_name = "";

// intial pseudo-count parameters (non-logged)
double expr_alpha = .005;
double fld_alpha = 1;
double bias_alpha = 1;
double mm_alpha = 1;

size_t bias_model_order = 3;

// fragment length parameters
size_t def_fl_max = 800;
size_t def_fl_mean = 200;
size_t def_fl_stddev = 80;
size_t def_fl_kernel_n = 4;
double def_fl_kernel_p = 0.5;

// option parameters
bool edit_detect = false;
bool error_model = true;
bool bias_correct = true;
bool calc_covar = false;
bool output_align_prob = false;
bool output_align_samp = false;
bool output_running_rounds = false;
bool output_running_reads = false;
size_t num_threads = 2;
size_t num_neighbors = 0;
size_t library_size = 0;

// directional parameters
Direction direction = BOTH;

bool running = true;

// used for multiple rounds of EM
bool first_round = true;
bool last_round = true;
bool batch_mode = false;
bool online_additional = false;
bool both = false;
size_t remaining_rounds = 0;

bool spark_pre = false;

typedef boost::unordered_map<string, double> AlphaMap;
AlphaMap* expr_alpha_map = NULL;

AlphaMap* parse_priors(string in_file) {
  ifstream ifs(in_file.c_str());
  if (!ifs.is_open()) {
    logger.severe("Unable to open input priors file '%s'.", in_file.c_str());
  }
  AlphaMap* alphas = new AlphaMap();

  string line;

  while(ifs.good()) {
    getline(ifs,line);

    size_t idx = line.find_first_of("\t ");
    if (idx!=string::npos) {
      string name = line.substr(0,idx);
      string val = line.substr(idx+1);
      (*alphas)[name] = atof(val.c_str());
    }
  }
  return alphas;
};

bool parse_options(int ac, char ** av) {

  size_t additional_online = 0;
  size_t additional_batch = 0;
  
  po::options_description standard("Standard Options");
  standard.add_options()
  ("help,h", "produce help message")
  ("output-dir,o", po::value<string>(&output_dir)->default_value(output_dir),
   "write all output files to this directory")
#ifdef PROTO
  ("preprocess,D", "run preprocess script for eXpressD")
#endif
  ("frag-len-mean,m", po::value<size_t>(&def_fl_mean)->default_value(def_fl_mean),
   "prior estimate for average fragment length")
  ("frag-len-stddev,s",
   po::value<size_t>(&def_fl_stddev)->default_value(def_fl_stddev),
   "prior estimate for fragment length std deviation")
  ("haplotype-file,H",
   po::value<string>(&haplotype_file_name)->default_value(haplotype_file_name),
   "path to a file containing haplotype pairs")
  ("additional-batch,B",
   po::value<size_t>(&additional_batch)->default_value(additional_batch),
   "number of additional batch EM rounds after initial online round")
  ("additional-online,O",
   po::value<size_t>(&additional_online)->default_value(additional_online),
   "number of additional online EM rounds after initial online round")
  ("max-read-len,L",
   po::value<size_t>(&max_read_len)->default_value(max_read_len),
   "maximum allowed length of a read")
  ("output-align-prob",
   "output alignments (sam/bam) with probabilistic assignments")
  ("output-align-samp",
   "output alignments (sam/bam) with sampled assignments")
  ("fr-stranded",
   "accept only forward->reverse alignments (second-stranded protocols)")
  ("rf-stranded",
   "accept only reverse->forward alignments (first-stranded protocols)")
  ("f-stranded",
   "accept only forward single-end alignments (second-stranded protocols)")
  ("r-stranded",
   "accept only reverse single-end alignments (first-stranded protocols)")
  ("no-update-check", "disables automatic check for update via web")
  ("logtostderr", "prints all logging messages to stderr")
  ;
  
  po::options_description advanced("Advanced Options");
  advanced.add_options()
  ("forget-param,f", po::value<double>(&ff_param)->default_value(ff_param),
   "sets the 'forgetting factor' parameter (0.5 < c <= 1)")
  ("library-size", po::value<size_t>(&library_size),
   "specifies library size for FPKM instead of calculating from alignments")
  ("max-indel-size",
   po::value<size_t>(&max_indel_size)->default_value(max_indel_size),
   "sets the maximum allowed indel size, affecting geometric indel prior")
  ("calc-covar", "calculate and output covariance matrix")
  ("expr-alpha", po::value<double>(&expr_alpha)->default_value(expr_alpha),
   "sets the strength of the prior, per bp")
  ("stop-at", po::value<size_t>(&stop_at)->default_value(stop_at),
   "sets the number of fragments to process, disabled with 0")
  ("burn-out", po::value<size_t>(&burn_out)->default_value(burn_out),
   "sets number of fragments after which to stop updating auxiliary parameters")
  ("no-bias-correct", "disables bias correction")
  ("no-error-model", "disables error modelling")
  ("aux-param-file",
   po::value<string>(&param_file_name)->default_value(param_file_name),
   "path to file containing auxiliary parameters to use instead of learning")
  ;

  string prior_file = "";

  po::options_description hidden("Experimental/Debug Options");
  hidden.add_options()
  ("num-threads,p", po::value<size_t>(&num_threads)->default_value(num_threads),
   "number of threads (>= 2)")
  ("edit-detect","")
  ("single-round", "")
  ("output-running-rounds", "")
  ("output-running-reads", "")
  ("batch-mode","")
  ("both","")
  ("prior-params", po::value<string>(&prior_file)->default_value(""), "")
  ("sam-file", po::value<string>(&in_map_file_names)->default_value(""), "")
  ("fasta-file", po::value<string>(&fasta_file_name)->default_value(""), "")
  ("num-neighbors", po::value<size_t>(&num_neighbors)->default_value(0), "")
  ("bias-model-order",
   po::value<size_t>(&bias_model_order)->default_value(bias_model_order),
   "sets the order of the Markov chain used to model sequence bias")
  ;

  po::positional_options_description positional;
  positional.add("fasta-file",1).add("sam-file",1);

  po::options_description cmdline_options;
  cmdline_options.add(standard).add(advanced).add(hidden);

  bool error = false;
  po::variables_map vm;
  try {
    po::store(po::command_line_parser(ac, av).options(cmdline_options)
              .positional(positional).run(), vm);
  } catch (po::error& e) {
    logger.info("Command-Line Argument Error: %s.", e.what());
    error = true;
  }
  po::notify(vm);

  if (ff_param > 1.0 || ff_param < 0.5) {
    logger.info("Command-Line Argument Error: forget-param/f option must be "
                "between 0.5 and 1.0.");
    error= true;
  }

  if (fasta_file_name == "") {
    logger.info("Command-Line Argument Error: target sequence fasta file "
                "required.");
    error = true;
  }

  if (error || vm.count("help")) {
    cerr << "express v" << PACKAGE_VERSION << endl
         << "-----------------------------\n"
         << "File Usage:  express [options] <target_seqs.fa> <hits.(sam/bam)>\n"
         << "Piped Usage: bowtie [options] -S <index> <reads.fq> | express "
         << "[options] <target_seqs.fa>\n\n"
         << "Required arguments:\n"
         << " <target_seqs.fa>     target sequence file in fasta format\n"
         << " <hits.(sam/bam)>     read alignment file in SAM or BAM format\n\n"
         << standard
         << advanced;
    return 1;
  }

  if (param_file_name.size()) {
    burn_in = 0;
    burn_out = 0;
    burned_out = true;
  }
  
  size_t stranded_count = 0;
  if (vm.count("fr-stranded")) {
    direction = FR;
    stranded_count++;
  }
  if (vm.count("rf-stranded")) {
    direction = RF;
    stranded_count++;
  }
  if (vm.count("f-stranded")) {
    direction = F;
    stranded_count++;
  }
  if (vm.count("r-stranded")) {
    direction = R;
    stranded_count++;
  }
  if (stranded_count > 1) {
    logger.severe("Multiple strandedness flags cannot be specified in the same "
                  "run.");
  }
  if (vm.count("logtostderr")) {
    logger.info_out(&cerr);
  }
  
  edit_detect = vm.count("edit-detect");
  calc_covar = vm.count("calc-covar");
  bias_correct = !(vm.count("no-bias-correct"));
  error_model = !(vm.count("no-error-model"));
  output_align_prob = vm.count("output-align-prob");
  output_align_samp = vm.count("output-align-samp");
  output_running_rounds = vm.count("output-running-rounds");
  output_running_reads = vm.count("output-running-reads");
  batch_mode = vm.count("batch-mode");
  both = vm.count("both");
  remaining_rounds = max(additional_online, additional_batch);
  spark_pre = vm.count("preprocess");

  if (batch_mode) {
    ff_param = 1;
  }
  
  if (additional_online > 0 && additional_batch > 0) {
    logger.severe("Cannot add both online and batch rounds.");
  } else if (additional_online > 0) {
    online_additional = true;
  }
  
  if (output_align_prob && output_align_samp) {
    logger.severe("Cannot output both alignment probabilties and sampled "
                  "alignments.");
  }
  if ((output_align_prob || output_align_samp) && remaining_rounds == 0) {
    logger.warn("It is recommended that at least one additional round "
                "be used when outputting alignment probabilities or sampled "
                "alignments. Use the '-B' or '-O' option to enable.");
  }
  
  // We have 1 processing thread and 1 parsing thread always, so we should not
  // count these as additional threads.
  if (num_threads < 2) {
    num_threads = 0;
  }
  num_threads -= 2;
  if (num_threads > 0) {
    num_threads -= edit_detect;
  }
  if (remaining_rounds && in_map_file_names == "") {
    logger.severe("Cannot process multiple rounds from streaming input.");
  }
  if (remaining_rounds) {
    last_round = false;
  }
  if (prior_file != "") {
    expr_alpha_map = parse_priors(prior_file);
  }

#ifndef WIN32
  if (!vm.count("no-update-check")) {
    check_version(PACKAGE_VERSION);
  }
#endif

  return 0;
}

void output_results(Librarian& libs, size_t tot_counts, int n=-1) {
  char buff[500];
  string dir = output_dir;
  if (n >= 0) {
    sprintf(buff, "%s/x_%d", output_dir.c_str(), n);
    logger.info("Writing results to %s.", buff);
    dir = string(buff);
    try {
      fs::create_directories(dir);
    } catch (fs::filesystem_error& e) {
      logger.severe(e.what());
    }
  }
  libs[0].targ_table->output_results(dir, tot_counts, last_round&calc_covar,
                                     last_round&edit_detect);

  for (size_t l = 0; l < libs.size(); l++) {
    if (libs.size() > 1) {
      sprintf(buff, "%s/params.%d.xprs", dir.c_str(), (int)l+1);
    } else {
      sprintf(buff, "%s/params.xprs", dir.c_str());
    }
    ofstream paramfile(buff);
    (libs[l].fld)->append_output(paramfile, "Fragment");
    if (libs[l].mismatch_table) {
      (libs[l].mismatch_table)->append_output(paramfile);
    }
    if (libs[l].bias_table) {
      (libs[l].bias_table)->append_output(paramfile);
    }
    paramfile.close();
  }
}

void process_fragment(Fragment* frag_p) {
  Fragment& frag = *frag_p;
  const Library& lib = *frag.lib();

  // sort hits to avoid deadlock
  frag.sort_hits();
  double mass_n = frag.mass();

  assert(frag.num_hits());

  vector<double> masses(frag.num_hits(), 0);
  vector<double> variances(frag.num_hits(), 0);
  double total_likelihood = LOG_0;
  double total_mass = LOG_0;
  double total_variance = LOG_0;
  size_t num_solvable = 0;

  // set of locked targets
  boost::unordered_set<const Target*> targ_set;
  boost::unordered_set<const Target*> locked_set;

  // Update bundles and merge in first loop
  Bundle* bundle = frag.hits()[0]->target()->bundle();
  
  if (frag.num_hits() > 1) {
    // Calculate marginal likelihoods and lock targets.
    for (size_t i = 0; i < frag.num_hits(); ++i) {
      FragHit& m = *frag.hits()[i];
      Target* t = m.target();
      
      bundle = lib.targ_table->merge_bundles(bundle, t->bundle());
      t->bundle(bundle);
      
      if (locked_set.count(t) == 0) {
        t->lock();
        locked_set.insert(t);
      }
      targ_set = locked_set;
      // lock neighbors
      foreach (const Target* neighbor, *m.neighbors()) {
        if (locked_set.count(neighbor) == 0) {
          neighbor->lock();
          locked_set.insert(neighbor);
        }
      }
      m.params()->align_likelihood = t->align_likelihood(m);
      m.params()->full_likelihood = m.params()->align_likelihood +
                                    t->sample_likelihood(first_round,
                                                         m.neighbors());
      masses[i] = t->mass();
      variances[i] = t->mass_var();
      total_likelihood = log_add(total_likelihood, m.params()->full_likelihood);
      total_mass = log_add(total_mass, masses[i]);
      total_variance = log_add(total_variance, variances[i]);
      num_solvable += t->solvable();
      assert(!isnan(total_likelihood));
    }
  } else {
    FragHit& m = *frag.hits()[0];
    Target* t = m.target();
    t->lock();
    locked_set.insert(t);
    total_likelihood = 0;
    m.params()->align_likelihood = 0;
    m.params()->full_likelihood = 0;
    foreach (const Target* neighbor, *frag.hits()[0]->neighbors()) {
      if (targ_set.count(neighbor) == 0) {
        neighbor->lock();
        locked_set.insert(neighbor);
      }
    }
  }

  if (islzero(total_likelihood)){
    assert(expr_alpha_map);
    logger.warn("Fragment '%s' has 0 likelihood of originating from the "
                "transcriptome. Skipping...", frag.name().c_str());
    foreach (const Target* t, locked_set) {
      t->unlock();
    }
    return;
  }

  if (first_round) {
    bundle->incr_counts();
  }
  if (first_round || online_additional) {
    bundle->incr_mass(mass_n);
  }
  
  // normalize marginal likelihoods
  for (size_t i = 0; i < frag.num_hits(); ++i) {
    FragHit& m = *frag[i];
    Target* t  = m.target();
    
    double p = m.params()->full_likelihood-total_likelihood;
    m.params()->posterior = p;
    if (targ_set.size() > 1) {
      double v = log_add(variances[i] - 2*total_mass,
                  total_variance + 2*masses[i] - 4*total_mass);
      t->add_hit(m, v, mass_n);
    } else if (i == 0) {
      t->add_hit(m, LOG_0, mass_n);
    }

    // update parameters
    if (first_round) {
      double r = rand()/double(RAND_MAX);
      
      if (i == 0 || frag[i-1]->target_id() != t->id()) {
        t->incr_counts(targ_set.size() <= 1);
      }
      if (!t->solvable() && num_solvable == frag.num_hits()-1) {
        t->solvable(true);
      }
      if (edit_detect && lib.mismatch_table) {
        (lib.mismatch_table)->update(m, p, lib.mass_n);
      }
      if (!burned_out && r < sexp(p)) {
        if (lib.mismatch_table && !edit_detect) {
          (lib.mismatch_table)->update(m, LOG_1, lib.mass_n);
        }
        if (m.pair_status() == PAIRED) {
          (lib.fld)->add_val(m.length(), lib.mass_n);
        }
        if (lib.bias_table) {
          (lib.bias_table)->update_observed(m, lib.mass_n);
        }
      }
    }
    if (calc_covar && (last_round || online_additional)) {
      double var = 2*mass_n + p + log_sub(LOG_1, p);
      lib.targ_table->update_covar(m.target_id(), m.target_id(), var);
      for (size_t j = i+1; j < frag.num_hits(); ++j) {
        const FragHit& m2 = *frag.hits()[j];
        double p2 = m2.params()->full_likelihood-total_likelihood;
        if (sexp(p2) == 0) {
          continue;
        }
        double covar = 2*mass_n + p + p2;
        lib.targ_table->update_covar(m.target_id(), m2.target_id(), covar);
      }
    }
  }

  foreach (const Target* t, locked_set) {
    t->unlock();
  }
}

void proc_thread(ParseThreadSafety* pts) {
  while (true) {
    Fragment* frag = pts->proc_on.pop();
    if (!frag) {
      break;
    }
    process_fragment(frag);
    pts->proc_out.push(frag);
  }
}

size_t threaded_calc_abundances(Librarian& libs) {
  logger.info("Processing input fragment alignments...");
  boost::scoped_ptr<boost::thread> bias_update;

  size_t n = 1;
  size_t num_frags = 0;
  double mass_n = 0;

  // For log-scale output
  size_t i = 1;
  size_t j = 6;

  DirectionDetector dir_detector;
  Fragment* frag;
  
  while (true) {
    // Loop through libraries
    for (size_t l = 0; l < libs.size(); l++) {
      Library& lib = libs[l];
      libs.set_curr(l);
      MapParser& map_parser = *lib.map_parser;
      boost::mutex bu_mut;
      // Used to signal bias update thread
      running = true;
      ParseThreadSafety pts(max((int)num_threads,10));
      boost::thread parse(&MapParser::threaded_parse, &map_parser, &pts,
                          stop_at, num_neighbors);
      vector<boost::thread*> thread_pool;
      RobertsFilter frags_seen;

      burned_out = lib.n >= burn_out;
      while(true) {
        if (lib.n == burn_in) {
          bias_update.reset(new boost::thread(&TargetTable::asynch_bias_update,
                                              lib.targ_table, &bu_mut));
          if (lib.mismatch_table) {
            (lib.mismatch_table)->activate();
          }
        }
        if (lib.n == burn_out) {
          if (lib.mismatch_table) {
            (lib.mismatch_table)->fix();
          };
          burned_out = true;
        }
        // Start threads once aux parameters are burned out
        if (burned_out && num_threads && thread_pool.size() == 0) {
          lib.targ_table->enable_bundle_threadsafety();
          thread_pool = vector<boost::thread*>(num_threads);
          for (size_t k = 0; k < thread_pool.size(); k++) {
            thread_pool[k] = new boost::thread(proc_thread, &pts);
          }
        }

        // Pop next parsed fragment and set mass
        frag = pts.proc_in.pop();
        if (frag) {
          frag->mass(mass_n);
          dir_detector.add_fragment(frag);
        }

        // Test that we have not already seen this fragment
        if (frag && first_round && frags_seen.test_and_push(frag->name())) {
          logger.severe("Alignments are not properly sorted. Read '%s' has "
                        "alignments which are non-consecutive.",
                        frag->name().c_str());
        }

        // If multi-threaded and burned out, push to the processing queue
        if (num_threads && burned_out) {
          // If no more fragments, send stop signal (NULL) to processing threads
          if (!frag) {
            for (size_t k = 0; k < thread_pool.size(); ++k) {
              pts.proc_on.push(NULL);
            }
            break;
          }
          pts.proc_on.push(frag);
        } else {
          if (!frag) {
            break;
          }
          {
            // Block the bias update thread from updating the paramater tables
            // during processing. We don't need to do this during multi-threaded
            // processing since the parameters are burned out before we start
            // the threads.
            boost::unique_lock<boost::mutex> lock(bu_mut);
            process_fragment(frag);
            pts.proc_out.push(frag);
          }
        }

        // Output intermediate results, if necessary
        if (output_running_reads && n == i*pow(10.,(double)j)) {
          boost::unique_lock<boost::mutex> lock(bu_mut);
          output_results(libs, n, (int)n);
          if (i++ == 9) {
            i = 1;
            j++;
          }
        }
        num_frags++;

        // Output progress
        if (num_frags % 1000000 == 0) {
          logger.info("Fragments Processed (%s): %d\tNumber of Bundles: %d.",
                      lib.in_file_name.c_str(), num_frags,
                      lib.targ_table->num_bundles());
          dir_detector.report_if_improper_direction();
        }

        n++;
        lib.n++;
        mass_n += ff_param*log((double)n-1) - log(pow(n,ff_param) - 1);
        lib.mass_n += ff_param*log((double)lib.n-1) -
                      log(pow(lib.n,ff_param) - 1);
      }

      // Signal bias update thread to stop
      running = false;

      parse.join();
      foreach(boost::thread* t, thread_pool) {
        t->join();
      }

      lib.targ_table->disable_bundle_threadsafety();
      lib.targ_table->collapse_bundles();
      
      if (bias_update) {
        logger.info("Waiting for auxiliary parameter update to complete...");
        bias_update->join();
        bias_update.reset(NULL);
      }
    }

    if (online_additional && remaining_rounds--) {
      if (output_running_rounds) {
        output_results(libs, n, (int)remaining_rounds);
      }

      logger.info("%d remaining rounds.", remaining_rounds);
      first_round = false;
      last_round = (remaining_rounds==0 && !both);
      for (size_t l = 0; l < libs.size(); l++) {
        libs[l].map_parser->write_active(last_round);
        libs[l].map_parser->reset_reader();
      }
      num_frags = 0;
    } else {
      break;
    }
  }

  logger.info("COMPLETED: Processed %d mapped fragments, targets are in %d "
              "bundles.", num_frags, libs[0].targ_table->num_bundles());

  return num_frags;
}

int estimation_main() {
  
  if (output_dir != ".") {
    try {
      fs::create_directories(output_dir);
    } catch (fs::filesystem_error& e) {
      logger.info(e.what());
    }
  }
  
  if (!fs::exists(output_dir)) {
    logger.severe("Cannot create directory %s.", output_dir.c_str());
  }
  
  // Parse input file names and instantiate Libray structs.
  vector<string> file_names;
  char buff[999];
  strcpy(buff, in_map_file_names.c_str());
  char * pch = strtok (buff,",");
  while (pch != NULL) {
    file_names.push_back(pch);
    pch = strtok (NULL, ",");
  }
  if (file_names.size() == 0) {
    file_names.push_back("");
  }
  Librarian libs(file_names.size());
  for (size_t i = 0; i < file_names.size(); ++i) {
    char out_map_file_name[500] = "";
    if (output_align_prob) {
      sprintf(out_map_file_name, "%s/hits.%d.prob",
              output_dir.c_str(), (int)i+1);
    }
    if (output_align_samp) {
      sprintf(out_map_file_name, "%s/hits.%d.samp",
              output_dir.c_str(), (int)i+1);
    }
    
    libs[i].in_file_name = file_names[i];
    libs[i].out_file_name = out_map_file_name;
    libs[i].map_parser.reset(new MapParser(&libs[i], last_round));

    if (param_file_name.size()) {
      libs[i].fld.reset(new LengthDistribution(param_file_name, "Fragment"));
      libs[i].mismatch_table.reset((error_model) ?
                                              new MismatchTable(param_file_name)
                                              : NULL);
      libs[i].bias_table.reset((bias_correct) ? new BiasBoss(bias_model_order,
                                                         param_file_name):NULL);
    } else {
      libs[i].fld.reset(new LengthDistribution(fld_alpha, def_fl_max,
                                               def_fl_mean, def_fl_stddev,
                                               def_fl_kernel_n,
                                               def_fl_kernel_p));
      libs[i].mismatch_table.reset((error_model) ? new MismatchTable(mm_alpha)
                                                   :NULL);
      libs[i].bias_table.reset((bias_correct) ? new BiasBoss(bias_model_order,
                                                    bias_alpha):NULL);
    }
    if (i > 0 &&
        (libs[i].map_parser->targ_index() != libs[i-1].map_parser->targ_index()
         || libs[i].map_parser->targ_lengths() !=
         libs[i-1].map_parser->targ_lengths())) {
          logger.severe("Alignment file headers do not match for '%s' and '%s'.",
                        file_names[i-1].c_str(), file_names[i].c_str());
        }
  }
  
  boost::shared_ptr<TargetTable> targ_table(
                                  new TargetTable(fasta_file_name,
                                                  haplotype_file_name,
                                                  edit_detect,
                                                  param_file_name.size(),
                                                  expr_alpha, expr_alpha_map,
                                                  &libs));
  size_t max_target_length = 0;
  for(size_t tid=0; tid < targ_table->size(); tid++) {
    max_target_length = max(max_target_length,
                            targ_table->get_targ(tid)->length());
  }

  for (size_t i = 0; i < libs.size(); ++i) {
    libs[i].targ_table = targ_table;
    if (bias_correct) {
      libs[i].bias_table->copy_expectations(*(libs.curr_lib().bias_table));
    }
  }
  double num_targ = (double)targ_table->size();
  
  if (calc_covar && (double)SSIZE_MAX < num_targ*(num_targ+1)) {
    logger.warn("Your system is unable to represent large enough values for "
                "efficiently hashing target pairs. Covariance calculation will "
                "be disabled.");
    calc_covar = false;
  }
  
  if (batch_mode) {
    targ_table->round_reset();
  }
  
  size_t tot_counts = threaded_calc_abundances(libs);
  if (library_size) {
    tot_counts = library_size;
  }
  
  if (!burned_out && bias_correct && param_file_name == "") {
    logger.warn("Not enough fragments observed to accurately learn bias "
                "parameters. Either disable bias correction "
                "(--no-bias-correct) or provide a file containing auxiliary "
                "parameters (--aux-param-file).");
  }
  
  if (both) {
    remaining_rounds = 1;
    online_additional = false;
  }
  
  if (remaining_rounds) {
    targ_table->masses_to_counts();
  }
  
  targ_table->round_reset();
  ff_param = 1.0;

  first_round = false;
  
  while (!last_round) {
    if (output_running_rounds) {
      output_results(libs, tot_counts, (int)remaining_rounds);
    }
    remaining_rounds--;
    logger.info("\nRe-estimating counts with additional round of EM (%d "
                "remaining)...", remaining_rounds);
    last_round = (remaining_rounds == 0);
    for (size_t l = 0; l < libs.size(); l++) {
      libs[l].map_parser->write_active(last_round);
      libs[l].map_parser->reset_reader();
    }
    tot_counts = threaded_calc_abundances(libs);
    if (library_size) {
      tot_counts = library_size;
    }
    targ_table->round_reset();
  }
  
   logger.info("Writing results to file...");
  output_results(libs, tot_counts);
  logger.info("Done.");
  
  return 0;
}

#ifdef PROTO
inline string base64_encode(const string& to_encode) {
  using namespace boost::archive::iterators;
  typedef base64_from_binary<transform_width<string::const_iterator,6,8> > it_base64_t;
  unsigned int writePaddChars = (3-to_encode.length()%3)%3;
  string base64(it_base64_t(to_encode.begin()), it_base64_t(to_encode.end()));
  base64.append(writePaddChars,'=');
  return base64;
}

int preprocess_main() {
  try {
    fs::create_directories(output_dir);
  } catch (fs::filesystem_error& e) {
      logger.info(e.what());
  }
  
  if (!fs::exists(output_dir)) {
    logger.severe("Cannot create directory %s.", output_dir.c_str());
  }
  
  Librarian libs(1);
  Library& lib = libs[0];
  lib.in_file_name = in_map_file_names;
  lib.out_file_name = "";
  
  lib.map_parser.reset(new MapParser (&lib, false));
  lib.fld.reset(new LengthDistribution(0, 0, 0, 1, 2, 0));
  MarkovModel bias_model(3, 21, 21, 0);
  MismatchTable mismatch_table(0);
  lib.targ_table.reset(new TargetTable(fasta_file_name, "", 0, 0, 0.0, NULL,
                                       &libs));
  
  logger.info("Converting targets to Protocol Buffers...");
  fstream targ_out((output_dir + "/targets.pb").c_str(),
                   ios::out | ios::trunc);
  string out_buff;
  proto::Target target_proto;
  for (TargID id = 0; id < lib.targ_table->size(); ++id) {
    target_proto.Clear();
    Target& targ = *lib.targ_table->get_targ(id);
    target_proto.set_name(targ.name());
    target_proto.set_id((unsigned int)targ.id());
    target_proto.set_length((unsigned int)targ.length());
    target_proto.set_seq(targ.seq(0).serialize());
    
    target_proto.SerializeToString(&out_buff);
    targ_out << base64_encode(out_buff) << endl;
  }
  targ_out.close();
  
  logger.info("Converting fragment alignments to Protocol Buffers...");
  ostream frag_out(cout.rdbuf());
  
  size_t num_frags = 0;
  Fragment* frag;
  
  ParseThreadSafety pts(10);
  boost::thread parse(&MapParser::threaded_parse, lib.map_parser.get(), &pts,
                      stop_at, 0);
  RobertsFilter frags_seen;
  proto::Fragment frag_proto;
  while(true) {
    frag_proto.Clear();
    
    // Pop next parsed fragment and set mass
    frag = pts.proc_in.pop();
    
    if (!frag) {
      break;
    }
    
    // Test that we have not already seen this fragment
    if (frags_seen.test_and_push(frag->name())) {
      logger.severe("Alignments are not properly sorted. Read '%s' has "
                    "alignments which are non-consecutive.",
                    frag->name().c_str());
    }
    
    frag_proto.set_paired(frag->paired());
    
    vector<char> ref_left_mm_indices;
    vector<char> ref_left_mm_seq;
    vector<char> ref_left_mm_ref;
    vector<char> ref_right_mm_indices;
    vector<char> ref_right_mm_seq;
    vector<char> ref_right_mm_ref;
    bool ref_left_first = false;
    
    for (size_t i = 0; i < frag->num_hits(); ++i) {
      FragHit& fh = *(*frag)[i];
      proto::FragmentAlignment& align_proto = *frag_proto.add_alignments();
      align_proto.set_target_id((unsigned int)fh.target_id());
      
      vector<char> left_mm_indices;
      vector<char> left_mm_seq;
      vector<char> left_mm_ref;
      vector<char> right_mm_indices;
      vector<char> right_mm_seq;
      vector<char> right_mm_ref;
      
      mismatch_table.get_indices(fh, left_mm_indices, left_mm_seq, left_mm_ref,
                                 right_mm_indices, right_mm_seq, right_mm_ref);
      
      if (i == 0) {
        ref_left_mm_indices = left_mm_indices;
        ref_left_mm_seq = left_mm_seq;
        ref_left_mm_ref = left_mm_ref;
        ref_right_mm_indices = right_mm_indices;
        ref_right_mm_seq = right_mm_seq;
        ref_right_mm_ref = right_mm_ref;
      }
      
      ReadHit* read_l = fh.left_read();
      if (read_l) {
        if (i==0) { ref_left_first = read_l->first; }
        proto::ReadAlignment& read_proto = *align_proto.mutable_read_l();
        read_proto.set_first(read_l->first);
        read_proto.set_left_pos(read_l->left);
        read_proto.set_right_pos(read_l->right-1);
        if (i == 0 || read_l->first != ref_left_first ||
            left_mm_indices != ref_left_mm_indices ||
            left_mm_seq != ref_left_mm_seq ||
            left_mm_ref != ref_left_mm_ref) {
          read_proto.set_mismatch_indices(string(left_mm_indices.begin(),
                                                 left_mm_indices.end()));
          read_proto.set_mismatch_nucs(string(left_mm_seq.begin(),
                                              left_mm_seq.end()));
        }
      }
      
      ReadHit* read_r = fh.right_read();
      if (read_r) {
        if (i==0) { ref_left_first = !read_r->first; }
        proto::ReadAlignment& read_proto = *align_proto.mutable_read_r();
        read_proto.set_first(read_r->first);
        read_proto.set_left_pos(read_r->left);
        read_proto.set_right_pos(read_r->right-1);
        if (i == 0 || read_r->first == ref_left_first ||
            right_mm_indices != ref_right_mm_indices ||
            right_mm_seq != ref_right_mm_seq ||
            right_mm_ref != ref_right_mm_ref) {
          read_proto.set_mismatch_indices(string(right_mm_indices.begin(),
                                                 right_mm_indices.end()));
          read_proto.set_mismatch_nucs(string(right_mm_seq.begin(),
                                              right_mm_seq.end()));
        }
      }
    }
    frag_proto.SerializeToString(&out_buff);
    frag_out << base64_encode(out_buff) << endl;
    
    pts.proc_out.push(frag);
    
    num_frags++;
    
    // Output progress
    if (num_frags % 1000000 == 0) {
      logger.info("Fragments Processed: %d", num_frags);
    }
  }
  
  parse.join();
  
  return 0;
}
#endif


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

  srand((unsigned int)time(NULL));
  int parse_ret = parse_options(argc,argv);
  if (parse_ret) {
    return parse_ret;
  }
  
#ifdef PROTO
  if (spark_pre) {
    return preprocess_main();
  }
#endif
  
  return estimation_main();
}