# Head to the "Integration of CnaStruct in other tools" document in
# CnaStruct webpage for integration of CnaStruct with analysis methods. 
# If you wish to use it alone, the following example should aid you in 
# the matter. Anyway, do not hesitate to contact me for further queries.


# running this example gives new generated data each time

library(CnaGen) # http://web.bioinformatics.cicbiogune.es/cnagen/ to download
library(CnaStruct)

# here, we generate some example data:
# 2 normal diploid region of 100 and 200 SNPs
# 2 copy number alterations (copy numbers 1 and 3) of 50 SNPs, twice
# 3 somatic LOH regions (copy numbers 1, 2 and 3) of 20 SNPs
cnagen(
  lengths = list(c(100,200), c(50), c(20), NULL,NULL,NULL),
  times = c(1, 2, 1, 0,0,0),
  copynumbers = list(list(c(1),c(2),c(3)),NULL,NULL),
  purities = 0.5, noiseLRR = 0.25
)
# lrr and baf variables are loaded from the Rdata file generated by CnaGen
load("data/p0.5_n1_metainfo.Rdata")

# bps will contain the output of the segmentation, namely the breakpoints

######### #########
bps = breakpoints(lrr,baf, maxseg=10, maxk=length(lrr), beta=0.5, homthr=0.9)
######### #########
# the best segmentation found will contain as many as maxseg segments
# maxk is the maximum segment length allowed
#   (best performance if not greater than few hundreds)
# beta controls the relative importance of BAF with respect to LRR
# homthr sets the BAF limit to consider a SNP homozygous

# we plot the data and the breakpoints to check for consistence
par(mfrow=c(2,1))
plot(lrr)
abline(v=regions$start, col="blue", lwd=2)
abline(v=bps, col="red")
plot(baf)
abline(v=regions$start, col="blue", lwd=2)
abline(v=bps, col="red")