erg.txt

mfold

//calculate the energy of stacked base pairs
int erg1(int i, int j, int ip, int jp, structure *ct, datatable *data)
{

  int energy;

  if ((i==(ct->numofbases))||(j==((ct->numofbases)+1))) {
    //this is not allowed because n and n+1 are not cavalently attached
    energy = infinity;
  }
  else {
    energy = data->stack[(ct->numseq[i])][(ct->numseq[j])]
      [(ct->numseq[ip])][(ct->numseq[jp])]+data->eparam[1];
  }
  return energy;
}


//calculate the energy of a bulge/internal loop
//where i is paired to j; ip is paired to jp; ip > i; j > jp
int erg2(int i, int j, int ip, int jp, structure *ct, datatable *data,
	 int a, int b)
{
  int energy = 0, size, size1, size2, loginc, lopsid, energy2; //tlink, count, key, e[4]
  /* size, size1, size2 = size of a loop
     energy = energy calculated
     loginc = the value of a log used in large hairpin loops
  */
  if (((i<=(ct->numofbases))&&(ip>(ct->numofbases)))||((
							jp<=(ct->numofbases))&&(j>(ct->numofbases)))) {
    //A loop cannot contain the ends of the sequence
    energy=infinity;
    return energy;
  }


  size1 = ip-i-1;
  size2 = j - jp - 1;

  if ((a!=0)||(b!=0)) {
    if ((a==1)||(b==1)) return infinity; //the loop contains a nuc that
    //should be double stranded
    else if ((a==5)) {
      //the loop is actually between two strands (ie: intermolecular)

      if (size2>1) {//free energy is that of two terminal mismatches
	//and the intermolecular initiation
	energy = data->init + data->tstack[ct->numseq[i]][ct->numseq[j]]
	  [ct->numseq[i+1]][ct->numseq[j-1]] +
	  data->tstack[ct->numseq[jp]][ct->numseq[ip]]
	  [ct->numseq[jp+1]][ct->numseq[ip-1]];
      }
      else if (size2==1) {//find the best terminal mismatch and terminal
	//stack free energies combination

	energy = data->init + data->tstack[ct->numseq[i]][ct->numseq[j]]
	  [ct->numseq[i+1]][ct->numseq[j-1]] +
	  erg4 (jp, ip, ip-1, 2, ct, data, false)+penalty(jp, ip, ct, data);
	energy2 = data->init + data->tstack[ct->numseq[jp]][ct->numseq[ip]]
	  [ct->numseq[jp+1]][ct->numseq[ip-1]] +
	  erg4 (i, j, i+1, 1, ct, data, false)+penalty(i, j, ct, data);
	energy = min (energy, energy2);
	//if ((ct->numseq[i+1]!=5)&&(ct->numseq[ip-1]!=5)) {
	//now consider if coaxial stacking is better:
	energy2 = data->init + data->tstackcoax[ct->numseq[jp]]
	  [ct->numseq[ip]][ct->numseq[jp+1]][ct->numseq[ip-1]]
	  + data->coaxstack[ct->numseq[jp+1]][ct->numseq[ip-1]]
	  [ct->numseq[j]][ct->numseq[i]]+penalty(i, j, ct, data)+penalty(jp, ip, ct, data);
	energy = min(energy, energy2);
	energy2 = data->init + data->tstackcoax[ct->numseq[jp]]
	  [ct->numseq[ip]][ct->numseq[j-1]][ct->numseq[ip-1]]
	  + data->coaxstack[ct->numseq[j-1]][ct->numseq[ip-1]]
	  [ct->numseq[j]][ct->numseq[i]]+penalty(i, j, ct, data)+penalty(jp, ip, ct, data);
	energy = min(energy, energy2);
	//}
      }
      else if (size2==0) {//just have dangling ends or flush stacking
	energy = data->init + erg4 (jp, ip, ip-1, 2, ct, data, false) +
	  erg4 (i, j, i+1, 1, ct, data, false)+penalty(i, j, ct, data)+penalty(jp, ip, ct, data);
	energy2 = data->init + data->coax[ct->numseq[ip]][ct->numseq[jp]]
	  [ct->numseq[j]][ct->numseq[i]]+penalty(i, j, ct, data)+penalty(jp, ip, ct, data);
	energy = min(energy, energy2);
      }


      return energy;
    }
    else if (b==5) {
      //the loop is actually between two strands (ie: intermolecular)

      if (size1>1) {//free energy is that of two terminal mismatches
	//and the intermolecular initiation
	energy = data->init + data->tstack[ct->numseq[i]][ct->numseq[j]]
	  [ct->numseq[i+1]][ct->numseq[j-1]] +
	  data->tstack[ct->numseq[jp]][ct->numseq[ip]]
	  [ct->numseq[jp+1]][ct->numseq[ip-1]];
      }
      else if (size1==1) {//find the best terminal mismatch and terminal
	//stack free energies combination

	energy = data->init + data->tstack[ct->numseq[i]][ct->numseq[j]]
	  [ct->numseq[i+1]][ct->numseq[j-1]] +
	  erg4 (ip, jp, jp+1, 1, ct, data, false)+penalty(ip, jp, ct, data);
	energy2 = data->init + data->tstack[ct->numseq[jp]][ct->numseq[ip]]
	  [ct->numseq[jp+1]][ct->numseq[ip-1]] +
	  erg4 (i, j, j-1, 2, ct, data, false)+penalty(i, j, ct, data);

	energy = min (energy, energy2);
	//if ((ct->numseq[i+1]!=5)&&(ct->numseq[ip-1]!=5)) {
	//now consider if coaxial stacking is better:
	energy2 = data->init + data->tstackcoax[ct->numseq[i]]
	  [ct->numseq[j]][ct->numseq[i+1]][ct->numseq[j-1]]
	  + data->coaxstack[ct->numseq[i+1]][ct->numseq[j-1]]
	  [ct->numseq[ip]][ct->numseq[jp]]+penalty(i, j, ct, data)+penalty(jp, ip, ct, data);
	energy = min(energy, energy2);
	energy2 = data->init + data->tstackcoax[ct->numseq[i]]
	  [ct->numseq[j]][ct->numseq[ip-1]][ct->numseq[j-1]]
	  + data->coaxstack[ct->numseq[ip-1]][ct->numseq[j-1]]
	  [ct->numseq[ip]][ct->numseq[jp]]+penalty(i, j, ct, data)+penalty(jp, ip, ct, data);
	energy = min(energy, energy2);
	//}
      }
      else if (size1==0) {//just have dangling ends or flush stacking
	energy = data->init + erg4 (jp, ip, jp+1, 1, ct, data, false) +
	  erg4 (i, j, j-1, 2, ct, data, false)+penalty(i, j, ct, data)+penalty(jp, ip, ct, data);
	energy2 = data->init + data->coax[ct->numseq[j]][ct->numseq[i]]
	  [ct->numseq[ip]][ct->numseq[j]]+penalty(i, j, ct, data)+penalty(jp, ip, ct, data);
	energy = min(energy, energy2);
      }


      return energy;

    }
  }


  //a typical internal or bulge loop:
  size1 = ip-i-1;
  size2 = j - jp - 1;
  if (size1==0||size2==0) {//bulge loop


    size = size1+size2;
    if (size==1) {
      energy = data->stack[ct->numseq[i]][ct->numseq[j]]
	[ct->numseq[ip]][ct->numseq[jp]]
	+ data->bulge[size] + data->eparam[2];
    }
    else if (size>30) {

      loginc = int((data->prelog)*log(double ((size)/30.0)));
      energy = data->bulge[30] + loginc + data->eparam[2];
      energy = energy + penalty(i, j, ct, data) + penalty(jp, ip, ct, data);

    }
    else {
      energy = data->bulge[size] + data->eparam[2];
      energy = energy + penalty(i, j, ct, data) + penalty(jp, ip, ct, data);
    }
  }
  else {//internal loop
    size = size1 + size2;
    lopsid = abs(size1-size2);

    if (size>30) {

      loginc = int((data->prelog)*log((double ((size))/30.0)));
      if ((size1==1||size2==1)&&data->gail) {
	energy = data->tstki[ct->numseq[i]][ct->numseq[j]][1][1] +
	  data->tstki[ct->numseq[jp]][ct->numseq[ip]][1][1] +
	  data->inter[30] + loginc + data->eparam[3] +
	  min(data->maxpen, (lopsid*
			    data->poppen[min(2, min(size1, size2))]));

      }

      else {
	energy = data->tstki[ct->numseq[i]][ct->numseq[j]]
	  [ct->numseq[i+1]][ct->numseq[j-1]] +
	  data->tstki[ct->numseq[jp]][ct->numseq[ip]]
	  [ct->numseq[jp+1]][ct->numseq[ip-1]] +
	  data->inter[30] + loginc + data->eparam[3] +
	  min(data->maxpen, (lopsid*
			    data->poppen[min(2, min(size1, size2))]));
      }
    }
    else if ((size1==2)&&(size2==2))//2x2 internal loop
      energy = data->iloop22[ct->numseq[i]][ct->numseq[ip]]
	[ct->numseq[j]][ct->numseq[jp]]
	[ct->numseq[i+1]][ct->numseq[i+2]]
	[ct->numseq[j-1]][ct->numseq[j-2]];


    else if ((size1==1)&&(size2==2)) {//2x1 internal loop
      energy = data->iloop21[ct->numseq[i]][ct->numseq[j]][ct->numseq[i+1]]
	[ct->numseq[j-1]][ct->numseq[jp+1]][ct->numseq[ip]][ct->numseq[jp]];


    }
    else if ((size1==2)&&(size2==1)) {//1x2 internal loop
      energy = data->iloop21[ct->numseq[jp]][ct->numseq[ip]][ct->numseq[jp+1]]
	[ct->numseq[ip-1]][ct->numseq[i+1]][ct->numseq[j]][ct->numseq[i]];
      //cout << "2x1 internal loop energy = "<< energy <<"\n";
    }

    else if (size==2) //a single mismatch
      //energy = data->stack[ct->numseq[i]][ct->numseq[j]]
      //	[ct->numseq[i+1]][ct->numseq[j-1]] +
      //	data->stack[ct->numseq[jp]][ct->numseq[ip]]
      //	[ct->numseq[jp+1]][ct->numseq[ip-1]];
      energy = data->iloop11[ct->numseq[i]][ct->numseq[i+1]][ct->numseq[ip]]
	[ct->numseq[j]][ct->numseq[j-1]][ct->numseq[jp]];
    else if ((size1==1||size2==1)&&data->gail) { //this loop is lopsided
      //note, we treat this case as if we had a loop composed of all As
      //if and only if the gail rule is set to 1 in miscloop.dat
      energy = data->tstki[ct->numseq[i]][ct->numseq[j]][1][1] +
	data->tstki[ct->numseq[jp]][ct->numseq[ip]][1][1] +
	data->inter[size] + data->eparam[3] +
	min(data->maxpen, (lopsid*
			  data->poppen[min(2, min(size1, size2))]));
    }
    /*else if (size2==1) { //this loop is lopsided - one side has a terminal
    //mismatch and a dangle - the other side just a terminal mismatch

    energy = data->stack[ct->numseq[i]][ct->numseq[j]]
    [ct->numseq[i+1]][ct->numseq[j-1]] +
    data->dangle[ct->numseq[jp]][ct->numseq[ip]]
    [ct->numseq[ip-1]][2]+penalty(ip, jp, ct, data);
    energy = min(energy, data->stack[ct->numseq[jp]][ct->numseq[ip]]
    [ct->numseq[jp+1]][ct->numseq[ip-1]]+data->dangle[ct->numseq[i]]
    [ct->numseq[j]][ct->numseq[i+1]][1])+penalty(i, j, ct, data);
    energy = energy + data->inter[size] + data->eparam[3] +
    min(data->maxpen, (lopsid*
    data->poppen[min(2, min(size1, size2))]));

    } */
    else {
      //debug:
      //if (i==54&&j==96&&ip==57&&jp==89)
      //{
      //	data->poppen[1] = data->poppen[2];
      //}


      energy = data->tstki[ct->numseq[i]][ct->numseq[j]][ct->numseq[i+1]][ct->numseq[j-1]] +
	data->tstki[ct->numseq[jp]][ct->numseq[ip]][ct->numseq[jp+1]][ct->numseq[ip-1]] +
	data->inter[size] + data->eparam[3] +
	min(data->maxpen, (lopsid*data->poppen[min(2, min(size1, size2))]));
    }
  }
  /////
  //      cout << "erg2: "<< energy<<"\n";
  //            cout << "i: " <<i<<"\n";
  //      cout << "j: " <<j<<"\n";
  //      cout << "ip: "<<ip<<"\n";
  //      cout << "jp: "<<jp<<"\n";
  return energy;
}


//calculate the energy of a hairpin loop:
int erg3(int i, int j, structure *ct, datatable *data, int dbl)
{
  int energy, size, loginc, tlink, count, key, k;
  /* size, size1, size2 = size of a loop
     energy = energy calculated
     loginc = the value of a log used in large hairpin loops
  */


  if (dbl==1) return infinity; //the loop contains a base that should be
  //double stranded

  else if (dbl==5) {//intermolecular interaction
    //intermolecular "hairpin" free energy is that of intermolecular
    //	initiation plus the stacked mismatch

    energy = data->init + data->tstack[ct->numseq[i]][ct->numseq[j]]
      [ct->numseq[i+1]][ct->numseq[j-1]];
    return energy;
  }


  if ((i<=(ct->numofbases))&&(j>(ct->numofbases))) {
    //A hairpin cannot contain the ends of the sequence
    energy = infinity;
    return energy;
  }

  size = j-i-1;



  if (size>30) {
    //cout << "erg3:  i = "<<i<<"   j = "<<j<<"   "<<log(double ((size)/30.0))<<"\n";
    loginc = int((data->prelog)*log((double ((size))/30.0)));



    energy = data->tstkh[ct->numseq[i]][ct->numseq[j]]
      [ct->numseq[i+1]][ct->numseq[j-1]]
      + data->hairpin[30]+loginc+data->eparam[4];
  }
  else if (size<3) {
    energy = data->hairpin[size] + data->eparam[4];
    if (ct->numseq[i]==4||ct->numseq[j]==4) energy = energy+6;
  }
  else if (size==4) {
    tlink = 0;
    key = (ct->numseq[j])*3125 + (ct->numseq[i+4])*625 +
      (ct->numseq[i+3])*125 + (ct->numseq[i+2])*25+(ct->numseq[i+1])*5+(ct->numseq[i]);
    for (count=1; count<=data->numoftloops&&tlink==0; count++) {
      if (key==data->tloop[count][0]) tlink = data->tloop[count][1];
    }
    energy = data->tstkh[ct->numseq[i]][ct->numseq[j]]
      [ct->numseq[i+1]][ct->numseq[j-1]]
      + data->hairpin[size] + data->eparam[4] + tlink;
  }
  else if (size==3) {
    tlink = 0;
    key = (ct->numseq[j])*625 +
      (ct->numseq[i+3])*125 + (ct->numseq[i+2])*25+(ct->numseq[i+1])*5+(ct->numseq[i]);
    for (count=1; count<=data->numoftriloops&&tlink==0; count++) {
      if (key==data->triloop[count][0]) tlink = data->triloop[count][1];
    }
    energy = data->tstkh[ct->numseq[i]][ct->numseq[j]]
      [ct->numseq[i+1]][ct->numseq[j-1]];
    energy =	data->hairpin[size] + data->eparam[4] + tlink
      +penalty(i, j, ct, data);
  }

  else {
    energy = data->tstkh[ct->numseq[i]][ct->numseq[j]]
      [ct->numseq[i+1]][ct->numseq[j-1]]
      + data->hairpin[size] + data->eparam[4];
  }
  /////
  //      cout << "erg3: "<< energy<<"\n";
  //            cout << "i: " <<i<<"\n";
  //      cout << "j: " <<j<<"\n";



  //check for GU closeure preceded by GG
  if (ct->numseq[i]==3&&ct->numseq[j]==4) {
    if ((i>2&&i<ct->numofbases)||(i>ct->numofbases+2))
      if (ct->numseq[i-1]==3&&ct->numseq[i-2]==3) {

	energy = energy + data->gubonus;
	//if (ct->numseq[i+1]==4&&ct->numseq[j-1]==4)
	//	energy = energy - data->uubonus;
	//if (ct->numseq[i+1]==3&&ct->numseq[j-1]==1)
	//	energy = energy - data->uubonus;


      }
  }

  //check for a poly-c loop
  tlink = 1;
  for (k=1; (k<=size)&&(tlink==1); k++) {
    if (ct->numseq[i+k] != 2) tlink = 0;
  }
  if (tlink==1) {  //this is a poly c loop so penalize
    if (size==3) energy = energy + data->c3;
    else energy = energy + data->cint + size*data->cslope;
  }


  return energy;
}

int erg4(int i, int j, int ip, int jp, structure *ct, datatable *data, bool lfce)
{
  int energy;
  //dangling base
  // jp = 1 => 3' dangle
  // jp = 2 => 5' dangle

  if (lfce) return infinity; //stacked nuc should be double stranded

  if (ip==5) return 0; //dangling nuc is an intermolecular linker

  energy = data->dangle[ct->numseq[i]][ct->numseq[j]][ct->numseq[ip]][jp];
  return energy;
}

//this function calculates whether a terminal pair i, j requires the end penalty
int penalty(int i, int j, structure* ct, datatable *data) {

  /*if ((ct->numseq[i]==1)||(ct->numseq[j]==1)) {
  //this is an A-U pair
  return data->auend;
  }

  else if ((ct->numseq[i]==4)&&(ct->numseq[j]==3)) {
  //this is a G-U pair with 3' terminal U
  return data->auend;
  }
  else if ((ct->numseq[i]==4)&&(ct->numseq[j]==3)) {
  //this is a G-U pair with 3' terminal G
  return data->auend;
  } */
  if (ct->numseq[i]==4||ct->numseq[j]==4)
    return data->auend;
  else return 0; //no end penalty


}

//this function calculates whether a terminal pair i, j requires the end penalty
int penalty2(int i, int j, datatable *data) {


  if (i==4||j==4)
    return data->auend;
  else return 0; //no end penalty


}