# Statistical calculations with sets of genes

The following piece of code executes 20 million times each time the program is called, so I need a way to make this code as optimized as possible.

int WilcoxinRST::work(GeneSet originalGeneSet, vector<string> randomGenes) {
vector<string> geneIDs;
vector<bool> isOriginal;
vector<int> rank;
vector<double> value;
vector<int> score;
int genesPerSet = originalGeneSet.geneCount();
unsigned int totalGenes, tempScore;
/**
* Fill the first half of the vectors with original gene set data
*/
totalGenes = genesPerSet * 2;
for (int i = 0; i < genesPerSet; i++) {
geneIDs.push_back(originalGeneSet.getMemberGeneAt(i));
isOriginal.push_back(true);
value.push_back(fabs(expressionLevels.getValue(geneIDs[i], statType)));
}
/**
* Fill the second half with random data
*/
for (unsigned int i = genesPerSet; i < totalGenes; i++) {
geneIDs.push_back(randomGenes.at(i - genesPerSet));
isOriginal.push_back(false);
value.push_back(fabs(expressionLevels.getValue(geneIDs[i], statType)));
}
totalGenes = geneIDs.size();
/**
* calculate the scores
*/
if (statType == Fold_Change || statType == T_Statistic
|| statType == Z_Statistic) {
//      Higher value is a winner
for (unsigned int i = 0; i < totalGenes; i++) {
tempScore = 0;
if (!isOriginal[i]) {
for (int j = 0; j < genesPerSet; j++) {
if (value.at(i) > value.at(j)) {
tempScore++;
}
}

} else {
for (unsigned int j = genesPerSet; j < totalGenes; j++) {
if (value.at(i) > value.at(j)) {
tempScore++;
}
}

}

score.push_back(tempScore);
}

} else if (statType == FDR_PValue || statType == PValue) {
// Lower value is a winner
for (unsigned int i = 0; i < totalGenes; i++) {
tempScore = 0;
if (!isOriginal[i]) {
for (int j = 0; j < genesPerSet; j++) {
if (value.at(i) < value.at(j)) {
tempScore++;
}
}

} else {
for (unsigned int j = genesPerSet; j < totalGenes; j++) {
if (value.at(i) < value.at(j)) {
tempScore++;
}
}

}

score.push_back(tempScore);
}

} else {
cout << endl << "ERROR. Statistic type not defined." << endl;
}

/**
* calculate Ua, Ub and U
*/
int U_Original = 0, U_Random = 0, U_Final;
for (int j = 0; j < genesPerSet; j++) {
U_Original += score[j];
}
for (unsigned int j = genesPerSet; j < totalGenes; j++) {
U_Random += score[j];
}
U_Final = (U_Original < U_Random) ? U_Original : U_Random;

/**
* calculate z
*/
double Zn, Zd, Z;
Zn = U_Final - ((genesPerSet * genesPerSet) / 2);
Zd = sqrt(
(double) (((genesPerSet * genesPerSet
* (genesPerSet + genesPerSet + 1)))) / 12.0);
Z = Zn / Zd;

/**
* Return 0/1/2
* 2: p value < 0.01
* 1: 0.01 < p value < 0.05
* 0: p value > 0.05
*/
if (fabs(Z) > 2.303)
return 2;
else if (fabs(Z) > 1.605)
return 1;
else
return 0;
}

• Your best course of action is to use a profiler. Other classes e.g. GeneSet need looking at too. – parkydr Mar 28 '13 at 8:43

## 3 Answers

Your code have a complexity of O(N*N) [genesPerSet=N]. But using the fact that the order of the values is irrelevant for you we can order it in O(N•log(N)) and compute the “scores” in O(N). (With potentially can be thousands time faster).

Also, we have a total of N*N comparisons. Then U_Original + U_Random = N*N, meaning we don’t need to compute U_Random. Also your statistic Zn= Umin-N*N/2;[Umix=min(U_Original,U_Random)], when you only abs(Zn/Zd) is symmetric around N*N/2. We need only one algorithm.

1.- A first thing can be to take arguments by (const) reference :

int WilcoxinRST::work(const GeneSet &originalGeneSet, const vector<string> &randomGenes)


2.- You fill into vector geneIDs; but dont use it ? Why?

3.- You can iterate only 2 times.

4.- You keep the signal values (probe intensitat?) together in one vector and use another vector to signal what is each item – simple keep in two vector.

5.- You don’t need the score vector, only the summa.

6.- Why 20 millions times? I gess your are computing some “statistic” stability or BStrap. Probably you use the same originalGeneSet many time. I think you can in another question post the code that call this function to spare in making each time the vector of values and sorting.

7.- If you aproximatly know the espected size of the vectors you can reserve() it, add new item with push_back(), and access with .at() (somewhere write a try/catsh). But here we know exactly the numer of item, and we can set it in the cinstructor and then use [] with is the faster way.

Here is first the new O(N•log(N)) code.

What follows is a cleanup of your code but still O(N*N), fast but only by a constant factor.

Then the same code but mixed with yours original code and with more comments.

Please, debugg this and tell me how was.

#include<vector>
#include<algorithm>

int WilcoxinRST::work(const GeneSet &originalGeneSet , const vector<string>& randomGenes)
{
size_t genesPerSet = originalGeneSet.geneCount();
std::vector<double> valueOri(genesPerSet), valueRnd(genesPerSet);
/**
* Fill the valueOri vector with original gene set data, and valueRnd with random data
*/
for (size_t i = 0; i < genesPerSet; i++)
{
valueOri[i] = std::fabs(expressionLevels.getValue( originalGeneSet.getMemberGeneAt(i) , statType ));
valueRnd[i] = std::fabs(expressionLevels.getValue( randomGenes.at(i)                  , statType ));
}
std::sort(valueOri.begin(),valueOri.end());
std::sort(valueRnd.begin(),valueRnd.end());

/**
* calculate the scores Ua, Ub and U
*/
long U_Original=0 ;

if (statType == Fold_Change || statType == T_Statistic  || statType == Z_Statistic
statType == FDR_PValue  || statType == PValue )
{
//      Higher value is a winner
size_t j=0;
for (size_t i = 0; i < genesPerSet /*totalGenes*/; i++)      // i   -  2x
{
while(valueOri[i]  > valueRnd[j]) ++j ;
U_Original += j;
}
} else { cout << endl << "ERROR. Statistic type not defined." << endl;  }

/**
* calculate z
*/
double Zn, Zd, Z;
Zn = U_Original - ((genesPerSet * genesPerSet) / 2);
Zd = std::sqrt( (double) (((genesPerSet * genesPerSet* (genesPerSet + genesPerSet + 1)))) / 12.0);
Z = Zn / Zd;

/**
* Return 0/1/2
* 2: p value < 0.01
* 1: 0.01 < p value < 0.05
* 0: p value > 0.05
*/
if (std::fabs(Z) > 2.303)  return 2;
else if (std::fabs(Z) > 1.605)  return 1;
else                            return 0;
}


What follows is a cleanup of your code but still O(N*N), fast but only by a constant factor.

#include<vector>
using namespace std;
class GeneSet ;
class WilcoxinRST;

int WilcoxinRST::work(const GeneSet &originalGeneSet , const vector<string>& randomGenes)
{
size_t genesPerSet = originalGeneSet.geneCount();
vector<double> valueOri(genesPerSet), valueRnd(genesPerSet);
/**
* Fill the valueOri vector with original gene set data, and valueRnd with random data
*/
for (size_t i = 0; i < genesPerSet; i++)
{
valueOri[i]   =  fabs(expressionLevels.getValue( originalGeneSet.getMemberGeneAt(i) , statType ));
valueRnd[i]   =  fabs(expressionLevels.getValue( randomGenes.at(i)                  , statType ));
}
/**
* calculate the scores Ua, Ub and U
*/
long U_Original = 0, U_Random = 0, U_Final;

if (statType == Fold_Change || statType == T_Statistic  || statType == Z_Statistic)
{
//      Higher value is a winner
for (size_t i = 0; i < genesPerSet /*totalGenes*/; i++)      // i   -  2x
{   for (size_t j = 0; j < genesPerSet; j++)
{   U_Random  +=  (valueRnd[i]  > valueOri[j]);// i en 2 set=Rnd, j in 1 set=Ori. count how many Ori are less than this Rnd
U_Original+=  (valueOri[i]  > valueRnd[j]);// i in 1 set=Ori, j in 2 set=Rnd, count how many Rnd are less than this Ori
}
}
} else
if (statType == FDR_PValue || statType == PValue)
{
// Lower value is a winner
for (size_t i = 0; i < genesPerSet; i++)
{
for (size_t j = 0; j < genesPerSet; j++)
{   U_Random  +=  (valueRnd[i]  < valueOri[j]);// i en 2 set=Rnd, j in 1 set=Ori. count how many Ori are > than this Rnd
U_Original+=  (valueOri[i]  < valueRnd[j]);// i in 1 set=Ori, j in 2 set=Rnd, count how many Rnd are > than this Ori
}
}
} else { cout << endl << "ERROR. Statistic type not defined." << endl;  }

U_Final = (U_Original < U_Random) ? U_Original : U_Random;

/**
* calculate z
*/
double Zn, Zd, Z;
Zn = U_Final - ((genesPerSet * genesPerSet) / 2);
Zd = sqrt(
(double) (((genesPerSet * genesPerSet
* (genesPerSet + genesPerSet + 1)))) / 12.0);
Z = Zn / Zd;

/**
* Return 0/1/2
* 2: p value < 0.01
* 1: 0.01 < p value < 0.05
* 0: p value > 0.05
*/
if (fabs(Z) > 2.303)       return 2;
else if (fabs(Z) > 1.605)       return 1;
else                            return 0;
}


the same code but mixed with yours original code and with more comments.

int WilcoxinRST::work(const GeneSet &originalGeneSet , const vector<string>& randomGenes)
{
size_t genesPerSet = originalGeneSet.geneCount();
unsigned int totalGenes, tempScore;
totalGenes = genesPerSet * 2;

//vector<string> geneIDs;
//vector<bool>   isOriginal;
//vector<int>    rank;
vector<double> valueOri(genesPerSet), valueRnd(genesPerSet);
//vector<int>    score;
/**
* Fill the first half of the vectors with original gene set data
*/

for (size_t i = 0; i < genesPerSet; i++)
{
//geneIDs.push_back( originalGeneSet.getMemberGeneAt(i)  );
//isOriginal.push_back(true);

valueOri[i]   =  fabs(expressionLevels.getValue( originalGeneSet.getMemberGeneAt(i) , statType ));
valueRnd[i]   =  fabs(expressionLevels.getValue( randomGenes.at(i)                  , statType ));
}
/**
* Fill the second half with random data
*/
//for (unsigned int i = genesPerSet; i < totalGenes; i++) {
//  geneIDs.push_back(randomGenes.at(i - genesPerSet));
//  isOriginal.push_back(false);
//  value.push_back(fabs(expressionLevels.getValue(geneIDs[i], statType)));
//}
//totalGenes = geneIDs.size();
/**
* calculate the scores
*/
/**
* calculate Ua, Ub and U
*/
long U_Original = 0, U_Random = 0, U_Final;
//for (int j = 0; j < genesPerSet; j++) // j in 1 set=Ori. count how many Ori are less than this Rnd
//{
//  U_Original += score[j];
//}
//for (unsigned int j = genesPerSet; j < totalGenes; j++) // j in 2 set=Rnd, count how many Rnd are less than this Ori
//{
//  U_Random += score[j];
//}

if (statType == Fold_Change || statType == T_Statistic  || statType == Z_Statistic)
{
//      Higher value is a winner
for (size_t i = 0; i < genesPerSet /*totalGenes*/; i++)      // i   -  2x
{   //tempScore = 0;
//if (!isOriginal[i])  // i en 2 set=Rnd, j in 1 set=Ori. count how many Ori are less than this Rnd
for (size_t j = 0; j < genesPerSet; j++)
{   U_Random  +=  (valueRnd[i]  > valueOri[j]);// i en 2 set=Rnd, j in 1 set=Ori. count how many Ori are less than this Rnd
U_Original+=  (valueOri[i]  > valueRnd[j]);// i in 1 set=Ori, j in 2 set=Rnd, count how many Rnd are less than this Ori
}
//} else
//{
//  for (unsigned int j = genesPerSet; j < totalGenes; j++)  // i in 1 set=Ori, j in 2 set=Rnd, count how many Rnd are less than this Ori
//  {   if (value.at(i) > value.at(j)) {    tempScore++;        }
//  }

//}
//score.push_back(tempScore);
}

} else
if (statType == FDR_PValue || statType == PValue)
{
// Lower value is a winner
for (size_t i = 0; i < genesPerSet; i++)
{
for (size_t j = 0; j < genesPerSet; j++)
{   U_Random  +=  (valueRnd[i]  < valueOri[j]);// i en 2 set=Rnd, j in 1 set=Ori. count how many Ori are > than this Rnd
U_Original+=  (valueOri[i]  < valueRnd[j]);// i in 1 set=Ori, j in 2 set=Rnd, count how many Rnd are > than this Ori
}
//} else
//{
//  for (unsigned int j = genesPerSet; j < totalGenes; j++)  // i in 1 set=Ori, j in 2 set=Rnd, count how many Rnd are less than this Ori
//  {   if (value.at(i) > value.at(j)) {    tempScore++;        }
//  }

//}
//score.push_back(tempScore);
}

//for (unsigned int i = 0; i < totalGenes; i++)
//{   tempScore = 0;
//  if (!isOriginal[i])
//  {   for (int j = 0; j < genesPerSet; j++) {
//          if (value.at(i) < value.at(j)) {  // Rnd i < Ori j increm U_Random
//              tempScore++;
//          }
//      }

//  } else {
//      for (unsigned int j = genesPerSet; j < totalGenes; j++) { // Ori i < Rnd j. Increm U_Original
//          if (value.at(i) < value.at(j)) {
//              tempScore++;
//          }
//      }

//  }

//  score.push_back(tempScore);
//}

} else { cout << endl << "ERROR. Statistic type not defined." << endl;  }

U_Final = (U_Original < U_Random) ? U_Original : U_Random;

/**
* calculate z
*/
double Zn, Zd, Z;
Zn = U_Final - ((genesPerSet * genesPerSet) / 2);
Zd = sqrt(
(double) (((genesPerSet * genesPerSet
* (genesPerSet + genesPerSet + 1)))) / 12.0);
Z = Zn / Zd;

/**
* Return 0/1/2
* 2: p value < 0.01
* 1: 0.01 < p value < 0.05
* 0: p value > 0.05
*/
if (fabs(Z) > 2.303)
return 2;
else if (fabs(Z) > 1.605)
return 1;
else
return 0;
}

• stop using using namespace std; stackoverflow.com/a/1453605/14065 – Martin York Mar 29 '13 at 16:18
• @LokiAstari OK. Actualized. More suggestion, please! – qPCR4vir Mar 29 '13 at 19:05
• @Pranjal. I forget one initialization. Now edited. – qPCR4vir Apr 3 '13 at 20:03

It's somewhat difficult without more context, but there are definitely a few things I can see:

int WilcoxinRST::work(GeneSet originalGeneSet, vector<string> randomGenes)


This copies randomGenes every time it is called. You should pass by const &:

int WilcoxinRST::work(GeneSet originalGeneSet, const vector<string>& randomGenes)


With some of your vectors, you already know how much space they will need up front - so initialize it - this will save reallocations later.

vector<string> geneIDs(genesPerSet);
vector<double> value(totalGenes);


Don't use vector<bool>. It's not a vector, and it doesn't hold bools (to paraphrase from Scott Myers). Use a std::bitset instead. This is also stack allocated, and so may be (slightly) faster, although be careful if it gets too large - you may run out of stack space:

std::bitset<genesPerSet * 2> isOriginal;


rank seems to be unused, get rid of it.

at is range checked and will throw an exception if it is out of bounds. Utilize operator[] if you're really worried about speed (although this will likely be a tiny speed-up).

These are probably the most obvious things. There's not really enough context for any other suggestions.

• A template's argument has to be a compile time constant. Your declaration for isOriginal will not compile. – Happy Green Kid Naps Mar 28 '13 at 14:44

I took many Ideas from your posts and finally the code runs in 40 seconds from 3 hours before all this. Here is the code that I am using now

/*
* WilcoxinRST.cpp
*
*  Created on: 2013-04-07
*      Author: Pranjal
*/

#include "WilcoxinRST.h"
#include <iostream>
#include <math.h>
#include <stdio.h>
#include "data/ExpressionLevels.h"
#include "data/GeneSetDB.h"
#include <string.h>
#include <vector>
#include <algorithm>

std::WilcoxinRST::WilcoxinRST(ExpressionLevels* exp, GeneSetDB* gsdb,
StatType stat, int maxGenesPerGeneSet) {
this->expressionLevels = *exp;
this->maxGenesPerGeneSet = maxGenesPerGeneSet;
this->geneSetDB = *gsdb;
this->statType = stat;
this->opFile.open("results.txt");
this->opFile << "Gene Set ID\t#Genes\tp<0.01\tp<0.05\tp>0.05" << endl;
srand(time(NULL));
go();
this->opFile.close();
}

void std::WilcoxinRST::go() {
int x = 0, p005 = 0, p001 = 0, pRem = 0, geneCount = 0, genesWithoutValues;
GeneSet originalGeneSet;
vector<string> randomGenes;
vector<string> originalGenes;
vector<string> originalGenesRefined;
unsigned int totalGeneSets = geneSetDB.totalGeneSets();

// Runs for each gene set
for (unsigned int i = 0; i < totalGeneSets; i++) {
p001 = 0;
p005 = 0;
pRem = 0;
genesWithoutValues = 0;
// get the original gene set
originalGeneSet = geneSetDB.getGeneSetNumber(i);
originalGenes = originalGeneSet.getMemberGenes();

// refine original gene set based on the availability of expression values.
for (unsigned int k = 0; k < originalGenes.size(); ++k) {
if (expressionLevels.getValue(originalGenes[k], statType)
!= -999.999) {
// cout << expressionLevels.getValue(originalGenes[i], statType);
originalGenesRefined.push_back(originalGenes[k]);
} else
genesWithoutValues++;
}

originalGeneSet.setMemberGenes(originalGenesRefined);
originalGenesRefined.clear();
geneCount = originalGeneSet.geneCount();
// check if geneCount > max limit;
if (geneCount == 4/*>= this->maxGenesPerGeneSet*/) {
//geneCount = originalGeneSet.geneCount();
vector<double> valueOriginal(geneCount);
vector<double> valueRandom(geneCount);
vector<string> originalGenes = originalGeneSet.getMemberGenes();

cout << endl << originalGeneSet.getGeneSetName() << endl;
// Fill the valueOriginal vector with expression values
for (int j = 0; j < geneCount; j++) {
valueOriginal[j] = fabs(
expressionLevels.getValue(originalGenes[j], statType));
}

// Sort the valueOriginal vector.
std::sort(valueOriginal.begin(), valueOriginal.end());

// iterate for 1000 random gene sets.
for (int j = 0; j < 1000; j++) {
//randomGenes = expressionLevels.getRandomGenes(geneCount);
// x = work3(originalGeneSet, randomGenes);
valueRandom = expressionLevels.getRandomValues(geneCount,
statType, true);

std::sort(valueRandom.begin(), valueRandom.end());
x = work3(valueOriginal, valueRandom);
switch (x) {
case 1: {
p005++;
break;
}
case 2: {
p001++;
break;
}
case 0: {
pRem++;
break;
}
}
}

this->opFile << originalGeneSet.getGeneSetId() << "\t"
<< originalGeneSet.geneCount() << "\t" << p001 << "\t"
<< p005 << "\t" << pRem << "\t" << endl;
cout << endl << i + 1 << "\t" << expressionLevels.printTime();
cout.flush();
break;
}
this->opFile.flush();
}
}

int std::WilcoxinRST::work3(vector<double> originalGenes,
vector<double> randomGenes) {

/*
* @param geneSetSize contains the size of the original gene set
* @param uOriginal stores the final U for the original gene set
* @param uRandom stores the final U for the random gene set
* @param originalTemp will be set such that
*/
int geneSetSize = originalGenes.size();
int uOriginal = 0;
int uRandom = 0;
int originalTemp = 0;
int randomTemp = 0;

for (int i = 0; i < geneSetSize; ++i) {
originalTemp = 0;
while ((originalGenes[i] > randomGenes[randomTemp])
&& randomTemp < geneSetSize) {
randomTemp++;
}
while ((randomGenes[i] > originalGenes[originalTemp])
&& originalTemp < geneSetSize) {
originalTemp++;
}
uOriginal += randomTemp;
uRandom += originalTemp;
}

int uFinal = (uOriginal < uRandom) ? uOriginal : uRandom;

/**
* calculate z
*/
double Zn, Zd, Z;
Zn = (double) uFinal
- (((double) geneSetSize * (double) geneSetSize) / 2.0);
Zd = sqrt(
(double) (((geneSetSize * geneSetSize
* (geneSetSize + geneSetSize + 1)))) / 12.0);
Z = Zn / Zd;

//  cout << endl << "U Original:\t" << uOriginal << "\tU Random:\t" << uRandom
//          << "\tSum = " << uOriginal + uRandom << "\t Z = " << Z << endl;
/**
* Return 0/1/2
* 2: p value < 0.01
* 1: 0.01 < p value < 0.05
* 0: p value > 0.05
*/

cout << Z << "\t";
cout.flush();
if (fabs(Z) > 2.303)
return 2;
else if (fabs(Z) > 1.605)
return 1;
else
return 0;
}

• You're still copying the vector arguments every time work3 is called. Change the argument types to std::vector<double> const & instead. Also, I'm not sure how this code can compile, since vector is not a type, it's std::vector. Maybe one of your includes does a using namespace std? In which case, as mentioned by Loki in another comment: stackoverflow.com/a/1453605/14065. – Darhuuk Apr 25 '17 at 16:48