Subscribed unsubscribe Subscribe Subscribe

CPSAM.org

computer, programming, statistics and more!

BWA

BWAアルゴリズムをRで実装してみる

require('stringr')

REF <- "CAGGTTGGTTGGTTGG"
reads <- c("AGGT", "GGTT", "TTCA")

wheel <- function(x){
  len <- str_length(x)
  L <- str_sub(x, 1, 1)
  R <- str_sub(x, 2, len)
  return(str_c(R, L))
}

myAline <- function(read, ref){

  tmp <- str_c(REF, "$")
  SA <- tmp

  for(i in 1:(str_length(tmp))-1){
    t <- wheel(SA[i])
    SA <- c(SA, t)
  }

  SA.trim <- (str_split_fixed(SA[str_order(SA)], "\\$", n=2))[,1]
  SA.add <- str_c(SA.trim, "$")

  hits <- list(seq(1:length(SA.add)))

  for(i in 1:str_length(read)){
    hit <- which(str_sub(SA.add[hits[[i]]], i, i)==str_sub(read, i, i))
    if(length(hit) < 1){
      break
    } else {
      hits[[i+1]] <- hits[[i]][hit]
    }
  }

  if(length(hit)<1){
    return(0)
  } else {
    hits.pos <- str_length(REF) - str_length(SA.add[hits[[str_length(read)]]]) + 2
    return(hits.pos)
  }

}

print(REF)

for(i in reads){
  print(c(i, myAline(i, REF)))
}
Remove all ads