###################################################
### chunk number 1: load.data
###################################################
#line 59 "vignettes/CNVtools/inst/doc/CNVtools-vignette.Rnw"
#source("../CNVtools.r"); dyn.load("../../src/CNVtools.so"); load("../../CNVtools/data/A112.RData")
library(CNVtools) 
data(A112)
head(A112)
raw.signal <- as.matrix(A112[, -c(1,2)])
dimnames(raw.signal)[[1]] <- A112$subject
  

mean.signal <- apply(raw.signal, MAR=1, FUN=mean)
pca.signal <- apply.pca(raw.signal)


pdf("fig/mean_pca_signal.pdf", width=10, height=5)
par(mfrow=c(1,2))
hist(mean.signal, breaks=50, main='Mean signal', cex.lab=1.3)
hist(pca.signal, breaks=50, main='First PCA signal', cex.lab=1.3)  
dev.off()


###################################################
### chunk number 2: model.select
###################################################
#line 92 "vignettes/CNVtools/inst/doc/CNVtools-vignette.Rnw"
batches <- factor(A112$cohort)
sample <- factor(A112$subject)
set.seed(0)
results <- CNVtest.select.model(signal=pca.signal, batch = batches, sample = sample, n.H0 = 3, method="BIC", v.ncomp = 1:5, v.model.component = rep('gaussian',5), v.model.mean = rep("~ strata(cn)",5), v.model.var = rep("~1", 5))
ncomp <- results$selected
pdf("fig/modelselect.pdf",width=5,height=5)
plot(-results$BIC, xlab="n comp", ylab="-BIC", type="b", lty=2, col="red", pch = '+')
dev.off()


###################################################
### chunk number 3: cluster.pca
###################################################
#line 130 "vignettes/CNVtools/inst/doc/CNVtools-vignette.Rnw"
ncomp <- 3
batches <- factor(A112$cohort)
sample <- factor(A112$subject)
fit.pca <- CNVtest.binary ( signal = pca.signal, sample = sample, batch = batches, ncomp = ncomp, n.H0=3, n.H1=0, model.var= '~ strata(cn)')
print(fit.pca$status.H0)

pdf("fig/pca-fit.pdf", width=10, height=5)
par(mfrow=c(1,2))
cnv.plot(fit.pca$posterior.H0, batch = '58C', main = 'Cohort 58C', breaks = 50, col = 'red')
cnv.plot(fit.pca$posterior.H0, batch = 'NBS', main = 'Cohort NBS', breaks = 50, col = 'red')
dev.off()  


###################################################
### chunk number 4: genotype.assignment
###################################################
#line 157 "vignettes/CNVtools/inst/doc/CNVtools-vignette.Rnw"
head(fit.pca$posterior.H0)


###################################################
### chunk number 5: ldf.improve
###################################################
#line 165 "vignettes/CNVtools/inst/doc/CNVtools-vignette.Rnw"
ncomp <- 3
pca.posterior <- as.matrix((fit.pca$posterior.H0)[, paste('P',seq(1:ncomp),sep='')])
dimnames(pca.posterior)[[1]] <- (fit.pca$posterior.H0)$subject
ldf.signal <- apply.ldf(raw.signal, pca.posterior)

pdf("fig/ldf_pca_signal.pdf", width=10, height=5)
par(mfrow=c(1,2))
hist(pca.signal, breaks=50, main='First PCA signal', cex.lab=1.3)  
hist(ldf.signal, breaks=50, main='LDF signal', cex.lab=1.3)  
dev.off()


###################################################
### chunk number 6: test.association
###################################################
#line 219 "vignettes/CNVtools/inst/doc/CNVtools-vignette.Rnw"
ncomp <- 3
trait <- ifelse( A112$cohort == '58C', 0, 1)
fit.ldf <- CNVtest.binary ( signal = ldf.signal, sample = sample, batch = batches, disease.status = trait, ncomp = ncomp, n.H0=3, n.H1=1, model.var = "~cn")
print(fit.ldf$status.H0)
print(fit.ldf$status.H1)

pdf("fig/ldf-fit.pdf", width=10, height=5)
par(mfrow=c(1,2))
cnv.plot(fit.ldf$posterior.H0, batch = '58C', main = 'Cohort 58C', breaks = 50, col = 'red')
cnv.plot(fit.ldf$posterior.H0, batch = 'NBS', main = 'Cohort NBS', breaks = 50, col = 'red')
dev.off()

LR.statistic <- -2*(fit.ldf$model.H0$lnL - fit.ldf$model.H1$lnL)
print(LR.statistic)


###################################################
### chunk number 7: test.association.allelic
###################################################
#line 250 "vignettes/CNVtools/inst/doc/CNVtools-vignette.Rnw"
fit.ldf <- CNVtest.binary ( signal = ldf.signal, sample = sample, batch = batches, disease.status = trait, ncomp = 3, n.H0=3, n.H1=1, model.disease = " ~ as.factor(cn)")
print(fit.ldf$status.H0)
print(fit.ldf$status.H1)
LR.statistic <- -2*(fit.ldf$model.H0$lnL - fit.ldf$model.H1$lnL)
print(LR.statistic)


###################################################
### chunk number 8: test.association.qt
###################################################
#line 265 "vignettes/CNVtools/inst/doc/CNVtools-vignette.Rnw"
batches <- rep("ALL",length(sample))
qt <- rnorm(length(sample), mean=9.0, sd=1.0)
fit.ldf <- CNVtest.qt(signal = ldf.signal, sample = sample, batch = batches, qt = qt, ncomp = ncomp, n.H0=3, n.H1=1, model.var = "~strata(cn)")
print(fit.ldf$status.H0)
print(fit.ldf$status.H1)
LR.statistic <- -2*(fit.ldf$model.H0$lnL - fit.ldf$model.H1$lnL)
print(LR.statistic)
pdf("fig/qt-fit.pdf", width=15, height=5)
qt.plot(fit.ldf) 
dev.off()