GenomeGraphs.R
In GenomeGraphs: Plotting genomic information from Ensembl

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

###################################################
### code chunk number 1: GenomeGraphs.Rnw:18-19
###################################################
options(width=50)


###################################################
### code chunk number 2: GenomeGraphs.Rnw:41-42
###################################################
library(GenomeGraphs)


###################################################
### code chunk number 3: GenomeGraphs.Rnw:45-51 (eval = FALSE)
###################################################
## require(biomaRt)
## require(grid)
## source("../../R/GenomeGraphs-classes.R")
## source("../../R/GenomeGraphs-methods.R")
## source("../../R/Overlay.R")
## source("../../R/GenomeGraphs.R")


###################################################
### code chunk number 4: GenomeGraphs.Rnw:68-69
###################################################
mart <- useMart("ensembl", dataset="hsapiens_gene_ensembl")


###################################################
### code chunk number 5: GenomeGraphs.Rnw:74-76
###################################################
gene <- makeGene(id = "ENSG00000095203", type="ensembl_gene_id", biomart = mart)
gdPlot(gene) 


###################################################
### code chunk number 6: GenomeGraphs.Rnw:83-85
###################################################
transcript <- makeTranscript(id = "ENSG00000095203", type="ensembl_gene_id", biomart = mart)
gdPlot(list(gene, transcript))


###################################################
### code chunk number 7: GenomeGraphs.Rnw:96-99
###################################################
plusStrand <- makeGeneRegion(chromosome = 19, start = 12050000, end = 12230000, strand = "+", biomart = mart)
genomeAxis <- makeGenomeAxis(add53 = TRUE)
gdPlot(list(genomeAxis, plusStrand))


###################################################
### code chunk number 8: GenomeGraphs.Rnw:105-109
###################################################
minStrand <- makeGeneRegion( chromosome = 19, start = 12050000, end = 12230000, strand = "-", biomart = mart)
ideogram <- makeIdeogram(chromosome = 19)
genomeAxis <- makeGenomeAxis(add53=TRUE, add35=TRUE)
gdPlot(list(ideogram, plusStrand, genomeAxis, minStrand))


###################################################
### code chunk number 9: GenomeGraphs.Rnw:123-146
###################################################
data("exampleData", package="GenomeGraphs")

minbase <- 180292097 
maxbase <- 180492096

genesplus <- makeGeneRegion(start = minbase, end = maxbase, 
                            strand = "+", chromosome = "3", biomart=mart)
genesmin <- makeGeneRegion(start = minbase, end = maxbase, 
                           strand = "-", chromosome = "3", biomart=mart)

seg <- makeSegmentation(segStart[[1]], segEnd[[1]], segments[[1]], 
                        dp = DisplayPars(color = "black", lwd=2,lty = "solid"))

cop <- makeGenericArray(intensity  = cn, probeStart = probestart, 
                        trackOverlay =  seg, dp = DisplayPars(size=3, color = "seagreen", type="dot"))

ideog <- makeIdeogram(chromosome = 3)
expres <- makeGenericArray(intensity = intensity, probeStart = exonProbePos, 
                           dp = DisplayPars(color="darkred", type="point"))
genomeAxis <- makeGenomeAxis(add53 = TRUE, add35=TRUE)

gdPlot(list(a=ideog,b=expres,c=cop,d=genesplus,e=genomeAxis,f=genesmin), 
       minBase = minbase, maxBase =maxbase, labelCex = 2)


###################################################
### code chunk number 10: GenomeGraphs.Rnw:154-175
###################################################
data("unrData", package="GenomeGraphs")

title <- makeTitle(text ="ENSG00000009307", color = "darkred")
exon <- makeExonArray(intensity = unrData, probeStart = unrPositions[,3], 
            probeEnd=unrPositions[,4], probeId = as.character(unrPositions[,1]), 
            nProbes = unrNProbes, dp = DisplayPars(color = "blue", mapColor = "dodgerblue2"), 
            displayProbesets=FALSE)
affyModel.model <- makeGeneModel(start = unrPositions[,3], end = unrPositions[,4])
affyModel <- makeAnnotationTrack(start = unrPositions[,3], end = unrPositions[,4],
                                 feature = "gene_model", group = "ENSG00000009307", 
                                 dp = DisplayPars(gene_model = "darkblue"))

gene <- makeGene(id = "ENSG00000009307", biomart = mart)


transcript <- makeTranscript( id ="ENSG00000009307" , biomart = mart)
legend <- makeLegend(c("affyModel","gene"), fill = c("darkgreen","orange"))
rOverlay <- makeRectangleOverlay(start = 115085100, end = 115086500, region = c(3,5),  
                                 dp = DisplayPars(alpha = .2, fill = "olivedrab1"))
gdPlot(list(title, exon, affyModel, gene, transcript, legend), 
       minBase = 115061061, maxBase=115102147, overlay = rOverlay)


###################################################
### code chunk number 11: GenomeGraphs.Rnw:183-206
###################################################
yeastMart <- useMart("ensembl", dataset = "scerevisiae_gene_ensembl")
minB <- 10000
maxB <- 20000

chrRoman <- as.character(as.roman(1))
grP <- makeGeneRegion(start = minB, end = maxB, strand = "+", 
                      chromosome = chrRoman, biomart = yeastMart)
grM <- makeGeneRegion(start = minB, end = maxB, strand = "-", 
                      chromosome = chrRoman, biomart = yeastMart)
gaxis <- makeGenomeAxis(add53 = TRUE, add35 = TRUE)
conserv <- yeastCons1[yeastCons1[,1] > minB & yeastCons1[,1] < maxB, ]

s1 <- makeSmoothing(x = lowess(conserv[,1], conserv[,2], f = .01)$x,
                    y = lowess(conserv[,1], conserv[,2], f = .01)$y, 
                    dp = DisplayPars(lwd = 3, color = "green"))
s2 <- makeSmoothing(x = lowess(conserv[,1], conserv[,2], f = .1)$x,
                    y = lowess(conserv[,1], conserv[,2], f = .1)$y, 
                    dp = DisplayPars(lwd = 3, color = "purple"))

consTrack <- makeBaseTrack(base = conserv[, 1], value = conserv[,2],
                           dp = DisplayPars(lwd=.2, ylim = c(0, 1.25), 
                           color = "darkblue"), trackOverlay = list(s1, s2))
gdPlot(list(grP, gaxis, grM, "conservation" = consTrack)) 


###################################################
### code chunk number 12: GenomeGraphs.Rnw:213-223
###################################################
plotGeneRegion <- function(chr = 1, minB = 9000, maxB = 13000, rot = 0, col = "green") {
  chrRoman <- as.character(as.roman(1:17)[chr])
  grP <- makeGeneRegion(start = minB, end = maxB, 
                        strand = "+", chromosome = chrRoman, biomart = yeastMart, 
                        dp = DisplayPars(plotId = TRUE, idRotation = rot, 
                        idColor = col))
  gaxis <- makeGenomeAxis( add53 = TRUE, add35 = TRUE, littleTicks = FALSE)
  gdPlot(list(grP, gaxis), minBase = minB, maxBase = maxB)
}
plotGeneRegion(col = "yellow", rot=90)


###################################################
### code chunk number 13: GenomeGraphs.Rnw:229-232
###################################################
gdPlot(makeGeneRegion(start = 9000, end = 15000, biomart = yeastMart,
                      strand = "-", chromosome = "I", 
                      dp = DisplayPars(plotId=TRUE)))


###################################################
### code chunk number 14: GenomeGraphs.Rnw:239-246
###################################################
ga <- makeGenomeAxis()
grF <- makeGeneRegion(start = 10000, end = 20000, chromosome = "I", strand = "+", biomart = yeastMart)
grR <- makeGeneRegion(start = 10000, end = 20000, chromosome = "I", strand = "-", biomart = yeastMart)
bt <- makeBaseTrack(base = yeastCons1[,1], value = yeastCons1[,2])
hr1 <- makeRectangleOverlay(start = 11000, end = 13000)
hr2 <- makeRectangleOverlay(start = 15900, end = 16500)
gdPlot(list(grF, ga, grR, bt), overlays = list(hr1, hr2))


###################################################
### code chunk number 15: GenomeGraphs.Rnw:253-257
###################################################
ro <- makeRectangleOverlay(start = 11000, end = 13000, region = c(1,3), 
          dp = DisplayPars(color = "green", alpha = .3))
to <- makeTextOverlay("here is some text", xpos = 15000, ypos = .95)
gdPlot(list(grF, ga, grR, bt), overlay = c(ro, to))


###################################################
### code chunk number 16: GenomeGraphs.Rnw:262-266
###################################################
roR <- makeRectangleOverlay(start = .1, end = .3, coords = "absolute",
                            dp = DisplayPars(fill = "grey", alpha = .2, lty = "dashed"),
                            region = c(.4,.7))
gdPlot(list(grF, ga, grR, bt), overlays = list(ro, roR))


###################################################
### code chunk number 17: GenomeGraphs.Rnw:302-328
###################################################
data("seqDataEx", package = "GenomeGraphs")
str = seqDataEx$david[,"strand"] == 1
biomart = useMart("ensembl", "scerevisiae_gene_ensembl")
pList = list("-" = makeGeneRegion(chromosome = "IV", start = 1300000, end = 1310000, 
                                   strand = "-", biomart = biomart, 
                                   dp = DisplayPars(plotId = TRUE, idRotation = 0, cex = .5)),
              makeGenomeAxis(dp = DisplayPars(size = 3)),
              "+" = makeGeneRegion(chromosome = "IV", start = 1300000, end = 1310000, 
                                   strand = "+", biomart = biomart, 
                                   dp = DisplayPars(plotId = TRUE, idRotation = 0, cex = .5)),
              "Nagalakshmi" = makeBaseTrack(base = seqDataEx$snyder[, "location"], value = seqDataEx$snyder[, "counts"], 
                                            dp = DisplayPars(lwd = .3, color = "darkblue", ylim = c(0,300))),
              "David +" = makeGenericArray(probeStart = seqDataEx$david[str, "location"], 
                                           intensity = seqDataEx$david[str, "expr", drop=FALSE], 
                                           dp = DisplayPars(pointSize = .5)),
              "David -" = makeGenericArray(probeStart = seqDataEx$david[!str, "location"], 
                                           intensity = seqDataEx$david[!str, "expr", drop=FALSE], 
                dp = DisplayPars(color = "darkgreen", pointSize = .5)),
              "Lee" = makeBaseTrack(base = seqDataEx$nislow[, "location"], 
                                    value = seqDataEx$nislow[, "evalue"], dp = DisplayPars(color="grey", lwd=.25)),
              "Conservation" = makeBaseTrack(base = seqDataEx$conservation[, "location"], 
                                             value = seqDataEx$conservation[, "score"], 
                                             dp = DisplayPars(color="gold4", lwd=.25)))

gdPlot(pList, minBase = 1301500, maxBase = 1302500, 
       overlay = makeRectangleOverlay(start = 1302105, end = 1302190, region = c(4,8),  dp = DisplayPars(alpha = .2)))


###################################################
### code chunk number 18: GenomeGraphs.Rnw:333-340
###################################################
setPar(pList$Lee@dp, "type", "h")
setPar(pList$Lee@dp, "color", "limegreen")
setPar(pList$Lee@dp, "lwd", 2)

gdPlot(pList, minBase = 1301500, maxBase = 1302500, 
       overlay = makeRectangleOverlay(start = 1302105, end = 1302190, region = c(4,8),  
       dp = DisplayPars(alpha = .2)))