<< "(" << rd2name[mtp.first] << ") "
<< "(" << rd1.getEID() << ") "
int main(int argc, char **argv){  if (! GetOptions(argc, argv)){    cerr << "Command line parsing failed" << endl;    printHelpText();    exit(1);  }  // set up some global variables  int MIN_OBS = 100;   // min. num. of observations required to change lib size  float NUM_SD = 3.00; // num. of standard deviations within which mate is good  int MIN_GOOD_CVG = 0; // good mate coverages below this will be reported  int MAX_SHORT_CVG = 2; // short mate coverages over this will be reported  int MAX_LONG_CVG = 2; // long mate coverages over this will be reported  int MAX_NORMAL_CVG = 3; // normal orientation read coverages > will be reptd.  int MAX_OUTIE_CVG = 3; // outie orientation read coverages > will be reptd.  int MAX_LINKING_CVG = 3; // linking mate coverage > will be reported  int MAX_SINGLEMATE_CVG = 3; // single mate coverage > will be reported  int MEA_CHANGE = 50;// libraries that change by less are considered unchanged  int SD_CHANGE = 10;// libraries that change by less are considered unchanged  bool recompute = false; // recompute library sizes  bool update = false;    // update library sizes in bank  bool feats = false;     // add features to bank  bool debug = false;     // write debugging info to STDERR  bool perinsert = false; // do per-insert output  if (globals.find("debug") != globals.end()){    debug = true;    cout << "Debugging mode: additional information written to STDERR" << endl;  }  if (globals.find("meanchange") != globals.end())      MEA_CHANGE = strtol(globals["meanchange"].c_str(), NULL, 10);  if (globals.find("sdchange") != globals.end())      SD_CHANGE = strtol(globals["sdchange"].c_str(), NULL, 10);  if (globals.find("recompute") != globals.end()){    recompute = true;    cout << "Will recompute library sizes" << endl;    cout << "Libraries whose mean changes by less than " << MEA_CHANGE	 << "\nand whose std. dev. changes by less than " << SD_CHANGE 	 << " are considered to not have changed" << endl;  } else    cout << "Will not recompute library sizes" << endl;  if (globals.find("update") != globals.end()){    update = true;    cout << "Will update libraries in bank" << endl;  } else    cout << "Will not update libraries in bank" << endl;  if (globals.find("feat") != globals.end()){    feats = true;    cout << "Will write features in bank" << endl;  } else     cout << "Will not write features in bank" << endl;  if (recompute){    if (globals.find("minobs") != globals.end())      MIN_OBS = strtol(globals["minobs"].c_str(), NULL, 10);    cout << "Minimum number of observations required to change library size: " 	 << MIN_OBS << endl;  }  if (globals.find("numsd") != globals.end())    NUM_SD = strtof(globals["numsd"].c_str(), NULL);  cout << "Mates within " << NUM_SD        << " standard deviations from mean considered good" << endl;  if (globals.find("goodcvg") != globals.end())    MIN_GOOD_CVG = strtol(globals["goodcvg"].c_str(), NULL, 10);  cout << "Reporting regions with less than " << MIN_GOOD_CVG       << " good mate coverage" << endl;  if (globals.find("shortcvg") != globals.end())    MAX_SHORT_CVG = strtol(globals["shortcvg"].c_str(), NULL, 10);  cout << "Reporting regions with more than " << MAX_SHORT_CVG       << " short mate coverage" << endl;  if (globals.find("longcvg") != globals.end())    MAX_LONG_CVG = strtol(globals["longcvg"].c_str(), NULL, 10);  cout << "Reporting regions with more than " << MAX_LONG_CVG       << " long mate coverage" << endl;  if (globals.find("normalcvg") != globals.end())    MAX_NORMAL_CVG = strtol(globals["normalcvg"].c_str(), NULL, 10);  cout << "Reporting regions with more than " << MAX_NORMAL_CVG        << " coverage by mates oriented the same way" << endl;  if (globals.find("outiecvg") != globals.end())    MAX_OUTIE_CVG = strtol(globals["outiecvg"].c_str(), NULL, 10);  cout << "Reporting regions with more than " << MAX_OUTIE_CVG       << " coverage by outie mates" << endl;  if (globals.find("outiecvg") != globals.end())    MAX_OUTIE_CVG = strtol(globals["outiecvg"].c_str(), NULL, 10);  cout << "Reporting regions with more than " << MAX_OUTIE_CVG       << " coverage by outie mates" << endl;  if (globals.find("linking") != globals.end())    MAX_LINKING_CVG = strtol(globals["linking"].c_str(), NULL, 10);  cout << "Reporting regions with more than " << MAX_LINKING_CVG       << " coverage by linking mates" << endl;  if (globals.find("singlemate") != globals.end())    MAX_SINGLEMATE_CVG = strtol(globals["singlemate"].c_str(), NULL, 10);  cout << "Reporting regions with more than " << MAX_SINGLEMATE_CVG       << " coverage by singleton mates" << endl;  cout << endl;  if (globals.find("perinsert") != globals.end()){    perinsert = true;    cout << "Per-insert information reported in file: " 	 << globals["perinsert"] << endl;  }  cout << endl;  // open necessary files  if (globals.find("bank") == globals.end()){ // no bank was specified    cerr << "A bank must be specified" << endl;    exit(1);  }  // prepare banks we'll need    Bank_t library_bank (Library_t::NCODE);  if (! library_bank.exists(globals["bank"])){    cerr << "No library account found in bank " << globals["bank"] << endl;    exit(1);  }  try {    if (update)      library_bank.open(globals["bank"], B_READ|B_WRITE);    else      library_bank.open(globals["bank"], B_READ);  } catch (Exception_t & e)    {      cerr << "Failed to open library account in bank " << globals["bank"]            << ": " << endl << e << endl;      exit(1);    }  BankStream_t contig_stream (Contig_t::NCODE);  if (! contig_stream.exists(globals["bank"])){    cerr << "No contig account found in bank " << globals["bank"] << endl;    exit(1);  }  try {    contig_stream.open(globals["bank"], B_READ);  } catch (Exception_t & e)    {      cerr << "Failed to open contig account in bank " << globals["bank"]            << ": " << endl << e << endl;      exit(1);    }  Bank_t read_bank (Read_t::NCODE);  if (! read_bank.exists(globals["bank"])){    cerr << "No read account found in bank " << globals["bank"] << endl;    exit(1);  }  try {    read_bank.open(globals["bank"], B_READ);  } catch (Exception_t & e)    {      cerr << "Failed to open read account in bank " << globals["bank"]            << ": " << endl << e << endl;      exit(1);    }  BankStream_t frag_bank (Fragment_t::NCODE);  if (! frag_bank.exists(globals["bank"])){    cerr << "No frag account found in bank " << globals["bank"] << endl;    exit(1);  }  try {    frag_bank.open(globals["bank"], B_READ);  } catch (Exception_t & e)    {      cerr << "Failed to open frag account in bank " << globals["bank"]            << ": " << endl << e << endl;      exit(1);    }  BankStream_t feat_stream (Feature_t::NCODE);  if (feats)    {      try {        if ( feat_stream.exists(globals["bank"]) ){          feat_stream.open(globals["bank"], B_READ|B_WRITE);        } else {          feat_stream.create(globals["bank"], B_READ|B_WRITE);        }      } catch (Exception_t & e)        {          cerr << "Failed to open feat account in bank " << globals["bank"]                << ": " << endl << e << endl;          exit(1);        }    }  ofstream insertFile;  if (perinsert){    insertFile.open(globals["perinsert"].c_str(), ofstream::out);    if (! insertFile.is_open()){      cerr << "Failed to open insert file " << globals["perinsert"] << endl;      exit(1);    }  }  // Stream through the contigs and retrieve map of read positions within ctgs.  Contig_t ctg;    hash_map<ID_t, ID_t, hash<ID_t>, equal_to<ID_t> > rd2ctg;     // map from read to contig  hash_map<ID_t, Range_t, hash<ID_t>, equal_to<ID_t> > rd2posn; // position in contig  hash_map<ID_t, Size_t, hash<ID_t>, equal_to<ID_t> > ctglen;   // length of contig  hash_map<ID_t, string, hash<ID_t>, equal_to<ID_t> > ctgname;   // name of contig  set<ID_t> ctgIDs;  while (contig_stream >> ctg)    for (vector<Tile_t>::iterator ti = ctg.getReadTiling().begin();          ti != ctg.getReadTiling().end(); ti++){      rd2ctg[ti->source] = ctg.getIID();      ctgIDs.insert(ctg.getIID());      Pos_t f, l; // coords of beginning/end of read      if (ti->range.getEnd() < ti->range.getBegin()) { // reverse	f = ti->offset + ti->range.getBegin() - ti->range.getEnd() + 	  ti->gaps.size();	l = ti->offset;      } else {	f = ti->offset;	l = ti->offset + ti->range.getEnd() - ti->range.getBegin() + 	  ti->gaps.size();      }      rd2posn[ti->source] =         Range_t(f, l);      ctglen[ctg.getIID()] = ctg.getLength();      ctgname[ctg.getIID()] = ctg.getEID();    }  contig_stream.close();  AnnotatedFragment frag;  pair<ID_t, ID_t> mtp;  list<AnnotatedFragment> mtl;  Read_t rd1, rd2;  set<ID_t> libIDs;  hash_map<ID_t, ID_t, hash<ID_t>, equal_to<ID_t> > rd2lib;  hash_map<ID_t, ID_t, hash<ID_t>, equal_to<ID_t> > rd2frg;  hash_map<ID_t, list<list<AnnotatedFragment>::iterator>, hash<ID_t>, equal_to<ID_t> > lib2frag;  hash_map<ID_t, list<list<AnnotatedFragment>::iterator>, hash<ID_t>, equal_to<ID_t> > ctg2frag;  hash_map<ID_t, string, hash<ID_t>, equal_to<ID_t> > rd2name;  while (frag_bank >> frag){    frag.status = MP_GOOD;    mtp = frag.getMatePair();    read_bank.fetch(mtp.first, rd1); // get the read    read_bank.fetch(mtp.second, rd2);    rd2name[mtp.first] = rd1.getEID();    rd2name[mtp.second] = rd2.getEID();    rd2frg[mtp.first] = frag.getIID();    rd2frg[mtp.second] = frag.getIID();    rd2lib[mtp.first] = frag.getLibrary();    rd2lib[mtp.second] = frag.getLibrary();    libIDs.insert(frag.getLibrary());    if (frag.getType() != Fragment_t::INSERT) {// we only handle end reads      //      cerr << "Type is " << frag.getType() << endl;      if (perinsert)	insertFile << "NOT END MATES: "		   << mtp.first 		   << "(" << rd2name[mtp.first] << ") " 		   << mtp.second 		   << "(" << rd2name[mtp.second] << ") "		   << "\n";      continue;    }    list<AnnotatedFragment>::iterator ti;     if (rd2ctg[mtp.first] == rd2ctg[mtp.second]){      // mate-pair between reads within the same contig      //      ctgIDs.insert(rd2ctg[mtp.first]);      ti = mtl.insert(mtl.end(), frag);      lib2frag[frag.getLibrary()].push_back(ti);      ctg2frag[rd2ctg[mtp.first]].push_back(ti);    } else {      if (rd2ctg[mtp.first] == 0 ||	  rd2ctg[mtp.second] == 0) { // mate should be in contig	frag.status = MP_SINGMATE;	ti = mtl.insert(mtl.end(), frag);	lib2frag[frag.getLibrary()].push_back(ti);	ctg2frag[rd2ctg[mtp.first]].push_back(ti);	if (perinsert)	  insertFile << "SINGLETON_MATE: " 		     << mtp.first 		     << "(" << rd2name[mtp.first] << ") " 		     << mtp.second 		     << "(" << rd2name[mtp.second] << ") "		     << " in separate contigs: "		     << rd2ctg[mtp.first] << ", "		     << rd2ctg[mtp.second] 		     << "\n";      } else {	frag.status = MP_LINKING;	ti = mtl.insert(mtl.end(), frag);	lib2frag[frag.getLibrary()].push_back(ti);	ctg2frag[rd2ctg[mtp.first]].push_back(ti);		if (perinsert)	  insertFile << "LINKING: " 		     << mtp.first 		     << "(" << rd2name[mtp.first] << ") " 		     << mtp.second 		     << "(" << rd2name[mtp.second] << ") "		     << " in separate contigs: "		     << rd2ctg[mtp.first] << ", "		     << rd2ctg[mtp.second] 		     << "\n";      }    }  }    // now we get the library information for each library  hash_map<ID_t, pair<Pos_t, SD_t>, hash<ID_t>, equal_to<ID_t> > lib2size;  hash_map<ID_t, string, hash<ID_t>, equal_to<ID_t> > lib2name;  Library_t lib;  for (set<ID_t>::iterator li = libIDs.begin(); li != libIDs.end(); li++){    library_bank.fetch(*li, lib);    lib2name[*li] = lib.getEID();    lib2size[*li] = pair<Pos_t, SD_t> (lib.getDistribution().mean, 				       lib.getDistribution().sd);  }    // now we're ready to figure out new library sizes  hash_map<ID_t, pair<Pos_t, SD_t>, hash<ID_t>, equal_to<ID_t> > NewLib2size;    for (set<ID_t>::iterator li = libIDs.begin(); li != libIDs.end(); li++){    // for each library    list<Pos_t> sizes;    pair<Pos_t, SD_t> libsize = lib2size[*li];        for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = lib2frag[*li].begin(); mi != lib2frag[*li].end(); mi++){      if ((*mi)->status != MP_GOOD)	continue;      // for each mate pair      Pos_t sz = mateLen(**mi, rd2posn[(*mi)->getMatePair().first], 			 rd2posn[(*mi)->getMatePair().second], 			 ctglen[rd2ctg[(*mi)->getMatePair().first]], 			 libsize.first + 3 * libsize.second);      if (recompute && sz != 0) {	sizes.push_back(sz);      }    } // if recomputing library sizes    // compute new sizes    pair<Pos_t, SD_t> newSize;    if (! recompute || sizes.size() < MIN_OBS ){      newSize = libsize;    } else {      newSize = getSz(sizes);      if ((abs (((int) newSize.first) - ((int) libsize.first)) > MEA_CHANGE) ||	  (abs (((int) newSize.second) - ((int) libsize.second)) > SD_CHANGE) ){ // if new size significantly different	cout <<  "Library " << lib2name[*li] << " recomputed as mean = " 	     << newSize.first << " SD = " << newSize.second << endl;	if (update){	  Library_t newlib;	  Distribution_t dist;	  dist.mean = newSize.first;	  dist.sd = newSize.second;	  newlib.setIID(*li);	  newlib.setEID(lib2name[*li]);	  newlib.setDistribution(dist);	  try {	    library_bank.replace(*li, newlib);	  }  	  catch (Exception_t & e)	    {	      cerr << "Failed to replace library " << *li << "(" 		   << lib2name[*li] << ") in bank " << globals["bank"] 		   << ": " << endl << e << endl;	      exit(1);	    }	} // if updating library in bank      } else {	cout << "Size of library " << lib2name[*li] 	     << " has not changed enough: MEA= " << newSize.first 	     << " SD= " << newSize.second << " versus MEA= "	     << libsize.first << " SD= " << libsize.second << endl;      }// if size actually changed    }     NewLib2size[*li] = newSize;    // recheck mate lengths and set status    for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = lib2frag[*li].begin(); mi != lib2frag[*li].end(); mi++){      // for each mate pair      // turn off end check       if ((*mi)->status != MP_GOOD){	if (perinsert){	  if ((*mi)->status == MP_OUTIE)	    insertFile << "OUTIE: ";	  if ((*mi)->status == MP_NORMAL)	    insertFile << "SAME: ";	  	  insertFile << ((*mi)->getMatePair()).first 		     << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 		     << ((*mi)->getMatePair()).second 		     << "(" << rd2name[((*mi)->getMatePair()).second] 		     << ") LIBRARY: "		     << *li << "(" << lib2name[*li] << ")"		     << "\n";	}	continue; // only interested in mates not already tagged      }      Pos_t sz = mateLen(**mi, rd2posn[((*mi)->getMatePair()).first], 			 rd2posn[((*mi)->getMatePair()).second], 			 ctglen[rd2ctg[((*mi)->getMatePair()).first]], 0);      if (debug && sz != 0)	cerr << "MATE " << *li << " " << sz << endl;            // sz == 0 implies mis-orientation      if (abs(sz - newSize.first) > newSize.second * NUM_SD) {	// farther than NUM_SD standard deviations from the mean	if (sz < newSize.first){	  (*mi)->status = MP_SHORT; // too short	  (*mi)->deviation = newSize.first - sz;	  if (perinsert)	    insertFile << "SHORT " 		       << ((*mi)->getMatePair()).first 		       << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 		       << ((*mi)->getMatePair()).second 		       << "(" << rd2name[((*mi)->getMatePair()).second] << ") "		       << "SIZE: " << sz << " EXPECTED: " << newSize.first 		       << "(" << newSize.second << ") LIBRARY: " 		       << *li << "(" << lib2name[*li] << ")"		       << "\n";			} else { 	  (*mi)->status = MP_LONG; // too long	  (*mi)->deviation = sz - newSize.first;	  if (perinsert)	    insertFile << "LONG " 		       << ((*mi)->getMatePair()).first 		       << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 		       << ((*mi)->getMatePair()).second 		       << "(" << rd2name[((*mi)->getMatePair()).second] << ") "		       << "SIZE: " << sz << " EXPECTED: " << newSize.first 		       << "(" << newSize.second << ") LIBRARY: "		       << *li << "(" << lib2name[*li] << ")"		       << "\n";	}      } else // long or short	if (perinsert)	  insertFile << "GOOD " 		     << ((*mi)->getMatePair()).first 		     << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 		     << ((*mi)->getMatePair()).second 		     << "(" << rd2name[((*mi)->getMatePair()).second] << ") "		     << "SIZE: " << sz << " EXPECTED: " << newSize.first 		     << "(" << newSize.second << ") LIBRARY: "		     << *li << "(" << lib2name[*li] << ")"		     << "\n";    } // for each mate pair  }// for each library  if (recompute)     cout << endl;      for (set<ID_t>::iterator ctg = ctgIDs.begin(); ctg != ctgIDs.end(); ctg++) {    list<pair<Pos_t, Pos_t> > ranges;  // ranges we are interested in    list<pair<Pos_t, Pos_t> > interest;    cout << ">Contig_" << *ctg << " " << ctgname[*ctg] << " " 	 << ctglen[*ctg] << " bases" << endl;    // find coverage by good mates    for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = ctg2frag[*ctg].begin(); mi != ctg2frag[*ctg].end(); mi++){      if ((*mi)->status == MP_GOOD){	Pos_t left, right;	left = rd2posn[((*mi)->getMatePair()).first].getEnd(); // use end of reads	right = rd2posn[((*mi)->getMatePair()).second].getEnd(); 	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));	if (debug)	  cerr << "Good mate " 	       << ((*mi)->getMatePair()).first 	       << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 	       << ((*mi)->getMatePair()).second 	       << "(" << rd2name[((*mi)->getMatePair()).second] << ") "	       << endl;      }    } // for each mate in the contig    getCvg(ranges, interest, ctglen[*ctg], MIN_GOOD_CVG, false);     // report interesting ranges    if (interest.size() > 0){      cout << "Regions of " << MIN_GOOD_CVG << "X" 	   << ((MIN_GOOD_CVG != 0) ? " or less": "") << " clone coverage" 	   <<endl;            for (list<pair<Pos_t, Pos_t> >::iterator ii = interest.begin(); 	   ii != interest.end(); ii++){	cout << ii->first << "\t" << ii->second << endl;	if (feats){	  Feature_t f;	  f.setComment("LOW_GOOD_CVG");	  f.setRange(Range_t(ii->first, ii->second));          f.setSource(make_pair(*ctg, Contig_t::NCODE));	  f.setType(Feature_t::COVERAGE);          feat_stream << f;	}      }      cout << endl;    }    ranges.clear();    // find coverage by short mates    for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = ctg2frag[*ctg].begin(); mi != ctg2frag[*ctg].end(); mi++){      if ((*mi)->status == MP_SHORT){	Pos_t left, right;	// use beginnings of reads	left = rd2posn[((*mi)->getMatePair()).first].getBegin(); 	right = rd2posn[((*mi)->getMatePair()).second].getBegin(); 	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));      }    } // for each mate in the contig    getCvg(ranges, interest, ctglen[*ctg], MAX_SHORT_CVG, true);     // report interesting ranges    if (interest.size() > 0){      cout << "Regions of " << MAX_SHORT_CVG << "X or more short mate coverage"	   <<endl;            for (list<pair<Pos_t, Pos_t> >::iterator ii = interest.begin(); 	   ii != interest.end(); ii++){	cout << ii->first << "\t" << ii->second << endl;	if (feats){	  Feature_t f;	  f.setComment("HIGH_SHORT_CVG");	  f.setRange(Range_t(ii->first, ii->second));          f.setSource(make_pair(*ctg, Contig_t::NCODE));	  f.setType(Feature_t::COVERAGE);          feat_stream << f;	}      }      cout << endl;    }    ranges.clear();    // find coverage by long mates    for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = ctg2frag[*ctg].begin(); mi != ctg2frag[*ctg].end(); mi++){      if ((*mi)->status == MP_LONG){	Pos_t left, right;	// use beginnings of reads	left = rd2posn[((*mi)->getMatePair()).first].getBegin(); 	right = rd2posn[((*mi)->getMatePair()).second].getBegin(); 	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));      }    } // for each mate in the contig    getCvg(ranges, interest, ctglen[*ctg], MAX_LONG_CVG, true);     // report interesting ranges    if (interest.size() > 0){      cout << "Regions of " << MAX_LONG_CVG << "X or more long mate coverage" 	   <<endl;            for (list<pair<Pos_t, Pos_t> >::iterator ii = interest.begin(); 	   ii != interest.end(); ii++){	cout << ii->first << "\t" << ii->second << endl;	if (feats){	  Feature_t f;	  f.setComment("HIGH_LONG_CVG");	  f.setRange(Range_t(ii->first, ii->second));          f.setSource(make_pair(*ctg, Contig_t::NCODE));	  f.setType(Feature_t::COVERAGE);          feat_stream << f;	}      }      cout << endl;    }        ranges.clear();    // find coverage by normal reads    for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = ctg2frag[*ctg].begin(); mi != ctg2frag[*ctg].end(); mi++){      if ((*mi)->status == MP_NORMAL){	Pos_t left, right;	left = rd2posn[((*mi)->getMatePair()).first].getBegin(); 	right = rd2posn[((*mi)->getMatePair()).first].getEnd(); 	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));	left = rd2posn[((*mi)->getMatePair()).second].getBegin(); 	right = rd2posn[((*mi)->getMatePair()).second].getEnd(); 	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));	if (debug)	  cerr << "Same orientation mates " 	       << ((*mi)->getMatePair()).first 	       << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 	       << ((*mi)->getMatePair()).second 	       << "(" << rd2name[((*mi)->getMatePair()).second] << ") "	       << endl;      }    } // for each mate in the contig    getCvg(ranges, interest, ctglen[*ctg], MAX_NORMAL_CVG, true);     // report interesting ranges    if (interest.size() > 0){      cout << "Regions of " << MAX_NORMAL_CVG 	   << "X or more same-orientation coverage" <<endl;            for (list<pair<Pos_t, Pos_t> >::iterator ii = interest.begin(); 	   ii != interest.end(); ii++){	cout << ii->first << "\t" << ii->second << endl;	if (feats){	  Feature_t f;	  f.setComment("HIGH_NORMAL_CVG");	  f.setRange(Range_t(ii->first, ii->second));          f.setSource(make_pair(*ctg, Contig_t::NCODE));	  f.setType(Feature_t::COVERAGE);          feat_stream << f;	}      }      cout << endl;    }    ranges.clear();    // find coverage by outie reads    for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = ctg2frag[*ctg].begin(); mi != ctg2frag[*ctg].end(); mi++){      if ((*mi)->status == MP_OUTIE){	Pos_t left, right;	//	if ((*mi)->getType() == Fragment_t::TRANSPOSON)	//	  continue;  // transposons are not errors	left = rd2posn[((*mi)->getMatePair()).first].getBegin(); 	right = rd2posn[((*mi)->getMatePair()).first].getEnd(); 	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));	left = rd2posn[((*mi)->getMatePair()).second].getBegin(); 	right = rd2posn[((*mi)->getMatePair()).second].getEnd(); 	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));	if (debug)	  cerr << "Outie mates " 	       << ((*mi)->getMatePair()).first 	       << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 	       << ((*mi)->getMatePair()).second 	       << "(" << rd2name[((*mi)->getMatePair()).second] << ") "	       << endl;      }    } // for each mate in the contig    getCvg(ranges, interest, ctglen[*ctg], MAX_OUTIE_CVG, true);     // report interesting ranges    if (interest.size() > 0){      cout << "Regions of " << MAX_OUTIE_CVG << "X or more outie coverage" 	   <<endl;            for (list<pair<Pos_t, Pos_t> >::iterator ii = interest.begin(); 	   ii != interest.end(); ii++){	cout << ii->first << "\t" << ii->second << endl;	if (feats){	  Feature_t f;	  f.setComment("HIGH_OUTIE_CVG");	  f.setRange(Range_t(ii->first, ii->second));          f.setSource(make_pair(*ctg, Contig_t::NCODE));	  f.setType(Feature_t::COVERAGE);          feat_stream << f;	}      }      cout << endl;    }    ranges.clear();    // find coverage by linking reads    for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = ctg2frag[*ctg].begin(); mi != ctg2frag[*ctg].end(); mi++){      if ((*mi)->status == MP_LINKING){	Pos_t left, right;	if (rd2ctg[((*mi)->getMatePair()).first] == *ctg) {	  left = rd2posn[((*mi)->getMatePair()).first].getBegin(); 	  right = rd2posn[((*mi)->getMatePair()).first].getEnd(); 	} else {	  left = rd2posn[((*mi)->getMatePair()).second].getBegin(); 	  right = rd2posn[((*mi)->getMatePair()).second].getEnd(); 	}	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));	if (debug)	  cerr << "Linking mates " 	       << ((*mi)->getMatePair()).first 	       << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 	       << ((*mi)->getMatePair()).second 	       << "(" << rd2name[((*mi)->getMatePair()).second] << ") "	       << endl;      }    } // for each mate in the contig    getCvg(ranges, interest, ctglen[*ctg], MAX_LINKING_CVG, true);     // report interesting ranges    if (interest.size() > 0){      cout << "Regions of " << MAX_LINKING_CVG << "X or more linking coverage" 	   <<endl;            for (list<pair<Pos_t, Pos_t> >::iterator ii = interest.begin(); 	   ii != interest.end(); ii++){	cout << ii->first << "\t" << ii->second << endl;	if (feats){	  Feature_t f;	  f.setComment("HIGH_LINKING_CVG");	  f.setRange(Range_t(ii->first, ii->second));          f.setSource(make_pair(*ctg, Contig_t::NCODE));	  f.setType(Feature_t::COVERAGE);          feat_stream << f;	}      }      cout << endl;    }    ranges.clear();    // find coverage by single mate reads    for (list<list<AnnotatedFragment>::iterator>::iterator 	   mi = ctg2frag[*ctg].begin(); mi != ctg2frag[*ctg].end(); mi++){      if ((*mi)->status == MP_SINGMATE){	Pos_t left, right;	if (rd2ctg[((*mi)->getMatePair()).first] == *ctg) {	  left = rd2posn[((*mi)->getMatePair()).first].getBegin(); 	  right = rd2posn[((*mi)->getMatePair()).first].getEnd(); 	} else {	  left = rd2posn[((*mi)->getMatePair()).second].getBegin(); 	  right = rd2posn[((*mi)->getMatePair()).second].getEnd(); 	}	if (left < right)	  ranges.push_back(pair<Pos_t, Pos_t>(left, right));	else	  ranges.push_back(pair<Pos_t, Pos_t>(right, left));	if (debug)	  cerr << "Singleton mates " 	       << ((*mi)->getMatePair()).first 	       << "(" << rd2name[((*mi)->getMatePair()).first] << ") " 	       << ((*mi)->getMatePair()).second 	       << "(" << rd2name[((*mi)->getMatePair()).second] << ") "	       << endl;      }    } // for each mate in the contig    getCvg(ranges, interest, ctglen[*ctg], MAX_SINGLEMATE_CVG, true);     // report interesting ranges    if (interest.size() > 0){      cout << "Regions of " << MAX_SINGLEMATE_CVG << "X or more singleton mate coverage" 	   <<endl;            for (list<pair<Pos_t, Pos_t> >::iterator ii = interest.begin(); 	   ii != interest.end(); ii++){	cout << ii->first << "\t" << ii->second << endl;	if (feats){	  Feature_t f;	  f.setComment("HIGH_SINGLEMATE_CVG");	  f.setRange(Range_t(ii->first, ii->second));          f.setSource(make_pair(*ctg, Contig_t::NCODE));	  f.setType(Feature_t::COVERAGE);          feat_stream << f;	}      }      cout << endl;    }    ranges.clear();  } // for each contig  frag_bank.close();  read_bank.close();  library_bank.close();  feat_stream.close();  if (perinsert)    insertFile.close();  return(0);} // main