### R code from vignette source 'vignettes/genoset/inst/doc/genoset.Rnw'

###################################################
### code chunk number 1: <options
###################################################
  options(width=90)


###################################################
### code chunk number 2: objectcreation
###################################################
library(genoset)
data(genoset)
genoset.ds = GenoSet( locData=locData.rd, foo=fake.lrr, pData=fake.pData, annotation="SNP6", universe="hg19" )

baf.ds = BAFSet( locData=locData.rd, lrr=fake.lrr, baf=fake.baf, pData=fake.pData, annotation="SNP6", universe="hg19" )

cn.ds = CNSet( locData=locData.rd, cn=fake.lrr, pData=fake.pData, annotation="SNP6", universe="hg19" )
cn.ds

rle.baf.ds = BAFSet( locData=locData.rd,
  lrr=DataFrame(K=Rle(c(rep(1.5,300),rep(2.3,700))),L=Rle( c(rep(3.2,700),rep(2.1,300)) ), M=Rle(rep(1.1,1000)),row.names=rownames(fake.lrr)),
  baf=DataFrame(K=Rle(c(rep(0.05,600),rep(0.5,400))),L=Rle( c(rep(0,700),rep(0.5,300)) ), M=Rle(rep(1,1000)),row.names=rownames(fake.baf)),
  pData=fake.pData,
  annotation="SNP6"
  )


###################################################
### code chunk number 3: <objectinspection>>
###################################################
slotNames(genoset.ds)
phenoData(genoset.ds)
pData(genoset.ds)
universe(genoset.ds)
annotation(genoset.ds)
locData(genoset.ds)


###################################################
### code chunk number 4: objectassaydata
###################################################
assayDataElementNames(baf.ds)
head( baf(baf.ds) )
head( lrr(baf.ds) )
head( cn(cn.ds) )
head( assayDataElement(genoset.ds,"foo") )


###################################################
### code chunk number 5: accessgenomeinfo
###################################################
chrNames(genoset.ds)
chrOrder(c("chr12","chr12","chrX","chr8","chr7","chrY"))
chrInfo(genoset.ds)
chrIndices(genoset.ds)
elementLengths(genoset.ds)

head(ranges(genoset.ds))
head(chr(genoset.ds))
head(start(genoset.ds))
head(end(genoset.ds))
head(pos(genoset.ds))
tail(pos(genoset.ds))
tail(genoPos(genoset.ds))


###################################################
### code chunk number 6: subsetbychr
###################################################
chr12.ds = cn.ds[ chrIndices(cn.ds,"chr12"), ]
chr12.ds


###################################################
### code chunk number 7: subsetbygene
###################################################
gene.rd = RangedData(ranges=IRanges(start=c(35e6,127e6),end=c(35.5e6,129e6),                       
                       names=c("HER2","CMYC")), space=c("chr17","chr8"))
gene.ds = cn.ds[ gene.rd, ]
gene.ds


###################################################
### code chunk number 8: subsetsamples
###################################################
cn.ds[1:4,1:2]
cn.ds[ c("p1","p2"), cn.ds$b %in% c("E","F")]


###################################################
### code chunk number 9: subsetassaydata
###################################################
all( cn.ds[ 1:4,1:2,"cn"] == cn(cn.ds)[1:4,1:2] )
all( cn.ds[ 1:4,1:2,"cn"] == assayDataElement(cn.ds,"cn")[1:4,1:2] )
cn.ds[ gene.rd, 1:2, "cn" ]


###################################################
### code chunk number 10: genomeorder
###################################################
chrOrder(names(genoset.ds))
isGenomeOrder(genoset.ds, strict=FALSE)
isGenomeOrder(genoset.ds, strict=TRUE)
genoset.ds = toGenomeOrder(genoset.ds, strict=TRUE)
isGenomeOrder(genoset.ds, strict=TRUE)


###################################################
### code chunk number 11: genelevel
###################################################
boundingIndices( c(127e6,127.5e6), c(127e6,128e6), start(chr12.ds) )
rangeSampleMeans(gene.rd, baf.ds, "lrr")


###################################################
### code chunk number 12: cbsdirect
###################################################
library(DNAcopy)
cbs.cna = CNA(cn(cn.ds), chr(cn.ds), pos(cn.ds) )
cbs.smoothed.CNA = smooth.CNA( cbs.cna )
cbs.segs = segment( cbs.cna )


###################################################
### code chunk number 13: runCBS
###################################################
assayDataElement(cn.ds,"cn.segs") = runCBS(cn(cn.ds),locData(cn.ds))


###################################################
### code chunk number 14: cbsmulticore (eval = FALSE)
###################################################
## library(multicore)
## assayDataElement(cn.ds,"cn.segs") = runCBS(cn(cn.ds),locData(cn.ds), n.cores=3)
## assayDataElement(cn.ds,"cn.segs")[1:5,1:3]
## segTable( assayDataElement(cn.ds,"cn.segs"), locData(cn.ds) )


###################################################
### code chunk number 15: GC (eval = FALSE)
###################################################
## library(BSgenome.Hsapiens.UCSC.hg19)
## cn.ds2 = cn.ds
## cn.ds2 = loadGC(cn.ds2,expand=1e6,bsgenome=Hsapiens)
## head(locData(cn.ds2)$gc)
## genoPlot(genoPos(cn.ds),locData(cn.ds2)$gc,locs=locData(cn.ds))
## cn(cn.ds2) = gcCorrect(cn(cn.ds2),locData(cn.ds2)$gc)


###################################################
### code chunk number 16: genoset.Rnw:336-337
###################################################
png(file="genoset-plotgenome.png",width=12,height=5,res=150,units="in")


###################################################
### code chunk number 17: plotgenome
###################################################
genoPlot(cn.ds, 1)
genoPlot(cn.ds, 1, element="cn.segs", add=TRUE, col="red")


###################################################
### code chunk number 18: genoset.Rnw:345-346
###################################################
invisible(dev.off())


###################################################
### code chunk number 19: genoset.Rnw:358-359
###################################################
png(file="genoset-plotchr.png",width=12,height=5,res=150,units="in")


###################################################
### code chunk number 20: plotchr
###################################################
genoPlot(cn.ds,1,chr="chr12")
genoPlot(cn.ds,1,chr="chr12",element="cn.segs",add=TRUE, col="red")


###################################################
### code chunk number 21: genoset.Rnw:367-368
###################################################
invisible(dev.off())


###################################################
### code chunk number 22: plotchrsimple (eval = FALSE)
###################################################
## chr12.ds = cn.ds[chr(cn.ds) == "chr12",]
## genoPlot(pos(chr12.ds),cn(chr12.ds)[,1],locs=locData(chr12.ds))  # Numeric data and location
## genoPlot(pos(chr12.ds),assayDataElement(chr12.ds,"cn.segs")[,1],add=TRUE, col="red") # Rle data and numeric position


###################################################
### code chunk number 23: mbafcutoff
###################################################
mbaf.data = baf2mbaf(baf(baf.ds),hom.cutoff = 0.90)
assayDataElement(baf.ds,"mbaf") = mbaf.data


###################################################
### code chunk number 24: mbaftorle
###################################################
mbaf.data = DataFrame( sapply(colnames( mbaf.data),
  function(x) { Rle( mbaf.data[,x] ) },
  USE.NAMES=TRUE, simplify=FALSE ) )
as.numeric(object.size( assayDataElement(baf.ds,"mbaf")))  / as.numeric( object.size(mbaf.data))


###################################################
### code chunk number 25: segment
###################################################
assayDataElement(baf.ds,"baf.segs") = runCBS( assayDataElement(baf.ds,"mbaf"), locData(baf.ds) )
assayDataElement(baf.ds,"lrr.segs") = runCBS( lrr(baf.ds), locData(baf.ds) )


###################################################
### code chunk number 26: genoset.Rnw:442-443
###################################################
png(file="genoset-plotlrr.png",width=7,height=4,res=150,units="in")


###################################################
### code chunk number 27: plotlrr
###################################################
genoPlot(baf.ds,1,chr="chr12",element="lrr",main="LRR of chr12")
genoPlot(baf.ds,1,chr="chr12",element="lrr.segs",add=TRUE,col="red")


###################################################
### code chunk number 28: genoset.Rnw:449-450
###################################################
invisible(dev.off())


###################################################
### code chunk number 29: genoset.Rnw:463-464
###################################################
png(file="genoset-plotbaf.png",width=7,height=8,res=150,units="in")


###################################################
### code chunk number 30: plotbaf
###################################################
par(mfrow=c(2,1))
genoPlot(baf.ds,1,chr="chr12",element="baf", main="BAF of chr12")
genoPlot(baf.ds,1,chr="chr12",element="mbaf", main="mBAF of chr12")
genoPlot(baf.ds,1,chr="chr12",element="baf.segs", add=TRUE,col="red")


###################################################
### code chunk number 31: genoset.Rnw:472-473
###################################################
invisible(dev.off())


###################################################
### code chunk number 32: crosssample
###################################################
gain.list = lapply(sampleNames(baf.ds),
  function(sample.name) {
    as.logical( assayDataElement(baf.ds,"lrr.segs")[,sample.name] > 0.3 )
})
gain.mat = do.call(cbind, gain.list)
gain.freq = rowMeans(gain.mat,na.rm=TRUE)


###################################################
### code chunk number 33: bigmatrix (eval = FALSE)
###################################################
## cn.ds = convertToBigMatrix(cn.ds)
## save(cn.ds,file="cn.RData")
## reloaded.cn.ds = readGenoSet("cn.RData")