different sequence length alignment
0
0
Entering edit mode
8.9 years ago
gizemtatar • 0

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 • 2.3k views
ADD COMMENT
0
Entering edit mode

Indices start at 1 in R...

ADD REPLY
0
Entering edit mode

I don't understand.

ADD REPLY
0
Entering edit mode

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.

ADD REPLY

Login before adding your answer.

Traffic: 1737 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6