Finding mismatches between two sequences using R

[chudar]

Hi,

I am new to R programming. I have two fasta files namely WT and MT that contains 3 protein sequences which is given below

WT:
>seq1
MGFTIIIJKLPADFDEMMJDHGASFWEIHDI
>seq2
IIJHGTIPIOIZTRWERQQMSCVYMMDMMLI
>seq3
KPOIHUTWMMKLIZGTIIWHSAKDFJGVAHD

MT:
>seq1
TGFTHIHJKLPADFDEMTJDHGASFWEIHDH
>seq2
HHJHGTIPIOIZTRWERQQMSCVYTMDTMLH
>seq3
KPOIHUTWTTKLIZGTHHWHSAKDFJGVAHD


WT contains set of fasta sequences without any mutated aminoacids while the MT contains same protein sequences of WT but some amino acid are mutated.

I want compare the two sequences and find out counts for particular mutations wherever M mutated to T, I mutated to H for every protein and list them like Seq1 MT_counts:xx;IH_counts:YY Seq2:MT_counts:xx; IH_counts:yy Seq3:MT_counts:xx IH_counts:yy

I have written a small function in R where I read the fasta sequencces using read.fasta from seqinr package and later using two for loops for iterating through all sequences and another for iterating through every letter of seq in WT and comparing with its counterpart in MT.

library(seqinr)
wt=read.fasta("C:/Users/tsekaran/Documents/sample_ref_protein.fasta")
mt=read.fasta("C:/Users/tsekaran/Documents/sample_mut_protein.fasta")
mismatch_finder=function(wt,mt)
{
for (i in 1:length(names(wt)))
{ 
  MT_count=0
  IH_count=0
  #for(j in 1:length(wt$seq1))
  for(j in 1:length(wt[[i]]))
  {  
    if(wt[j]=="m" && mt[j]=="t" )
    {
      MT_count=MT_count+1
    }
  else if(wt[j]=="i" && mt[j]=="h" )
    {
      IH_count=IH_count+1
    }
  }
  print(names(wt[i]))
  print(MT_count)
  print(IH_count)
}
}

But it gives me counts for MT and IH as 0 for all seq1 ,seq2 and seq3. I know there are some mutated aminoacid but I am surprised as the function reports the count as 0. Could some guide me where if there is any mistake in comparison. Kindly help me. Thanks in advance.

[dpryan]

If you ask people to help you with an error message then you need to post the error message. I see a number of problems with your code, from using non-existent variables ("wtf") to incorrect syntax ("wt$[j]").

[chudar]

If you ask people to help you with an error message then you need to post the error message. I see a number of problems with your code, from using non-existent variables ("wtf") to incorrect syntax ("wt$[j]").

Dear DpRyan,

Thank you very much for your concern. I have updated my R script and have edited my question also. When I ran my updated script it expected it give number of mismatches but it gave me only 0 for all my sequences. I dont whether my comparison code is corrrect. Kindly take a look and guide me please. Thanks in advance

[dpryan]

Grr, the site ate my initial reply!

I assume you meant to write this, since otherwise you're trying to compare a character and a list (wt[i] == "m"), which will always be false and thus yield counts of 0:

for (i in 1:length(names(wt)))
{ 
  MT_count=0
  IH_count=0
  #for(j in 1:length(wt$seq1))
  for(j in 1:length(wt[[i]]))
  {  
    if(wt[[i]][j]=="m" && mt[[i]][j]=="t" )
    {
      MT_count=MT_count+1
    }
  else if(wt[[i]][j]=="i" && mt[[i]][j]=="h" )
    {
      IH_count=IH_count+1
    }
  }
  print(names(wt[i]))
  print(MT_count)
  print(IH_count)
}

This could be done more succinctly (and probably with better performance) with:

for (i in 1:length(names(wt))) {
    print(names(wt)[i])
    print(length(intersect(which(wt[[i]] == "m"), which(mt[[i]] == "t"))))
    print(length(intersect(which(wt[[i]] == "i"), which(mt[[i]] == "h"))))
}

[chudar]

Hi Dpryan,

Thank you very much for code. it works fine. I would like to output the data into a data frame so that it looks like

ID     MT  IH
seq1   2    5
seq2   4    7
seq3   6    9

So I edited the code like below where I created an data frame with NA and used it inside the for loop

mismatch=function(wt,mt)
{
df=data.frame(ID=NA,MT=NA,IH=NA)
for (i in 1:length(names(wt))) {
  df$ID[i]=names(wt)[i]
  df$MT[i]= length(intersect(which(wt[[i]] == "m"), which(mt[[i]] == "t")))
  df$IH[i]=length(intersect(which(wt[[i]] == "i"), which(mt[[i]] == "h")))
}
 return (df)
}

but it gives me an error as follows

Error in `$<-.data.frame`(`*tmp*`, "ID", value = c("seq1", "seq2")) : 
  replacement has 2 rows, data has 1

I have no clue how am I ending up like this.

[dpryan]

Dataframes need to have identical length rows, so you can't add something to ID without also adding values to MT and IH. The actual R way to do this would be something like:

library(seqinr)
wt=read.fasta("C:/Users/tsekaran/Documents/sample_ref_protein.fasta")
mt=read.fasta("C:/Users/tsekaran/Documents/sample_mut_protein.fasta")

getDiffs <- function(wt, mt) {
    c(names(wt),
    length(intersect(which(wt == "m"), which(mt == "t"))),
    length(intersect(which(wt == "i"), which(mt == "h"))))
}

res <- as.data.frame(t(mapply(getDiffs, wt, mt)))
names(res) <- c("MT","IT")

[dpryan]

BTW, you could also just use rbind() if you want to keep the for loop.