Determine TP, TN, FP, FN stats

FusionDetection/Evaluator/.gitignore

@@ -0,0 +1,8 @@
+bin/
+CMakeCache*
+CMakeFiles
+CTestTestfile*
+cmake_install.cmake
+testout*
+Makefile
+*.a

FusionDetection/Evaluator/src/BedCompare.cpp

@@ -1,4 +1,5 @@
 #include "BedCompare.h"
+#include <stdio.h>
 
 bool isStrandKnown(string strand)
 {
@@ -106,7 +107,7 @@ int  getNewTruePos(const bedpe_t & tt, const bedpe_t & rr, const uint32_t & tt_p
 }
 
 
-int transcript_compare(Bedpe & res, Bedpe & truth, int resolution, int max_diff, double & sensitivity, double & precision, pseudo_counts_t & pct)
+int transcript_compare(Bedpe & res, Bedpe & truth, int resolution, int max_diff, double & sensitivity, double & precision, pseudo_counts_t & pct, evaluate_t & et)
 {
     int diff=resolution-1;
     if(diff<0)
@@ -131,6 +132,7 @@ int transcript_compare(Bedpe & res, Bedpe & truth, int resolution, int max_diff,
         {
              bedpe_t tt=truth.getBedpe(i);
              bedpe_t rr=res.getBedpe(j);
+         	   //cerr<<"check " + tt.name + " " + rr.name<<endl;
 
              if(tt.chr1==rr.chr1 && tt.chr2==rr.chr2 && tt.strand1==rr.strand1 && tt.strand2==rr.strand2)
              {
@@ -170,14 +172,25 @@ int transcript_compare(Bedpe & res, Bedpe & truth, int resolution, int max_diff,
 
     for(int i=0;i<found.size();i++)
     {
-        if(found[i]==1) fd++;
+        if(found[i]==1) {
+          fd++;
+        } else {
+          // False negatives
+          et.false_negatives.push_back(truth.getBedpe(i));
+        }
     }
 
     int ct=0; 
     for(int j=0;j<correct.size();j++)
     {
-        if(correct[j]==1)
+        if(correct[j]==1) {
             ct++;
+          // True positive
+          et.true_positives.push_back(res.getBedpe(j));
+        } else {
+          // False positives
+          et.false_positives.push_back(res.getBedpe(j));
+        }
     }
 
     //code checking 0 is not here
@@ -197,6 +210,9 @@ int BedpeCompare::compare(Bedpe & res, Bedpe & truth, Gene & g, evaluate_t & et,
 //trans 1. 0 0; 2. compare;
     int num_res_trans=res.size()+pct.t_t+pct.f_t;;
     int num_truth_trans=truth.size()+pct.truth_t;
+    printf ("Num res trans: %d\n", num_res_trans);
+    printf ("Num truth trans: %d\n", num_truth_trans);
+
 
     et.num_res_trans=num_res_trans;
     et.num_truth_trans=num_truth_trans;   
@@ -250,7 +266,7 @@ int BedpeCompare::compare(Bedpe & res, Bedpe & truth, Gene & g, evaluate_t & et,
     else
     {
         double sens_t, prec_t;
-        transcript_compare(res, truth, resolution, max_diff, sens_t, prec_t, pct);
+        transcript_compare(res, truth, resolution, max_diff, sens_t, prec_t, pct, et);
         double f1_t=f_score(sens_t,prec_t);
 
         et.sensitivity_t=my_db2string(sens_t);
@@ -270,16 +286,16 @@ int BedpeCompare::compare(Bedpe & res, Bedpe & truth, Gene & g, evaluate_t & et,
 //gene 1. id set pair; 2. 0 0; 3. compare;
 //1..
     Bed2SetPair bp;
-    vector<set_pair_t> resSP;
-    bp.getSetPair(res, resSP, g);
-    vector<set_pair_t> truthSP;
-    bp.getSetPair(truth, truthSP, g);
+    bp.getSetPair(res, et.resSP, g);
+    bp.getSetPair(truth, et.truthSP, g);
     SetPair sp;   
-    sp.merge(resSP);
-    sp.merge(truthSP);       
+    sp.merge(et.resSP);
+    sp.merge(et.truthSP);       
 //2..
-    int num_res_gene=resSP.size()+pct.t_g+pct.f_g;
-    int num_truth_gene=truthSP.size()+pct.truth_g; 
+    int num_res_gene=et.resSP.size()+pct.t_g+pct.f_g;
+    int num_truth_gene=et.truthSP.size()+pct.truth_g; 
+    printf ("Num res gene: %d\n", num_res_gene);
+    printf ("Num truth gene: %d\n", num_truth_gene);
 
 
     et.num_res_gene=num_res_gene;
@@ -333,7 +349,7 @@ int BedpeCompare::compare(Bedpe & res, Bedpe & truth, Gene & g, evaluate_t & et,
     }
     else //3
     {
-        int num_inter_gene=sp.numIntersectSet(resSP, truthSP);
+        int num_inter_gene=sp.numIntersectSet(et.resSP, et.truthSP, et);
         num_inter_gene+=pct.t_g;
         //evaluate;
         double sens_g=(num_inter_gene+0.0)/num_truth_gene;

FusionDetection/Evaluator/src/ExecuteByRule.cpp

@@ -32,6 +32,7 @@ int ExecuteByRule::execute(Bedpe & res, Bedpe & truth, char * outputfile, char *
       string tmpfile1=string(outputfile)+"."+std::to_string((long long int)i)+"."+"result.bedpe";
       string tmpfile2=string(outputfile)+"."+std::to_string((long long int)i)+"."+"truth.bedpe";
       string ruleString=rtv[i].rule;
+	printf("\nRule: %s\n", ruleString.c_str());
       Bedpe  outRes;
       Bedpe  outTruth;
       cbr.cut(res, outRes, ruleString, tmpfile1);

FusionDetection/Evaluator/src/MyTypes.h

@@ -91,6 +91,12 @@ typedef struct
 
 } bedpe_t;
 
+typedef struct
+{
+    vector<int> ids1;
+    vector<int> ids2;
+} set_pair_t;
+
 typedef struct
 {
     string name;
@@ -109,14 +115,16 @@ typedef struct
     string precision_g;
     string f_g;
 
+    vector<bedpe_t> false_negatives;
+    vector<bedpe_t> false_positives;
+    vector<bedpe_t> true_positives;
+    vector<set_pair_t> gene_true_positives;
+    vector<set_pair_t> gene_false_positives;
+    vector<set_pair_t> gene_false_negatives;
+    vector<set_pair_t> resSP;
+    vector<set_pair_t> truthSP;
 } evaluate_t;
 
-typedef struct
-{
-    vector<int> ids1;
-    vector<int> ids2;
-} set_pair_t;
-
 typedef struct
 {
     int t_t;

FusionDetection/Evaluator/src/Printer.cpp

@@ -1,9 +1,9 @@
 #include "Printer.h"
 
-int Printer::print(vector<evaluate_t> & evals, char * outputfile,  string gene_file, int base_resolution, int max_diff, string pseudo_counts, string version)
+int Printer::print(vector<evaluate_t> & evals, string outputfile,  string gene_file, int base_resolution, int max_diff, string pseudo_counts, string version, Gene g)
 {
     ofstream ofs;
-    ofs.open (outputfile, std::ofstream::out);
+    ofs.open ((char *)outputfile.c_str(), std::ofstream::out);
 
     ofs<<"#fusionToolEvaluator version "+version<<endl;
     ofs<<"#--gene-annotation-file="+gene_file+"; --base-resolution="+my_db2string(base_resolution)+"; --max-diff="+my_db2string(max_diff)+"; --pseudo-counts="+pseudo_counts<<endl;
@@ -22,6 +22,89 @@ int Printer::print(vector<evaluate_t> & evals, char * outputfile,  string gene_f
         ofs<<et.sensitivity_g<<"\t";
         ofs<<et.precision_g<<"\t";
         ofs<<et.f_g<<"\n";
+
+        ofstream ofs_fuse;
+        ofs_fuse.open ((char *)(outputfile + std::to_string(i) + ".fusions").c_str(), std::ofstream::out);
+
+        // Dump transcript fusion results (TP, FP, FN)
+        for (int i = 0; i < et.false_negatives.size(); i++)
+        {
+            std::string n = et.false_negatives[i].name;
+            n = n.replace(n.find("-"), 1, "\t");
+            ///std::replace(n.begin(), n.end(), "-", "\t");
+            ofs_fuse<<"T_FN\t" + n <<"\n";
+        }
+
+        for (int i = 0; i < et.false_positives.size(); i++)
+        {
+            std::string n = et.false_positives[i].name;
+            n = n.replace(n.find("-"), 1, "\t");
+            //std::replace(n.begin(), n.end(), '-', '\t');
+            ofs_fuse<<"T_FP\t" + n <<"\n";
+        }
+
+        for (int i = 0; i < et.true_positives.size(); i++)
+        {
+            std::string n = et.true_positives[i].name;
+            n = n.replace(n.find("-"), 1, "\t");
+            //std::replace(n.begin(), n.end(), '-', '\t');
+            ofs_fuse<<"T_TP\t" + n <<"\n";
+        }
+
+        for (int i = 0; i < et.gene_true_positives.size(); i++)
+        {
+          set_pair_t sp = et.gene_true_positives[i];
+          for (int j = 0; j < sp.ids1.size(); j++)
+          {
+            gene_t *x = g.getGene(sp.ids1[j]);
+            for (int k = 0; k < sp.ids2.size(); k++)
+            {
+              gene_t *y = g.getGene(sp.ids2[k]);
+              ofs_fuse<<"G_TP\t" + x->name2 + "\t" + y->name2<<"\n";
+            }
+          }
+        }
+
+        for (int i = 0; i < et.gene_false_positives.size(); i++)
+        {
+          set_pair_t sp = et.gene_false_positives[i];
+          for (int j = 0; j < sp.ids1.size(); j++)
+          {
+            gene_t *x = g.getGene(sp.ids1[j]);
+            for (int k = 0; k < sp.ids2.size(); k++)
+            {
+              gene_t *y = g.getGene(sp.ids2[k]);
+              ofs_fuse<<"G_FP\t" + x->name2 + "\t" + y->name2<<"\n";
+            }
+          }
+        }
+
+        for (int i = 0; i < et.gene_false_negatives.size(); i++)
+        {
+          set_pair_t sp = et.gene_false_negatives[i];
+          for (int j = 0; j < sp.ids1.size(); j++)
+          {
+            gene_t *x = g.getGene(sp.ids1[j]);
+            for (int k = 0; k < sp.ids2.size(); k++)
+            {
+              gene_t *y = g.getGene(sp.ids2[k]);
+              ofs_fuse<<"G_FN\t" + x->name2 + "\t" + y->name2<<"\n";
+            }
+          }
+        }
+
+        /*
+        ofs<<"resSP\n";
+                  ofs<<et.resSP[i].ids1+"-"+et.resSP[i].ids2<<"\n";
+              }
+
+        ofs<<"\n";
+        ofs<<"truthSP\n";
+        for (int i = 0; i < et.truthSP.size(); i++) {
+                  ofs<<et.truthSP[i].ids1+"-"+et.truthSP[i].ids2<<"\n";
+              }
+          */
+        ofs_fuse.close();
     }
     ofs.close();
 

FusionDetection/Evaluator/src/Printer.h

@@ -7,6 +7,7 @@
 #ifndef PRINTER_H_
 #define PRINTER_H_
 
+#include "Gene2.h"
 #include "MyTypes.h"
 #include "Util.h"
 
@@ -19,7 +20,7 @@ class Printer
 {
     public:
         Printer(){};
-        int print(vector<evaluate_t> & evals, char * outputfile,  string gene_file, int base_resolution, int max_diff, string pseudo_counts, string version);
+        int print(vector<evaluate_t> & evals, string outputfile,  string gene_file, int base_resolution, int max_diff, string pseudo_counts, string version, Gene g);
 };
 
 

FusionDetection/Evaluator/src/SetPair.cpp

@@ -67,22 +67,40 @@ int SetPair::merge(vector<set_pair_t> & spv)
 }
 
 //use after merge
-int SetPair::numIntersectSet(vector<set_pair_t> & spv1,vector<set_pair_t> & spv2)
+int SetPair::numIntersectSet(vector<set_pair_t> & res,vector<set_pair_t> & truth, evaluate_t & et)
 {
-    vector<int> intersect(spv1.size(),0);
-    for(int i=0;i<spv1.size();i++)
+    vector<int> found(truth.size(),0);
+    vector<int> correct(res.size(),0);
+
+    for(int i=0;i<res.size();i++)
     {
-        for(int j=0;j<spv2.size();j++)
+        for(int j=0;j<truth.size();j++)
         {
-            int inter1=numIntersect(spv1[i].ids1, spv2[j].ids1);
-            int inter2=numIntersect(spv1[i].ids2, spv2[j].ids2);
-            if(inter1>0 && inter2>0)
-                intersect[i]=1;
+            int inter1=numIntersect(res[i].ids1, truth[j].ids1);
+            int inter2=numIntersect(res[i].ids2, truth[j].ids2);
+            if(inter1>0 && inter2>0) {
+                found[j]=1;
+                correct[i]=1;
+                et.gene_true_positives.push_back(truth[j]);
+            }
         }
     }
+
     int num=0;
-    for(int i=0;i<intersect.size();i++)
-        num+=intersect[i];
+    for(int i=0;i<found.size();i++) {
+      if (found[i] == 1) {
+        num += 1;
+      } else {
+        et.gene_false_negatives.push_back(truth[i]);
+      }
+    }
+
+    for(int i=0;i<correct.size();i++) {
+      if (correct[i] == 0) {
+        et.gene_false_positives.push_back(res[i]);
+      }
+    }
+
     return num;
 }
 

FusionDetection/Evaluator/src/SetPair.h

@@ -20,7 +20,7 @@ class SetPair
     public:
        SetPair(){};
        int merge(vector<set_pair_t> & spv);
-       int numIntersectSet(vector<set_pair_t> & spv1,vector<set_pair_t> & spv2); 
+       int numIntersectSet(vector<set_pair_t> & spv1,vector<set_pair_t> & spv2, evaluate_t & et); 
 };
 
 

FusionDetection/Evaluator/src/main.cpp

@@ -181,7 +181,7 @@ int main(int argc, char * argv[])
     cerr<<"printing evaluation results..."<<endl;
 
     Printer per;
-    per.print(evaluates, (char*)output_file.c_str(), gene_file, base_resolution, max_diff, pseudo_counts, version);
+    per.print(evaluates, (char*)output_file.c_str(), gene_file, base_resolution, max_diff, pseudo_counts, version, g);
 
     return 0;
-}
+}