Fixed numerical discrepancies for consistency with versions <= 1.8; added normalization test

Author: iainrb

Fcr.cpp

@@ -169,13 +169,14 @@ void FcrWriter::write(Egt *egt, Manifest *manifest, ostream *outStream,
     else sampleName = gtc->sampleName;
     for (unsigned int j = 0; j < manifest->snps.size(); j++) {
       string snpName = manifest->snps[j].name;
-      double x_raw = gtc->xRawIntensity[j];
-      double y_raw = gtc->yRawIntensity[j];
+      unsigned short x_raw = gtc->xRawIntensity[j];
+      unsigned short y_raw = gtc->yRawIntensity[j];
       float score = gtc->scores[j];
       double x_norm;
       double y_norm;
       unsigned int norm_id = manifest->normIdMap[manifest->snps[j].normId];
-      gtc->normalizeIntensity(x_raw, y_raw, x_norm, y_norm, norm_id);
+      XFormClass xf = gtc->XForm[norm_id];
+      xf.normalize(x_raw, y_raw, x_norm, y_norm);
       // correction of negative intensities, for consistency with GenomeStudio
       if (x_norm < epsilon) { x_norm = 0.0; }
       if (y_norm < epsilon) { y_norm = 0.0; }

Gtc.cpp

@@ -34,15 +34,86 @@
 //
 #include "Gtc.h"
 #include <cstring>
+#include <cstdio>
 #include <iostream>
 #include <fstream>
 #include <sstream>
 #include <vector>
 
 using namespace std;
 
-XFormClass::XFormClass(void)
+XFormClass::XFormClass(int version, float xOffset, float yOffset, 
+                       float xScale, float yScale, float shear, float theta)
 {
+  this->version = version;
+  this->xOffset = xOffset;
+  this->yOffset = yOffset;
+  this->xScale = xScale;
+  this->yScale = yScale;
+  this->shear = shear;
+  this->theta = theta;  
+}
+
+void XFormClass::normalize(unsigned short x_raw, unsigned short y_raw, 
+                           double &x_norm, double &y_norm)
+{
+  // This is the intensity normalization calculation, according to Illumina
+  double tempx = x_raw - xOffset;
+  double tempy = y_raw - yOffset;
+  double cos_theta = cos(theta);
+  double sin_theta = sin(theta);
+  double tempx2 = cos_theta * tempx + sin_theta * tempy;
+  double tempy2 = -sin_theta * tempx + cos_theta * tempy;
+  double tempx3 = tempx2 - shear * tempy2;
+  double tempy3 = tempy2;
+  x_norm = tempx3 / xScale;
+  y_norm = tempy3 / yScale;
+}
+
+string XFormClass::toString(void)
+{
+  int bufSize = 1000;
+  char buffer[bufSize];
+  string v, xoff, yoff, xsca, ysca, shr, th;
+  int n;
+  const char *float_format = "%.6f";
+  try { // attempt to convert other data types to strings
+      n = snprintf(buffer, bufSize, "%d", version);
+      if (n >= bufSize) { throw 1; }
+      v = string(buffer, n);
+
+      n = snprintf(buffer, bufSize, float_format, xOffset);
+      if (n >= bufSize) { throw 1; }
+      xoff = string(buffer, n);
+
+      n = snprintf(buffer, bufSize, float_format, yOffset);
+      if (n >= bufSize) { throw 1; }
+      yoff = string(buffer, n);
+
+      n = snprintf(buffer, bufSize, float_format, xScale);
+      if (n >= bufSize) { throw 1; }
+      xsca = string(buffer, n);
+
+      n = snprintf(buffer, bufSize, float_format, yScale);
+      if (n >= bufSize) { throw 1; }
+      ysca = string(buffer, n);
+
+      n = snprintf(buffer, bufSize, float_format, shear);
+      if (n >= bufSize) { throw 1; }
+      shr = string(buffer, n);
+
+      n = snprintf(buffer, bufSize, float_format, theta);
+      if (n >= bufSize) { throw 1; }
+      th = string(buffer, n);
+  } catch (int param) {
+    cerr << "Error: Unable to convert variable from GTC XForm to "
+         << "string format. Bad input or corrupt GTC file?" << endl;
+    throw;
+  }
+  string xf = "Version: "+v+"\nXOffset: "+xoff+"\nYOffset: "+yoff\
+    +"\nXScale: "+xsca+"\nYScale: "+ysca+"\nShear: "+shr+"\nTheta: "+th;
+  return xf;
+
 }
 
 Gtc::Gtc(void) 
@@ -219,48 +290,38 @@ string Gtc::json_dump(void)
 	return s.str();
 }
 
-void Gtc::normalizeIntensity(double x_raw, double y_raw,
-                             double &x_norm, double &y_norm,
-                             unsigned int norm_id) {
-  // This is the normalization calculation, according to Illumina
-  XFormClass *XF = &(this->XForm[norm_id]);
-  double tempx = x_raw - XF->xOffset;
-  double tempy = y_raw - XF->yOffset;
-  double cos_theta = cos(XF->theta);
-  double sin_theta = sin(XF->theta);
-  double tempx2 = cos_theta * tempx + sin_theta * tempy;
-  double tempy2 = -sin_theta * tempx + cos_theta * tempy;
-  double tempx3 = tempx2 - XF->shear * tempy2;
-  double tempy3 = tempy2;
-  x_norm = tempx3 / XF->xScale;
-  y_norm = tempy3 / XF->yScale;
-}
-
 void Gtc::readXForm(ifstream &file, int offset)
 {
 	int arrayLen;
+        int total_reserved = 6; // 'reserved' floats after the xform fields
 	float reserved;
 
+        int version;
+        float xOffset;
+        float yOffset;
+        float xScale;
+        float yScale;
+        float shear;
+        float theta;
+
 	ios::pos_type pos = file.tellg();
 	file.seekg(offset);
 	file.read((char*)&arrayLen,4);
 
-	for (int n=0; n<arrayLen; n++) {
-		XFormClass X;
-		file.read((char*)&X.version, 4);
-		file.read((char*)&X.xOffset, 4);
-		file.read((char*)&X.yOffset, 4);
-		file.read((char*)&X.xScale, 4);
-		file.read((char*)&X.yScale, 4);
-		file.read((char*)&X.shear, 4);
-		file.read((char*)&X.theta, 4);
-		file.read((char*)&reserved, 4);
-		file.read((char*)&reserved, 4);
-		file.read((char*)&reserved, 4);
-		file.read((char*)&reserved, 4);
-		file.read((char*)&reserved, 4);
-		file.read((char*)&reserved, 4);
-		this->XForm.push_back(X);
+	for (int i=0; i<arrayLen; i++) {
+            file.read((char*)&version, 4);
+            file.read((char*)&xOffset, 4);
+            file.read((char*)&yOffset, 4);
+            file.read((char*)&xScale, 4);
+            file.read((char*)&yScale, 4);
+            file.read((char*)&shear, 4);
+            file.read((char*)&theta, 4);
+            XFormClass X = XFormClass(version, xOffset, yOffset, xScale, 
+                                      yScale, shear, theta);
+            this->XForm.push_back(X);
+            for (int j=0; j<total_reserved; j++) {
+              file.read((char*)&reserved, 4);
+            }
 	}
 	file.seekg(pos);
 }

Gtc.h

@@ -45,14 +45,21 @@ using namespace std;
 
 class XFormClass {
 public:
-	XFormClass();
-	int version;
-	float xOffset;
-	float yOffset;
-	float xScale;
-	float yScale;
-	float shear;
-	float theta;
+  XFormClass(int version, float xOffset, float yOffset, float xScale, 
+             float yScale, float shear, float theta);
+  int version;
+  float xOffset;
+  float yOffset;
+  float xScale;
+  float yScale;
+  float shear;
+  float theta;
+  
+  void normalize(unsigned short x_raw, unsigned short y_raw, 
+                 double &x_norm, double &y_norm);
+
+  string toString();
+
 };
 
 class BaseCallClass {
@@ -81,8 +88,6 @@ class Gtc {
 	double passRate(double cutoff);
 	double correctedPassRate(double cutoff);
 
-        void normalizeIntensity(double x_raw, double y_raw, double &x_norm, double &y_norm, unsigned int norm_id);
-
 	string errorMsg;
 	string filename;
 	int version;			// file version (expected to be 3

commands.cpp

@@ -218,11 +218,12 @@ void Commander::commandCreate(string infile, string outfile, bool normalize, str
       double xn;
       double yn;
       int idx = snp->index - 1;   // index is zero based in arrays, but starts from 1 in the map file
-      double x_raw = gtc->xRawIntensity[idx];
-      double y_raw = gtc->yRawIntensity[idx];
+      unsigned short x_raw = gtc->xRawIntensity[idx];
+      unsigned short y_raw = gtc->yRawIntensity[idx];
       unsigned int norm_id = manifest->normIdMap[snp->normId];
       if (normalize) {
-        gtc->normalizeIntensity(x_raw, y_raw, xn, yn, norm_id);
+        XFormClass xf = gtc->XForm[norm_id];
+        xf.normalize(x_raw, y_raw, xn, yn);
       } else {
 	xn = gtc->xRawIntensity[idx];
 	yn = gtc->yRawIntensity[idx];

test_simtools.h

@@ -232,7 +232,7 @@ class FcrTest : public TestBase
   }
 };
 
-class NormalizeTest : public TestBase
+class ManifestTest : public TestBase
 {
  public:
 
@@ -257,7 +257,7 @@ class NormalizeTest : public TestBase
     TS_TRACE("Finished manifest test");
   }
 
-  void testNormalize(void)
+  void testManifestNormalize(void)
   {
     // compare normalized output with reference file
     string infile = "data/mock.bpm.csv";
@@ -489,5 +489,44 @@ class Win2UnixTest : public TestBase
 
 };
 
+class XFormTest : public TestBase
+{
+  // tests the intensity normalization for writing .sim files
+
+ public:
+
+  void testXForm(void) {
+
+    int version= 1;
+    float xOffset = 150.0;
+    float yOffset = 90.0;
+    float xScale = 12000.0;
+    float yScale = 8000.0;
+    float shear = 0.02;
+    float theta = 0.01;
+
+    XFormClass *X;
+    TS_ASSERT_THROWS_NOTHING(X = new XFormClass(version, xOffset, yOffset, 
+                                                xScale, yScale, shear, theta)
+                             );
+    TS_TRACE("Created XFormClass object");
+
+    unsigned short x_raw = 400;
+    unsigned short y_raw = 200;
+    double x_norm;
+    double y_norm;
+    double x_norm_expected = 0.020744799066257; // 14 significant digits
+    double y_norm_expected = 0.013436817602487;
+    TS_ASSERT_THROWS_NOTHING(X->normalize(x_raw, y_raw, x_norm, y_norm));
+    TS_TRACE("Ran intensity normalization");
+
+    double epsilon = 1e-14;
+    TS_ASSERT(abs(x_norm - x_norm_expected) < epsilon);
+    TS_ASSERT(abs(y_norm - y_norm_expected) < epsilon);
+    TS_TRACE("Normalized intensities checked against expected values");
+  }
+
+};
+
 // Putting TestSuite classes in separate files appears not to work
 // See Cxx manual section 4.4; possible issue with compiler options