New option --cf-quartet to print sCF of all resampled quartets (requested by Benjamin Rosenzweig and Matt Hahn)

Author: bqminh

main/phyloanalysis.cpp

@@ -4095,6 +4095,42 @@ void assignBranchSupportNew(Params &params) {
     }
     out.close();
     cout << "Concordance factors per branch printed to " << filename << endl;
+    
+    if (params.print_cf_quartets) {
+        filename = prefix + ".cf.quartet";
+        out.open(filename);
+        out << "# Site concordance factor for all resampled quartets (with replacement)" << endl
+            << "# This file can be read in MS Excel or in R with command:" << endl
+            << "#   tab=read.table('" <<  filename << "',header=TRUE)" << endl
+            << "# Columns are tab-separated with following meaning:" << endl
+            << "#   ID: Branch ID" << endl
+            << "#   QuartID: Quartet ID" << endl
+            << "#   Seq1: ID of sequence 1 on 'left' side of the branch" << endl
+            << "#   Seq2: ID of sequence 2 on 'left' side of the branch" << endl
+            << "#   Seq3: ID of sequence 3 on 'right' side of the branch" << endl
+            << "#   Seq4: ID of sequence 4 on 'right' side of the branch" << endl
+            << "#   qCF: Fraction of concordant sites supporting quartet Seq1,Seq2|Seq3,Seq4 (=qCF_N/qN)" << endl
+            << "#   qCF_N: Number of concordant sites supporting quartet Seq1,Seq2|Seq3,Seq4" << endl
+            << "#   qDF1: Fraction of discordant sites supporting quartet Seq1,Seq3|Seq2,Seq4  (=qDF1_N/qN)" << endl
+            << "#   qDF1_N: Number of discordant sites supporting quartet Seq1,Seq3|Seq2,Seq4" << endl
+            << "#   qDF2: Fraction of discordant sites supporting quartet Seq1,Seq4|Seq2,Seq3 (=qDF2_N/qN)" << endl
+            << "#   qDF2_N: Number of discordant sites supporting quartet Seq1,Seq4|Seq2,Seq3" << endl
+            << "#   qN: Number of decisive sites with four taxa Seq1,Seq2,Seq3,Seq4 (=qCF_N+qDF1_N+qDF2_N)" << endl
+            << "ID\tQuartID\tSeq1\tSeq2\tSeq3\tSeq4\tqCF\tqCF_N\tqDF1\tqDF1_N\tqDF2\tqDF2_N\tqN" << endl;
+        for (brit = branches.begin(); brit != branches.end(); brit++) {
+            Neighbor *branch = brit->second->findNeighbor(brit->first);
+            int ID = brit->second->id;
+            for (int qid = 0; ; qid++) {
+                string qstr;
+                if (branch->attributes.find("q" + convertIntToString(qid)) == branch->attributes.end())
+                    break;
+                out << ID << '\t' << qid+1 << '\t' << branch->attributes["q" + convertIntToString(qid)] << endl;
+            }
+        }
+        out.close();
+        cout << "Site concordance factors for quartets printed to " << filename << endl;
+    }
+    
     if (!params.site_concordance_partition)
         return;
 

tree/discordance.cpp

@@ -189,6 +189,7 @@ void PhyloTree::computeSiteConcordance(Branch &branch, int nquartets, int *rstre
             }
         }
     }
+    Neighbor *nei = branch.second->findNeighbor(branch.first);
     for (i = 0; i < nquartets; i++) {
         int j;
         // get a random quartet
@@ -209,6 +210,16 @@ void PhyloTree::computeSiteConcordance(Branch &branch, int nquartets, int *rstre
             sDF1_N += support[1];
             sDF2_N += support[2];
         }
+        if (params->print_cf_quartets) {
+            // print sCF for each quartet
+            stringstream ss;
+            ss << quartet[0]+1 << '\t' << quartet[1]+1 << '\t' << quartet[2]+1 << '\t' << quartet[3]+1
+               << '\t' << (((double)support[0]) / sum) << '\t' << support[0]
+               << '\t' << (((double)support[1]) / sum) << '\t' << support[1]
+               << '\t' << (((double)support[2]) / sum) <<  '\t' << support[2]
+               << '\t' << sum;
+            nei->putAttr("q" + convertIntToString(i), ss.str());
+        }
     }
     sN = (double)sum_sites / nquartets;
     // rounding
@@ -218,7 +229,6 @@ void PhyloTree::computeSiteConcordance(Branch &branch, int nquartets, int *rstre
     sCF_N = round(sCF_N / nquartets * 100)/100;
     sDF1_N = round(sDF1_N / nquartets * 100)/100;
     sDF2_N = round(sDF2_N / nquartets * 100)/100;
-    Neighbor *nei = branch.second->findNeighbor(branch.first);
     PUT_ATTR(nei, sCF);
     PUT_ATTR(nei, sN);
     PUT_ATTR(nei, sDF1);

utils/tools.cpp

@@ -811,6 +811,7 @@ void parseArg(int argc, char *argv[], Params &params) {
     params.support_tag = NULL;
     params.site_concordance = 0;
     params.site_concordance_partition = false;
+    params.print_cf_quartets = false;
     params.print_df1_trees = false;
     params.internode_certainty = 0;
     params.tree_weight_file = NULL;
@@ -1586,6 +1587,10 @@ void parseArg(int argc, char *argv[], Params &params) {
                 params.site_concordance_partition = true;
                 continue;
             }
+            if (strcmp(argv[cnt], "--cf-quartet") == 0) {
+                params.print_cf_quartets = true;
+                continue;
+            }
             if (strcmp(argv[cnt], "--df-tree") == 0) {
                 params.print_df1_trees = true;
                 continue;
@@ -4584,6 +4589,7 @@ void usage_iqtree(char* argv[], bool full_command) {
     << "  -s FILE              Sequence alignment for --scf" << endl
     << "  -p FILE|DIR          Partition file or directory for --scf" << endl
     << "  --cf-verbose         Write CF per tree/locus to cf.stat_tree/_loci" << endl
+    << "  --cf-quartet         Write sCF for all resampled quartets to .cf.quartet" << endl
 
 #ifdef USE_LSD2
     << endl << "TIME TREE RECONSTRUCTION:" << endl

utils/tools.h

@@ -1268,6 +1268,9 @@ class Params {
      TRUE to print concordant sites per partition
      */
     bool site_concordance_partition;
+    
+    /** TRUE to print sCF for all sampled quartet */
+    bool print_cf_quartets;
 
     /** TRUE to print trees associated with discordance factor 1 (NNI-1 tree) */
     bool print_df1_trees;