#######################################################################
## ##
## Package: BatchMap ##
## ##
## File: group.R ##
## Contains: check.linkage, group, print.group ##
## ##
## Written by Gabriel Rodrigues Alves Margarido and Marcelo Mollinari ##
## copyright (c) 2007-9, Gabriel R A Margarido and Marcelo Mollinari ##
## ##
## First version: 11/07/2007 ##
## Last update: 01/14/2016 ##
## License: GNU General Public License version 2 (June, 1991) or later ##
## ##
#######################################################################
## Function to assign markers to linkage groups
##' Assign markers to linkage groups
##'
##' Identifies linkage groups of markers, using results from two-point
##' (pairwise) analysis and the \emph{transitive} property of linkage.
##'
##' If the arguments specifying thresholds used to group markers, i.e., minimum
##' LOD Score and maximum recombination fraction, are \code{NULL} (default),
##' the values used are those contained in object \code{input.seq}. If not
##' using \code{NULL}, the new values override the ones in object
##' \code{input.seq}.
##'
##' @aliases group
##' @param input.seq an object of class \code{sequence}.
##' @param LOD a (positive) real number used as minimum LOD score
##' (threshold) to declare linkage.
##' @param max.rf a real number (usually smaller than 0.5) used as
##' maximum recombination fraction to declare linkage.
##' @param verbose logical. If \code{TRUE}, current progress is shown;
##' if \code{FALSE}, no output is produced.
##' @return Returns an object of class \code{group}, which is a list
##' containing the following components: \item{data.name}{name of
##' the object of class \code{onemap} that contains the raw
##' data.} \item{twopt}{name of the object of class \code{rf.2ts}
##' used as input, i.e., containing information used to assign
##' markers to linkage groups.} \item{marnames}{marker names,
##' according to the input file.} \item{n.mar}{total number of
##' markers.} \item{LOD}{minimum LOD Score to declare linkage.}
##' \item{max.rf}{maximum recombination fraction to declare
##' linkage.} \item{n.groups}{number of linkage groups found.}
##' \item{groups}{number of the linkage group to which each marker
##' is assigned.}
##' @author Gabriel R A Margarido, \email{gramarga@@gmail.com} and
##' Marcelo Mollinari, \email{mmollina@@usp.br}
##' @seealso \code{\link[BatchMap]{rf.2pts}} and
##' \code{\link[BatchMap]{make.seq}}
##' @references Lincoln, S. E., Daly, M. J. and Lander, E. S. (1993)
##' Constructing genetic linkage maps with MAPMAKER/EXP Version
##' 3.0: a tutorial and reference manual. \emph{A Whitehead
##' Institute for Biomedical Research Technical Report}.
##' @keywords misc
##' @examples
##'
##' data(example.out)
##' twopts <- rf.2pts(example.out)
##'
##' all.data <- make.seq(twopts,"all")
##' link_gr <- group(all.data)
##' link_gr
##' print(link_gr, details=FALSE) #omit the names of the markers
##'
group <- function(input.seq, LOD=NULL, max.rf=NULL, verbose=TRUE)
{
## checking for correct object
if(!any(class(input.seq)=="sequence")) stop(deparse(substitute(input.seq))," is not an object of class 'sequence'")
## determining thresholds
if (is.null(LOD))
LOD <- get(input.seq$twopt, pos=1)$LOD
if (is.null(max.rf))
max.rf <- get(input.seq$twopt, pos=1)$max.rf
cl<-class(get(input.seq$data.name))[2]
geno<-get(input.seq$data.name)$geno[,input.seq$seq.num]
st<-get(input.seq$data.name)$segr.type.num
groups<-rep(0, length(input.seq$seq.num))
tp<-list(unlk=input.seq$seq.num)
i<-1
if(verbose) cat(" Selecting markers: \n")
while(length(tp$unlk) > 0)
{
g<-tp$unlk[1]
s<-tp$unlk
j<-1
tp<-check.linkage(i=g[j], s=s, cl=cl, geno=geno, st=st, max.rf=max.rf, LOD=LOD)
gt<-tp$lk
if(length(gt) > 0)
{
if(verbose) cat("\t group ", i,"\n\t ")
g<-c(g,gt)
while(!is.na(g[j])){
if(verbose)
{
cat(".")
if(j %% 60 == 0) cat("\n\t ")
}
tp<-check.linkage(i=g[j+1], s=tp$unlk, cl=cl, geno=geno, st=st, max.rf=max.rf, LOD=LOD)
gt<-tp$lk
g<-c(g,gt)
j<-j+1
}
if(verbose) cat("\n")
groups[g]<-i
i<-i+1
}
}
if(all(groups==0)) cat("\t No group found.\n")
## results
structure(list(data.name=input.seq$data.name, input.name=deparse(substitute(input.seq)),
twopt=input.seq$twopt, marnames=colnames(geno),
n.mar=length(input.seq$seq.num), seq.num=input.seq$seq.num, LOD=LOD, max.rf=max.rf,
n.groups=i-1, groups=groups), class = "group")
}
##' Print group result
##'
##' Generic print methods
##'
##' @param x The input object
##' @param detailed Print detailed output
##' @param ... Not used
##'
##' @method print group
print.group <-
function(x, detailed=TRUE,...) {
## checking for correct object
if(!any(class(x)=="group"))
stop(deparse(substitute(x))," is not an object of class 'group'")
cat(" This is an object of class 'group'\n")
cat(paste(" It was generated from the object \"", x$input.name,
"\"\n\n",sep=""))
## criteria
cat(" Criteria used to assign markers to groups:\n")
cat(" LOD =", x$LOD, ", Maximum recombination fraction =",
x$max.rf, "\n")
## printing summary
cat("\n No. markers: ", x$n.mar, "\n")
cat(" No. groups: ", x$n.groups, "\n")
cat(" No. linked markers: ", sum(!is.na(x$groups)), "\n")
cat(" No. unlinked markers: ", sum(x$groups==0), "\n")
if (detailed) {
## printing detailed results (markers in each linkage group)
cat("\n Printing groups:")
for (i in 1:x$n.groups) {
cat("\n Group", i, ":", length(which(x$groups==i)) , "markers\n ")
cat(x$marnames[x$seq.num[which(x$groups==i)]], "\n")
}
if (any(is.na(x$groups))) {
cat("\n Unlinked markers:", length(which(is.na(x$groups))) ," markers\n ")
cat(x$marnames[x$seq.num[which(is.na(x$groups))]], "\n")
}
}
}
##Checks if a marker i is linked with markers in a vector s
check.linkage<-function(i, s, cl, geno, st=NULL, max.rf, LOD)
{
s<-s[is.na(match(s,i))]
r<-est_rf_out(geno = geno[,c(i,s)], mrk = 1, seg_type = st[c(i,s)], nind = nrow(geno),verbose = FALSE)
sig<-apply(r[[1]], 2, function(x,y) min(x) <= y, y=max.rf) &
apply(r[[2]], 2, function(x,y) max(x) >= y, y=LOD)
return(list(lk=s[sig[-1]], unlk=s[!(sig[-1])]))
}
###end of file
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.