R/ripple_seq.R
In augusto-garcia/onemap: Construction of Genetic Maps in Experimental Crosses
Defines functions
ripple_seq_inbred
ripple_seq_outcross
ripple_seq
Documented in
ripple_seq
########################################################################
## #
## Package: onemap #
## #
## File: ripple_seq.R #
## Contains: ripple_seq #
## #
## Written by Gabriel Rodrigues Alves Margarido #
## copyright (c) 2009, Gabriel R A Margarido #
## #
## First version: 02/27/2009 #
## License: GNU General Public License version 2 (June, 1991) or later #
## #
########################################################################
## This function searches for alternative orders, by comparing all possible
## orders of subsets of markers
##' Compares and displays plausible alternative orders for a given linkage
##' group
##'
##' For a given sequence of ordered markers, computes the multipoint likelihood
##' of alternative orders, by shuffling subsets (windows) of markers within the
##' sequence. For each position of the window, all possible \eqn{(ws)!}{(ws)!}
##' orders are compared.
##'
##' Large values for the window size make computations very slow, specially if
##' there are many partially informative markers.
##'
##'@importFrom methods is
##'@importFrom utils head tail
##'
##'
##' @param input.seq an object of class \code{sequence} with a
##' predefined order.
##' @param ws an integer specifying the length of the window size
##' (defaults to 4).
##' @param ext.w an integer specifying how many markers should be
##' considered in the vicinity of the permuted window. If
##' \code{ext.w=NULL} all markers in the sequence are
##' considered. In this version, it is used only in backcross,
##' \eqn{F_2}{F_2} or RIL crosses.
##' @param LOD threshold for the LOD-Score, so that alternative orders
##' with LOD less then or equal to this threshold will be
##' displayed.
##' @param tol tolerance for the C routine, i.e., the value used to
##' evaluate convergence.
##' @param verbose A logical, if TRUE it output progress status
##' information.
##'
##' @return This function does not return any value; it just produces
##' text output to suggest alternative orders.
##' @author Gabriel R A Margarido, \email{gramarga@@gmail.com} and
##' Marcelo Mollinari, \email{mmollina@@usp.br}
##' @seealso \code{\link[onemap]{make_seq}},
##' \code{\link[onemap]{compare}}, \code{\link[onemap]{try_seq}}
##' and \code{\link[onemap]{order_seq}}.
##' @references Broman, K. W., Wu, H., Churchill, G., Sen, S., Yandell, B.
##' (2008) \emph{qtl: Tools for analyzing QTL experiments} R package version
##' 1.09-43
##'
##' Jiang, C. and Zeng, Z.-B. (1997). Mapping quantitative trait loci with
##' dominant and missing markers in various crosses from two inbred lines.
##' \emph{Genetica} 101: 47-58.
##'
##' Lander, E. S., Green, P., Abrahamson, J., Barlow, A., Daly, M. J., Lincoln,
##' S. E. and Newburg, L. (1987) MAPMAKER: An interactive computer package for
##' constructing primary genetic linkage maps of experimental and natural
##' populations. \emph{Genomics} 1: 174-181.
##'
##' Mollinari, M., Margarido, G. R. A., Vencovsky, R. and Garcia, A. A. F.
##' (2009) Evaluation of algorithms used to order markers on genetics maps.
##' \emph{Heredity} 103: 494-502.
##'
##' Wu, R., Ma, C.-X., Painter, I. and Zeng, Z.-B. (2002a) Simultaneous maximum
##' likelihood estimation of linkage and linkage phases in outcrossing species.
##' \emph{Theoretical Population Biology} 61: 349-363.
##'
##' Wu, R., Ma, C.-X., Wu, S. S. and Zeng, Z.-B. (2002b). Linkage mapping of
##' sex-specific differences. \emph{Genetical Research} 79: 85-96
##' @keywords utilities
##' @examples
##'
##' \donttest{
##' #Outcross example
##' data(onemap_example_out)
##' twopt <- rf_2pts(onemap_example_out)
##' markers <- make_seq(twopt,c(27,16,20,4,19,21,23,9,24,29))
##' markers.map <- map(markers)
##' ripple_seq(markers.map)
##'
##' #F2 example
##' data(onemap_example_f2)
##' twopt <- rf_2pts(onemap_example_f2)
##' all_mark <- make_seq(twopt,"all")
##' groups <- group(all_mark)
##' LG3 <- make_seq(groups,1)
##' LG3.ord <- order_seq(LG3, subset.search = "twopt", twopt.alg = "rcd", touchdown=TRUE)
##' LG3.ord
##' make_seq(LG3.ord) # get safe sequence
##' ord.1<-make_seq(LG3.ord,"force") # get forced sequence
##' ripple_seq(ord.1, ws=5)
##' }
##'
##'@export
ripple_seq<-function(input.seq, ws=4, ext.w=NULL, LOD=3, tol=10E-2, verbose=TRUE)
{
if(inherits(input.seq$data.name, "outcross") || inherits(input.seq$data.name, "f2"))
return(ripple_seq_outcross(input.seq=input.seq,
ws=ws,
LOD=LOD,
tol=tol, verbose=TRUE))
else
return(ripple_seq_inbred(input.seq=input.seq,
ws=ws,
ext.w=ext.w,
LOD=LOD,
tol=tol, verbose=TRUE))
}
## This function searches for alternative orders, by comparing all possible
## orders of subsets of markers (for outcrosses)
ripple_seq_outcross<-function(input.seq,ws=4,LOD=3,tol=10E-2, verbose=TRUE) {
## checking for correct objects
if(!inherits(input.seq,"sequence")) stop(deparse(substitute(input.seq))," is not an object of class 'sequence'")
if(ws < 2) stop("ws must be greater than or equal to 2")
if(ws > 5) warning("this operation may take a VERY long time\n\n")
flush.console()
len <- length(input.seq$seq.num)
## computations unnecessary in this case
if (len <= ws) stop("Length of sequence ", deparse(substitute(input.seq))," is smaller than ws. You can use the compare function instead")
## convert numerical linkage phases to strings, to facilitate rearrangements
link.phases <- matrix(NA,len,2)
link.phases[1,] <- rep(1,2)
for (i in 1:length(input.seq$seq.phases)) {
switch(EXPR=input.seq$seq.phases[i],
link.phases[i+1,] <- link.phases[i,]*c(1,1),
link.phases[i+1,] <- link.phases[i,]*c(1,-1),
link.phases[i+1,] <- link.phases[i,]*c(-1,1),
link.phases[i+1,] <- link.phases[i,]*c(-1,-1),
)
}
## allocate variables
rf.init <- rep(NA,len-1)
phase <- rep(NA,len-1)
tot <- prod(1:ws)
best.ord.phase <- matrix(NA,tot,len-1)
best.ord.like <- best.ord.LOD <- rep(-Inf,tot)
## gather two-point information
list.init <- phases(input.seq)
#### first position
if(verbose) cat(input.seq$seq.num[1:ws],"|",input.seq$seq.num[ws+1], "...", sep="-")
## create all possible alternative orders for the first subset
all.ord <- t(apply(perm_tot(head(input.seq$seq.num,ws)),1,function(x) c(x,tail(input.seq$seq.num,-ws))))
for(i in 1:nrow(all.ord)){
all.match <- match(all.ord[i,],input.seq$seq.num)
## rearrange linkage phases according to the current order
for(j in 1:ws) {
temp <- paste(as.character(link.phases[all.match[j],]*link.phases[all.match[j+1],]),collapse=".")
phase[j] <- switch(EXPR=temp,
'1.1' = 1,
'1.-1' = 2,
'-1.1' = 3,
'-1.-1' = 4
)
}
if(len > (ws+1)) phase[(ws+1):length(phase)] <- tail(input.seq$seq.phase,-ws)
## get initial values for recombination fractions
for(j in 1:(len-1)){
if(all.match[j] > all.match[j+1]) ind <- acum(all.match[j]-2)+all.match[j+1]
else ind <- acum(all.match[j+1]-2)+all.match[j]
if(length(which(list.init$phase.init[[ind]] == phase[j])) == 0) rf.init[j] <- 0.49 ## for safety reasons
else rf.init[j] <- list.init$rf.init[[ind]][which(list.init$phase.init[[ind]] == phase[j])]
}
## estimate parameters
final.map <- est_map_hmm_out(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
type=input.seq$data.name$segr.type.num[all.ord[i,]],
phase=phase,
rf.vec=rf.init,
verbose=FALSE,
tol=tol)
best.ord.phase[i,] <- phase
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(best.ord.LOD > -LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.phase <- best.ord.phase[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose) {
cat("\n Alternative orders:\n")
for(j in which.LOD) {
if(inherits(input.seq$data.name,"outcross"))
cat(" ",all.ord[j,1:(ws+1)],ifelse(len > (ws+1),"... : ",": "),LOD.print[j],"( linkage phases:",best.ord.phase[j,1:ws],ifelse(len > (ws+1),"... )\n",")\n"))
else
cat(" ",all.ord[j,1:(ws+1)],ifelse(len > (ws+1),"... : ",": "),LOD.print[j],"\n")
}
cat("\n")
}
}
else if(verbose) cat(" OK\n\n")
#### middle positions
if (len > (ws+1)) {
for (p in 2:(len-ws)) {
if(verbose) cat("...", input.seq$seq.num[p-1], "|", input.seq$seq.num[p:(p+ws-1)],"|", input.seq$seq.num[p+ws],"...", sep="-")
## create all possible alternative orders for the first subset
all.ord <- t(apply(perm_tot(input.seq$seq.num[p:(p+ws-1)]),1,function(x) c(head(input.seq$seq.num,p-1),x,tail(input.seq$seq.num,-p-ws+1))))
for(i in 1:nrow(all.ord)){
all.match <- match(all.ord[i,],input.seq$seq.num)
## rearrange linkage phases according to the current order
if(p > 2) phase[1:(p-2)] <- head(input.seq$seq.phase,p-2)
for(j in (p-1):(p+ws-1)) {
temp <- paste(as.character(link.phases[all.match[j],]*link.phases[all.match[j+1],]),collapse=".")
phase[j] <- switch(EXPR=temp,
'1.1' = 1,
'1.-1' = 2,
'-1.1' = 3,
'-1.-1' = 4
)
}
if(p < (len-ws)) phase[(p+ws):length(phase)] <- tail(input.seq$seq.phase,len-p-ws)
## get initial values for recombination fractions
for(j in 1:(len-1)){
if(all.match[j] > all.match[j+1]) ind <- acum(all.match[j]-2)+all.match[j+1]
else ind <- acum(all.match[j+1]-2)+all.match[j]
if(length(which(list.init$phase.init[[ind]] == phase[j])) == 0) rf.init[j] <- 0.49 ## for safety reasons
else rf.init[j] <- list.init$rf.init[[ind]][which(list.init$phase.init[[ind]] == phase[j])]
}
## estimate parameters
final.map <- est_map_hmm_out(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
type=input.seq$data.name$segr.type.num[all.ord[i,]],
phase=phase,
rf.vec=rf.init,
verbose=FALSE,
tol=tol)
best.ord.phase[i,] <- phase
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(best.ord.LOD > -LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.phase <- best.ord.phase[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose){
cat("\n Alternative orders:\n")
for(j in which.LOD) {
if(inherits(input.seq$data.name,"outcross"))
cat(ifelse(p>2," ..."," "),all.ord[j,(p-1):(p+ws)],ifelse((p+ws)<len,"... : ",": "),LOD.print[j],"( linkage phases:",ifelse(p>2,"...","\b"),best.ord.phase[j,(p-1):(p+ws-1)],ifelse((p+ws)<len,"... )\n",")\n"))
else
cat(ifelse(p>2," ..."," "),all.ord[j,(p-1):(p+ws)],ifelse((p+ws)<len,"... : ",": "),LOD.print[j],"\n")
}
cat("\n")
}
}
else if(verbose) cat(" OK\n\n")
}
}
###### last position
if(verbose) cat(input.seq$seq.num[len-ws], "|" , tail(input.seq$seq.num,ws), sep="-")
## create all possible alternative orders for the first subset
all.ord <- t(apply(perm_tot(tail(input.seq$seq.num,ws)),1,function(x) c(head(input.seq$seq.num,-ws),x)))
for(i in 1:nrow(all.ord)){
all.match <- match(all.ord[i,],input.seq$seq.num)
## rearrange linkage phases according to the current order
if(len > (ws+1)) phase[1:(len-ws-1)] <- head(input.seq$seq.phase,-ws)
for(j in (len-ws):(len-1)) {
temp <- paste(as.character(link.phases[all.match[j],]*link.phases[all.match[j+1],]),collapse=".")
phase[j] <- switch(EXPR=temp,
'1.1' = 1,
'1.-1' = 2,
'-1.1' = 3,
'-1.-1' = 4
)
}
## get initial values for recombination fractions
for(j in 1:(len-1)){
if(all.match[j] > all.match[j+1]) ind <- acum(all.match[j]-2)+all.match[j+1]
else ind <- acum(all.match[j+1]-2)+all.match[j]
if(length(which(list.init$phase.init[[ind]] == phase[j])) == 0) rf.init[j] <- 0.49 ## for safety reasons
else rf.init[j] <- list.init$rf.init[[ind]][which(list.init$phase.init[[ind]] == phase[j])]
}
## estimate parameters
final.map <- est_map_hmm_out(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
type=input.seq$data.name$segr.type.num[all.ord[i,]],
phase=phase,
rf.vec=rf.init,
verbose=FALSE,
tol=tol)
best.ord.phase[i,] <- phase
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(best.ord.LOD > -LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.phase <- best.ord.phase[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose) {
cat("\n Alternative orders:\n")
for(j in which.LOD) {
if(inherits(input.seq$data.name,"outcross"))
cat(ifelse(len > (ws+1)," ..."," "),all.ord[j,(len-ws):len],": ",LOD.print[j],"( linkage phases:",ifelse(len > (ws+1),"...","\b"),best.ord.phase[j,(len-ws):(len-1)],")\n")
else
cat(ifelse(len > (ws+1)," ..."," "),all.ord[j,(len-ws):len],": ",LOD.print[j],"\n")
}
cat("\n")
}
}
else if(verbose) cat(" OK\n\n")
}
## This function searches for alternative orders, by comparing all possible
## orders of subsets of markers (for crosses derived from inbred lines)
ripple_seq_inbred<-function(input.seq, ws=4, ext.w=NULL, LOD=3, tol=10E-2, verbose=TRUE)
{
## checking for correct objects
if(!inherits(input.seq,"sequence"))
stop(deparse(substitute(input.seq))," is not an object of class 'sequence'")
if(ws < 2)
stop("ws must be greater than or equal to 2")
if(ws > 5)
warning("This operation may take a VERY long time\n\n")
flush.console()
len <- length(input.seq$seq.num)
## computations unnecessary in this case
if (len <= ws)
stop("Length of sequence ", deparse(substitute(input.seq))," is smaller than ws. You can use the 'compare' function instead")
## allocate variables
rf.init <- rep(NA,len-1)
tot <- prod(1:ws)
best.ord.like <- best.ord.LOD <- rep(-Inf,tot)
## first position
if(verbose) cat(input.seq$seq.num[1:ws],"|",input.seq$seq.num[ws+1], "...", sep="-")
## create all possible alternative orders for the first subset
if(is.null(ext.w) || ext.w >= (length(input.seq$seq.num)-ws))
all.ord <- t(apply(perm_tot(head(input.seq$seq.num,ws)),1,function(x) c(x,tail(input.seq$seq.num,-ws))))
else
all.ord <- t(apply(perm_tot(head(input.seq$seq.num,ws)),1,function(x) c(x,input.seq$seq.num[(ws+1):(ws+1+ext.w)])))
for(i in 1:nrow(all.ord)){
## estimate parameters
seq.temp<-make_seq(input.seq$twopt, arg=all.ord[i,])
seq.temp$twopt<-input.seq$twopt
rf.temp<-get_vec_rf_in(seq.temp, acum=FALSE)
final.map<-est_map_hmm_bc(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
rf.vec=rf.temp,
verbose=FALSE,
tol=tol)
if(inherits(input.seq$data.name, c("riself", "risib"))){
crosstype <- ifelse(inherits(input.seq$data.name, "riself"), "riself", "risib")
final.map$rf<-adjust_rf_ril(final.map$rf,
type=crosstype,
expand = FALSE)
}
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(abs(best.ord.LOD) < LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose){
cat("\n Alternative orders:\n")
for(j in which.LOD)
cat(" ",all.ord[j,1:(ws)], ifelse(len > (ws+1),"... : ",": "), LOD.print[j],"\n")
cat("\n")
}
}
else
if(verbose) cat(" OK\n\n")
## middle positions
if (len > (ws+1)) {
for (p in 2:(len-ws)) {
if(verbose) cat("...", input.seq$seq.num[p-1], "|", input.seq$seq.num[p:(p+ws-1)],"|", input.seq$seq.num[p+ws],"...", sep="-")
## create all possible alternative orders for the first subset
if(is.null(ext.w) || ext.w >= (length(input.seq$seq.num)-ws))
all.ord <- t(apply(perm_tot(input.seq$seq.num[p:(p+ws-1)]),1,function(x) c(head(input.seq$seq.num,p-1),x,tail(input.seq$seq.num,-p-ws+1))))
else{
x0<-(p-ext.w):(p-1)
x1<-(p+ws):((p+ws)+ext.w)
all.ord <- t(apply(perm_tot(input.seq$seq.num[p:(p+ws-1)]),1,function(x) c(input.seq$seq.num[x0[x0>0]],x,input.seq$seq.num[x1[x1 <= len]])))
}
for(i in 1:nrow(all.ord)){
## estimate parameters
seq.temp<-make_seq(input.seq$twopt, arg=all.ord[i,])
seq.temp$twopt<-input.seq$twopt
rf.temp<-get_vec_rf_in(seq.temp, acum=FALSE)
final.map<-est_map_hmm_bc(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
rf.vec=rf.temp,
verbose=FALSE,
tol=tol)
if(inherits(input.seq$data.name, c("riself", "risib"))){
crosstype <- ifelse(inherits(input.seq$data.name, "riself"), "riself", "risib")
final.map$rf<-adjust_rf_ril(final.map$rf,
type=crosstype,
expand = FALSE)
}
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(best.ord.LOD > -LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD, digits=2, nsmall=2)
if(verbose) {
cat("\n Alternative orders:\n")
for(j in which.LOD)
{
fin<-which(names(all.ord[j,])=="M")
fin
cat(ifelse(p>2," ..."," "),all.ord[j,(fin-ws+1):fin],ifelse((p+ws)<len,"... : ",": "),LOD.print[j],"\n")
}
cat("\n")
}
}
else if(verbose) cat(" OK\n\n")
}
}
## last position
if(verbose) cat(input.seq$seq.num[len-ws], "|" , tail(input.seq$seq.num,ws), sep="-")
## create all possible alternative orders for the first subset
if(is.null(ext.w) || ext.w >= (length(input.seq$seq.num)-ws))
all.ord <- t(apply(perm_tot(tail(input.seq$seq.num,ws)),1,function(x) c(head(input.seq$seq.num,-ws),x)))
else
all.ord <- t(apply(perm_tot(tail(input.seq$seq.num,ws)),1,function(x) c(input.seq$seq.num[(len-ws-ext.w+1):(len-ws)],x)))
for(i in 1:nrow(all.ord)){
## estimate parameters
seq.temp<-make_seq(input.seq$twopt, arg=all.ord[i,])
seq.temp$twopt<-input.seq$twopt
rf.temp<-get_vec_rf_in(seq.temp, acum=FALSE)
final.map<-est_map_hmm_bc(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
rf.vec=rf.temp,
verbose=FALSE,
tol=tol)
if(inherits(input.seq$data.name, c("riself", "risib"))){
crosstype <- ifelse(inherits(input.seq$data.name, "riself"), "riself", "risib")
final.map$rf<-adjust_rf_ril(final.map$rf,
type=crosstype,
expand = FALSE)
}
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(abs(best.ord.LOD) < LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose) {
cat("\n Alternative orders:\n")
for(j in which.LOD)
cat(" ",all.ord[j,(ncol(all.ord)-ws+1):ncol(all.ord)],ifelse(len > (ws+1),"... : ",": "),LOD.print[j],"\n")
cat("\n")
}
}
else
if(verbose) cat(" OK\n\n")
}
# end of file
augusto-garcia/onemap documentation built on Nov. 30, 2022, 9:13 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions ripple_seq_inbred ripple_seq_outcross ripple_seq
Documented in ripple_seq
########################################################################
## #
## Package: onemap #
## #
## File: ripple_seq.R #
## Contains: ripple_seq #
## #
## Written by Gabriel Rodrigues Alves Margarido #
## copyright (c) 2009, Gabriel R A Margarido #
## #
## First version: 02/27/2009 #
## License: GNU General Public License version 2 (June, 1991) or later #
## #
########################################################################
## This function searches for alternative orders, by comparing all possible
## orders of subsets of markers
##' Compares and displays plausible alternative orders for a given linkage
##' group
##'
##' For a given sequence of ordered markers, computes the multipoint likelihood
##' of alternative orders, by shuffling subsets (windows) of markers within the
##' sequence. For each position of the window, all possible \eqn{(ws)!}{(ws)!}
##' orders are compared.
##'
##' Large values for the window size make computations very slow, specially if
##' there are many partially informative markers.
##'
##'@importFrom methods is
##'@importFrom utils head tail
##'
##'
##' @param input.seq an object of class \code{sequence} with a
##' predefined order.
##' @param ws an integer specifying the length of the window size
##' (defaults to 4).
##' @param ext.w an integer specifying how many markers should be
##' considered in the vicinity of the permuted window. If
##' \code{ext.w=NULL} all markers in the sequence are
##' considered. In this version, it is used only in backcross,
##' \eqn{F_2}{F_2} or RIL crosses.
##' @param LOD threshold for the LOD-Score, so that alternative orders
##' with LOD less then or equal to this threshold will be
##' displayed.
##' @param tol tolerance for the C routine, i.e., the value used to
##' evaluate convergence.
##' @param verbose A logical, if TRUE it output progress status
##' information.
##'
##' @return This function does not return any value; it just produces
##' text output to suggest alternative orders.
##' @author Gabriel R A Margarido, \email{gramarga@@gmail.com} and
##' Marcelo Mollinari, \email{mmollina@@usp.br}
##' @seealso \code{\link[onemap]{make_seq}},
##' \code{\link[onemap]{compare}}, \code{\link[onemap]{try_seq}}
##' and \code{\link[onemap]{order_seq}}.
##' @references Broman, K. W., Wu, H., Churchill, G., Sen, S., Yandell, B.
##' (2008) \emph{qtl: Tools for analyzing QTL experiments} R package version
##' 1.09-43
##'
##' Jiang, C. and Zeng, Z.-B. (1997). Mapping quantitative trait loci with
##' dominant and missing markers in various crosses from two inbred lines.
##' \emph{Genetica} 101: 47-58.
##'
##' Lander, E. S., Green, P., Abrahamson, J., Barlow, A., Daly, M. J., Lincoln,
##' S. E. and Newburg, L. (1987) MAPMAKER: An interactive computer package for
##' constructing primary genetic linkage maps of experimental and natural
##' populations. \emph{Genomics} 1: 174-181.
##'
##' Mollinari, M., Margarido, G. R. A., Vencovsky, R. and Garcia, A. A. F.
##' (2009) Evaluation of algorithms used to order markers on genetics maps.
##' \emph{Heredity} 103: 494-502.
##'
##' Wu, R., Ma, C.-X., Painter, I. and Zeng, Z.-B. (2002a) Simultaneous maximum
##' likelihood estimation of linkage and linkage phases in outcrossing species.
##' \emph{Theoretical Population Biology} 61: 349-363.
##'
##' Wu, R., Ma, C.-X., Wu, S. S. and Zeng, Z.-B. (2002b). Linkage mapping of
##' sex-specific differences. \emph{Genetical Research} 79: 85-96
##' @keywords utilities
##' @examples
##'
##' \donttest{
##' #Outcross example
##' data(onemap_example_out)
##' twopt <- rf_2pts(onemap_example_out)
##' markers <- make_seq(twopt,c(27,16,20,4,19,21,23,9,24,29))
##' markers.map <- map(markers)
##' ripple_seq(markers.map)
##'
##' #F2 example
##' data(onemap_example_f2)
##' twopt <- rf_2pts(onemap_example_f2)
##' all_mark <- make_seq(twopt,"all")
##' groups <- group(all_mark)
##' LG3 <- make_seq(groups,1)
##' LG3.ord <- order_seq(LG3, subset.search = "twopt", twopt.alg = "rcd", touchdown=TRUE)
##' LG3.ord
##' make_seq(LG3.ord) # get safe sequence
##' ord.1<-make_seq(LG3.ord,"force") # get forced sequence
##' ripple_seq(ord.1, ws=5)
##' }
##'
##'@export
ripple_seq<-function(input.seq, ws=4, ext.w=NULL, LOD=3, tol=10E-2, verbose=TRUE)
{
if(inherits(input.seq$data.name, "outcross") || inherits(input.seq$data.name, "f2"))
return(ripple_seq_outcross(input.seq=input.seq,
ws=ws,
LOD=LOD,
tol=tol, verbose=TRUE))
else
return(ripple_seq_inbred(input.seq=input.seq,
ws=ws,
ext.w=ext.w,
LOD=LOD,
tol=tol, verbose=TRUE))
}
## This function searches for alternative orders, by comparing all possible
## orders of subsets of markers (for outcrosses)
ripple_seq_outcross<-function(input.seq,ws=4,LOD=3,tol=10E-2, verbose=TRUE) {
## checking for correct objects
if(!inherits(input.seq,"sequence")) stop(deparse(substitute(input.seq))," is not an object of class 'sequence'")
if(ws < 2) stop("ws must be greater than or equal to 2")
if(ws > 5) warning("this operation may take a VERY long time\n\n")
flush.console()
len <- length(input.seq$seq.num)
## computations unnecessary in this case
if (len <= ws) stop("Length of sequence ", deparse(substitute(input.seq))," is smaller than ws. You can use the compare function instead")
## convert numerical linkage phases to strings, to facilitate rearrangements
link.phases <- matrix(NA,len,2)
link.phases[1,] <- rep(1,2)
for (i in 1:length(input.seq$seq.phases)) {
switch(EXPR=input.seq$seq.phases[i],
link.phases[i+1,] <- link.phases[i,]*c(1,1),
link.phases[i+1,] <- link.phases[i,]*c(1,-1),
link.phases[i+1,] <- link.phases[i,]*c(-1,1),
link.phases[i+1,] <- link.phases[i,]*c(-1,-1),
)
}
## allocate variables
rf.init <- rep(NA,len-1)
phase <- rep(NA,len-1)
tot <- prod(1:ws)
best.ord.phase <- matrix(NA,tot,len-1)
best.ord.like <- best.ord.LOD <- rep(-Inf,tot)
## gather two-point information
list.init <- phases(input.seq)
#### first position
if(verbose) cat(input.seq$seq.num[1:ws],"|",input.seq$seq.num[ws+1], "...", sep="-")
## create all possible alternative orders for the first subset
all.ord <- t(apply(perm_tot(head(input.seq$seq.num,ws)),1,function(x) c(x,tail(input.seq$seq.num,-ws))))
for(i in 1:nrow(all.ord)){
all.match <- match(all.ord[i,],input.seq$seq.num)
## rearrange linkage phases according to the current order
for(j in 1:ws) {
temp <- paste(as.character(link.phases[all.match[j],]*link.phases[all.match[j+1],]),collapse=".")
phase[j] <- switch(EXPR=temp,
'1.1' = 1,
'1.-1' = 2,
'-1.1' = 3,
'-1.-1' = 4
)
}
if(len > (ws+1)) phase[(ws+1):length(phase)] <- tail(input.seq$seq.phase,-ws)
## get initial values for recombination fractions
for(j in 1:(len-1)){
if(all.match[j] > all.match[j+1]) ind <- acum(all.match[j]-2)+all.match[j+1]
else ind <- acum(all.match[j+1]-2)+all.match[j]
if(length(which(list.init$phase.init[[ind]] == phase[j])) == 0) rf.init[j] <- 0.49 ## for safety reasons
else rf.init[j] <- list.init$rf.init[[ind]][which(list.init$phase.init[[ind]] == phase[j])]
}
## estimate parameters
final.map <- est_map_hmm_out(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
type=input.seq$data.name$segr.type.num[all.ord[i,]],
phase=phase,
rf.vec=rf.init,
verbose=FALSE,
tol=tol)
best.ord.phase[i,] <- phase
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(best.ord.LOD > -LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.phase <- best.ord.phase[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose) {
cat("\n Alternative orders:\n")
for(j in which.LOD) {
if(inherits(input.seq$data.name,"outcross"))
cat(" ",all.ord[j,1:(ws+1)],ifelse(len > (ws+1),"... : ",": "),LOD.print[j],"( linkage phases:",best.ord.phase[j,1:ws],ifelse(len > (ws+1),"... )\n",")\n"))
else
cat(" ",all.ord[j,1:(ws+1)],ifelse(len > (ws+1),"... : ",": "),LOD.print[j],"\n")
}
cat("\n")
}
}
else if(verbose) cat(" OK\n\n")
#### middle positions
if (len > (ws+1)) {
for (p in 2:(len-ws)) {
if(verbose) cat("...", input.seq$seq.num[p-1], "|", input.seq$seq.num[p:(p+ws-1)],"|", input.seq$seq.num[p+ws],"...", sep="-")
## create all possible alternative orders for the first subset
all.ord <- t(apply(perm_tot(input.seq$seq.num[p:(p+ws-1)]),1,function(x) c(head(input.seq$seq.num,p-1),x,tail(input.seq$seq.num,-p-ws+1))))
for(i in 1:nrow(all.ord)){
all.match <- match(all.ord[i,],input.seq$seq.num)
## rearrange linkage phases according to the current order
if(p > 2) phase[1:(p-2)] <- head(input.seq$seq.phase,p-2)
for(j in (p-1):(p+ws-1)) {
temp <- paste(as.character(link.phases[all.match[j],]*link.phases[all.match[j+1],]),collapse=".")
phase[j] <- switch(EXPR=temp,
'1.1' = 1,
'1.-1' = 2,
'-1.1' = 3,
'-1.-1' = 4
)
}
if(p < (len-ws)) phase[(p+ws):length(phase)] <- tail(input.seq$seq.phase,len-p-ws)
## get initial values for recombination fractions
for(j in 1:(len-1)){
if(all.match[j] > all.match[j+1]) ind <- acum(all.match[j]-2)+all.match[j+1]
else ind <- acum(all.match[j+1]-2)+all.match[j]
if(length(which(list.init$phase.init[[ind]] == phase[j])) == 0) rf.init[j] <- 0.49 ## for safety reasons
else rf.init[j] <- list.init$rf.init[[ind]][which(list.init$phase.init[[ind]] == phase[j])]
}
## estimate parameters
final.map <- est_map_hmm_out(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
type=input.seq$data.name$segr.type.num[all.ord[i,]],
phase=phase,
rf.vec=rf.init,
verbose=FALSE,
tol=tol)
best.ord.phase[i,] <- phase
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(best.ord.LOD > -LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.phase <- best.ord.phase[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose){
cat("\n Alternative orders:\n")
for(j in which.LOD) {
if(inherits(input.seq$data.name,"outcross"))
cat(ifelse(p>2," ..."," "),all.ord[j,(p-1):(p+ws)],ifelse((p+ws)<len,"... : ",": "),LOD.print[j],"( linkage phases:",ifelse(p>2,"...","\b"),best.ord.phase[j,(p-1):(p+ws-1)],ifelse((p+ws)<len,"... )\n",")\n"))
else
cat(ifelse(p>2," ..."," "),all.ord[j,(p-1):(p+ws)],ifelse((p+ws)<len,"... : ",": "),LOD.print[j],"\n")
}
cat("\n")
}
}
else if(verbose) cat(" OK\n\n")
}
}
###### last position
if(verbose) cat(input.seq$seq.num[len-ws], "|" , tail(input.seq$seq.num,ws), sep="-")
## create all possible alternative orders for the first subset
all.ord <- t(apply(perm_tot(tail(input.seq$seq.num,ws)),1,function(x) c(head(input.seq$seq.num,-ws),x)))
for(i in 1:nrow(all.ord)){
all.match <- match(all.ord[i,],input.seq$seq.num)
## rearrange linkage phases according to the current order
if(len > (ws+1)) phase[1:(len-ws-1)] <- head(input.seq$seq.phase,-ws)
for(j in (len-ws):(len-1)) {
temp <- paste(as.character(link.phases[all.match[j],]*link.phases[all.match[j+1],]),collapse=".")
phase[j] <- switch(EXPR=temp,
'1.1' = 1,
'1.-1' = 2,
'-1.1' = 3,
'-1.-1' = 4
)
}
## get initial values for recombination fractions
for(j in 1:(len-1)){
if(all.match[j] > all.match[j+1]) ind <- acum(all.match[j]-2)+all.match[j+1]
else ind <- acum(all.match[j+1]-2)+all.match[j]
if(length(which(list.init$phase.init[[ind]] == phase[j])) == 0) rf.init[j] <- 0.49 ## for safety reasons
else rf.init[j] <- list.init$rf.init[[ind]][which(list.init$phase.init[[ind]] == phase[j])]
}
## estimate parameters
final.map <- est_map_hmm_out(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
type=input.seq$data.name$segr.type.num[all.ord[i,]],
phase=phase,
rf.vec=rf.init,
verbose=FALSE,
tol=tol)
best.ord.phase[i,] <- phase
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(best.ord.LOD > -LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.phase <- best.ord.phase[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose) {
cat("\n Alternative orders:\n")
for(j in which.LOD) {
if(inherits(input.seq$data.name,"outcross"))
cat(ifelse(len > (ws+1)," ..."," "),all.ord[j,(len-ws):len],": ",LOD.print[j],"( linkage phases:",ifelse(len > (ws+1),"...","\b"),best.ord.phase[j,(len-ws):(len-1)],")\n")
else
cat(ifelse(len > (ws+1)," ..."," "),all.ord[j,(len-ws):len],": ",LOD.print[j],"\n")
}
cat("\n")
}
}
else if(verbose) cat(" OK\n\n")
}
## This function searches for alternative orders, by comparing all possible
## orders of subsets of markers (for crosses derived from inbred lines)
ripple_seq_inbred<-function(input.seq, ws=4, ext.w=NULL, LOD=3, tol=10E-2, verbose=TRUE)
{
## checking for correct objects
if(!inherits(input.seq,"sequence"))
stop(deparse(substitute(input.seq))," is not an object of class 'sequence'")
if(ws < 2)
stop("ws must be greater than or equal to 2")
if(ws > 5)
warning("This operation may take a VERY long time\n\n")
flush.console()
len <- length(input.seq$seq.num)
## computations unnecessary in this case
if (len <= ws)
stop("Length of sequence ", deparse(substitute(input.seq))," is smaller than ws. You can use the 'compare' function instead")
## allocate variables
rf.init <- rep(NA,len-1)
tot <- prod(1:ws)
best.ord.like <- best.ord.LOD <- rep(-Inf,tot)
## first position
if(verbose) cat(input.seq$seq.num[1:ws],"|",input.seq$seq.num[ws+1], "...", sep="-")
## create all possible alternative orders for the first subset
if(is.null(ext.w) || ext.w >= (length(input.seq$seq.num)-ws))
all.ord <- t(apply(perm_tot(head(input.seq$seq.num,ws)),1,function(x) c(x,tail(input.seq$seq.num,-ws))))
else
all.ord <- t(apply(perm_tot(head(input.seq$seq.num,ws)),1,function(x) c(x,input.seq$seq.num[(ws+1):(ws+1+ext.w)])))
for(i in 1:nrow(all.ord)){
## estimate parameters
seq.temp<-make_seq(input.seq$twopt, arg=all.ord[i,])
seq.temp$twopt<-input.seq$twopt
rf.temp<-get_vec_rf_in(seq.temp, acum=FALSE)
final.map<-est_map_hmm_bc(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
rf.vec=rf.temp,
verbose=FALSE,
tol=tol)
if(inherits(input.seq$data.name, c("riself", "risib"))){
crosstype <- ifelse(inherits(input.seq$data.name, "riself"), "riself", "risib")
final.map$rf<-adjust_rf_ril(final.map$rf,
type=crosstype,
expand = FALSE)
}
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(abs(best.ord.LOD) < LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose){
cat("\n Alternative orders:\n")
for(j in which.LOD)
cat(" ",all.ord[j,1:(ws)], ifelse(len > (ws+1),"... : ",": "), LOD.print[j],"\n")
cat("\n")
}
}
else
if(verbose) cat(" OK\n\n")
## middle positions
if (len > (ws+1)) {
for (p in 2:(len-ws)) {
if(verbose) cat("...", input.seq$seq.num[p-1], "|", input.seq$seq.num[p:(p+ws-1)],"|", input.seq$seq.num[p+ws],"...", sep="-")
## create all possible alternative orders for the first subset
if(is.null(ext.w) || ext.w >= (length(input.seq$seq.num)-ws))
all.ord <- t(apply(perm_tot(input.seq$seq.num[p:(p+ws-1)]),1,function(x) c(head(input.seq$seq.num,p-1),x,tail(input.seq$seq.num,-p-ws+1))))
else{
x0<-(p-ext.w):(p-1)
x1<-(p+ws):((p+ws)+ext.w)
all.ord <- t(apply(perm_tot(input.seq$seq.num[p:(p+ws-1)]),1,function(x) c(input.seq$seq.num[x0[x0>0]],x,input.seq$seq.num[x1[x1 <= len]])))
}
for(i in 1:nrow(all.ord)){
## estimate parameters
seq.temp<-make_seq(input.seq$twopt, arg=all.ord[i,])
seq.temp$twopt<-input.seq$twopt
rf.temp<-get_vec_rf_in(seq.temp, acum=FALSE)
final.map<-est_map_hmm_bc(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
rf.vec=rf.temp,
verbose=FALSE,
tol=tol)
if(inherits(input.seq$data.name, c("riself", "risib"))){
crosstype <- ifelse(inherits(input.seq$data.name, "riself"), "riself", "risib")
final.map$rf<-adjust_rf_ril(final.map$rf,
type=crosstype,
expand = FALSE)
}
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(best.ord.LOD > -LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD, digits=2, nsmall=2)
if(verbose) {
cat("\n Alternative orders:\n")
for(j in which.LOD)
{
fin<-which(names(all.ord[j,])=="M")
fin
cat(ifelse(p>2," ..."," "),all.ord[j,(fin-ws+1):fin],ifelse((p+ws)<len,"... : ",": "),LOD.print[j],"\n")
}
cat("\n")
}
}
else if(verbose) cat(" OK\n\n")
}
}
## last position
if(verbose) cat(input.seq$seq.num[len-ws], "|" , tail(input.seq$seq.num,ws), sep="-")
## create all possible alternative orders for the first subset
if(is.null(ext.w) || ext.w >= (length(input.seq$seq.num)-ws))
all.ord <- t(apply(perm_tot(tail(input.seq$seq.num,ws)),1,function(x) c(head(input.seq$seq.num,-ws),x)))
else
all.ord <- t(apply(perm_tot(tail(input.seq$seq.num,ws)),1,function(x) c(input.seq$seq.num[(len-ws-ext.w+1):(len-ws)],x)))
for(i in 1:nrow(all.ord)){
## estimate parameters
seq.temp<-make_seq(input.seq$twopt, arg=all.ord[i,])
seq.temp$twopt<-input.seq$twopt
rf.temp<-get_vec_rf_in(seq.temp, acum=FALSE)
final.map<-est_map_hmm_bc(geno=t(input.seq$data.name$geno[,all.ord[i,]]),
error=input.seq$data.name$error[all.ord[i,] +
rep(c(0:(input.seq$data.name$n.ind-1))*input.seq$data.name$n.mar,
each=length(all.ord[i,])),],
rf.vec=rf.temp,
verbose=FALSE,
tol=tol)
if(inherits(input.seq$data.name, c("riself", "risib"))){
crosstype <- ifelse(inherits(input.seq$data.name, "riself"), "riself", "risib")
final.map$rf<-adjust_rf_ril(final.map$rf,
type=crosstype,
expand = FALSE)
}
best.ord.like[i] <- final.map$loglike
}
## calculate LOD-Scores for alternative orders
best.ord.LOD <- round((best.ord.like-max(best.ord.like))/log(10),2)
## which orders will be printed
which.LOD <- which(abs(best.ord.LOD) < LOD)
if(length(which.LOD) > 1) {
## if any order to print, sort by LOD-Score
order.print <- order(best.ord.LOD,decreasing=TRUE)
all.ord <- all.ord[order.print,]
best.ord.LOD <- best.ord.LOD[order.print]
## display results
which.LOD <- which(best.ord.LOD > -LOD)
LOD.print <- format(best.ord.LOD,digits=2,nsmall=2)
if(verbose) {
cat("\n Alternative orders:\n")
for(j in which.LOD)
cat(" ",all.ord[j,(ncol(all.ord)-ws+1):ncol(all.ord)],ifelse(len > (ws+1),"... : ",": "),LOD.print[j],"\n")
cat("\n")
}
}
else
if(verbose) cat(" OK\n\n")
}
# end of file
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.