Lets say I have a DNA sequence with an ambiguous base, N, where N can represent any base (its a flex position).
dna.seq <- 'ATGCN'
I want a vector of every possible DNA sequence this could represent. It would look like:
c('ATGCA','ATGCT','ATGCG','ATGCC')
The solution needs to account for dna sequences with multiple N characters as well, which will create many more potential DNA sequences.
CJ from data.table can help you here:
library(data.table)
dna.seq <- 'ATGCN'
# split into components
l = tstrsplit(dna.seq, '', fixed = TRUE)
# replace N with all possibilities
all_bases = c('A', 'T', 'C', 'G')
l = lapply(l, function(x) if (x == 'N') all_bases else x)
# use CJ and reduce to strings:
Reduce(paste0, do.call(CJ, l))
# [1] "ATGCA" "ATGCC" "ATGCG" "ATGCT"
Flexibility to handle multiple N:
dna.seq <- 'ATNCN'
Reduce(paste0, do.call(CJ, l))
# [1] "ATACA" "ATACC" "ATACG" "ATACT" "ATCCA" "ATCCC" "ATCCG" "ATCCT"
# [9] "ATGCA" "ATGCC" "ATGCG" "ATGCT" "ATTCA" "ATTCC" "ATTCG" "ATTCT"
If you wanted to drop the data.table dependency you could replace tstrsplit with t(strsplit()) and CJ with expand.grid; you'll just be sacrificing computational speed.
dna.seq <- 'ATNGCN'
dna.seq.copy = dna.seq
while(grepl("N", dna.seq.copy[1])){
dna.seq.copy = as.vector(sapply(c("A", "C", "T", "G"), function(x) sub("N", x, dna.seq.copy)))
}
dna.seq.copy
# [1] "ATAGCA" "ATCGCA" "ATTGCA" "ATGGCA" "ATAGCC" "ATCGCC" "ATTGCC" "ATGGCC" "ATAGCT" "ATCGCT" "ATTGCT" "ATGGCT" "ATAGCG" "ATCGCG" "ATTGCG"
#[16] "ATGGCG"
If you love us? You can donate to us via Paypal or buy me a coffee so we can maintain and grow! Thank you!
Donate Us With