######################################################
### chunk number 0: Get data
######################################################
library(CNVassoc)
data(NeveData)
intensities<-NeveData$data
pheno<-NeveData$pheno

######################################################
# Calling takes about 25 min. Use 'data(NeveCalled)' 
# to get results from chunks 1 to 5
######################################################
data(NeveCalled) 


######################################################
### chunk number 1: Class of aCGH data 
######################################################
library(CGHcall)
Neve <- cghRaw(intensities)


######################################################
### chunk number 2: Preprocessing
######################################################
cghdata <- preprocess(Neve, maxmiss=30, nchrom=22)


######################################################
### chunk number 3: Normalization
######################################################
norm.cghdata <- normalize(cghdata, method="median", smoothOutliers=TRUE)


######################################################
### chunk number 4: Segmentation
######################################################
seg.cghdata <- segmentData(norm.cghdata, method="DNAcopy")


######################################################
### chunk number 5: Calling
######################################################
NeveCalled <- CGHcall(seg.cghdata)


######################################################
### chunk number 6: Get posterior probabilities 
######################################################
probs<-getProbs(NeveCalled)



######################################################
# Blocking takes about 3 min. Use 'data(NeveRegions)' 
# to get results from chunk 7
######################################################
data(NeveRegions)


######################################################
### chunk number 7: Blocking/Regions
######################################################
library(CGHregions)
NeveRegions <- CGHregions(NeveCalled)

# probabilities 
probsRegions<-getProbsRegions(probs, NeveRegions, intensities) 


######################################################
### chunk number 8: Association Analysis
######################################################
pvals<-multiCNVassoc(probsRegions, formula="pheno~CNV", model="mult", num.copies=0:2, cnv.tol=0.01)
pvalsBH<-getPvalBH(pvals)
# Table 6 from Gonzalez et al, (2009). BMC Bioinformatics
cumsum(table(cut(pvalsBH[,1],c(-Inf,1e-5,1e-4,1e-3,1e-2,0.05))))