Needleman_Wunsch_Alignment.cpp

#include <iostream>
#include<vector>

using namespace std;

pair<string,string> NWA(string S,string T, int g, int mis, int match)
{
    int m=S.size(), n=T.size(), i, j;
    vector <vector <pair<int,int>>> F(m+1, vector <pair<int,int>>(n+1));
    for(i=0;i<n+1;i++) //n+1 columns fir T
    {
        F[0][i].first=-i*g;
        F[0][i].second = -1; // 1 means UP, -1 left
    }
    for(i=0;i<m+1;i++) //n+1 rows fir S
    {
        F[i][0].first=-i*g;
        F[i][0].second=1;//left
    }

    int tptr=0,sptr=0;
    //MATCH, MIS->> 2
    int Sgap,Tgap; //Make a visualiser
    for(i=1;i<m+1;i++){
        for(j=1;j<n+1;j++){
            Sgap=F[i][j-1].first-g;
            Tgap=F[i-1][j].first-g;
            if(Tgap>Sgap)///case when GAP IN T (1) PUT AT THE JTH INDEX Hence S aligned upto i-1 and new xi in s has none but a gap to align to is better scoring than a gap in S (2)
            {
                if(S[i-1]==T[j-1]){
                    if(Tgap>F[i-1][j-1].first+match){
                        F[i][j].first=Tgap;
                        F[i][j].second=1;
                    }
                    else{
                        F[i][j].first=F[i-1][j-1].first+match;
                        F[i][j].second=2;
                    }
                }
                else{
                    if(Tgap>F[i-1][j-1].first-mis){
                        F[i][j].first=Tgap;
                        F[i][j].second=1;
                    }
                    else{
                        F[i][j].first=F[i-1][j-1].first-mis;
                        F[i][j].second=2;
                    }
                }
            }
            else{ ///SGAP>TGAP
                if(S[i-1]==T[j-1]){
                    if(Sgap>F[i-1][j-1].first+match){
                        F[i][j].first=Sgap;
                        F[i][j].second=-1;
                    }
                    else{
                        F[i][j].first=F[i-1][j-1].first+match;
                        F[i][j].second=2;
                    }
                }
                else{
                    if(Sgap>F[i-1][j-1].first-mis){
                        F[i][j].first=Sgap;
                        F[i][j].second=-1;
                    }
                    else{
                        F[i][j].first=F[i-1][j-1].first-mis;
                        F[i][j].second=2;
                    }
                }
            }
        }
    }

    //rebuilding the sequences for matching S and T
    // vector <string> S_al (m);
    // vector <string> T_al (n);

    string S_al="";
    string T_al ="";
    i=m;j=n;
    while(i>0 || j>0)
    {
        if(F[i][j].second==2)
        {
            S_al=S[i-1]+S_al;
            T_al=T[j-1]+T_al;
            i--;
            j--;
        }
        else if(F[i][j].second==-1)
        {
            S_al="-"+S_al;
            T_al=T[j-1]+T_al;
            j--;
        }
        else if(F[i][j].second==1)
        {
            S_al=S[i-1]+S_al;
            T_al="-"+T_al;
            i--;
        }
    }
    return {S_al,T_al};
}

int main()
{
    cout<<NWA("AAGC", "AGT", 2,1,1).first<<" "<<NWA("AAGC", "AGT", 2,1,1).second<<"\n";
}