This is a test version of Biostars. For the public version, visit https://www.biostars.org.
different sequence length alignment

I'm creating for global alignment in this R code.But this code is working same length amino acid sequence.If I want to different length sequence, this program gives Error: subscript out of bounds. Can you help this topic?

#scoring matrix

s1string <- strsplit(seq1, "")[[1]]
s2string <- strsplit(seq2, "")[[1]]
n=length(s1string) #seq1
m=length(s2string) #seq2


s=matrix(data = 0, nrow = m, ncol =n, dimnames = list(s2string,s1string))

for (i in seq_len(m))
    for (j in seq_len(n))
    {
        if (s1string[j]==s2string[i])
        s[i,j]<- 1
        else
        s[i,j]<- -1
    }

seq1string=strsplit(seq1, "")[[1]]
seq2string=strsplit(seq2, "")[[1]]


n1 <- length(seq1string)
n2 <- length(seq2string)

#GLOBALALIGNMENT

fmatrix <-template<- matrix(0, nrow = n2+1 , ncol = n1+1)

d = -1  #this is gap penalty



for(i in 0 : nchar(seq2))
{
    fmatrix[i+1,1] = d * I }

for(j in 0 : nchar(seq1))
{
    fmatrix[1,j+1] = d * j
}

for(i in seq_len(n2))
    for(j in  seq_len(n1))
    {
        match  =  fmatrix[i,j] + sbmx[i,j]
        delete = fmatrix[i,j+1] + d
        insert = fmatrix[i+1,j] + d

        cho = c(fmatrix[i,j] + sbmx[i,j],fmatrix[i,j+1] + d,fmatrix[i+1,j] + d)

        fmatrix[i+1,j+1] = max(match,delete,insert)
        template[i+1,j+1] <- which.max(cho)
    }

colnames(fmatrix) = strsplit( paste(" " , seq1, sep=""), "")[[1]];
rownames(fmatrix) = strsplit( paste(" " , seq2, sep=""), "")[[1]];

#traceback

AlignmentA = ""
AlignmentB = "" I = nchar(seq1)+1
j = nchar(seq2)+1

while(i >0 && j>0)
{
   if( template[i,j]==1 )
   {
        AlignmentA=paste(substr(seq1,j-1,j-1),AlignmentA, sep="")
        AlignmentB=paste(substr(seq2,i-1,i-1),AlignmentB, sep="")
        j=j-1
        i=i-1
    }
    else if (template[i,j]==3 )
    {

        AlignmentA=paste(substr(seq1,j-1,j-1),AlignmentA, sep="")
        AlignmentB=paste("-",AlignmentB, sep="")
        j=j-1
    }
    else (template[i,j]==2)
    {
        AlignmentA=paste("-",AlignmentA, sep="")
        AlignmentB=paste(substr(seq2,i-1,i-1),AlignmentB, sep="")
        i=i-1
    }    

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

Indices start at 1 in R...

Ah, I see you're adding 1 to things, never mind that. Anyway, it'd be helpful if you mentioned exactly where you're getting the error.

0 answers

No answers yet.

Log in to answer this question.