jayoung
commented
5 years ago PS no need to hurry to implement this for my benefit - I made my own functions that (I think) do what I want. Not elegant, and I'm sure not robust to the user doing unintended things, but they work at least on my own alignments.
Still might be a useful thing to implement in Biostrings (although maybe I should be using some other package for this?)
if you're curious:
## myAln is a DNAStringSet with a gapped alignment
getCodons <- function(myAln) {
seqs <- as.character(myAln)
len <- width(myAln)[1]
starts <- seq(from=1, to=len, by=3)
ends <- starts + 2
myViews <- lapply(myAln, function(x) {
Views(x, starts, ends)
})
myCodons <- lapply(myViews, function(x) {
as.character(DNAStringSet(x))
})
myCodons
}
## myCodons is a simple character vector, each item is a codon (like one of the items in a list generated by getCodons)
translateCodons <- function(myCodons, unknownCodonTranslatesTo="-") {
## make new genetic code
gapCodon <- "-"
names(gapCodon) <- "---"
my_GENETIC_CODE <- c(GENETIC_CODE, gapCodon)
## translate the codons
pep <- my_GENETIC_CODE[myCodons]
## check for codons that were not possible to translate, e.g. frameshift codons
if (sum(is.na(pep))>0) {
cat("\nwarning - there were codons I could not translate. Using this character", unknownCodonTranslatesTo, "\n\n")
pep[ which(is.na(pep)) ] <- unknownCodonTranslatesTo
}
## prep for output
pep <- paste(pep, collapse="")
return(pep)
}
## wrap those functions together into one:
translateGappedAln <- function(myAln, unknownCodonTranslatesTo="-") {
myCodons <- getCodons(myAln)
myAAaln <- AAStringSet(unlist(lapply(myCodons, translateCodons, unknownCodonTranslatesTo=unknownCodonTranslatesTo)))
return(myAAaln)
}
## tests:
test1 <- DNAStringSet( c("ATGGCTGCGCGGGGC", "ATGGCTGGGCGGGGC") )
test2 <- DNAStringSet( c("ATGGCTGCGCGGGGC", "ATGGCTG-GCGGGGC") )
translateGappedAln(test1)
translateGappedAln(test2)
translateGappedAln(test2, unknownCodonTranslatesTo="X")
hpages
sdalin
ahl27
hi there,
I've been reading in some multiple sequence alignments, as DNAStringSet objects. They're nucleotide alignments that encode proteins. I've been playing with using 'translate', but it looks like it's not set up to deal with gap characters.
I might be missing some nice alternative way to do this, but if not, I guess I'd like to suggest enhancement to better deal with in-frame nucleotide alignments of coding sequences. I think the code below will show you what I mean, but if it's not clear please let me know.
thanks!
Janet
Dr. Janet Young
Malik lab http://research.fhcrc.org/malik/en.html
Division of Basic Sciences Fred Hutchinson Cancer Research Center 1100 Fairview Avenue N., A2-025, P.O. Box 19024, Seattle, WA 98109-1024, USA.
tel: (206) 667 4512 email: jayoung ...at... fredhutch.org