fix bug set minblen=-1 in relative date case, and do not remove outgroup

Author: Thu-Hien To

bin/lsd2_mac

@@ -46,7 +46,6 @@ int lsd::buildTimeTree( int argc, char** argv, InputOutputStream *inputOutput)
     if (io->inOutgroup){
         extrait_outgroup(io, opt);
     }
-    ifstream gr(opt->rate.c_str());
     *(io->outTree1)<<"#NEXUS\n";
     *(io->outTree2)<<"#NEXUS\n";
     bool constraintConsistent=true;
@@ -100,7 +99,7 @@ int lsd::buildTimeTree( int argc, char** argv, InputOutputStream *inputOutput)
         double br=0;
         if (opt->givenRate[0]){
             string line;
-            if( getline(gr, line)) {
+            if( getline(*(io->inRate), line)) {
                 vector<double> all_rates = read_double_from_line(line);
                 opt->rho = all_rates[0];
                 if (all_rates.size() > 1){
@@ -254,7 +253,6 @@ int lsd::buildTimeTree( int argc, char** argv, InputOutputStream *inputOutput)
     *(io->outTree1)<<"End;\n";
     *(io->outTree2)<<"End;\n";
     cout<<"\nTOTAL ELAPSED TIME: "<<elapsed_time<<" seconds"<<endl;
-    gr.close();
     if (!inputOutput)
         delete io;
     return EXIT_SUCCESS;

bin/lsd2_unix

@@ -26,6 +26,9 @@ class InputOutputStream {
     /** input partition stream */
     istream *inPartition;
 
+    /** input rate stream */
+    istream *inRate;
+    
     /** output result stream */
     ostream *outResult;
 
@@ -47,7 +50,7 @@ class InputOutputStream {
      @param outgroup outgroup string
      @param date date string
      */
-    InputOutputStream(string tree, string outgroup, string date,string partition);
+    InputOutputStream(string tree, string outgroup, string date,string rate,string partition);
 
     /** destructor */
     virtual ~InputOutputStream();
@@ -66,16 +69,22 @@ class InputOutputStream {
     virtual void setOutgroup(string str);
 
     /**
-     set the content of the outgroup stream
+     set the content of the date stream
      @param str a string
      */
     virtual void setDate(string str);
 
     /**
-     set the content of the outgroup stream
+     set the content of the partition stream
      @param str a string
      */
     virtual void setPartition(string str);
+    
+    /**
+     set the content of the rate stream
+     @param str a string
+     */
+    virtual void setRate(string str);
 };
 
 /**

src/lsd.cpp

@@ -25,7 +25,7 @@ Pr* getOptions( int argc, char** argv )
 
 Pr* getCommandLine( int argc, char** argv)
 {
-    const string VERSION="v1.7";
+    const string VERSION="v1.7.1";
     Pr* opt = new Pr();
     int c;
     string s;
@@ -362,7 +362,7 @@ Pr* getInterface()
 
 void printInterface(ostream& in, Pr* opt)
 {
-    const string VERSION = "v1.7";
+    const string VERSION = "v1.7.1";
 
     in<<"\nLEAST-SQUARE METHODS TO ESTIMATE RATES AND DATES - "<<VERSION<<" \n\n";
     in<<"\nInput files:\n";
@@ -505,7 +505,7 @@ void printHelp( void )
     const string BOLD = "\033[00;01m";
     const string LINE = "\033[00;04m";
     const string FLAT = "\033[00;00m";
-    const string VERSION = "v1.7";
+    const string VERSION = "v1.7.1";
 
     cout<<BOLD<<"LSD: LEAST-SQUARES METHODS TO ESTIMATE RATES AND DATES - "<<VERSION<<"\n\n";
     cout<<BOLD<<"DESCRIPTION\n"

src/lsd.h

@@ -418,7 +418,8 @@ void extrait_outgroup(InputOutputStream *io, Pr* pr){
     list<string> outgroups = getOutgroup(*(io->inOutgroup), pr->fnOutgroup);
     ostringstream w;
     for (int y=0;y<pr->nbData;y++){
-        cout<<"Removing outgroups of tree "<<y+1<<" ...\n";
+        if (pr->keepOutgroup) cout<<"Reroot the tree "<<y+1<<" using given outgroups ...\n";
+        else cout<<"Removing outgroups of tree "<<y+1<<" ...\n";
         Node** nodes = tree2data(*(io->inTree),pr,s);
         if (!pr->rooted) {
             nodes=unrooted2rootedS(pr, nodes, s);
@@ -432,23 +433,22 @@ void extrait_outgroup(InputOutputStream *io, Pr* pr){
             Node** nodes_new = cloneLeaves(pr,nodes,0);
             int p_r=reroot_rootedtree(r, pr, nodes, nodes_new);
             computeSuc_polytomy(pr, nodes_new);
-            if (pr->keepOutgroup) {
-                w << newick(0,0,pr,nodes_new,nbTips);
+            if (keepBelow) {
+                w << newick(r, r, pr, nodes_new,nbTips);
             }
             else{
-                if (keepBelow) {
-                    w << newick(r, r, pr, nodes_new,nbTips);
-                }
-                else{
-                    w << newick(p_r, p_r,pr, nodes_new,nbTips).c_str();
-                }
+                w << newick(p_r, p_r,pr, nodes_new,nbTips).c_str();
             }
-            for (int i=0;i<=pr->nbBranches;i++) delete nodes_new[i];
-            delete[] nodes_new;
             if ((nbTips+outgroups.size()) != (pr->nbBranches+1 - pr->nbINodes)){
                 cerr<<"Error: The outgroups do not form a monophyletic in the input tree."<<endl;
                 exit(EXIT_FAILURE);
             }
+            if (pr->keepOutgroup) {
+                w.str("");
+                w << newick(0,0,pr,nodes_new,nbTips);
+            }
+            for (int i=0;i<=pr->nbBranches;i++) delete nodes_new[i];
+            delete[] nodes_new;
         }
         else{
             cout<<"The outgroups are not in the tree "<<y+1<<endl;

src/options.cpp

@@ -21,21 +21,24 @@ InputOutputStream::InputOutputStream () {
     inTree = nullptr;
     inOutgroup = nullptr;
     inDate = nullptr;
+    inRate = nullptr;
     inPartition = nullptr;
     outResult = nullptr;
     outTree1 = nullptr;
     outTree2 = nullptr;
     outTree3 = nullptr;
 }
 
-InputOutputStream::InputOutputStream(string tree, string outgroup, string date, string partition) {
+InputOutputStream::InputOutputStream(string tree, string outgroup, string date, string rate, string partition) {
     inTree = nullptr;
     inOutgroup = nullptr;
     inDate = nullptr;
+    inRate = nullptr;
     inPartition = nullptr;
     setTree(tree);
     setOutgroup(outgroup);
     setDate(date);
+    setRate(rate);
     setPartition(partition);
     outResult = new ostringstream;
     outTree1 = new ostringstream;
@@ -60,6 +63,10 @@ InputOutputStream::~InputOutputStream() {
         delete inPartition;
         inPartition = nullptr;
     }
+    if (inRate) {
+        delete inRate;
+        inRate = nullptr;
+    }
     if (outResult) {
         delete outResult;
         outResult = nullptr;
@@ -108,6 +115,14 @@ void InputOutputStream::setPartition(string str) {
     inPartition = new istringstream(str);
 }
 
+void InputOutputStream::setRate(string str) {
+    if (str.empty())
+        return;
+    if (inRate)
+        delete inRate;
+    inRate = new istringstream(str);
+}
+
 InputOutputFile::InputOutputFile(Pr *opt) : InputOutputStream() {
     treeIsFile = true;
     // open the tree file
@@ -138,7 +153,7 @@ InputOutputFile::InputOutputFile(Pr *opt) : InputOutputStream() {
         }
     }
 
-    // open date file
+    // open partition file
     if (opt->partitionFile != "") {
         ifstream *part_file = new ifstream(opt->partitionFile);
         inPartition = part_file;
@@ -148,6 +163,17 @@ InputOutputFile::InputOutputFile(Pr *opt) : InputOutputStream() {
         }
     }
 
+    // open given rate file
+    if (opt->rate != "") {
+        ifstream *rate_file = new ifstream(opt->rate);
+        inRate = rate_file;
+        if (!rate_file->is_open()) {
+            cerr << "Error: cannot open rate file " << opt->rate << endl;
+            exit(EXIT_FAILURE);
+        }
+    }
+
+    
     // open the result file
     ofstream *result_file = new ofstream(opt->outFile);
     outResult = result_file;
@@ -191,6 +217,9 @@ InputOutputFile::~InputOutputFile() {
     if (inPartition) {
         ((ifstream*)inPartition)->close();
     }
+    if (inRate) {
+        ((ifstream*)inRate)->close();
+    }
     if (outResult) {
         ((ofstream*)outResult)->close();
     }
@@ -2835,6 +2864,7 @@ void imposeMinBlen(ostream& file,Pr* pr, Node** nodes, double median_rate,bool m
             }
         }
         else {
+            minblen = m;
             if (pr->minblenL < 0){
                 cout<<"Minimum branch length of time scaled tree (settable via option -u and -U): "<<m<<endl;
                 file<<"Minimum branch length of time scaled tree (settable via option -u and -U): "<<m<<"\n";