---

TranslationModelMinimizable.cpp

Go to the documentation of this file.

#include "TranslationModelMinimizable.h"
#include <iostream>
#include <fstream>
#include <math.h>

const double TranslationModelMinimizable::TINY = 1.0e-16;

const int TranslationModelMinimizable::MAX_RNA_LENGTH = 4000;

const int TranslationModelMinimizable::MAX_RNA_MOLS = 2000;

const int TranslationModelMinimizable::MAX_STRING_LENGTH = 100;

TranslationModelMinimizable::TranslationModelMinimizable(LocateInputFileDirectory *locate, const char* infilename, int ribosome_length)
{
  rib_size = ribosome_length;

  // open file of genes and RNA and protein data
        ifstream input;
        input.open(infilename);

//  ifstream input(infilename,ios::nocreate,filebuf::sh_read);
  if (!input)
  {
          std::cout << "1 Error opening file " << infilename << std::endl;
    return;
  }

  input >> num_genes;  // number of genes to be used
  input >> num_data_sets;  // number of sets of fold change data
                           // (1 less than number of expt conditions)

  fm = new double*[num_data_sets];
  fp = new double*[num_data_sets];
  fmerr = new double*[num_data_sets];
  fperr = new double*[num_data_sets];

  for (int i=0; i<num_data_sets; i++)
  {
    fm[i] = new double[num_genes];
    fp[i] = new double[num_genes];
    fmerr[i] = new double[num_genes];
    fperr[i] = new double[num_genes];
  }

  geneticcode = new GeneticCode;

  char** seqfilename;
  seqfilename = new char*[MAX_RNA_MOLS]; // hold filenames of seq. files

  // get experimental data from file
  for (int i=0; i<num_genes; i++)
  {
    seqfilename[i] = new char[MAX_STRING_LENGTH];
    input >> seqfilename[i];

    for (int j=0; j<num_data_sets; j++)
    {
      input >> fm[j][i];
      input >> fmerr[j][i];
    }
    for (int j=0; j<num_data_sets; j++)
    {
      input >> fp[j][i];
      input >> fperr[j][i];
    }
  }

  RNAseq = new Codon*[num_genes];

  // for now, set all efficiencies to 1
  // These will correspond to enum Codon
  efficiency = new double[geneticcode->GetNumCodons()];
  for (int i=0; i<geneticcode->GetNumCodons(); i++)
    efficiency[i]=1;

  Codon* temp_seq_array;
  temp_seq_array = new Codon[MAX_RNA_LENGTH+1];

  length = new int[num_genes];

  // begin opening sequence files and reading
  for (int i=0; i<num_genes; i++)
  {
    // open a sequence file
    ifstream inputseq;
        char fileName[1000];

        sprintf(fileName, "%s%s", locate->Get("Sequences/"), seqfilename[i]);
        inputseq.open(fileName);

// cout << "opened inputseq file " << seqfilename[i] << endl;

    if (!inputseq)
    {
                std::cout << "2 Error opening file " << seqfilename[i] << std::endl;
      return;
    }

    eatwhite(inputseq);

    // RNA sequence to a temporary array
    int RNAlength=0;  // length of RNA being read in
    char codonstring[6];
    Codon tempcodon;
    Anticodon temptRNA = uuu;

    while (temptRNA!=stop)
    {
      inputseq.get(codonstring, 4);
      tempcodon = geneticcode->read_seq(codonstring);
      temp_seq_array[RNAlength]=tempcodon;

      temptRNA = geneticcode->get_anticodon(tempcodon);

      RNAlength++;
    }

    // transfer seqs to appropriate sequence array
    length[i]=RNAlength;
// cout << seqfilename[i] << '\t' << length[i] << endl;
    RNAseq[i] = new Codon[length[i]];

    for (int j=0; j<RNAlength; j++)
      RNAseq[i][j]=temp_seq_array[j];
  }
  // end opening sequence files and reading

  delete [] temp_seq_array;
  for (int i=0; i<num_genes; i++)
    delete [] seqfilename[i];
  delete [] seqfilename;  

  // compute number of residuals
  int num_resid=num_data_sets*num_genes;
  SetNResiduals(num_resid);

  // compute number of fit parameters
  num_params=0;
  num_params+=num_data_sets;  // ribosome concs., where Rn:=1
  num_params+=(geneticcode->GetNumAnticodons()-1)*(num_data_sets+1);  // tRNA elong. rates
  num_params+=num_genes-1;  // init. rates, where kinit0:=1
  num_params+=1;  // term. rate

/* fit parameters ordered as follows:
kelong[0]:  0,...,NUM_ANTICODONS-2
R0:  NUM_ANTICODONS-1
kelong[1]:  NUM_ANTICODONS,...,2*NUM_ANTICODONS-2
R1:  2*NUM_ANTICODONS-1
...
kelong[num_data_sets]:  num_data_sets*NUM_ANTICODONS,...,
                                (num_data_sets+1)*NUM_ANTICODONS-2
kterm:  (num_data_sets+1)*NUM_ANTICODONS-1
kinit[1,...,num_genes-1]:  (num_data_sets+1)*NUM_ANTICODONS,...,
                             (num_data_sets+1)*NUM_ANTICODONS+num_genes-2
*/

  // allocate space for things that will be computed
  Allocate(nResiduals);
}

TranslationModelMinimizable::~TranslationModelMinimizable()
{
  delete geneticcode;
  delete [] efficiency;
  delete [] length;

  for (int i=0; i<num_genes; i++)
    delete [] RNAseq[i];
  delete [] RNAseq;

  for (int i=0; i<num_data_sets; i++)
  {
    delete [] fm[i];
    delete [] fp[i];
    delete [] fmerr[i];
    delete [] fperr[i];
  }

  delete [] fm;
  delete [] fp;
  delete [] fmerr;
  delete [] fperr;
}

double TranslationModelMinimizable::EntropyShift(double T)
{
  return 0.0;
}

double TranslationModelMinimizable::F(double *parameters, double T)
{
  double objFuncVal = 0.0;
  objFuncVal = ComputeResiduals(parameters);
  return objFuncVal;
}

double TranslationModelMinimizable::F0(double *parameters, double T)
{
  return 0.0;
}

int TranslationModelMinimizable::GetNParameters()
{
  return num_params;
}

double TranslationModelMinimizable::ComputeResiduals(double *parameters)
{
  double* k = new double[MAX_RNA_LENGTH];  // holds translation rates
                          // k[0]=init rate
                          // k[1]=rate for 1st codon after AUG
                          // k[length-1]=term rate

  // compute prot prod rates for reference state
  double* P0 = new double[num_genes];
  for (int i=0; i<num_genes; i++)
  {
    // find initiation rates
    if (i==0)
      k[0]=1*parameters[geneticcode->GetNumAnticodons()-1];  // kinit[0]*R0
    else
      k[0]=parameters[((num_data_sets+1)*geneticcode->GetNumAnticodons()-1+i)]
                * parameters[geneticcode->GetNumAnticodons()-1];  // kinit[i]*R0


    // find elongation rates
    for (int j=1; j<length[i]-1; j++)
    {
      int tRNA = int(geneticcode->get_anticodon(RNAseq[i][j]));
      k[j]=parameters[tRNA];
    }
    // find termination rate
    k[length[i]-1]=parameters[(num_data_sets+1)*geneticcode->GetNumAnticodons()-1];

    // compute prot prod rate
    P0[i]=compute_current(k,length[i],rib_size);
// cout << i << '\t' << P0[i] << endl;
/*
if (i==0)
{
  cout << "state 0" << endl;
  for (int q=0; q<=length[i]-1; q++)
    cout << k[q] << endl;
}
*/
  }

  for (int i=0; i<num_data_sets; i++)
  {
    for (int j=0; j<num_genes; j++)
    {
      // work with gene j in dataset i

      // find initiation rate
      double R;  // ribosome conc
      if (i==num_data_sets-1)
        R=1;
      else
        R=parameters[(i+2)*geneticcode->GetNumAnticodons()-1];
      double kinit;
      if (j==0)
        kinit=1;
      else
        kinit=parameters[((num_data_sets+1)*geneticcode->GetNumAnticodons()-1+j)];
      k[0]=kinit*R;

      // find elongation rates
      for (int q=1; q<length[j]-1; q++)
      {
        int tRNA = int(geneticcode->get_anticodon(RNAseq[j][q]));
        k[q]=parameters[(i+1)*geneticcode->GetNumAnticodons()+tRNA];
      }
      // find termination rate
      k[length[j]-1]=parameters[(num_data_sets+1)*geneticcode->GetNumAnticodons()-1];

      // compute prot prod rate
      double Pij=compute_current(k,length[j],rib_size);
// cout << j << '\t' << Pij << endl;
/*
if (j==0)
{
  cout << "state " << i+1 << endl;
  for (int q=0; q<=length[j]-1; q++)
    cout << k[q] << endl;
}
*/
      // compute residual
      double fpij=Pij/P0[j]*fm[i][j];

      // compute residual for ln(fp)
      residuals[i*num_genes+j]=(log(fpij)-log(fp[i][j]))/sqrt(pow(fperr[i][j]/fp[i][j],2)+pow(fmerr[i][j]/fm[i][j],2));
    }
  }

  delete [] k;
  delete [] P0;

  return SumOfSquares();
}


/* fit parameters ordered as follows:
kelong[0]:  0,...,NUM_ANTICODONS-2
R0:  NUM_ANTICODONS-1
kelong[1]:  NUM_ANTICODONS,...,2*NUM_ANTICODONS-2
R1:  2*NUM_ANTICODONS-1
...
kelong[num_data_sets]:  num_data_sets*NUM_ANTICODONS,...,
                                (num_data_sets+1)*NUM_ANTICODONS-2
kterm:  (num_data_sets+1)*NUM_ANTICODONS-1
kinit[1,...,num_genes-1]:  (num_data_sets+1)*NUM_ANTICODONS,...,
                             (num_data_sets+1)*NUM_ANTICODONS+num_genes-2
*/

// N is length of RNA
// L is ribosome size
double TranslationModelMinimizable::compute_current(double* k, int N, int L)
{
/*
cout << endl << endl;
for (int i=0; i<N; i++)
  cout << k[i] << endl;
cout << endl << endl;
*/
  // find upper and lower bounds on current
  double kmin=k[0];
  double* kinv = new double[N];
  double* KL = new double[N];
  for (int i=0; i<N; i++)
  {
    if (k[i]<kmin)
      kmin=k[i];
    kinv[i]=1/k[i];
  }
  double Jmin=kmin/((1+sqrt((double)L))*(1+sqrt((double)L)));
  for (int i=L-1; i<N; i++)
  {
    double sum=0;
    for (int j=0; j<L; j++)
      sum+=kinv[i-j];
    KL[i]=1/sum;
  }
  double Jmax=KL[L-1];
  for (int i=L; i<N; i++)
    if (KL[i]<Jmax)
      Jmax=KL[i];

  // allocate space for ribosome densities
  double* n = new double[N];

  // find current
  double Jub=Jmax;
  double Jlb=Jmin;
//cout << "Jub=" << Jub << " Jlb=" << Jlb << endl;
  // start bisection
  int done=0;
  while (Jub-Jlb>TINY && done==0)
  {
    double Jnew=(Jlb+Jub)/2;

    update_dens(n, k, Jnew, N-1, L);

    int bisected=0;
    // check for negative ribosome densities
    int i=N-1;
    while (i>=0 && bisected==0)
    {
      if (n[i]<0)  // negative rib dens--Jnew is too large
      {
        Jub=Jnew;
//cout << "neg n; new Jub=" << Jub << endl;
        bisected=1;
      }
      i--;
    }
    if (bisected==0)
    {
      if (n[0]>k[0])
      {
        Jub=Jnew;
//cout << "pred kinit too big; new Jub=" << Jub << endl;
        bisected=1;
      }
      else if (n[0]<k[0])
      {
        Jlb=Jnew;
//cout << "pred kinit too low; new Jlb=" << Jlb << endl;
        bisected=1;
      }
      else
      {
        done=1;
        Jlb=Jnew;
        Jub=Jnew;
//cout << "kinit fine; J="  << Jlb << endl;
      }
    }
  }
  // end bisecting current

  delete [] n;
  delete [] kinv;
  delete [] KL;

  return ((Jlb+Jub)/2);
}

// N is coordinate of final ribosome density, so array length is N+1
void TranslationModelMinimizable::update_dens(double* n, double* k, double J, int N, 
int L)
{
  for (int i=N; i>=N-L+1; i--)
    n[i]=J/k[i];
  for (int i=N-L; i>=1; i--)
  {
    double sum=0;
    for (int j=1; j<=L; j++)
      sum+=n[i+j];
    n[i]=(J/k[i])*(1-sum+n[i+L])/(1-sum);
  }

  // store predicted kinit in n[0]
  double sum=0;
  for (int j=1; j<=L; j++)
    sum+=n[j];
  n[0]=J/(1-sum);
}

// from Stroustrup's book
ifstream& TranslationModelMinimizable::eatwhite(ifstream& is)
{
  char c;
  while (is.get(c))
  {
    if (!isspace(c)) // is c a whitespace character?
    {
      is.putback(c);  // put c back into the input buffer
      break;
    }
  }
  return is;
}

---