### 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(SNPchip)
library(GenomicRanges)
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: 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 8: loadTrioSetListExample
###################################################
data(trioSetListExample)


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


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


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


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


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


###################################################
### code chunk number 14: computeBayesFactor
###################################################
gr <- MAP(object=trioSetList, ranges=md.segs2, nupdates=5)
show(gr)


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


###################################################
### code chunk number 16: triofig
###################################################
denovo.range <- GRanges("chr22", IRanges(20.8e6, 21.4e6))
i <- subjectHits(findOverlaps(denovo.range, gr))
gr.denovo <- gr[i, ]
gr.denovo <- gr.denovo[state(gr.denovo)=="221", ]
library(lattice)
library(foreach)
trioSet <- trioSetList[[2]][, match(sampleNames(gr.denovo), sampleNames(trioSetList))]
mindist(trioSet) <- md[[2]][, match(sampleNames(gr.denovo), sampleNames(trioSetList)), drop=FALSE]
fig <- xyplotTrio(gr.denovo, trioSet, frame=500e3,
	   ylab="log R ratio and BAFs",
	   xlab="physical position (Mb)",
	   panel=xypanelTrio,
	   scales=list(x="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)))),
	   lrr.segments=lrr.segs,
	   md.segments=md.segs,
	   layout=c(1, 4),
	   col.hom="grey50",
	   col.het="grey50",
	   col.np="grey20",
	   state.cex=0.8,
	   cex=0.3,
	   ylim=c(-3, 1.5),
	   par.strip.text=list(lines=0.8, cex=0.6),
	   key=list(text=list(c(expression(log[2]("R ratios")), expression("B allele freqencies")),
		    col=c("black", "blue")), columns=2))


###################################################
### code chunk number 17: displayFigure
###################################################
pdf("triofig.pdf", width=10, height=6)
print(fig)
dev.off()


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