#Use help(package=CGHcall)
library(CGHcall)
library(CGHbase)

#CONSIDER CLEANING R'S MEMORY FIRST USING:
#rm(list=ls());gc();

#Enter the folder where the input file is

setwd("C:\\VUData\\CGHdata\\NoWaves")

#Name of the input file
inputfile <- "Corrected_Dutch_tumor.txt" 

NoBpEnd <- TRUE #TRUE if Basepair end info is absent (3 data info columns only)
#NoBpEnd <- FALSE

#CGHcall functions. Use ?function name (eg ?segmentData) to see the options etc.
dfraw <- read.table(inputfile,header=TRUE,sep="\t",fill=TRUE)
if(NoBpEnd) dfraw <- data.frame(dfraw[,1:3],bpend=dfraw[,3]+60,dfraw[,-(1:3)])
dfraw[,1] <- sapply(1:nrow(dfraw),function(x){paste("clone",x,sep="")}) #if clone names contain missings/errors then try this
raw <- cghRaw(dfraw)
rm(dfraw);gc();

prep <- preprocess(raw, maxmiss = 30, nchrom = 23)
rm(raw);gc();

nor <-  normalize(prep,method = "median", cellularity = 1, smoothOutliers = TRUE)  
rm(prep);gc();

#For noisy profiles undo.SD may be increased to enforce less segments. 
#Advise is to consider the segment-plots before calling  
seg <-  segmentData(nor, method = "DNAcopy", nperm=1000,undo.splits="sdundo",min.width=5, undo.SD=1.5)
rm(nor);gc();
save(seg, file="seg.Rdata")
#load("seg.Rdata")
#plot the first segmented profile    
plot(seg[,1])

#load("seg.Rdata")
segnorm <- postsegnormalize(seg,inter=c(-0.1,0.1))
rm(seg);gc();

listcalls <- CGHcall(segnorm,nclass=4,robustsig="yes")
save(listcalls, file="listcalls.Rdata")
#load("listcalls.Rdata")

calls <- ExpandCGHcall(listcalls,segnorm, divide=5, memeff=FALSE)

#IN CASE OF R MEMORY PROBLEMS USE (see also ?ExpandCGHcall); possibly increase divide.  Results can be combined using 'combine'.
#calls <- ExpandCGHcall(listcalls,segnorm, divide=5, memeff=TRUE)

rm(segnorm);gc();
save(calls, file="calls.Rdata")   
plot(calls)

#saves a PostScript file with plots in your working directory (similar for segmentation plots; replace calls by seg)
postscript("CallPlot.ps")
plot(calls)
dev.off()

#pdf plots are more difficult to browse than ps
#pdf("CallPlot.pdf")
#plot(calls)
#dev.off()


############## Export segmented, hard calls (-1,0,1,2) or soft calls (probabilities) as .txt  ############

# Instructions: Indicate which output files you want. Then run everything below (select and copy into R console). 
# You may change the output file names below. Files are store in your working directory (where the input file is as well)

HardCalls <- TRUE
SoftCalls <- TRUE
Segmented <- FALSE
Normalized <- FALSE

Clone <- featureNames(calls)
Chromo <- chromosomes(calls)
BPstart <- bpstart(calls)
BPend <- bpend(calls)
Calls <- calls(calls)
Probloss <- probloss(calls)
Probnorm <- probnorm(calls)
Probgain <- probgain(calls)
colnam <- c(colnames(Probloss),colnames(Probnorm),colnames(Probgain))

#Here, the call probabilities are ordered by sample rather than by aberration state
if(is.null(probamp(calls))) 
{
    allprob <- c()
    ncl <- ncol(Probnorm)
    for(i in 1:ncl)
    {
        Probsall <- cbind(Probloss[,i],Probnorm[,i],Probgain[,i])
        colnames(Probsall) <- c(colnam[i],colnam[ncl+i],colnam[2*ncl+i])
        allprob <- cbind(allprob,Probsall)
    }
} else {
    Probamp <- probamp(calls)
    colnam <- c(colnam,colnames(Probamp))
    allprob <- c()
    ncl <- ncol(Probnorm)
    for(i in 1:ncl)
        {   
            Probsall <- cbind(Probloss[,i],Probnorm[,i],Probgain[,i],Probamp[,i])
            colnames(Probsall) <- c(colnam[i],colnam[ncl+i],colnam[2*ncl+i],colnam[3*ncl+i])
            allprob <- cbind(allprob,Probsall)
        }
}
allprob2 <- round(allprob,3)
probs <- data.frame(Clone,Chromo,BPstart,BPend,allprob2)
exportdata <- data.frame(Clone,Chromo,BPstart,BPend,Calls)
exportdata_seg <- data.frame(Clone,Chromo,BPstart,BPend,segmented(calls))
exportdata_norm <- data.frame(Clone,Chromo,BPstart,BPend,copynumber(calls))

if(HardCalls) {write.table(exportdata, file="CallsTxtFile.txt", sep="\t", quote=F,row.names=F)}
if(SoftCalls) {
    write.table(probs, file="ProbsTxtFile.txt", sep="\t", quote=F,row.names=F)
    #save(probs,file="ProbsRFile.Rdata")
    }
if(Segmented) {write.table(exportdata_seg, file="SegTxtFile.txt", sep="\t", quote=F,row.names=F)}
if(Normalized){write.table(exportdata_norm, file="NormTxtFile.txt", sep="\t", quote=F,row.names=F)}