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)
}
Indices start at 1 in R...
I don't understand.
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.