.gitignore

@@ -1,3 +1,6 @@
+# Emacs temporary files
+*~
+
 # Compiled Object files
 *.slo
 *.lo
@@ -11,3 +14,19 @@
 *.lai
 *.la
 *.a
+
+# Temporary test directories
+temp*
+
+# Files generated by cxxtest
+test_simtools.cpp
+runner
+runner.cpp
+
+# Compiled executables
+sim
+simtools
+g2i
+gtc
+gtc_process
+normalize_manifest

Egt.cpp

@@ -0,0 +1,287 @@
+
+// Egt.cpp
+//
+// Author: Iain Bancarz <ib5@sanger.ac.uk>
+//
+// Copyright (c) 2014 Genome Research Ltd.
+//
+// Redistribution and use in source and binary forms, with or without 
+// modification, are permitted provided that the following conditions are met:
+// 1. Redistributions of source code must retain the above copyright notice, 
+// this list of conditions and the following disclaimer.
+// 2. Redistributions in binary form must reproduce the above copyright 
+// notice, this list of conditions and the following disclaimer in the 
+// documentation and/or other materials provided with the distribution.
+// 3. Neither the name of Genome Research Ltd nor the names of the 
+// contributors may be used to endorse or promote products derived from 
+// software without specific prior written permission.
+//
+// THIS SOFTWARE IS PROVIDED BY THE AUTHOR ``AS IS'' AND ANY EXPRESS OR 
+// IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 
+// OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. 
+// IN NO EVENT SHALL GENOME RESEARCH LTD. BE LIABLE FOR ANY DIRECT, INDIRECT, 
+// INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, 
+// BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF 
+// USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 
+// THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 
+// (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF 
+// THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+//
+
+/*
+ * Egt is a class to parse Illumina EGT cluster files.
+ *
+ * This is a port of the EGT Python class in zCall to C++.
+ * There is no public documentation for the EGT binary file format.
+ * (However, the zCall EGT class is claimed to originate from Illumina.)
+ *
+ * Repository for zCall:  https://github.com/wtsi-npg/zCall
+ * 
+ * Native EGT format stores coordinates as (R, Theta) polar coordinates
+ * Python EGT class converts to Cartesian for storage in memory
+ * This class stores as polar
+ * (More efficient than converting and storing as Cartesian, since both are 
+needed for FCR output)
+ *
+ * IMPORTANT: Illumina defines its own polar coordinates:
+ * - Theta is rescaled, such that 1 unit = pi/2 radians
+ * - R is equal to x+y, *not* sqrt(x^2 + y^2)
+ * - R is a "Manhattan distance", not Euclidean distance
+ */
+
+#include <string> 
+#include <iostream>
+#include <sstream>
+#include "Egt.h"
+
+using namespace std;
+
+Egt::Egt(bool verbose)
+{
+  this->verbose = verbose;
+  GENOTYPES_PER_SNP = 3; 
+  PARAMS_PER_SNP = 12;
+  NUMERIC_BYTES = 4;
+  ENTRIES_IN_RECORD = 30;
+  BYTES_IN_RECORD = NUMERIC_BYTES * ENTRIES_IN_RECORD;
+  ENTRIES_TO_USE = 15;
+}
+
+void Egt::open(string filename)
+{
+  this->filename = filename;
+  ifstream file;
+
+  file.open(filename.c_str());
+  if (!file) {
+    cerr << "Can't open file: " << filename << endl << flush;
+    exit(1);
+  }
+  // read header data
+  readHeader(file);
+  readPreface(file);
+  if (verbose) {
+    printHeader();
+    printPreface();
+  }
+  // read cluster data
+  // expected cluster positions, in polar coordinates (R, Theta) for (AA,AB,BB)
+  // Order of params is: devRAA, devRAB, devRBB, meanRAA, meanRAB, meanRBB, devThetaAA, devThetaAB, devThetaBB, meanThetaAA, meanThetaAB, meanThetaBB
+  // Can't use a multidimensional array because we don't know snpTotal at compile time; instead use a pseudo-multidimensional array such that the (i,j)th value is (i*WIDTH + j)
+  counts = new int[GENOTYPES_PER_SNP*snpTotal]; // nAA, nAB, nBB
+  params = new float[PARAMS_PER_SNP*snpTotal];
+  char *block = new char[BYTES_IN_RECORD];
+  for (int i=0; i<snpTotal; i++) { 
+    file.read(block, BYTES_IN_RECORD);
+    int *ints = bytesToInts(block, 0, GENOTYPES_PER_SNP);
+    float *floats = bytesToFloats(block, GENOTYPES_PER_SNP, ENTRIES_IN_RECORD);
+    for (int j=0;j<GENOTYPES_PER_SNP;j++)
+      counts[i*GENOTYPES_PER_SNP + j] = ints[j];
+    for (int j=0;j<PARAMS_PER_SNP;j++)
+      params[i*PARAMS_PER_SNP + j] = floats[j];
+    delete floats;
+    delete ints;
+  }
+  delete block;
+  snpNames = new string[snpTotal];
+  readSNPNames(file, snpNames);
+  file.close();
+}
+
+void Egt::open(char *filename)
+{
+  string f = filename;
+  open(f);
+}
+
+
+int* Egt::bytesToInts(char block[], int start, int end) {
+  // convert a section of a byte array into ints
+  // start, end indices refer to positions in the array of ints (not bytes)
+  int *results = new int[end - start];
+  numericConverter converter;
+  for (int i=start;i<end;i++) {
+    for (int j=0;j<NUMERIC_BYTES;j++) {
+      converter.ncChar[j] = block[i*NUMERIC_BYTES + j];
+    }
+    results[i-start] = converter.ncInt;
+  }
+  return results;
+}
+
+float* Egt::bytesToFloats(char block[], int start, int end) {
+  // convert a section of a byte array into floats
+  // start, end indices refer to positions in the array of floats (not bytes)
+  float *results = new float[end - start];
+  numericConverter converter;
+  for (int i=start;i<end;i++) {
+    for (int j=0;j<NUMERIC_BYTES;j++) {
+      converter.ncChar[j] = block[i*NUMERIC_BYTES + j];
+    }
+    results[i-start] = converter.ncFloat;
+  }
+  return results;
+}
+
+void Egt::getClusters(long index, float snp_params[]) {
+  // return all (theta, R) params for given position in the SNP manifest
+  // includes s.d. and mean of (r, theta) for AA, AB, BB
+  int start = index*PARAMS_PER_SNP;
+  for (int j=0; j < PARAMS_PER_SNP; j++) {
+    snp_params[j] = this->params[start + j];
+  }
+}
+
+void Egt::getMeanR(long index, float means[]) {
+  // find mean polar radius for AA, AB, BB at given index
+  int start = index*PARAMS_PER_SNP + GENOTYPES_PER_SNP;
+  for (int j=0; j < GENOTYPES_PER_SNP; j++) {
+    means[j] = this->params[start + j];
+  }
+}
+
+void Egt::getMeanTheta(long index, float means[]) {
+  // find mean polar angle for AA, AB, BB at given index
+  int start = index*PARAMS_PER_SNP + 3*GENOTYPES_PER_SNP;
+  for (int j=0; j < GENOTYPES_PER_SNP; j++) {
+    means[j] = this->params[start + j];
+  }
+}
+
+numericConverter Egt::getNextConverter(ifstream &file) {
+  // convenience method to read the next few bytes into a union
+  // can then use the union for numeric conversion
+  char * buffer;
+  buffer = new char[NUMERIC_BYTES];
+  file.read(buffer, NUMERIC_BYTES);
+  numericConverter converter;
+  for (int i=0; i<NUMERIC_BYTES; i++) {
+    converter.ncChar[i] = buffer[i];
+  }
+  delete buffer;
+  return converter;
+}
+
+float Egt::readFloat(ifstream &file) {
+  float result;
+  numericConverter converter = getNextConverter(file);
+  result = converter.ncFloat;
+  return result;
+
+}
+
+void Egt::readHeader(ifstream &file) 
+{
+  // populate instance variables with header values
+  file.seekg(0); // set read position to zero, if not already there
+  fileVersion = readInteger(file);
+  gcVersion = readString(file, "GC version");
+  clusterVersion = readString(file, "Cluster version");
+  callVersion = readString(file, "Call version");
+  normalizationVersion = readString(file, "Normalization version");
+  dateCreated = readString(file, "Date created");
+  mode = file.get();
+  manifest = readString(file, "Manifest name");
+}
+
+int Egt::readInteger(ifstream &file)
+{
+  int result;
+  numericConverter converter = getNextConverter(file);
+  result = converter.ncInt;
+  return result; 
+}
+
+void Egt::readPreface(ifstream &file) {
+  // read the 'preface' from the body of an EGT file
+  // assumes file is positioned at start of the body
+  dataVersion = readInteger(file);
+  opa = readString(file, "OPA");
+  snpTotal = readInteger(file);
+}
+
+void Egt::readSNPNames(ifstream &file, string names[]) {
+  // read SNP names from an EGT file
+  // assumes file is positioned at end of cluster (mean, sd) data
+  int pos = file.tellg();
+  file.seekg(pos + 13 * snpTotal); // skip SNP quality scores
+  for (int i=0;i<snpTotal;i++) {
+    // skip genotype scores
+    // length of strings is unknown, so must read each one and discard it
+    readString(file, "genotype score");
+  }
+  for (int i=0;i<snpTotal;i++) {
+    stringstream sstream;
+    sstream << "SNP name at index " << i;
+    names[i] = readString(file, sstream.str());
+  }
+}
+
+string Egt::readString(ifstream &file, string name) {
+  // EGT string format is as follows: 
+  // - First byte is a *signed* char encoding the string length
+  // - Subsequent bytes contain the string
+  // Total bytes read is (length encoded in first byte)+1 -- at most 128
+  // Second argument is an (optional) identifier used for logging
+  int length = int(file.get()); // get a single byte
+  if (not file.good()) {
+    throw("Cannot read length from EGT file, file state is not good");
+  } else if (length < 0) {
+    throw("Illegal string length in EGT file");
+  } else if (length == 0) {
+    cerr << "Warning: String '" << name << "' has zero length." << endl;
+    return "";
+  }
+  char * buffer;
+  buffer = new char[length+1];
+  buffer[length] = '\0'; // ensure buffer ends with a null character
+  file.read(buffer, length);
+  if (not file.good()) {
+    stringstream sstream;
+    sstream << "Cannot read string '" << name << \
+      "' from EGT file, file state is not good";
+    throw(sstream.str());
+  }
+  string result = string(buffer);
+  delete buffer;
+  return result;
+}
+
+void Egt::printHeader() {
+  // convenience method to print file header
+  cout << "FILE_VERSION " << fileVersion << endl;
+  cout << "GC_VERSION " << gcVersion << endl;
+  cout << "CLUSTER_VERSION " << clusterVersion << endl;
+  cout << "CALL_VERSION " << callVersion << endl;
+  cout << "NORMALIZATION_VERSION " << normalizationVersion << endl;
+  cout << "DATE_CREATED " << dateCreated << endl;
+  cout << "MODE " << (int) mode << endl;
+  cout << "MANIFEST " << manifest << endl;
+}
+
+void Egt::printPreface() {
+  // convenience method to print preface of "main" data section
+  cout << "DATA_VERSION " << dataVersion << endl;
+  cout << "OPA " << opa << endl;
+  cout << "TOTAL_SNPS " << snpTotal << endl;
+}