inst/doc/genoset.R

In genoset: A RangedSummarizedExperiment with methods for copy number analysis

## ----style, echo = FALSE, results='hide'--------------------------------------
BiocStyle::markdown()

## ----objectcreation-----------------------------------------------------------
library(genoset)

sample.names = LETTERS[11:13]
probe.names = paste("p", 1:1000, sep="")
num.samples = length(sample.names)
num.probes  = length(probe.names)

locs = GRanges(
  ranges= IRanges(
      start=c(seq(from=125e6,by=3e4,length=400), seq(from=1,length=400,by=3.25e4),
              seq(from=30e6,length=200,by=3e4)),width=1,names=probe.names),
  seqnames=factor(c(rep("chr8",400), rep("chr12",400),rep("chr17",200)),levels=c("chr8","chr12","chr17")))

fake.cn = matrix(c(
  c(rnorm(200,.4,0.05),rnorm(200,.2,0.05),rnorm(200,0.23,0.05),rnorm(200,.15,0.05),rnorm(200,2.,0.05)), 
  c(rnorm(200,0,0.05), rnorm(200,3,0.05),rnorm(200,14,0.05),rnorm(200,.1,0.05),rnorm(200,-0.05,0.05)),
  c(rnorm(200,.1,0.05),rnorm(200,1,0.05),rnorm(200,-.5,0.05),rnorm(200,3,0.05),rnorm(200,3,0.05))
  ),
  nrow=num.probes,ncol=num.samples,dimnames=list(probe.names,sample.names))
fake.pData=data.frame(matrix(LETTERS[1:15],nrow=3,ncol=5,dimnames=list(sample.names,letters[1:5])))

gs = GenoSet( rowRanges=locs, assays=list(cn=fake.cn), colData=fake.pData )
gs

rle.ds = GenoSet( rowRanges=locs,
                 assays=list(cn = fake.cn,
                             cn.segments=RleDataFrame(
                                 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.cn))),
                             colData=fake.pData
  )


## ----objectassaydata----------------------------------------------------------
names(assays(rle.ds))
head( rle.ds[,,"cn"] )
head( rle.ds[,,"cn.segments"] )

## ----accessgenomeinfo---------------------------------------------------------
head( rowRanges(gs) )
chrNames(gs)
chrOrder(c("chr12","chr12","chrX","chr8","chr7","chrY"))
chrInfo(gs)
chrIndices(gs)

head(chr(gs))
head(start(gs))
head(end(gs))
head(pos(gs))
head(genoPos(gs))

## ----genomeorder--------------------------------------------------------------
chrOrder(chrNames(gs))
gs = toGenomeOrder(gs, strict=TRUE)
isGenomeOrder(gs, strict=TRUE)

## ----subsetbychr--------------------------------------------------------------
chr12.ds = gs[ chrIndices(gs,"chr12"), ]
dim(chr12.ds)
chrIndices(chr12.ds)

## ----subsetbygene-------------------------------------------------------------
gene.gr = GRanges(ranges=IRanges(start=c(35e6,127e6),end=c(35.5e6,129e6),                       
                       names=c("HER2","CMYC")), seqnames=c("chr17","chr8"))
gene.ds = gs[ gene.gr, ]
dim(gene.ds)
chrIndices(gene.ds)

## ----subsetsamples------------------------------------------------------------
dim(gs[1:4,1:2])

## ----subsetassaydata----------------------------------------------------------
all( gs[ 1:4,1:2,"cn"] == assay(gs,"cn")[1:4,1:2] )

## ----GC, eval=FALSE-----------------------------------------------------------
#  library(BSgenome.Hsapiens.UCSC.hg19)
#  gc = rnorm(nrow(gs))
#  gs[,,"cn"] = gcCorrect(gs[,,"cn"],gc)

## ----cbsdirect----------------------------------------------------------------
library(DNAcopy)
cbs.cna = CNA(gs[,,"cn"], chr(gs), pos(gs) )
cbs.smoothed.CNA = smooth.CNA( cbs.cna )
cbs.segs = segment( cbs.cna )

## ----runCBS-------------------------------------------------------------------
gs[,,"cn.segs"] = runCBS(gs[,,"cn"],rowRanges(gs))

## ----cbsmulticore,eval=FALSE--------------------------------------------------
#  library(parallel)
#  gs[,,"cn.segs"] = runCBS(gs[, , "cn"],rowRanges(gs), n.cores=3)
#  gs[,,"cn.segs"][1:5,1:3]

## ----segments-----------------------------------------------------------------
head( gs[,,"cn.segs"] )

segs = segTable( gs[,2,"cn.segs"], rowRanges(gs) )
list.of.segs = segTable( gs[,,"cn.segs"], rowRanges(gs) )
rbind.list.of.segs = segTable( gs[,,"cn.segs"], rowRanges(gs), stack=TRUE )
two.kinds.of.segs = segPairTable( gs[,2,"cn.segs"], gs[,3,"cn.segs"], rowRanges(gs) )

rle = segs2Rle( segs, rowRanges(gs) )
rle.df = segs2RleDataFrame( list.of.segs, rowRanges(gs) )

bounds = matrix( c(1,3,4,6,7,10),ncol=2,byrow=TRUE)
cn = c(1,3,2)
rle = bounds2Rle( bounds, cn, 10 )

## ----plotgenome, echo=TRUE----------------------------------------------------
genoPlot(gs, gs[ , 1, "cn"])
genoPlot(gs, gs[ , 1, "cn.segs"], add=TRUE, col="red")

## ----plotchr, echo=TRUE-------------------------------------------------------
genoPlot(gs,gs[,1,"cn"],chr="chr12")
genoPlot(gs,gs[,1,"cn.segs"],chr="chr12",add=TRUE, col="red")

## ----plotchrsimple, eval=FALSE------------------------------------------------
#  chr12.ds = gs[chr(gs) == "chr12",]
#  genoPlot(pos(chr12.ds),chr12.ds[,1,"cn"],locs=rowRanges(chr12.ds))  # Numeric data and location
#  genoPlot(pos(chr12.ds),chr12.ds[,1,"cn.segs"],add=TRUE, col="red") # Rle data and numeric position

## ----mbafcutoff---------------------------------------------------------------
fake.baf  = sample(c(0,0.5,1), length(probe.names), replace=TRUE) + rnorm(length(probe.names),0,0.01)
fake.baf[ fake.baf > 1 ] = fake.baf[ fake.baf > 1 ] - 1
fake.baf[ fake.baf < 0 ] = fake.baf[ fake.baf < 0 ] + 1
hets = fake.baf < 0.75 & fake.baf > 0.25
fake.baf[ 101:200 ][ hets[101:200] ] = fake.baf[ 101:200 ][ hets[101:200] ] + rep(c(-0.2,0.2),50)[hets[101:200]]
fake.baf = matrix(fake.baf,nrow=num.probes,ncol=num.samples,dimnames=list(probe.names,sample.names))

baf.ds = GenoSet( rowRanges=locs, assays=list(lrr=fake.cn, baf=fake.baf), colData=fake.pData )
baf.ds[, , "mbaf"] = baf2mbaf(baf.ds[, , "baf"], hom.cutoff = 0.90)

## ----mbaftorle----------------------------------------------------------------
mbaf.data = RleDataFrame( sapply(colnames( baf.ds),
  function(x) { Rle( baf.ds[,x, "mbaf"] ) },
  USE.NAMES=TRUE, simplify=FALSE ) )
as.numeric(object.size( baf.ds[, ,"mbaf"]))  / as.numeric( object.size(mbaf.data))

## ----segment------------------------------------------------------------------
baf.ds[,,"baf.segs"] = runCBS( baf.ds[, ,"mbaf"], rowRanges(baf.ds) )
baf.ds[,,"lrr.segs"] = runCBS( baf.ds[, , "lrr"], rowRanges(baf.ds) )

## ----plotlrr------------------------------------------------------------------
genoPlot(baf.ds,baf.ds[,1,"lrr"],chr="chr12",main="LRR of chr12")
genoPlot(baf.ds,baf.ds[,1,"lrr.segs"],chr="chr12",add=TRUE,col="red")

## ----plotbaf, echo=TRUE-------------------------------------------------------
par(mfrow=c(2,1))
genoPlot(baf.ds,baf.ds[,1,"baf"],chr="chr12", main="BAF of chr12")
genoPlot(baf.ds,baf.ds[,1,"mbaf"],chr="chr12", main="mBAF of chr12")
genoPlot(baf.ds,baf.ds[,1,"baf.segs"],chr="chr12", add=TRUE,col="red")

## ----crosssample--------------------------------------------------------------
gain.list = lapply(colnames(baf.ds),
  function(sample.name) {
    as.logical( baf.ds[, sample.name, "lrr.segs"] > 0.3 )
})
gain.mat = do.call(cbind, gain.list)
gain.freq = rowMeans(gain.mat,na.rm=TRUE)