###################################################
### chunk number 1: 
###################################################
require(MANOR)


###################################################
### chunk number 2: 
###################################################
SNR.FUN <- function(arrayCGH, var.FG, var.BG, snr.thr) {
  which(arrayCGH$arrayValues[[var.FG]] < arrayCGH$arrayValues[[var.BG]]*snr.thr)
}
SNR.char <- "B"
SNR.label <- "Low signal to noise ratio"
SNR.flag <- to.flag(SNR.FUN, SNR.char, args=alist(var.FG="REF_F_MEAN", var.BG="REF_B_MEAN", snr.thr=3))


###################################################
### chunk number 3: 
###################################################
global.spatial.FUN <- function(arrayCGH, var)
  {
    if (!is.null(arrayCGH$arrayValues$Flag))
        arrayCGH$arrayValues$LogRatio[which(arrayCGH$arrayValues$Flag!="")] <- NA
##     Trend <- arrayTrend(arrayCGH, var, span=0.03, degree=1, iterations=3, family="symmetric")
    Trend <- arrayTrend(arrayCGH, var, span=0.03, degree=1, iterations=3)
    arrayCGH$arrayValues[[var]] <- Trend$arrayValues[[var]]-Trend$arrayValues$Trend
    arrayCGH
  }
global.spatial.flag <- to.flag(global.spatial.FUN, args=alist(var="LogRatio"))


###################################################
### chunk number 4: 
###################################################
chromosome.FUN <- function(arrayCGH, var) {
  var.rep <- arrayCGH$id.rep
  w <- which(!is.na(match(as.character(arrayCGH$cloneValues[[var]]), c("X", "Y"))))
  l <- arrayCGH$cloneValues[w, var.rep]
  which(!is.na(match(arrayCGH$arrayValues[[var.rep]], as.character(l))))
}

chromosome.char <- "X"
chromosome.label <- "Sexual chromosome"
chromosome.flag <- to.flag(chromosome.FUN, chromosome.char, type="temp.flag", args=alist(var="Chromosome"), label=chromosome.label)


###################################################
### chunk number 5: 
###################################################
args(to.flag)


###################################################
### chunk number 6: 
###################################################
args(flag.arrayCGH)


###################################################
### chunk number 7: 
###################################################
args(flag.summary.arrayCGH)


###################################################
### chunk number 8: 
###################################################
args(flag.summary.default)


###################################################
### chunk number 9: 
###################################################
pct.spot.FUN <- function(arrayCGH, var) {
  100*sum(!is.na(arrayCGH$arrayValues[[var]]))/dim(arrayCGH$arrayValues)[1]
}
pct.spot.name <- "SPOT_PCT"
pct.spot.label <- "Proportion of spots after normalization"
pct.spot.qscore <- to.qscore(pct.spot.FUN, name=pct.spot.name, args=alist(var="LogRatioNorm"), label=pct.spot.label)


###################################################
### chunk number 10: 
###################################################
args(to.qscore)


###################################################
### chunk number 11: 
###################################################
args(qscore.arrayCGH)


###################################################
### chunk number 12: 
###################################################
args(qscore.summary.arrayCGH)


###################################################
### chunk number 13: 
###################################################
data(spatial)


###################################################
### chunk number 14: 
###################################################
data(spatial)

## edge: example of array with local spatial effects
arrayPlot(edge, "LogRatio", main="Local spatial effects", zlim=c(-1,1), mediancenter=TRUE, bar="h")


###################################################
### chunk number 15: 
###################################################
data(spatial)
arrayPlot(gradient, "LogRatio", main="Spatial gradient" , zlim=c(-2,2), mediancenter=TRUE, bar="h")


###################################################
### chunk number 16: 
###################################################
dev.off()


###################################################
### chunk number 17: 
###################################################
data(spatial)
genome.plot(edge.norm, chrLim="LimitChr", cex=1)


###################################################
### chunk number 18: 
###################################################
data(spatial)
edge.norm$cloneValues$ZoneGNL <- as.factor(edge.norm$cloneValues$ZoneGNL)
genome.plot(edge.norm, col.var="ZoneGNL", chrLim="LimitChr", cex=1)


###################################################
### chunk number 19: 
###################################################
data(spatial)
report.plot(edge.norm, chrLim="LimitChr", zlim=c(-1,1), cex=1)


###################################################
### chunk number 20: 
###################################################
dev.off()


###################################################
### chunk number 21: 
###################################################
dir.in <- system.file("data", package="MANOR")

## import from 'spot' files
spot.names <- c("LogRatio", "RefFore", "RefBack", "DapiFore", "DapiBack", "SpotFlag", "ScaledLogRatio")
clone.names <- c("PosOrder", "Chromosome")
edge <- import(paste(dir.in, "/edge.txt", sep=""), type="spot",
spot.names=spot.names, clone.names=clone.names, add.lines=TRUE)


###################################################
### chunk number 22: 
###################################################
data(flags)
data(spatial)

## local.spatial.flag$args <- alist(var="ScaledLogRatio", by.var=NULL, nk=5, prop=0.25, thr=0.15, beta=1, family="symmetric")
local.spatial.flag$args <- alist(var="ScaledLogRatio", by.var=NULL, nk=5, prop=0.25, thr=0.15, beta=1, family="gaussian")
flag.list <- list(spatial=local.spatial.flag, spot=spot.corr.flag, ref.snr=ref.snr.flag, dapi.snr=dapi.snr.flag, rep=rep.flag, unique=unique.flag)

edge.norm <- norm(edge, flag.list=flag.list, FUN=median, na.rm=TRUE)
edge.norm <- sort(edge.norm, position.var="PosOrder")


###################################################
### chunk number 23: 
###################################################
report.plot(edge.norm, chrLim="LimitChr", zlim=c(-1,1), cex=1)


###################################################
### chunk number 24: 
###################################################
dev.off()


###################################################
### chunk number 25: 
###################################################
##DNA copy number assessment: GLAD
profileCGH <- as.profileCGH(edge.norm$cloneValues)

profileCGH <- daglad(profileCGH, smoothfunc="lawsglad", lkern="Exponential", model="Gaussian", qlambda=0.999,  bandwidth=10, base=FALSE, round=2, lambdabreak=6, lambdaclusterGen=20, param=c(d=6), alpha=0.001, msize=5, method="centroid", nmin=1, nmax=8, amplicon=1, deletion=-5, deltaN=0.10,  forceGL=c(-0.15,0.15), nbsigma=3, MinBkpWeight=0.35, verbose=FALSE)

edge.norm$cloneValues <- as.data.frame(profileCGH)
edge.norm$cloneValues$ZoneGNL <- as.factor(edge.norm$cloneValues$ZoneGNL)

data(qscores)
## list of relevant quality scores
qscore.list <- list(smoothness=smoothness.qscore,
                    var.replicate=var.replicate.qscore,
                    dynamics=dynamics.qscore)
edge.norm$quality <- qscore.summary.arrayCGH(edge.norm, qscore.list)
edge.norm$quality


###################################################
### chunk number 26: 
###################################################
html.report(edge.norm, dir.out=".", array.name="an array with local bias", chrLim="LimitChr", light=FALSE, pch=20, zlim=c(-2,2), file.name="edge")


###################################################
### chunk number 27: 
###################################################
## import from 'gpr' files
spot.names <- c("Clone", "FLAG", "TEST_B_MEAN", "REF_B_MEAN", "TEST_F_MEAN", "REF_F_MEAN", "ChromosomeArm")
clone.names <- c("Clone", "Chromosome", "Position", "Validation")

ac <- import(paste(dir.in, "/gradient.gpr", sep=""), type="gpr", spot.names=spot.names, clone.names=clone.names, sep="\t", comment.char="@", add.lines=TRUE)

## compute log-ratio
ac$arrayValues$F1 <- log(ac$arrayValues[["TEST_F_MEAN"]], 2)
ac$arrayValues$F2 <- log(ac$arrayValues[["REF_F_MEAN"]], 2)
ac$arrayValues$B1 <- log(ac$arrayValues[["TEST_B_MEAN"]], 2)
ac$arrayValues$B2 <- log(ac$arrayValues[["REF_B_MEAN"]], 2)

Ratio <- (ac$arrayValues[["TEST_F_MEAN"]]-ac$arrayValues[["TEST_B_MEAN"]])/
    (ac$arrayValues[["REF_F_MEAN"]]-ac$arrayValues[["REF_B_MEAN"]])
Ratio[(Ratio<=0)|(abs(Ratio)==Inf)] <- NA
ac$arrayValues$LogRatio <- log(Ratio, 2)
gradient <- ac


###################################################
### chunk number 28: 
###################################################
data(spatial)
data(flags)

flag.list <- list(local.spatial=local.spatial.flag, spot=spot.flag, SNR=SNR.flag, global.spatial=global.spatial.flag, val.mark=val.mark.flag, position=position.flag, unique=unique.flag, amplicon=amplicon.flag, chromosome=chromosome.flag, replicate=replicate.flag)

gradient.norm <- norm(gradient, flag.list=flag.list, FUN=median, na.rm=TRUE)
gradient.norm <- sort(gradient.norm)


###################################################
### chunk number 29: 
###################################################
genome.plot(gradient.norm, chrLim="LimitChr", cex=1)


###################################################
### chunk number 30: 
###################################################
##DNA copy number assessment: GLAD
profileCGH <- as.profileCGH(gradient.norm$cloneValues)

profileCGH <- daglad(profileCGH, smoothfunc="lawsglad", lkern="Exponential", model="Gaussian", qlambda=0.999,  bandwidth=10, base=FALSE, round=2, lambdabreak=6, lambdaclusterGen=20, param=c(d=6), alpha=0.001, msize=5, method="centroid", nmin=1, nmax=8, amplicon=1, deletion=-5, deltaN=0.10,  forceGL=c(-0.15,0.15), nbsigma=3, MinBkpWeight=0.35, verbose=FALSE)

gradient.norm$cloneValues <- as.data.frame(profileCGH)
gradient.norm$cloneValues$ZoneGNL <- as.factor(gradient.norm$cloneValues$ZoneGNL)

data(qscores)
## list of relevant quality scores
qscore.list <- list(smoothness=smoothness.qscore, var.replicate=var.replicate.qscore, dynamics=dynamics.qscore)
gradient.norm$quality <- qscore.summary.arrayCGH(gradient.norm, qscore.list)
gradient.norm$quality


###################################################
### chunk number 31: 
###################################################
html.report(gradient.norm, dir.out=".", array.name="an array with spatial gradient", chrLim="LimitChr", light=FALSE, pch=20, zlim=c(-2,2), file.name="gradient")


###################################################
### chunk number 32: 
###################################################
sessionInfo()