#
# Load package
#

library(CNVassoc)


#
# Load data
#

data(dataMLPA)
head(dataMLPA)


#
# Density plot
#

par(mfrow=c(2,2),mar=c(3,4,3,1))
hist(dataMLPA$Gene1,main="gene 1 signal histogram",xlab="",ylab="frequency")
hist(dataMLPA$Gene2,main="gene 2 signal histogram",xlab="",ylab="frequency")
par(xaxs="i")
plot(density(dataMLPA$Gene1),main="gene 1 signal density function",xlab="",ylab="density")
plot(density(dataMLPA$Gene2),main="gene 2 signal density function",xlab="",ylab="density")



#
# Inferring CNV status (calling)
#

CNV.1<-cnv(x=dataMLPA$Gene1, threshold.0=0.06, num.class=3, mix.method = "mixdist")
CNV.1

CNV.2<-cnv(x=dataMLPA$Gene2, threshold.0=0.01, mix.method = "mixdist")
CNV.2

probs.2<-attr(CNV.2,"probabilities")
CNV.2probs<-cnv(probs.2)
CNV.2probs


plot(CNV.1, case.control=dataMLPA$casco, main="gen 1")
dev.off()
pdf("./figures/fig2b.pdf")
plot(CNV.2, case.control=dataMLPA$casco, main="gen 2")

plot(CNV.2probs, case.control=dataMLPA$casco)


#
# Measuring uncertainty
#

getQualityScore(CNV.1,type="class")
getQualityScore(CNV.2,type="class")

getQualityScore(CNV.1,type="CNVtools")
getQualityScore(CNV.2,type="CNVtools")


#
# Assesing association
# 

# 1. case-control

model1mul<-CNVassoc(casco~CNV.1,data=dataMLPA,model="mul")
model1mul
summary(model1mul)
CNVtest(model1mul,type="Wald")
CNVtest(model1mul,type="LRT")
coef(summary(model1mul,ref=2))

model2mul<-CNVassoc(casco~CNV.2,data=dataMLPA,model="mul")
model2mul
summary(model1mul)
CNVtest(model2mul,type="Wald")
CNVtest(model2mul,type="LRT")
coef(summary(model2mul,ref=2))

model2add<-CNVassoc(casco~CNV.2, data=dataMLPA, model="add")
summary(model2add)
anova(model2mul,model2add)



# 2. quantitative trait and adjust for covariate

mod<-CNVassoc(quanti ~ CNV.2 + cov, data = dataMLPA, model = 'add', emsteps=10)
mod