smoothing for long read plotting

OpenGene · Nov 20, 2017 · 708c858 · 708c858
1 parent 605e833
commit 708c858
Show file tree

Hide file tree

Showing 2 changed files with 68 additions and 35 deletions.
diff --git a/src/stats.cpp b/src/stats.cpp
@@ -330,7 +330,21 @@ string Stats::list2string(double* list, int size, long* coords) {
     stringstream ss;
     for(int i=0; i<size; i++) {
         // coords is 1,2,3,...
-        ss << list[coords[i]-1];
+        long start = 0;
+        if(i>0)
+            start = coords[i-1];
+        long end = coords[i];
+
+        double total = 0.0;
+        for(int k=start; k<end; k++)
+            total += list[k];
+
+        // get average
+        if(end == start)
+            ss << "0.0";
+        else
+            ss << total / (end - start);
+        //ss << list[coords[i]-1];
         if(i < size-1)
             ss << ",";
     }
@@ -352,6 +366,10 @@ void Stats::reportHtml(ofstream& ofs, string filteringType, string readName) {
     reportHtmlContents(ofs, filteringType, readName);
 }
 
+bool Stats::isLongRead() {
+    return mCycles > 300;
+}
+
 void Stats::reportHtmlQuality(ofstream& ofs, string filteringType, string readName) {
 
     // quality
@@ -369,26 +387,33 @@ void Stats::reportHtmlQuality(ofstream& ofs, string filteringType, string readNa
 
     long *x = new long[mCycles];
     int total = 0;
-    const int fullSampling = 300;
-    for(int i=0; i<mCycles && i<fullSampling; i++){
-        x[total] = i+1;
-        total++;
-    }
-    // down sampling if it's too long
-    if(mCycles>fullSampling) {
-        double pos = fullSampling;
-        while(true){
-            pos *= 1.05;
-            if(pos >= mCycles)
-                break;
-            x[total] = (int)pos;
+    if(!isLongRead()) {
+        for(int i=0; i<mCycles; i++){
+            x[total] = i+1;
             total++;
         }
-        // make sure lsat one is contained
-        if(x[total-1] != mCycles){
-            x[total] = mCycles;
+    } else {
+        const int fullSampling = 40;
+        for(int i=0; i<fullSampling && i<mCycles; i++){
+            x[total] = i+1;
             total++;
         }
+        // down sampling if it's too long
+        if(mCycles>fullSampling) {
+            double pos = fullSampling;
+            while(true){
+                pos *= 1.05;
+                if(pos >= mCycles)
+                    break;
+                x[total] = (int)pos;
+                total++;
+            }
+            // make sure lsat one is contained
+            if(x[total-1] != mCycles){
+                x[total] = mCycles;
+                total++;
+            }
+        }
     }
     // four bases
     for (int b = 0; b<5; b++) {
@@ -404,7 +429,7 @@ void Stats::reportHtmlQuality(ofstream& ofs, string filteringType, string readNa
     json_str += "];\n";
     json_str += "var layout={title:'" + title + "', xaxis:{title:'cycles'";
     // use log plot if it's too long
-    if(mCycles>fullSampling) {
+    if(isLongRead()) {
         json_str += ",type:'log'";
     }
     json_str += "}, yaxis:{title:'quality'}};\n";
@@ -433,26 +458,33 @@ void Stats::reportHtmlContents(ofstream& ofs, string filteringType, string readN
 
     long *x = new long[mCycles];
     int total = 0;
-    const int fullSampling = 300;
-    for(int i=0; i<mCycles && i<fullSampling; i++){
-        x[total] = i+1;
-        total++;
-    }
-    // down sampling if it's too long
-    if(mCycles>fullSampling) {
-        double pos = fullSampling;
-        while(true){
-            pos *= 1.05;
-            if(pos >= mCycles)
-                break;
-            x[total] = (int)pos;
+    if(!isLongRead()) {
+        for(int i=0; i<mCycles; i++){
+            x[total] = i+1;
             total++;
         }
-        // make sure lsat one is contained
-        if(x[total-1] != mCycles){
-            x[total] = mCycles;
+    } else {
+        const int fullSampling = 40;
+        for(int i=0; i<fullSampling && i<mCycles; i++){
+            x[total] = i+1;
             total++;
         }
+        // down sampling if it's too long
+        if(mCycles>fullSampling) {
+            double pos = fullSampling;
+            while(true){
+                pos *= 1.05;
+                if(pos >= mCycles)
+                    break;
+                x[total] = (int)pos;
+                total++;
+            }
+            // make sure lsat one is contained
+            if(x[total-1] != mCycles){
+                x[total] = mCycles;
+                total++;
+            }
+        }
     }
     // four bases
     for (int b = 0; b<6; b++) {
@@ -468,7 +500,7 @@ void Stats::reportHtmlContents(ofstream& ofs, string filteringType, string readN
     json_str += "];\n";
     json_str += "var layout={title:'" + title + "', xaxis:{title:'cycles'";
     // use log plot if it's too long
-    if(mCycles>300) {
+    if(isLongRead()) {
         json_str += ",type:'log'";
     }
     json_str += "}, yaxis:{title:'base content ratios'}};\n";

diff --git a/src/stats.h b/src/stats.h
@@ -32,6 +32,7 @@ class Stats{
     void reportHtml(ofstream& ofs, string filteringType, string readName);
     void reportHtmlQuality(ofstream& ofs, string filteringType, string readName);
     void reportHtmlContents(ofstream& ofs, string filteringType, string readName);
+    bool isLongRead();
 
 public:
     static string list2string(double* list, int size);