Add --cycle parameter to set max number for read balancing iterations

include/rpkm_calculation.h

@@ -34,7 +34,7 @@ vector <double> get_sum_density_by_all_intervals (const vector <vector <double>
 void update_isoforms_density_to_average_for_isoform (vector <vector <double> > & density_array, vector <vector <double> > & unique_density_array, bool keep_unique);
 void print_array (const vector <double> & intput_array, const string & title, string tab_spacer = "   ");
 void adjust_isoforms_density_by_coef (vector <vector <double> > & weight_array, const vector <double> & original_densities_by_isoforms,  const vector <double> & new_densities_by_isoforms);
-int run_cycle (vector <vector <double> > & weight_array, double & res_sum, vector <vector <double> > & unique_density_array, bool keep_unique);
+int run_cycle (vector <vector <double> > & weight_array, double & res_sum, vector <vector <double> > & unique_density_array, bool keep_unique, int max_cycles);
 double sum_all (const vector <vector <double> > dens_matrix);
 void subtract_matrix (vector <vector <double> > & first, const vector <vector <double> > & second);
 void calculate_totReads_density (const vector<vector<double> > & weight_array, std::map <string, Isoform> & iso_map,

src/main.cpp

@@ -46,6 +46,7 @@ bool verify_params (cxxopts::Options params){
         cerr << "  exclude: " << params["exclude"].as<std::string>() << endl;
         cerr << "  threads: " << params["threads"].as<int>() << endl;
         cerr << "  keepUnique: " << params["keepUnique"].as<bool>() << endl;
+        cerr << "  cycles: " << params["cycles"].as<int>() << endl;
         cerr << "  dutp: " << params["dutp"].as<bool>() << endl;
     } catch (...){
         cerr << "Some of the parameters have a mistake" << endl;
@@ -67,6 +68,10 @@ bool verify_params (cxxopts::Options params){
         cerr << "minReadLen cannot be less than 0" << endl;
         return false;
     }
+    if (params["cycles"].as<int>() < 0){
+        cerr << "cycles cannot be less than 0" << endl;
+        return false;
+    }    
     if (params["threads"].as<int>() < 1){
         cerr << "threads cannot be less than 1" << endl;
         return false;
@@ -102,6 +107,8 @@ int main(int argc, char **argv) {
              "Set the number of threads")
             ("u,keepUnique", "flag to fix unique reads for the specific isoromf interval", cxxopts::value<bool>(),
              "Set this flag if you want prevent distributing the isoform unique reads among other isoforms")
+            ("s,cycles", "Max cycles for read balancing", cxxopts::value<int>()->default_value("2000"),
+             "Set the maximum number of cycles used for read balancing")             
             ("d,dutp", "set dutp enabled", cxxopts::value<bool>(),
              "Set this dutp flag if strand specific analysys should be made");
 

src/rpkm_calculation.cpp

@@ -178,7 +178,7 @@ void subtract_matrix (vector <vector <double> > & first, const vector <vector <d
 }
 
 
-int run_cycle (vector <vector <double> > & density_array, double & res_sum, vector <vector <double> > & unique_density_array, bool keep_unique){
+int run_cycle (vector <vector <double> > & density_array, double & res_sum, vector <vector <double> > & unique_density_array, bool keep_unique, int max_cycles){
     double cutoff = 10e-9; // TODO put it in separate configuration file
     int cycles = 0;
     vector <vector <double> > tmp_matrix;
@@ -188,7 +188,7 @@ int run_cycle (vector <vector <double> > & density_array, double & res_sum, vect
 
 //    print_array (original_sum_dens, "Original density sums");
 
-    for (int i = 0; i < 2000; i++){
+    for (int i = 0; i < max_cycles; i++){
         // Update original density array with average values for isoforms
 //        cout << endl << "Set average by row" << endl;
         update_isoforms_density_to_average_for_isoform (density_array, unique_density_array, keep_unique);

src/thread.cpp

@@ -34,6 +34,7 @@ void process (   vector < std::map <string, multimap <long, GffRecordPtr> >::ite
 
     int min_length = params["minIntLen"].as<int>();
     bool keep_unique = params["keepUnique"].as<bool>();
+    int max_cycles = params["cycles"].as<int>();
     int min_read_segment_length = params["minReadLen"].as<int>();;
     bool dUTP = params["dutp"].as<bool>();
     string bam_full_path_name = params["bam"].as<string>();
@@ -409,7 +410,7 @@ void process (   vector < std::map <string, multimap <long, GffRecordPtr> >::ite
 
 //            cerr << "[" << thread_number << "] " << "Started to run cycles" << endl;
             double res_sum;
-            int cycles = run_cycle(weight_array, res_sum, unique_weight_array, keep_unique);
+            int cycles = run_cycle(weight_array, res_sum, unique_weight_array, keep_unique, max_cycles);
 //            cerr << "[" << thread_number << "] " << "Finished to run cycles : " << cycles << endl;
 //            cerr << "[" << thread_number << "] " << "Result sum : " << res_sum << endl;
             print_weight_array(weight_array, correspondence_map, min_weight, "Final density array");