### R code from vignette source 'Clonality.Rnw'

###################################################
### code chunk number 1: read data
###################################################
 library(Clonality)
set.seed(100)
chrom<-paste("chr",rep(c(1:22),each=100),"p",sep="")
chrom[nchar(chrom)==5]<-paste("chr0",substr(chrom[nchar(chrom)==5] ,4,5),sep="")
maploc<- rep(c(1:100),22)
data<-NULL
for (pt in 1:9)  #first 9 patients have independent tumors
{
tumor1<-tumor2<- NULL
mean1<- rnorm(22) 
mean2<- rnorm(22)
for (chr in 1:22)
{ 
  r<-runif(2) 
if (r[1]<=0.5) tumor1<-c(tumor1,rep(0,100))   
  else if   (r[1]>0.7)  tumor1<-c(tumor1,rep(mean1[chr],100))
  else  { i<-sort(sample(1:100,2))
        tumor1<-c(tumor1,mean1[chr]*c(rep(0,  i[1]),rep(1, i[2]-i[1]), rep(0,  100-i[2])))
        }
if (r[2]<=0.5) tumor2<-c(tumor2,rep(0,100))
  else if   (r[2]>0.7)  tumor2<-c(tumor2,rep(mean2[chr],100))
  else   {i<-sort(sample(1:100,2))
       tumor2<-c(tumor2,mean2[chr]*c(rep(0,  i[1]),rep(1, i[2]-i[1]), rep(0,  100-i[2])))
         }
}
data<-cbind(data,tumor1,tumor2)
}

#last patient has identical profiles
tumor1<- NULL
mean1<- rnorm(22) 
for (chr in 1:22)
{ 
  r<-runif(1) 
if (r<=0.4) tumor1<-c(tumor1,rep(0,100))   
  else if   (r>0.6)  tumor1<-c(tumor1,rep(mean1[chr],100))
  else  { i<-sort(sample(1:100,2))
        tumor1<-c(tumor1,mean1[chr]*c(rep(0,  i[1]),rep(1, i[2]-i[1]), rep(0,  100-i[2])))
        }

}
data<-cbind(data,tumor1,tumor1)

data<-data+matrix(rnorm( 44000,mean=0,sd=0.4) ,nrow=2200,ncol=20)


samnms<-paste("pt",rep(1:10,each=2),rep(1:2,10),sep=".")



###################################################
### code chunk number 2: Clonality.Rnw:104-105
###################################################
dim(data)


###################################################
### code chunk number 3: CNA
###################################################

dataCNA<-CNA(data,chrom=chrom,maploc=maploc,sampleid=samnms)
as.matrix(dataCNA)[1:5,1:10]



###################################################
### code chunk number 4: averaging
###################################################
dataAve<- ave.adj.probes(dataCNA,2)


###################################################
### code chunk number 5: Clonality.Rnw:130-131
###################################################
ptlist<- paste("pt",rep(1:10,each=2),sep=".")


###################################################
### code chunk number 6: main clonality analysis
###################################################
results<-clonality.analysis(dataCNA, ptlist,  pfreq = NULL, refdata = NULL, nmad = 1,  reference = TRUE, allpairs = TRUE)


###################################################
### code chunk number 7: Clonality.Rnw:140-141
###################################################
results$LR


###################################################
### code chunk number 8: genomewidePlots
###################################################
genomewidePlots(results$OneStepSeg, results$ChromClass,ptlist , c("pt.10.1", "pt.10.2"),results$LR,  plot.as.in.analysis = TRUE) 


###################################################
### code chunk number 9: Clonality.Rnw:154-155
###################################################
chromosomePlots(results$OneStepSeg, ptlist,ptname="pt.10",nmad=1)


###################################################
### code chunk number 10: Clonality.Rnw:160-161
###################################################
histogramPlot(results$LR[,4], results$refLR[,4])


###################################################
### code chunk number 11: LOH data generation
###################################################
set.seed(25)
LOHtable<-cbind(1:20,matrix(sample(c(0,1,2),20*20,replace=TRUE),20))
LOHtable[,3]<-LOHtable[,2]
LOHtable[1,3]<-0


###################################################
### code chunk number 12: LOH data analysis
###################################################
LOHtable[,1:5]
LOHclonality(LOHtable,rep(1:10,each=2),pfreq=NULL,noloh=0,loh1=1,loh2=2)


###################################################
### code chunk number 13: LR estimation
###################################################
data(lcis)
n<-nrow(lcis)
summary(lcis$probi)
table(lcis$TK47IDC.TK47LCIS1 )
lcis$probi[lcis$TK47IDC.TK47LCIS1==1]


###################################################
### code chunk number 14: LR test
###################################################
x1<-x2<-rep(0,n)
x1[lcis$TK47IDC.TK47LCIS1==1]<-x2[lcis$TK47IDC.TK47LCIS1==1]<-1
x1[lcis$TK47IDC.TK47LCIS1==2]<-1
set.seed(1) 


###################################################
### code chunk number 15: LR test2
###################################################
SNVtest(x1,x2,lcis$probi)


###################################################
### code chunk number 16: REM estimation
###################################################
data(lcis)
mod <- mutation.rem(lcis)


###################################################
### code chunk number 17: REM estimation2
###################################################
print(mod)


###################################################
### code chunk number 18: REM proba
###################################################
data(lcis)
mod <- mutation.rem(lcis, proba=TRUE)



###################################################
### code chunk number 19: REM proba2
###################################################
print(mod)


###################################################
### code chunk number 20: REM proba new case
###################################################
# generate a case with 30 mutations
# probabilities of each observed mutation
set.seed(159)
pi <- runif(30,0.001,0.13)
# mutation 1=shared or 2=private
newpair <- cbind(pi,rbinom(30,1,1-pi^2)+1)
# generate the matrix of likelihood values
new.likmat <- grid.lik(xigrid=c(0, seq(0.0005, 0.9995, by=0.001)), 
                       as.matrix(newpair[,c(-1)]), newpair[,1])
# probability of being clonal using the model previoulsy estimated
proba <- mutation.proba(c(mod$mu, mod$sigma, mod$pi), t(as.matrix(new.likmat)) )


###################################################
### code chunk number 21: REM proba new case2
###################################################
print(proba)


###################################################
### code chunk number 22: sessionInfo
###################################################
sessionInfo()