### R code from vignette source 'ClassifySequences.Rnw' ################################################### ### code chunk number 1: ClassifySequences.Rnw:47-52 ################################################### options(continue=" ") options(width=80) options(SweaveHooks=list(fig=function() par(mar=c(4.1, 4.1, 0.3, 0.1)))) set.seed(123) ################################################### ### code chunk number 2: startup ################################################### library(DECIPHER) ################################################### ### code chunk number 3: expr1 (eval = FALSE) ################################################### ## # specify the path to your file of training sequences: ## seqs_path <- "<>" ## # read the sequences into memory ## seqs <- readDNAStringSet(seqs_path) ## # NOTE: use readRNAStringSet for RNA sequences ## ## # (optionally) specify a path to the taxid file: ## rank_path <- "<>" ## taxid <- read.table(rank_path, ## header=FALSE, ## col.names=c('Index', 'Name', 'Parent', 'Level', 'Rank'), ## sep="*", # asterisks delimited ## quote="", # preserve quotes ## stringsAsFactors=FALSE) ## # OR, if no taxid text file exists, use: ## #taxid <- NULL ################################################### ### code chunk number 4: expr2 (eval = FALSE) ################################################### ## # if they exist, remove any gaps in the sequences: ## seqs <- RemoveGaps(seqs) ################################################### ### code chunk number 5: expr3 (eval = FALSE) ################################################### ## # obtain the taxonomic assignments ## groups <- names(seqs) # sequence names ## # assume the taxonomy begins with 'Root;' ## groups <- gsub("(.*)(Root;)", "\\2", groups) # extract the group label ## groupCounts <- table(groups) ## u_groups <- names(groupCounts) # unique groups ## length(u_groups) # number of groups ################################################### ### code chunk number 6: expr4 (eval = FALSE) ################################################### ## maxGroupSize <- 10 # max sequences per label (>= 1) ## ## remove <- logical(length(seqs)) ## for (i in which(groupCounts > maxGroupSize)) { ## index <- which(groups==u_groups[i]) ## keep <- sample(length(index), ## maxGroupSize) ## remove[index[-keep]] <- TRUE ## } ## sum(remove) # number of sequences eliminated ################################################### ### code chunk number 7: expr5 (eval = FALSE) ################################################### ## maxIterations <- 3 # must be >= 1 ## allowGroupRemoval <- FALSE ## ## probSeqsPrev <- integer() # suspected problem sequences from prior iteration ## for (i in seq_len(maxIterations)) { ## cat("Training iteration: ", i, "\n", sep="") ## ## # train the classifier ## trainingSet <- LearnTaxa(seqs[!remove], ## names(seqs)[!remove], ## taxid) ## ## # look for problem sequences ## probSeqs <- trainingSet$problemSequences$Index ## if (length(probSeqs)==0) { ## cat("No problem sequences remaining.\n") ## break ## } else if (length(probSeqs)==length(probSeqsPrev) && ## all(probSeqsPrev==probSeqs)) { ## cat("Iterations converged.\n") ## break ## } ## if (i==maxIterations) ## break ## probSeqsPrev <- probSeqs ## ## # remove any problem sequences ## index <- which(!remove)[probSeqs] ## remove[index] <- TRUE # remove all problem sequences ## if (!allowGroupRemoval) { ## # replace any removed groups ## missing <- !(u_groups %in% groups[!remove]) ## missing <- u_groups[missing] ## if (length(missing) > 0) { ## index <- index[groups[index] %in% missing] ## remove[index] <- FALSE # don't remove ## } ## } ## } ## sum(remove) # total number of sequences eliminated ## length(probSeqs) # number of remaining problem sequences ################################################### ### code chunk number 8: expr6 ################################################### data("TrainingSet_16S") trainingSet <- TrainingSet_16S ################################################### ### code chunk number 9: expr7 ################################################### trainingSet ################################################### ### code chunk number 10: expr8 ################################################### getOption("SweaveHooks")[["fig"]]() plot(trainingSet) ################################################### ### code chunk number 11: expr9 ################################################### fas <- "<>" # OR use the example 16S sequences: fas <- system.file("extdata", "Bacteria_175seqs.fas", package="DECIPHER") # read the sequences into memory test <- readDNAStringSet(fas) # NOTE: use readRNAStringSet for RNA sequences ################################################### ### code chunk number 12: expr10 ################################################### # if they exist, remove any gaps in the sequences: test <- RemoveGaps(test) test ################################################### ### code chunk number 13: expr11 ################################################### ids <- IdTaxa(test, trainingSet, type="extended", strand="top", threshold=60, processors=1) ################################################### ### code chunk number 14: expr12 ################################################### ids ################################################### ### code chunk number 15: expr12 ################################################### ids[1:5] # summary results for the first 5 sequences ids[[1]] # results for the first sequence ################################################### ### code chunk number 16: expr13 ################################################### getOption("SweaveHooks")[["fig"]]() plot(ids, trainingSet) ################################################### ### code chunk number 17: expr14 ################################################### output <- sapply(ids, function (id) { paste(id$taxon, " (", round(id$confidence, digits=1), "%)", sep="", collapse="; ") }) tail(output) #WriteLines(output, "<>") ################################################### ### code chunk number 18: expr15 (eval = FALSE) ################################################### ## set.seed(123) # choose a whole number as the random seed ## # then classify sequences with IdTaxa (not shown) ## set.seed(NULL) # return to the original state by unsetting the seed ################################################### ### code chunk number 19: sessinfo ################################################### toLatex(sessionInfo(), locale=FALSE)