## Load data
library('derfinder')
## Annotate regions, first two regions only
regions <- genomeRegions$regions[1:2]
annotatedRegions <- annotateRegions(regions = regions,
genomicState = genomicState$fullGenome, minoverlap = 1)
## Find nearest annotation with bumphunter::matchGenes()
library('bumphunter')
library('TxDb.Hsapiens.UCSC.hg19.knownGene')
genes <- annotateTranscripts(txdb = TxDb.Hsapiens.UCSC.hg19.knownGene)
nearestAnnotation <- matchGenes(x = regions, subject = genes)
## Obtain fullCov object
fullCov <- list('21'=genomeDataRaw$coverage)
## Assign chr lengths using hg19 information
library('GenomicRanges')
data(hg19Ideogram, package = 'biovizBase', envir = environment())
seqlengths(regions) <- seqlengths(hg19Ideogram)[names(seqlengths(regions))]
## Get the region coverage
regionCov <- getRegionCoverage(fullCov=fullCov, regions=regions)
## Make plots for the regions
plotRegionCoverage(regions=regions, regionCoverage=regionCov,
groupInfo=genomeInfo$pop, nearestAnnotation=nearestAnnotation,
annotatedRegions=annotatedRegions, whichRegions=1:2)
## Re-make plots with transcript information
plotRegionCoverage(regions=regions, regionCoverage=regionCov,
groupInfo=genomeInfo$pop, nearestAnnotation=nearestAnnotation,
annotatedRegions=annotatedRegions, whichRegions=1:2,
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene)
## If you prefer, you can save the plots to a pdf file
pdf('ders.pdf', h = 6, w = 9)
plotRegionCoverage(regions=regions, regionCoverage=regionCov,
groupInfo=genomeInfo$pop, nearestAnnotation=nearestAnnotation,
annotatedRegions=annotatedRegions, whichRegions=1:2,
txdb = TxDb.Hsapiens.UCSC.hg19.knownGene, ask = FALSE)
dev.off()
Run the code above in your browser using DataLab