[R] function to compute consensus DNA sequence by plurality?

[Jean lobry]

Kim,

Is is what you want?

tmp <- readLines(textConnection(
"TGCATACACCGACAACATCCTCGACGACTACACCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG
AGCATACACCGACAACATCCTCGATGACTACTGCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG
AGCATACACCGACAACATCCTCGATGACTACTGCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG
CGCCTACACCAACGATGTCCTGGACGACTTCTGCTACTACG"))
library(seqinr)
mat <- matrix(sapply(tmp, s2c), nrow = length(tmp), byrow = TRUE)
bma <- function(nucl){
	nucl <- tolower(nucl)
	iupac <- sapply(amb(), amb)
	iupac$u <- NULL # no RNA
	names(iupac)[unlist(lapply(iupac, setequal, nucl))]
}
res <- apply(mat, 2, bma)
toupper(c2s(res))
[1] "HGCMTACACCRACRAYRTCCTSGAYGACTWCWSCTACTACG"

Best,

Jean
-- 
Jean R. Lobry            (lobry at biomserv.univ-lyon1.fr)
Laboratoire BBE-CNRS-UMR-5558, Univ. C. Bernard - LYON I,
43 Bd 11/11/1918, F-69622 VILLEURBANNE CEDEX, FRANCE
allo  : +33 472 43 27 56     fax    : +33 472 43 13 88
http://pbil.univ-lyon1.fr/members/lobry/