inst/doc/GeneBreak.R

In GeneBreak: Gene Break Detection

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

###################################################
### code chunk number 1: loadingPackage
###################################################
library(GeneBreak)


###################################################
### code chunk number 2: settingOptions
###################################################
options("GeneBreak::verbose"=NA)
options(width=75)


###################################################
### code chunk number 3: loadingCopynumberData
###################################################
library(CGHcall)
data( "copynumber.data.chr20" )


###################################################
### code chunk number 4: displayCopynumberData
###################################################
copynumber.data.chr20


###################################################
### code chunk number 5: getBreakpoints
###################################################
breakpoints <- getBreakpoints( data = copynumber.data.chr20 )
breakpoints


###################################################
### code chunk number 6: getBreakpointsAlternative (eval = FALSE)
###################################################
## library(CGHcall)
## cgh <- copynumber.data.chr20
## 
## segmented <- data.frame( Chromosome=chromosomes(cgh), Start=bpstart(cgh),
##  End=bpend(cgh), FeatureName=rownames(cgh), segmented(cgh))
## called <- data.frame( Chromosome=chromosomes(cgh), Start=bpstart(cgh),
##  End=bpend(cgh), FeatureName=rownames(cgh), calls(cgh))
## 
## breakpoints <- getBreakpoints( data = segmented, data2 = called )


###################################################
### code chunk number 7: bpFilter
###################################################
breakpointsFiltered <- bpFilter( breakpoints, filter = "CNA-ass" )
breakpointsFiltered


###################################################
### code chunk number 8: loadingAnnotation
###################################################
data( "ens.gene.ann.hg18" )


###################################################
### code chunk number 9: showHeadGeneAnnotation
###################################################
head( ens.gene.ann.hg18 )


###################################################
### code chunk number 10: addGeneAnnotation
###################################################
breakpointsAnnotated <- addGeneAnnotation( breakpointsFiltered, ens.gene.ann.hg18 )


###################################################
### code chunk number 11: showFeatures_PCMTD2
###################################################
featuresPerGene ( breakpointsAnnotated , geneName = "PCMTD2" )


###################################################
### code chunk number 12: showGeneAssociatedFeatureInfo
###################################################
geneFeatures <- geneInfo( breakpointsAnnotated )
head( geneFeatures[ , 
  c("Gene", "Chromosome", "Start", "End", "featureTotal", 
    "featureNames", "remarks") ] )


###################################################
### code chunk number 13: bpGenes
###################################################
breakpointGenes <- bpGenes( breakpointsAnnotated )


###################################################
### code chunk number 14: headBreakpointGenes
###################################################
result_BreakpointGenes <- geneInfo ( breakpointGenes )
head( result_BreakpointGenes[ which ( result_BreakpointGenes$sampleCount > 0 ) , 
  c( "Gene", "Chromosome", "Start", "End", "featureTotal", "nrOfBreakLocations",
     "sampleCount", "sampleNamesWithBreakpoints") ] )


###################################################
### code chunk number 15: bpStats
###################################################
breakpointStatistics <- bpStats( breakpointGenes, 
  level = "gene", method = "Gilbert" )


###################################################
### code chunk number 16: recurrentGenes
###################################################
head( recurrentGenes( breakpointStatistics ) )


###################################################
### code chunk number 17: bpStats
###################################################
breakpointStatistics <- bpStats( 
  breakpointStatistics, level = "feature", method = "BH" )


###################################################
### code chunk number 18: showStatsObject
###################################################
breakpointStatistics


###################################################
### code chunk number 19: recurrentGenes
###################################################
head( featureInfo( breakpointStatistics ) )


###################################################
### code chunk number 20: GeneBreak.Rnw:252-253
###################################################
pdf("bpPlot.png", width=10)


###################################################
### code chunk number 21: bpPlot
###################################################
bpPlot( breakpointStatistics, fdr.threshold = 0.1 )


###################################################
### code chunk number 22: GeneBreak.Rnw:258-259
###################################################
dev.off()


###################################################
### code chunk number 23: createAnnotationExample (eval = FALSE)
###################################################
## # gene annotations obtained via Biomart. 
## # HUGO gene names (HGNC symbol), Ensembl_ID and chromosomal location
## 
## # Used (and most) recent releases:
## # HG18: release54
## # HG19: release75
## # HG38: release80 (date: 150629)
## 
## library(biomaRt)
## 
## ensembl54 = useMart( 
##   host = 'may2009.archive.ensembl.org', 
##   biomart = 'ENSEMBL_MART_ENSEMBL', 
##   dataset = "hsapiens_gene_ensembl" 
## )
## ensembl75 = useMart( 
##   host = 'feb2014.archive.ensembl.org', 
##   biomart = 'ENSEMBL_MART_ENSEMBL', 
##   dataset = "hsapiens_gene_ensembl" 
## )
## ensembl80 = useMart( 
##   "ensembl", 
##   dataset = "hsapiens_gene_ensembl" 
## )
## 
## createAnnotationFile <- function( biomartVersion ) {
##   biomart_result <- getBM( 
##     attributes =  c( 
##       "hgnc_symbol", "ensembl_gene_id", "chromosome_name",  
##       "start_position", "end_position", "band", "strand" 
##     ), 
##     mart = biomartVersion
##   )
## 
##   biomart_result[ ,3] <- as.vector( biomart_result[ ,3] )
##   idx_x <- biomart_result$chromosome_name == "X"
##   idx_y <- biomart_result$chromosome_name == "Y"
##   biomart_result$chromosome_name[ idx_x ] <- "23"
##   biomart_result$chromosome_name[ idx_y ] <- "24"
##   
##   biomart_genes <- biomart_result[ which(biomart_result[ ,1] != "" & 
##     biomart_result[ ,3] %in% c(1:24)) , ]
##   colnames(biomart_genes)[1:5] <- c("Gene","EnsID","Chromosome","Start","End")
##   
##   cat( 
##     c( "Biomart version:", biomartVersion@host, 
##           "including:", dim(biomart_genes)[1], "genes\n"
##     )
##   ) 
##   
##   return( biomart_genes )
## }
## 
## ens.gene.ann.hg18 <- createAnnotationFile( ensembl54 )
## ens.gene.ann.hg19 <- createAnnotationFile( ensembl75 )
## ens.gene.ann.hg38 <- createAnnotationFile( ensembl80 )
## 


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