Algorytm naiwny

- liczymy ile jest jest poszczególnych aminokwasów obdarzonych ładunkiem:

- podstawiamy do wzoru:



Kod programu:

#include <iostream>
#include <cmath>                //biblioteka potrzebna do obliczeń matematycznych
#include <fstream>                //biblioteka do odczytu i zapisu plików
#include <string>                    //biblioteka obsługująca stringi - ciągi znaków

using namespace std;

int main(int argc, char *argv[])
{
      std::string protein;
          ifstream aa;
       aa.open("aa.txt");        //pobieramy dane z pliku zrodlowego czyli musimy mieć taki plik w folderze w którym wykonujemy
       aa>>protein;                //program, powinien on zawierać sekwencję zapisaną w kodzie jednoliterowym (duże litery)
           aa.close();
        
 int dlugosc;                       //obliczamy długość naszego białka - korzystamy z wbudowanej funkcji length
   dlugosc = protein.length();
 
 
   char Asp = 'D';
   char Glu = 'E';
   char Cys = 'C';
   char Tyr = 'Y';
   char His = 'H';
   char Lys = 'K';
   char Arg = 'R';

int IloscAsp = 0;
int IloscGlu = 0;
int IloscCys = 0;
int IloscTyr = 0;
int IloscHis = 0;
int IloscLys = 0;
int IloscArg = 0;

int i=0;


for ( i = 0; i <= protein.length() - 1; ++i)              //  sprawdzamy całe białko aminokwas po aminokwasie
                                                                            //w poszukiwaniu tych obdarzonych ładunkiem i je liczymy
    {
              if (protein[i] == Asp)
                 ++IloscAsp;

              if (protein[i] == Glu)
                 ++IloscGlu;

              if (protein[i] == Cys)
                 ++IloscCys;

              if (protein[i] == Tyr)
                 ++IloscTyr;

              if (protein[i] == His)
                 ++IloscHis;

              if (protein[i] == Lys)
                 ++IloscLys;

              if (protein[i] == Arg)
                 ++IloscArg;
    }
   
    double NQ = 0.0; //całkowit ładunek białka przy określonym 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;           //punkt startowy

for(;;)                //nieskończona pętla
{

// jak widać używamy wartości z Wikipedii, ponieważ dają one całkiem w przybliżeniu średni wynik dla wartości z innych źródeł
//  jeśli ktoś chce może dowolnie zmienić podane niżej wartości

    QN1=-1/(1+pow(10,(3.65-pH)));                                        
    QN2=-IloscAsp/(1+pow(10,(3.9-pH)));           
    QN3=-IloscGlu/(1+pow(10,(4.07-pH)));           
    QN4=-IloscCys/(1+pow(10,(8.18-pH)));           
    QN5=-IloscTyr/(1+pow(10,(10.46-pH)));        
    QP1=IloscHis/(1+pow(10,(pH-6.04)));            
    QP2=1/(1+pow(10,(pH-8.2)));                
    QP3=IloscLys/(1+pow(10,(pH-10.54)));           
    QP4=IloscArg/(1+pow(10,(pH-12.48)));            

    NQ=QN1+QN2+QN3+QN4+QN5+QP1+QP2+QP3+QP4;

  //odkomentuj linię poniżej, jeśli chcesz zobaczyć jak program liczy krok po kroku
  //      cout<<"NQ="<<NQ<<"\tprzy pH ="<<pH<<"\ti:" <<i++<<endl;      
      

if (pH>=14.0)
   {                                                 //ten komunikat nie ma prawa się pojawić
  cout<<"Cos poszlo nie tak - pH powyzej 14"<<endl;  //jeśli tak się stało skontaktuj się z
   break;                                            //producentem [zachowaj wartości, które
   }                                                 //spowodowały wyjście poza zasięg programu]  :-)

if (NQ<=0)                            // spełnienie tego warunku rozwiązuje nasze zadanie i możemy opuściś pętle
    break;

 pH+=0.01;                            // jeśli nie to zwiększamy pH i liczymy od nowa

 }
 
 ofstream outfile;                //zapisujemy szczegółowe wyniki do pliku outfile.txt
 outfile.open("outfile.txt");          
   outfile << "Dlugosc wynosi: "<<dlugosc<<endl;
   outfile << "Ilosc Asp = "<<IloscAsp<<endl;  
   outfile << "Ilosc Glu = "<<IloscGlu<<endl;
   outfile << "Ilosc Cys = "<<IloscCys<<endl;
   outfile << "Ilosc Tyr = "<<IloscTyr<<endl;
   outfile << "Ilosc His = "<<IloscHis<<endl;
   outfile << "Ilosc Lys = "<<IloscLys<<endl;
   outfile << "Ilosc Arg = "<<IloscArg<<endl;
   outfile << "Punkt izoelektryczny podanego białka wynosi: "<<pH<<endl;
 outfile.close();
   cout << "Punkt izoelektryczny podanego bialka wynosi: "<<pH<<endl;

    system("PAUSE");
    return EXIT_SUCCESS;
}

Poniżej możesz pobrać skompilowny program gotowy do użytku Download: Punkt izoelektryczny.exe

Obliczanie pI z wykorzystaniem bisekcji.

Jak łatwo zauważyć powyższe równanie będzie liczone wiele razy (sądząc po średniej wartości pI wszystkich białek około 650 razy przy dokładności 0.01) w związku z czym naiwny algorytm  będzie dość czasochłonny. Przy liczeniu pojedyńczych wartości pI nie ma to większego znaczenia jednak w innym przypadku algorytm taki jest niewydajny. W związku z czym, aby przyspieszyć obliczanie pI zastosujemy tzw. bisekcję.

Bisekcja to naukowe określenie jednej z najprostszych i zarazem najskuteczniejszej techniki obliczania kłopotliwych lub skomplikowanych równań.

Działanie bisekcji wytłumaczymy na naszym przykładzie:
1) wyobraźmy sobie przestrzeń możliwych rozwiązań w postaci odcinka o określonej długości (u nas będzie to odcinek o długości 14),
2) następnie dzielimy odcinek na połowę i sprawdzamy wynik naszego równania w punkcje środkowym. Jeśli NQ>0 to nasze rozwiązanie leży między 0 a 7, jeśli NQ<=0 to rozwiązanie leży po drugiej stronie
3) postępuj jak w punkcie 2 do momentu gdy osiągniesz zamierzoną dokładność

Podejście takie do obliczania punktu izoelektrycznego zaproponował niejaki pan Tabb już w 2001 tyle, że czytając artykuł można odnieść wrażenie, że nawet nie wiedział, że to co proponuje to tak naprawdę bisekcja i algorytm który wymyślił jest dosyć znany. Nie mniej jednak artykuł ten jest godny polecenia jako doskonałe źródło teoretycznych podstaw dla obliczania punktu izoelektrycznego (bo praktyki to tam nie ma wcale).

Tabb DL (2001) An algorithm for isoelectric point estimation.

Dla naszego przykładu w języku C++ wygląda to mniej więcej tak:

#include <iostream>
#include <cmath>                //biblioteka potrzebna do obliczeń matematycznych
#include <fstream>                //biblioteka do odczytu i zapisu plików
#include <string>                    //biblioteka obsługująca stringi - ciągi znaków

using namespace std;

int main(int argc, char *argv[])
{
      std::string protein;
          ifstream aa;
       aa.open("aa.txt");        //pobieramy dane z pliku zrodlowego czyli musimy mieć taki plik w folderze w którym wykonujemy
       aa>>protein;                //program, powinien on zawierać sekwencję zapisaną w kodzie jednoliterowym (duże litery)
           aa.close();
       
 int dlugosc;                       //obliczamy długość naszego białka - korzystamy z wbudowanej funkcji length
   dlugosc = protein.length();
 
 
   char Asp = 'D';
   char Glu = 'E';
   char Cys = 'C';
   char Tyr = 'Y';
   char His = 'H';
   char Lys = 'K';
   char Arg = 'R';

int IloscAsp = 0;
int IloscGlu = 0;
int IloscCys = 0;
int IloscTyr = 0;
int IloscHis = 0;
int IloscLys = 0;
int IloscArg = 0;

int i=0;


for ( i = 0; i <= protein.length() - 1; ++i)              //  sprawdzamy całe białko aminokwas po aminokwasie
                                                                            //w poszukiwaniu tych obdarzonych ładunkiem i je liczymy
    {
              if (protein[i] == Asp)
                 ++IloscAsp;

              if (protein[i] == Glu)
                 ++IloscGlu;

              if (protein[i] == Cys)
                 ++IloscCys;

              if (protein[i] == Tyr)
                 ++IloscTyr;

              if (protein[i] == His)
                 ++IloscHis;

              if (protein[i] == Lys)
                 ++IloscLys;

              if (protein[i] == Arg)
                 ++IloscArg;
    }
  
    double NQ = 0.0; //całkowit ładunek białka przy określonym 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 = 6.5;             //punkt środkowy pI = 6.5 - teoretycznie powinien wynosić 7,
                                 //ale białka mają pI średnio w pH=6.5, więc zwiększamy
    double pHprev = 0.0;         //prawdopodobieństwo, że trafimy od razu tam gdzie trzeba
    double pHnext = 14.0;        //0-14 zakres pH - obszar poszukiwań
    double X = 0.0;
    double E = 0.01;             //epsilon - dokładność [pI = pH ± E]
    double temp = 0.0;
      
cout<<endl<<endl;
for(;;)                //nieskończona pętla
{

// jak widać używamy wartości z Wikipedii, ponieważ dają one całkiem w przybliżeniu średni wynik dla wartości z innych źródeł
//  jeśli ktoś chce może dowolnie zmienić zmienić odane niżej wartości

    QN1=-1/(1+pow(10,(3.65-pH)));                                       
    QN2=-IloscAsp/(1+pow(10,(3.9-pH)));          
    QN3=-IloscGlu/(1+pow(10,(4.07-pH)));          
    QN4=-IloscCys/(1+pow(10,(8.18-pH)));          
    QN5=-IloscTyr/(1+pow(10,(10.46-pH)));       
    QP1=IloscHis/(1+pow(10,(pH-6.04)));           
    QP2=1/(1+pow(10,(pH-8.2)));               
    QP3=IloscLys/(1+pow(10,(pH-10.54)));          
    QP4=IloscArg/(1+pow(10,(pH-12.48)));           

    NQ=QN1+QN2+QN3+QN4+QN5+QP1+QP2+QP3+QP4;

  //odkomentuj linię poniżej, jeśli chcesz zobaczyć jak program liczy krok po kroku
  //      cout<<"NQ="<<NQ<<"\tprzy pH ="<<pH<<"\ti:" <<i++<<endl;     
     

if (pH>=14.0)
   {                                                 //ten komunikat nie ma prawa się pojawić
  cout<<"Cos poszlo nie tak - pH powyzej 14"<<endl;  //jeśli tak się stało skontaktuj się z
   break;                                            //producentem [zachowaj wartości, które
   }                                                 //spowodowały wyjście poza zasięg programu]  :-)

//%%%%%%%%%%%%%%%%%%%%%%%%%   BISEKCJA   %%%%%%%%%%%%%%%%%%%%%%%%

    if(NQ<0)              //wyszliśmy poza zakres, więc nowy punkt pH musi być mniejszy    
     {                   
        temp = pH;
        pH = pH-((pH-pHprev)/2);
        pHnext = temp;
        cout<<"pH: "<<pH<<", \tpHnext: "<<pHnext<<endl;
    }
    else                  //jeszcze jesteśmy na plusie, więc użyliśmy za małego pH, zwiększamy je
    {                     // odpowiednio
        temp = pH;
        pH = pH + ((pHnext-pH)/2);
        pHprev = temp;
        cout<<"pH: "<<pH<<",\t pHprev: "<<pHprev<<endl;
 
    }

       if ((pH-pHprev<E)&&(pHnext-pH<E)) // warunek wyjścia z pętli - osiągnięcie zamierzonej dokładności
          break;
   
//podsumowanie: zastosowanie bisekcji pozwoliło skrócić liczbę iteracji średnio do 10!!!!

 }
 
 ofstream outfile;                //zapisujemy szczegółowe wyniki do pliku outfile.txt
 outfile.open("outfile.txt");         
   outfile << "Dlugosc wynosi: "<<dlugosc<<endl;
   outfile << "Ilosc Asp = "<<IloscAsp<<endl; 
   outfile << "Ilosc Glu = "<<IloscGlu<<endl;
   outfile << "Ilosc Cys = "<<IloscCys<<endl;
   outfile << "Ilosc Tyr = "<<IloscTyr<<endl;
   outfile << "Ilosc His = "<<IloscHis<<endl;
   outfile << "Ilosc Lys = "<<IloscLys<<endl;
   outfile << "Ilosc Arg = "<<IloscArg<<endl;
   outfile << "Punkt izoelektryczny podanego białka wynosi: "<<pH<<endl;
 outfile.close();

   cout << "\n\nPunkt izoelektryczny podanego bialka wynosi: "<<pH<<endl;

    system("PAUSE");
    return EXIT_SUCCESS;
}

Czy  opłaciło się możesz przekonać się uruchamiając program. Jak widać tym razem program liczy główne równanie około 10 razy (zamiast średnio 650).

Poniżej możesz pobrać skompilowny program gotowy do użytku Download: pI_bisekcja.exe

Usprawnianie szybkości działania ciąg dalszy...

Choć teoretycznie można usprawnić działanie naszego programu narzucając odpowiednie warunki tak, aby wybór kolejnego punktu był bardziej trafny to nie przyniesie to zbyt wielkiego zysku, ponieważ sprawdzanie kolejnych warunków w dużym stopniu zniweczy to co zyskamy. Istnieje jeszcze jeden sposób na osiągnięcie znaczącej poprawy szybkości programu (10-krotne). W tym przypadku zmiana równoznaczna będzie ze "strzałem w dziesiątkę", a tego nie zapewni nam żaden algorytm. Więc jak to można osiągnąć?

Otóż jest to dosyć proste, jedyne do musimy zrobić to przyjrzeć się problemowi od innej strony.

Spójrzmy na główne równanie. Dosyć skomplikowane,  mnożymy, dzielimy, a nawet co gorsza potęgujemy (w wykładniku znajdują się ujemne liczby zmiennoprzecinkowe). Ale zastanówmy się czy wszystkie te operacje są konieczne. Otóż jeśli zadowolimy się określoną dokładnością  to możemy zauważyć, że jest z góry określona liczba możliwych rozwiązań mianownika cząstkowych składowych równania głównego (np. dla 0.01 wynosi ona 1400 x 9). Jeśli obliczymy mianowniki dla każdego pH buforu z określoną dokładnością i zapiszemy je w tablicy to łatwo pozbędziemy się potęgowania (przygotowanie takiej tablicy wymaga napisania prościutkiego programiku, który czytelnik powinien dosyć łatwo napisać na podstawie poprzednich przykładów). Teraz jedyne co będziemy musieli zrobić to dzielić wartość odpowiadającą liczbie danego aminokwasu przez odpowiednią wartość pobraną z tablicy. Ponadto ładunek C i N końca białka jest stały w określonym pH tak, więc w tych dwóch przypadkach możemy uzyskać dokładne wartości raz na zawsze i zrezygnować z ich liczenia całkowicie.

Na pierwszy rzut oka może się wydawać, że w ten sposób pozbawimy algorytm dokładności należy jednak pamiętać, że większa precyzja nie jest konieczna, ponieważ wynik i tak jest obarczony sporym błędem z powodów określonych powyżej. W tym momencie bardziej krytyczne staną się wartości pK jakich używamy niż wstawianie do równania wykładnika różniącego się o tysięczne części.

Poniżej możesz pobrać skompilowny program gotowy do użytku (folder sources zawiera żródła programu) Download: pI_optymal.exe

góra strony