### R code from vignette source 'vignettes/MinimumDistance/inst/doc/MinimumDistance.Rnw'

###################################################
### code chunk number 1: setup
###################################################
options(prompt="R> ", continue=" ", device=pdf, width=65)


###################################################
### code chunk number 2: enableLD (eval = FALSE)
###################################################
## library(snow)
## library(doSNOW)
## library(oligoClasses)
## cl <- makeCluster(22, type="SOCK")
## registerDoSNOW(cl)
## parStatus()


###################################################
### code chunk number 3: registerBackend
###################################################
library(oligoClasses)
library(MinimumDistance)
foreach::registerDoSEQ()


###################################################
### code chunk number 4: filenames
###################################################
path <- system.file("extdata", package="MinimumDistance")
fnames <- list.files(path, pattern=".txt")
fnames


###################################################
### code chunk number 5: pedigreeInfo
###################################################
pedigreeInfo <- data.frame(F="F.txt", M="M.txt", O="O.txt")


###################################################
### code chunk number 6: pedigreeConstructor
###################################################
ped <- Pedigree(pedigreeInfo)
ped


###################################################
### code chunk number 7: callDenovoSegments
###################################################
map.segs <- callDenovoSegments(path=path,
			       pedigreeData=ped,
			       cdfname="human610quadv1b",
			       chromosome=1,
			       segmentParents=FALSE,
			       genome="hg18")


###################################################
### code chunk number 8: constructSampleSheet
###################################################
library(human610quadv1bCrlmm)
path <- system.file("extdata", package="MinimumDistance")
load(file.path(path, "pedigreeInfo.rda"))
load(file.path(path, "sample.sheet.rda"))
load(file.path(path, "logRratio.rda"))
load(file.path(path, "baf.rda"))
stopifnot(colnames(logRratio) %in% allNames(Pedigree(pedigreeInfo)))
nms <- paste("NA",substr(sample.sheet$Sample.Name, 6, 10),sep="")
trioSetList <- TrioSetList(lrr=logRratio, ## must provide row.names
			   baf=baf,
			   pedigree=Pedigree(pedigreeInfo),
			   sample.sheet=sample.sheet,
			   row.names=nms,
			   cdfname="human610quadv1bCrlmm",
			   genome="hg18")
show(trioSetList)


###################################################
### code chunk number 9: loadTrioSetListExample
###################################################
data(trioSetListExample)


###################################################
### code chunk number 10: computeMinimumDistance
###################################################
md <- calculateMindist(lrr(trioSetList))


###################################################
### code chunk number 11: segmentMinimumDistance
###################################################
md.segs <- segment2(object=trioSetList, md=md)


###################################################
### code chunk number 12: showMd.Segs
###################################################
head(md.segs)


###################################################
### code chunk number 13: segmentLRR
###################################################
lrr.segs <- segment2(trioSetList, segmentParents=TRUE)


###################################################
### code chunk number 14: narrow
###################################################
mads.md <- mad2(md, byrow=FALSE)
md.segs2 <- narrow(md.segs, lrr.segs, thr=0.75, mad.minimumdistance=mads.md, fD=featureData(trioSetList), genome="hg18")


###################################################
### code chunk number 15: computeBayesFactor
###################################################
##trace(MinimumDistance:::joint4, browser)
map.segs <- computeBayesFactor(object=trioSetList, ranges=md.segs2)
show(map.segs)


###################################################
### code chunk number 16: computeBayesFactorOneRange
###################################################
computeBayesFactor(trioSetList, md.segs2[1, ])


###################################################
### code chunk number 17: triofig
###################################################
rd.denovoDel <- map.segs[state(map.segs) == 221 & numberProbes(map.segs) > 5, ]
rd <- rd.denovoDel[1, ]
CHR <- match(as.character(chromosome(rd)), paste("chr",chromosome(trioSetList), sep=""))
S <- match(sampleNames(rd), sampleNames(trioSetList))
trioSet <- trioSetList[[CHR]][, S]
mindist(trioSet) <- md[[CHR]][, S, drop=FALSE]
frame <- 100e3
fig <- xyplotTrio(rd=rd,
		  object=trioSet,
		  frame=frame,
		  ylab="log R ratio",
		  xlab="physical position (Mb)",
		  panel=xypanelTrio,
		  lrr.segments=lrr.segs,
		  md.segments=md.segs,
		  col.hom="grey50",
		  col.het="grey50",
		  col.np="grey20",
		  cex=0.3,
		  layout=c(1, 4),
		  ylim=c(-3, 1.5),
		  scales=list(cex=0.7, x=list(relation="same"),
		  y=list(alternating=1,
		  at=c(-1, 0, log2(3/2), log2(4/2)),
		  labels=expression(-1, 0, log[2](3/2), log[2](4/2)))),
		  par.strip.text=list(cex=0.7),
		  state.col="black",
		  segment.col="black",
		  state.cex=0.8,
		  key=list(text=list(c(expression(log[2]("R ratios")), expression("B allele freqencies")),
			   col=c("grey50", "blue")), columns=2))


###################################################
### code chunk number 18: displayFigure
###################################################
pdf("MinimumDistance-displayFigure.pdf")
print(fig)
dev.off()


###################################################
### code chunk number 19: sessionInfo
###################################################
toLatex(sessionInfo())