released/pIRS.old/stat_illumina_reads/gcvgstater/main.cpp

   1 #include <iostream>
   2 #include <string>
   3 #include <vector>
   4 #include <map>
   5 #include <iterator>
   6 #include <cstdlib>
   7 #include "gzstream.h"
   8 #include "self_util.h"
   9 #include "stat_soap_coverage.h"
  10
  11 using namespace std;
  12
  13 void usage()
  14 {
  15     cerr << "Description:" << endl
  16         << "    it is a program to stat about the gc_depth and base coverage ratio and base_depth with *.soap.coverage  result." << endl
  17         << "Program:" << endl
  18         << "Name:           gc_coverage_bias" << endl
  19         << "Compile Date:   2011-08-03" << endl
  20         << "Author:         ganjun" << endl
  21         << "Version:        1.01" << endl
  22         << "Contact:        ganjun@genomics.org,cn" << endl;
  23
  24     cerr << "\n\nUsage:\tgc_coverage_bias [options] <*.soap.coverage>" << endl
  25         << "Option: -r <string>     a reference sequence file about FA format " << endl
  26         << "        -o <string>     the the prefix about output file" << endl
  27         << "        -w <string>     the window length[such as:100,200,300] " << endl
  28         << "        --gcdump        output the gc ratio in the window length " << endl
  29         << "        --depwindump    output the avg depth in the window length " <<    endl;
  30
  31     cerr << "\n\nUsage:" << endl
  32         << "\t./gc_coverage_bias -r test.fa -o test -w 40,50 [--gcdump --depwindump] test.depth" << endl;
  33     exit(EXIT_FAILURE);
  34 }
  35
  36 /*
  37 map<string, string> DealReference(string str_ref_file_name);
  38 map<string, vector<vector<uint64_t, uint64_t, uint64_t> >, uint64_t> DealSoapCoverage()
  39 */
  40
  41
  42 int main(int argc, char* argv[])
  43 {
  44     if(argc < 8)
  45     {
  46         usage();
  47     }
  48
  49     string str_argv;
  50     string str_ref_file_name;
  51     string str_width;
  52     string str_output_prefix;
  53     bool b_gcdump = false;
  54     bool b_depwindump = false;
  55     vector<string> vec_soap_file_name;
  56     vector<string> vec_width;
  57
  58     for(int i=1; i<argc; ++i)
  59     {
  60         str_argv = argv[i];
  61
  62         if(str_argv.compare("-r") == 0)
  63         {
  64             ++i;
  65             str_ref_file_name = argv[i];
  66             continue;
  67         }
  68
  69         if(str_argv.compare("-o") == 0)
  70         {
  71             ++i;
  72             str_output_prefix = argv[i];
  73             continue;
  74         }
  75
  76         if(str_argv.compare("-w") == 0)
  77         {
  78             ++i;
  79             str_width = argv[i];
  80             continue;
  81         }
  82
  83         if(str_argv.compare("--gcdump") == 0)
  84         {
  85             b_gcdump = true;
  86             continue;
  87         }
  88
  89         if(str_argv.compare("--depwindump") == 0)
  90         {
  91             b_depwindump = true;
  92             continue;
  93         }
  94
  95         vec_soap_file_name.push_back(argv[i]);
  96     }
  97
  98     if(str_ref_file_name.empty())
  99     {
 100         usage();
 101     }else
 102     {
 103         igzstream in(str_ref_file_name.c_str());
 104
 105         if(!in)
 106         {
 107             cerr << "can't open the reference file " << str_ref_file_name << ", please check!" << endl;
 108             exit(EXIT_FAILURE);
 109         }
 110
 111         in.close();
 112     }
 113
 114     if(str_output_prefix.empty())
 115     {
 116         usage();
 117     }
 118
 119     if(str_width.empty())
 120     {
 121         usage();
 122     }
 123     else
 124     {
 125         vec_width = splitString(str_width, ",");
 126     }
 127
 128     //cerr << vec_soap_file_name.size() << endl;
 129     for(vector<string>::iterator it = vec_soap_file_name.begin(); it != vec_soap_file_name.end(); ++it)
 130     {
 131         igzstream in((*it).c_str());
 132
 133         if(!in)
 134         {
 135             cerr << "can't open the reference file " << *it << ", please check!" << endl;
 136             exit(EXIT_FAILURE);
 137         }
 138
 139         in.close();
 140     }
 141
 142     stat_soap_coverage test(str_ref_file_name, str_output_prefix, vec_soap_file_name, vec_width, b_gcdump, b_depwindump);
 143 }
 144
 145 /*
 146 map<string, string> DealReference(string str_ref_file_name)
 147 {
 148         igzstream in(str_ref_file.name);
 149         string line;
 150         map<string, string> map_reference_base;
 151         string keyname;
 152         string sequence;
 153
 154         while(getline(in, line))
 155         {
 156                 TrimLeft(line);
 157                 TrimRight(line);
 158
 159                 if(line[0] == '>')
 160                 {
 161                          if(sequence.length() != 0)
 162                         {
 163                                  map_reference_base[keyname] = sequence;
 164                         }
 165                         int index;
 166                         if(((index = line.find(" ")) == string::npos) || ((index = line.find("\t")) == string::npos) || ((index = line.find("\n")) == string::npos))
 167                         {
 168                                 keyname = line.substr(1, index);
 169                         }
 170                         sequence.clear();
 171                         continue;
 172                 }
 173
 174                 sequence += line;
 175         }
 176
 177          if(sequence.length() != 0)
 178         {
 179                  map_reference_base[keyname] = sequence;
 180         }
 181
 182         in.close();
 183 }
 184
 185 */