Question: traceback algorithm global alignment in R code
0
gravatar for gizemtatar
5.3 years ago by
gizemtatar0
Turkey
gizemtatar0 wrote:

  My needleman-wunsch global alignment scoring matrix :       

         G C A T  G C U
     0 -1 -2 -3 -4 -5 -6 -7
G -1  1  0 -1 -2 -3 -4 -5
A -2  0  0  1  0 -1 -2 -3
T -3 -1 -1  0  2  1  0 -1
T -4 -2 -2 -1  1  1  0 -1
A -5 -3 -3 -1  0  0  0 -1
C -6 -4 -2 -2 -1 -1  1  0
A -7 -5 -3 -1 -2 -2  0  0

I want traceback algorithm and find best alignment sequence.

I try  for this code.

 Astring <- strsplit(A, "")[[1]]

    Bstring <- strsplit(B, "")[[1]]

    n=length(Astring) #seq1

    m=length(Bstring) #seq2



    s=matrix(data = 0, nrow = m, ncol =n, dimnames = list(Bstring,Astring))



    for (i in seq_len(m))

    for (j in seq_len(n))

    {

        if (Astring[j]==Bstring[i])

        s[i,j]<- 1

        else

        s[i,j]<- -1

    }



    fmatrix = matrix(0, nrow = (nchar(A)+1) , ncol = nchar(B)+1)

    d = -1  #this is gap penalty

    for(i in 0 : nchar(B)){

        fmatrix[i+1,1] = d * i  #populates initial row with gap penalty

    }

    for(j in 0 : nchar(A)){

        fmatrix[1,j+1] = d * j

    }

    for(i in 1 : nchar(B)){

        for(j in 1 : nchar(A)) {

             #get me sccore for the pair of aa or nt

            match  =  fmatrix[i,j] + s[i,j]

            delete = fmatrix[i,j+1] + d

            insert = fmatrix[i+1,j] + d

            fmatrix[i+1,j+1] = max(match,delete,insert)

        }

    }

    colnames(fmatrix) = strsplit( paste(" " , A, sep=""), "")[[1]];

    rownames(fmatrix) = strsplit( paste(" " , B, sep=""), "")[[1]];

    #return(fmatrix)


AlignmentA = ""

AlignmentB = ""


i = nchar(B) +1

j = nchar(A) +1

maxscore=numeric(length=0)



    while(i > 1 && j > 1)

    {

        CurrentScore= fmatrix[i,j]

        ScoreDiag= fmatrix[i-1,j-1]

        ScoreLeft=fmatrix[i,j-1]

        ScoreUp=fmatrix[i-1,j]

        maxscore= max(CurrentScore,ScoreDiag,ScoreLeft,ScoreUp)



        if( CurrentScore <= ScoreDiag  ) {

            AlignmentA=paste(substr(A,j-1,j-1),AlignmentA, sep="")

            AlignmentB=paste(substr(B,i-1,i-1),AlignmentB, sep="")

            j=j-1

            i=i-1

        }

        else if (ScoreLeft >= CurrentScore) {

            AlignmentA=paste(substr(A,j-1,j-1),AlignmentA, sep="")

            AlignmentB=paste("-",AlignmentB, sep="")



            j=j-1

        } else if (CurrentScore>=ScoreUp+d) {

            AlignmentA=paste("-",AlignmentA, sep="")

            AlignmentB=paste(substr(B,i-1,i-1),AlignmentB, sep="")



            i=i-1

             }

        ans <- matrix(c(al1=paste(AlignmentA, collapse=""), al2=paste(AlignmentB, collapse="")))

        print(ans)

But I don't find best alignment. The best alignments :

      GCATG-CU      GCA-TGCU      GCAT-GCU
      G-ATTACA      G-ATTACA      G-ATTACA

Can you help this topic please? 

 

 

alignment sequence R • 3.3k views
ADD COMMENTlink modified 11 months ago by parkkh11880 • written 5.3 years ago by gizemtatar0
1

Why reinvent the wheel when pairwiseAlignment() is a well tested function from Biostrings that does all of this for you anyway?

ADD REPLYlink modified 11 months ago • written 5.3 years ago by RamRS30k

I need manual R code.

ADD REPLYlink written 5.3 years ago by gizemtatar0
1

For an assignment?

ADD REPLYlink written 5.3 years ago by Brice Sarver3.5k
0
gravatar for parkkh1188
11 months ago by
parkkh11880
parkkh11880 wrote:

fmatrix = matrix(0,nrow = (nchar(A)+1), ncol = nchar(B)+1)

Switch to

fmatrix = matrix(0,nrow = (nchar(B)+1), ncol = nchar(A)+1)

Switch nchar(A) to nchar(B), nchar(B) to nchar(A)

ADD COMMENTlink written 11 months ago by parkkh11880
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 2176 users visited in the last hour