Initial Commit

#util.h#

+#ifndef UTIL_H
+#define UTIL_H
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cstdlib>
+#include <sys/stat.h>
+#include <sys/types.h>
+#include <dirent.h>
+#include <errno.h>
+#include <sstream>
+#include <time.h>
+
+using namespace std;
+
+class util
+{
+ private:
+  
+ public:
+  util();
+  vector<string> load_names(string dir);
+  void extract_paths(string full_path, string& p, string& b);
+  int rm_mkdir(string path, string dirname);
+  int rm(string path, string dirname);
+  int rm(string full_path);
+  void open_ifile(ifstream& instream, string infile, string var_name);
+  void open_ofile(ofstream& outstream, string outfile, string var_name);
+  void open_ofile_ne(ofstream& outstream, string outfile, string var_name);
+  void open_ofile_app(ofstream& outstream, string outfile, string var_name);
+  void get_tl(string infile,int& tl);
+  void get_rss(string infile,int& rss);
+  string int_to_string(int i);
+  void write_out_time(ofstream& logstream);
+};
+
+util::util()
+{
+}
+
+void util::extract_paths(string full_path, string& p, string& b)
+{
+  if(full_path[full_path.length()-1]=='/')
+    {
+      full_path = full_path.substr(0,full_path.length()-1);
+    }
+
+  p=b="";
+  bool slash_seen=false;
+  for(int i=full_path.length()-1; i>=0; i--)
+    {
+      if(!slash_seen && full_path[i]!='/')
+	{
+	  b=full_path[i]+b;
+	}
+      else if(full_path[i]=='/' || slash_seen) 
+	{
+	  slash_seen=true;
+	  p=full_path[i]+p;
+	}
+    }
+}
+
+int util::rm_mkdir(string path, string dirname)
+{
+  string full_path = path+"/"+dirname;
+  string outcommand = "rm -r "+full_path;
+  if(system(NULL)) puts ("Ok");
+  else
+    {
+      exit(EXIT_FAILURE);
+    }
+  int system_return = system(outcommand.c_str());
+  //cout<<"rm -r "<<path<<" - system_return="<<system_return<<endl;
+  outcommand = "mkdir "+full_path;
+  system_return = system(outcommand.c_str());
+  //cout<<"mkdir "<<path<<" - system_return="<<system_return<<endl;  
+  return system_return;
+}
+
+
+int util::rm(string path, string dirname)
+{
+  string full_path = path+"/"+dirname;
+  string outcommand = "rm -r "+full_path;
+  if(system(NULL)) puts ("Ok");
+  else
+    {
+      exit(EXIT_FAILURE);
+    }
+  int system_return = system(outcommand.c_str());
+  return system_return;
+}
+
+int util::rm(string full_path)
+{
+  string outcommand = "rm "+full_path;
+  if(system(NULL)) puts ("Ok");
+  else
+    {
+      exit(EXIT_FAILURE);
+    }
+  int system_return = system(outcommand.c_str());
+  return system_return;
+}
+
+vector<string> util::load_names(string dir)
+{
+  vector<string> folders;
+
+  DIR *dp;
+  struct dirent *dirp;
+  if((dp = opendir(dir.c_str())) == NULL) 
+    {
+      cout<<"Failed to open path to dir="<<dir<<endl;
+      return folders;
+    }
+
+  while ((dirp = readdir(dp)) != NULL) 
+    {
+      string temp = string(dirp->d_name);
+      if(temp.length()>5)
+	{
+	  folders.push_back(string(dirp->d_name));
+	}
+      else
+	{
+	  //cout<<"Rejecting directory named "<<temp<<endl;
+	}
+    }
+  closedir(dp);
+
+  return folders;
+}
+
+void util::open_ifile(ifstream& instream, string infile, string var_name)
+{
+  instream.open(infile.c_str());
+  if(instream.fail())
+    {
+      cout<<"Failed to open ifstream to "<<var_name<<" = "<<infile<<endl;
+      exit(1);
+    }
+}
+
+void util::open_ofile(ofstream& outstream, string outfile, string var_name)
+{
+  outstream.open(outfile.c_str());
+  if(outstream.fail())
+    {
+      cout<<"Failed to open ofstream to "<<var_name<<" = "<<outfile<<endl;
+      exit(1);
+    }
+}
+
+void util::open_ofile_ne(ofstream& outstream, string outfile, string var_name)
+{
+  outstream.open(outfile.c_str());
+  while(outstream.fail())
+    {
+      cout<<"util::open_ofile_ne(ofstream&,string,string)"<<endl;
+      cout<<"Failed to open stream to outfile = "<<outfile<<endl;
+      cout<<"Does the directory exist? Please create it otherwise and then enter any int to continue: ";
+      int wfd;
+      cin>>wfd;
+      outstream.open(outfile.c_str());
+    }
+  cout<<"util::open_ofile_ne(ofstream&,string,string): Successfully opened ofstream to "<<outfile<<endl;
+}
+
+void util::open_ofile_app(ofstream& outstream, string outfile, string var_name)
+{
+  outstream.open(outfile.c_str(), ios::out | ios::app);
+  if(outstream.fail())
+    {
+      cout<<"Failed to open ofstream to "<<var_name<<" = "<<outfile<<endl;
+      exit(1);
+    }
+}
+
+void util::get_tl(string infile,int& tl)
+{
+  ifstream instream;
+  open_ifile(instream,infile,"infile");
+  int i;
+  instream>>i>>i>>tl;
+  instream.close();
+}
+
+void util::get_rss(string infile,int& rss)
+{
+  ifstream instream;
+  open_ifile(instream,infile,"infile");
+  instream>>rss;  
+  instream.close();
+}
+
+string util::int_to_string(int i)
+{
+  stringstream ss;
+  ss<<i;
+  string result;
+  ss>>result;
+  return result;
+}
+
+void util::write_out_time(ofstream& logstream)
+{
+  time_t current_time;
+  current_time=time(NULL);
+  struct tm* timeinfo;
+  timeinfo=localtime(&current_time);
+  logstream<<"Date of run: "<<asctime(timeinfo)<<endl;
+}
+
+#endif /*UTIL_H*/

Coptim1_3Dmat.h

+#include <iostream>
+#include <cstdlib>
+#include <vector>
+#include "Coptim1_layer.h"
+
+using namespace std;
+
+class Coptim1_3Dmat
+{
+ private:
+  vector<Coptim1_layer> mat;
+  int width;
+  int * hit_all;
+ public:
+  Coptim1_3Dmat(int number_of_genomes, int number_of_stacks);
+  void increment_element(int x, int y, int z);
+  int get_element(int x, int y, int z);
+  void increment_allhits(int index);
+  int get_allhits(int index);   
+};
+
+Coptim1_3Dmat::Coptim1_3Dmat(int number_of_genomes, int number_of_stacks)
+{
+  width = number_of_genomes;
+  hit_all = new int[number_of_stacks];
+
+  for(int i=0; i<number_of_stacks; i++)
+    {
+      Coptim1_layer temp(width);
+      mat.push_back(temp);
+      hit_all[i]=0;
+    }
+  cout<<"In Coptim1_3Dmat::Coptim1_3Dmat(int,int). Correlation array initialized to 0."<<endl;
+  cout<<"Size of 3Dmat (number of stacks) = "<<mat.size()<<endl;
+}
+
+void Coptim1_3Dmat::increment_element(int x, int y, int z)
+{
+  mat.at(z).increment_element(x,y);
+}
+
+int Coptim1_3Dmat::get_element(int x, int y, int z)
+{
+  return mat.at(z).get_element(x,y);
+}
+
+void Coptim1_3Dmat::increment_allhits(int index)
+{
+  if(index>=0 && index<mat.size())
+  {
+    hit_all[index]++;
+  }
+  else
+    {
+      cout<<"Invalid index in Coptim1_3Dmat::increment_allhits(int). "<<index<<endl;
+      exit(1);
+    }
+}
+
+int Coptim1_3Dmat::get_allhits(int index)
+{
+  if(index>=0 && index<mat.size())
+    {
+      return hit_all[index];
+    }
+  else
+    {
+      cout<<"Invalid index in Coptim1_3Dmat::get_allhits(int.): "<<index<<endl;
+      exit(1);
+      return 0;
+    }
+}
+		

Coptim1_layer.h

+#include <iostream>
+#include <cstdlib>
+
+using namespace std;
+
+class Coptim1_layer
+{
+ private:
+  int* array;
+  int* diagonal;
+  int array_size;
+  int diagonal_size;
+
+  int get_index(int x, int y);
+
+ public:
+  Coptim1_layer(int number_of_genomes);
+  void increment_element(int x, int y);
+  int get_element(int x, int y);
+  
+};
+
+Coptim1_layer::Coptim1_layer(int number_of_genomes)
+{
+  //array_size = (((number_of_genomes*number_of_genomes)-number_of_genomes)/2.0);
+  array_size = get_index(number_of_genomes-2,number_of_genomes-1)+1;
+  diagonal_size = number_of_genomes;
+  array = new int[array_size];
+  diagonal = new int[diagonal_size];
+
+  for(int i=0; i<array_size; i++)
+    {
+      array[i]=0;
+    }
+  for(int i=0; i<diagonal_size; i++)
+    {
+      diagonal[i]=0;
+    }
+}
+
+void Coptim1_layer::increment_element(int x, int y)
+{
+  if(x<0 || x>diagonal_size || y<0 || y>diagonal_size)
+    {
+      cout<<"x or y out of bounds in Coptim1_layer::increment_layer(int,int): "<<x<<" "<<y<<endl;
+      exit(1);
+    }
+  else if(x==y)
+    {
+      diagonal[x]++;
+    }
+  else
+    {
+      int index = get_index(x,y);
+      if(index>=array_size || index<0)
+	{
+	  cout<<"Invalid index size: "<<index<<" x,y="<<x<<","<<y<<endl;
+	  exit(1);
+	}
+      array[index]++;
+    }
+}
+
+int Coptim1_layer::get_index(int x, int y)
+{
+  if(x>y)
+    {
+      return (((x*x)+x)/2.0)+y-x;
+    }
+  else if(y>x)
+    {
+      return (((y*y)+y)/2.0)+x-y;
+    }
+  else
+    {
+      cout<<"Weird values of x and y in Coptim_layer::get_index(int,int): "<<x<<" "<<y<<endl;
+      exit(1);
+      return -1;
+    }
+
+  exit(1); //additional precaution
+  return -1; //to stop the compiler from complaining
+}
+
+int Coptim1_layer::get_element(int x, int y)
+{
+  if(x<0 || x>diagonal_size || y<0 || y>diagonal_size)
+    {
+      cout<<"Out of bounds x or y in Coptim1_layer::get_element(int,int): "<<x<<" "<<y<<endl;
+      exit(1);
+    }
+  else if(x==y)
+    {
+      return diagonal[x];
+    }
+  else
+    {
+      int index = get_index(x,y);
+      if(index<0 || index>=array_size)
+	{
+	  cout<<"Invalid value for index in Coptim1_layer::get_element: "<<index<<" x,y="<<x<<","<<y<<endl;
+	  exit(1);
+	}
+      return array[index];
+    }
+}

amino_acids.h

+#ifndef AMINO_ACIDS_H
+#define AMINO_ACIDS_H
+
+#include <iostream>
+#include <cstdlib>
+
+using namespace std;
+
+class amino_acids
+{
+ private:
+ public:
+  amino_acids();
+  char get_aa(int i);
+  int get_aa(char c);
+};
+
+amino_acids::amino_acids()
+{
+}
+
+char amino_acids::get_aa(int i)
+{
+  if(i==0)
+    {
+      return 'A';
+    }
+  else if(i==1)
+    {
+      return 'C';
+    }
+  else if(i==2)
+    {
+      return 'D';
+    }
+  else if(i==3)
+    {
+      return 'E';
+    }
+  else if(i==4)
+    {
+      return 'F';
+    }
+  else if(i==5)
+    {
+      return 'G';
+    }
+  else if(i==6)
+    {
+      return 'H';
+    }
+  else if(i==7)
+    {
+      return 'I';
+    }
+  else if(i==8)
+    {
+      return 'K';
+    }
+  else if(i==9)
+    {
+      return 'L';
+    }
+  else if(i==10)
+    {
+      return 'M';
+    }
+  else if(i==11)
+    {
+      return 'N';
+    }
+  else if(i==12)
+    {
+      return 'P';
+    }
+  else if(i==13)
+    {
+      return 'Q';
+    }
+  else if(i==14)
+    {
+      return 'R';
+    }
+  else if(i==15)
+    {
+      return 'S';
+    }
+  else if(i==16)
+    {
+      return 'T';
+    }
+  else if(i==17)
+    {
+      return 'V';
+    }
+  else if(i==18)
+    {
+      return 'W';
+    }
+  else if(i==19)
+    {
+      return 'Y';
+    }
+  else
+    {
+      cout<<"Invalid amino acid."<<endl;
+      exit(1);
+    }
+}
+
+int amino_acids::get_aa(char c)
+{
+  if(c=='A')
+    {
+      return 0;
+    }
+  else if(c=='C')
+    {
+      return 1;
+    }
+  else if(c=='D')
+    {
+      return 2;
+    }
+  else if(c=='E')
+    {
+      return 3;
+    }
+  else if(c=='F')
+    {
+      return 4;
+    }
+  else if(c=='G')
+    {
+      return 5;
+    }
+  else if(c=='H')
+    {
+      return 6;
+    }
+  else if(c=='I')
+    {
+      return 7;
+    }
+  else if(c=='K')
+    {
+      return 8;
+    }
+  else if(c=='L')
+    {
+      return 9;
+    }
+  else if(c=='M')
+    {
+      return 10;
+    }
+  else if(c=='N')
+    {
+      return 11;
+    }
+  else if(c=='P')
+    {
+      return 12;
+    }
+  else if(c=='Q')
+    {
+      return 13;
+    }
+  else if(c=='R')
+    {
+      return 14;
+    }
+  else if(c=='S')
+    {
+      return 15;
+    }
+  else if(c=='T')
+    {
+      return 16;
+    }
+  else if(c=='V')
+    {
+      return 17;
+    }
+  else if(c=='W')
+    {
+      return 18;
+    }
+  else if(c=='Y')
+    {
+      return 19;
+    }
+  else
+    {
+      return -1;
+    }
+}
+
+#endif

archaea_run_no_filtering.sh

+#!/bin/sh
+# $1 = full path to directory containing FAA directory
+# $2 = tag length
+
+./mif $1 $2
+#./sttag -p $1 -s B -f 1
+#./tmerg $1 $2
+#./filt -p $1 -f 10
+#./cm -p $1 -f -1 -o -1
+#./mdist $1
+#./fh $1

bitter.h

+#ifndef BITTER_H
+#define BITTER_H
+
+#include <iostream>
+#include <fstream>
+#include <vector>
+#include <cstdlib>
+#include <cmath>
+
+using namespace std;
+
+//const int MAX_KMERS_SIZE = 3000;
+
+class bitter
+{
+ private:  
+  int id1;
+  int id2; //id1 always < id2;
+  double min;
+  double numerator;
+  int denominator;
+  double rms_numerator;
+  double match_score_cutoff; //from prior SlopeTree run
+  double con_msc;
+  double hgt_msc;
+  vector<double> pids_id1;
+  vector<double> pids_id2;
+  double distance; //from prior SlopeTree run.
+  int nn1; //id1's nearest neighbor gnum
+  int nn2; //id2's nearest neighbor gnum
+  void set_nns(int n1, int n2);
+  int get_nn1();
+  int get_nn2();
+  
+ public:
+  bitter();