Hi all,
Here is the code in c++ and i am new to coding world hence i need little help in converting this code into R
#include <iostream>
#include <cmath>
#include <fstream>
#include <string>
using namespace std;
int main(int argc, char *argv[])
{
std::string protein;
ifstream aa;
aa.open("aa.txt"); //we are getting the data from the file (we assume that we have aa.txt file)
aa>>protein; //the sequence should be in one letter code (uppercase letters)
aa.close();
int ProtLength; //now we are getting protein length
ProtLength = protein.length();
char Asp = 'D';
char Glu = 'E';
char Cys = 'C';
char Tyr = 'Y';
char His = 'H';
char Lys = 'K';
char Arg = 'R';
int AspNumber = 0;
int GluNumber = 0;
int CysNumber = 0;
int TyrNumber = 0;
int HisNumber = 0;
int LysNumber = 0;
int ArgNumber = 0;
int i=0;
for ( i = 0; i <= protein.length() - 1; ++i) // we are looking for charged amino acids
{
if (protein[i] == Asp)
++AspNumber;
if (protein[i] == Glu)
++GluNumber;
if (protein[i] == Cys)
++CysNumber;
if (protein[i] == Tyr)
++TyrNumber;
if (protein[i] == His)
++HisNumber;
if (protein[i] == Lys)
++LysNumber;
if (protein[i] == Arg)
++ArgNumber;
}
double NQ = 0.0; //net charge in given pH
double QN1=0; //C-terminal charge
double QN2=0; //D charge
double QN3=0; //E charge
double QN4=0; //C charge
double QN5=0; //Y charge
double QP1=0; //H charge
double QP2=0; //NH2 charge
double QP3=0; //K charge
double QP4=0; //R charge
double pH = 0.0;
for(;;) //the infinite loop
{
// we are using pK values form Wikipedia as they give quite good approximation
// if you want you can change it
QN1=-1/(1+pow(10,(3.65-pH)));
QN2=-AspNumber/(1+pow(10,(3.9-pH)));
QN3=-GluNumber/(1+pow(10,(4.07-pH)));
QN4=-CysNumber/(1+pow(10,(8.18-pH)));
QN5=-TyrNumber/(1+pow(10,(10.46-pH)));
QP1=HisNumber/(1+pow(10,(pH-6.04)));
QP2=1/(1+pow(10,(pH-8.2)));
QP3=LysNumber/(1+pow(10,(pH-10.54)));
QP4=ArgNumber/(1+pow(10,(pH-12.48)));
NQ=QN1+QN2+QN3+QN4+QN5+QP1+QP2+QP3+QP4;
//if you want to see how the program works step by step uncomment below line
// cout<<"NQ="<<NQ<<"\tat pH ="<<pH<<"\ti:" <<i++<<endl;
if (pH>=14.0)
{ //you should never see this message
cout<<"Something is wrong - pH is higher than 14"<<endl; //
break;
}
if (NQ<=0) // if this condition is true we can stop calculate
break;
pH+=0.01; // if not increase pH
}
ofstream outfile; //we are writing results to outfile.txt
outfile.open("outfile.txt");
outfile << "Protein length: "<<ProtLength<<endl;
outfile << "Number of Asp = "<<AspNumber<<endl;
outfile << "Number of Glu = "<<GluNumber<<endl;
outfile << "Number of Cys = "<<CysNumber<<endl;
outfile << "Number of Tyr = "<<TyrNumber<<endl;
outfile << "Number of His = "<<HisNumber<<endl;
outfile << "Number of Lys = "<<LysNumber<<endl;
outfile << "Number of Arg = "<<ArgNumber<<endl;
outfile << "Protein isoelectric point: "<<pH<<endl;
outfile.close();
cout << "Protein isoelectric point: "<<pH<<endl;
system("PAUSE");
return EXIT_SUCCESS;
}