1.1 Install derfinderPlot

R is an open-source statistical environment which can be easily modified to enhance its functionality via packages. derfinderPlot is a R package available via the Bioconductor repository for packages. R can be installed on any operating system from CRAN after which you can install derfinderPlot by using the following commands in your R session:

install.packages("BiocManager")
BiocManager::install("derfinderPlot")

## Check that you have a valid Bioconductor installation
BiocManager::valid()

1.2 Required knowledge

derfinderPlot is based on many other packages and in particular in those that have implemented the infrastructure needed for dealing with RNA-seq data. A derfinderPlot user is not expected to deal with those packages directly but will need to be familiar with derfinder and for some plots with ggbio.

If you are asking yourself the question “Where do I start using Bioconductor?” you might be interested in this blog post.

1.3 Asking for help

As package developers, we try to explain clearly how to use our packages and in which order to use the functions. But R and Bioconductor have a steep learning curve so it is critical to learn where to ask for help. The blog post quoted above mentions some but we would like to highlight the Bioconductor support site as the main resource for getting help: remember to use the derfinder or derfinderPlot tags and check the older posts. Other alternatives are available such as creating GitHub issues and tweeting. However, please note that if you want to receive help you should adhere to the posting guidelines. It is particularly critical that you provide a small reproducible example and your session information so package developers can track down the source of the error.

1.4 Citing derfinderPlot

We hope that derfinderPlot will be useful for your research. Please use the following information to cite the package and the overall approach. Thank you!

## Citation info
citation('derfinderPlot')

## 
## Collado-Torres L, Jaffe AE, Leek JT (2017). _derfinderPlot: Plotting
## functions for derfinder_. doi: 10.18129/B9.bioc.derfinderPlot (URL:
## https://doi.org/10.18129/B9.bioc.derfinderPlot),
## https://github.com/leekgroup/derfinderPlot - R package version
## 1.18.1, <URL: http://www.bioconductor.org/packages/derfinderPlot>.
## 
## Collado-Torres L, Nellore A, Frazee AC, Wilks C, Love MI, Langmead
## B, Irizarry RA, Leek JT, Jaffe AE (2017). "Flexible expressed region
## analysis for RNA-seq with derfinder." _Nucl. Acids Res._. doi:
## 10.1093/nar/gkw852 (URL: https://doi.org/10.1093/nar/gkw852), <URL:
## http://nar.oxfordjournals.org/content/early/2016/09/29/nar.gkw852>.
## 
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.

## Load libraries
suppressPackageStartupMessages(library('derfinder'))
library('derfinderData')
library('derfinderPlot')

## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

3.1 Analyze data

For this example, we created a small table with the relevant phenotype data for 12 samples: 6 from fetal samples and 6 from adult samples. We chose at random a brain region, in this case the primary auditory cortex (core) and for the example we will only look at data from chromosome 21. Other variables include the age in years and the gender. The data is shown below.

library('knitr')
## Get pheno table
pheno <- subset(brainspanPheno, structure_acronym == 'A1C')

## Display the main information
p <- pheno[, -which(colnames(pheno) %in% c('structure_acronym',
    'structure_name', 'file'))]
rownames(p) <- NULL
kable(p, format = 'html', row.names = TRUE)

	gender	lab	Age	group
1	M	HSB114.A1C	-0.5192308	fetal
2	M	HSB103.A1C	-0.5192308	fetal
3	M	HSB178.A1C	-0.4615385	fetal
4	M	HSB154.A1C	-0.4615385	fetal
5	F	HSB150.A1C	-0.5384615	fetal
6	F	HSB149.A1C	-0.5192308	fetal
7	F	HSB130.A1C	21.0000000	adult
8	M	HSB136.A1C	23.0000000	adult
9	F	HSB126.A1C	30.0000000	adult
10	M	HSB145.A1C	36.0000000	adult
11	M	HSB123.A1C	37.0000000	adult
12	F	HSB135.A1C	40.0000000	adult

We can load the data from derfinderData (Collado-Torres, Jaffe, and Leek, 2019) by first identifying the paths to the BigWig files with derfinder::rawFiles() and then loading the data with derfinder::fullCoverage().

## Determine the files to use and fix the names
files <- rawFiles(system.file('extdata', 'A1C', package = 'derfinderData'),
    samplepatt = 'bw', fileterm = NULL)
names(files) <- gsub('.bw', '', names(files))

## Load the data from disk
system.time(fullCov <- fullCoverage(files = files, chrs = 'chr21'))

Alternatively, since the BigWig files are publicly available from BrainSpan (see here), we can extract the relevant coverage data using derfinder::fullCoverage(). Note that as of rtracklayer 1.25.16 BigWig files are not supported on Windows: you can find the fullCov object inside derfinderData to follow the examples.

## Determine the files to use and fix the names
files <- pheno$file
names(files) <- gsub('.A1C', '', pheno$lab)

## Load the data from the web
system.time(fullCov <- fullCoverage(files = files, chrs = 'chr21'))

Once we have the base-level coverage data for all 12 samples, we can construct the models. In this case, we want to find differences between fetal and adult samples while adjusting for gender and a proxy of the library size.

## Get some idea of the library sizes
sampleDepths <- sampleDepth(collapseFullCoverage(fullCov), 1)

## 2019-05-22 23:03:47 sampleDepth: Calculating sample quantiles

## 2019-05-22 23:03:47 sampleDepth: Calculating sample adjustments

## Define models
models <- makeModels(sampleDepths, testvars = pheno$group,
    adjustvars = pheno[, c('gender')])

Next, we can find candidate differentially expressed regions (DERs) using as input the segments of the genome where at least one sample has coverage greater than 3. In this particular example, we chose a low theoretical F-statistic cutoff and used 20 permutations.

## Filter coverage
filteredCov <- lapply(fullCov, filterData, cutoff = 3)

## 2019-05-22 23:03:49 filterData: originally there were 48129895 rows, now there are 90023 rows. Meaning that 99.81 percent was filtered.

## Perform differential expression analysis
suppressPackageStartupMessages(library('bumphunter'))
system.time(results <- analyzeChr(chr = 'chr21', filteredCov$chr21,
    models, groupInfo = pheno$group, writeOutput = FALSE,
    cutoffFstat = 5e-02, nPermute = 20, seeds = 20140923 + seq_len(20)))

## 2019-05-22 23:03:50 analyzeChr: Pre-processing the coverage data

## 2019-05-22 23:03:52 analyzeChr: Calculating statistics

## 2019-05-22 23:03:52 calculateStats: calculating the F-statistics

## 2019-05-22 23:03:52 analyzeChr: Calculating pvalues

## 2019-05-22 23:03:52 analyzeChr: Using the following theoretical cutoff for the F-statistics 5.31765507157871

## 2019-05-22 23:03:52 calculatePvalues: identifying data segments

## 2019-05-22 23:03:52 findRegions: segmenting information

## 2019-05-22 23:03:52 findRegions: identifying candidate regions

## 2019-05-22 23:03:52 findRegions: identifying region clusters

## 2019-05-22 23:03:52 calculatePvalues: calculating F-statistics for permutation 1 and seed 20140924

## 2019-05-22 23:03:53 findRegions: segmenting information

## 2019-05-22 23:03:53 findRegions: identifying candidate regions

## 2019-05-22 23:03:53 calculatePvalues: calculating F-statistics for permutation 2 and seed 20140925

## 2019-05-22 23:03:54 findRegions: segmenting information

## 2019-05-22 23:03:54 findRegions: identifying candidate regions

## 2019-05-22 23:03:54 calculatePvalues: calculating F-statistics for permutation 3 and seed 20140926

## 2019-05-22 23:03:54 findRegions: segmenting information

## 2019-05-22 23:03:54 findRegions: identifying candidate regions

## 2019-05-22 23:03:54 calculatePvalues: calculating F-statistics for permutation 4 and seed 20140927

## 2019-05-22 23:03:54 findRegions: segmenting information

## 2019-05-22 23:03:54 findRegions: identifying candidate regions

## 2019-05-22 23:03:54 calculatePvalues: calculating F-statistics for permutation 5 and seed 20140928

## 2019-05-22 23:03:54 findRegions: segmenting information

## 2019-05-22 23:03:54 findRegions: identifying candidate regions

## 2019-05-22 23:03:54 calculatePvalues: calculating F-statistics for permutation 6 and seed 20140929

## 2019-05-22 23:03:54 findRegions: segmenting information

## 2019-05-22 23:03:54 findRegions: identifying candidate regions

## 2019-05-22 23:03:54 calculatePvalues: calculating F-statistics for permutation 7 and seed 20140930

## 2019-05-22 23:03:55 findRegions: segmenting information

## 2019-05-22 23:03:55 findRegions: identifying candidate regions

## 2019-05-22 23:03:55 calculatePvalues: calculating F-statistics for permutation 8 and seed 20140931

## 2019-05-22 23:03:55 findRegions: segmenting information

## 2019-05-22 23:03:55 findRegions: identifying candidate regions

## 2019-05-22 23:03:55 calculatePvalues: calculating F-statistics for permutation 9 and seed 20140932

## 2019-05-22 23:03:55 findRegions: segmenting information

## 2019-05-22 23:03:55 findRegions: identifying candidate regions

## 2019-05-22 23:03:55 calculatePvalues: calculating F-statistics for permutation 10 and seed 20140933

## 2019-05-22 23:03:55 findRegions: segmenting information

## 2019-05-22 23:03:55 findRegions: identifying candidate regions

## 2019-05-22 23:03:55 calculatePvalues: calculating F-statistics for permutation 11 and seed 20140934

## 2019-05-22 23:03:55 findRegions: segmenting information

## 2019-05-22 23:03:55 findRegions: identifying candidate regions

## 2019-05-22 23:03:56 calculatePvalues: calculating F-statistics for permutation 12 and seed 20140935

## 2019-05-22 23:03:56 findRegions: segmenting information

## 2019-05-22 23:03:56 findRegions: identifying candidate regions

## 2019-05-22 23:03:56 calculatePvalues: calculating F-statistics for permutation 13 and seed 20140936

## 2019-05-22 23:03:56 findRegions: segmenting information

## 2019-05-22 23:03:56 findRegions: identifying candidate regions

## 2019-05-22 23:03:56 calculatePvalues: calculating F-statistics for permutation 14 and seed 20140937

## 2019-05-22 23:03:56 findRegions: segmenting information

## 2019-05-22 23:03:56 findRegions: identifying candidate regions

## 2019-05-22 23:03:56 calculatePvalues: calculating F-statistics for permutation 15 and seed 20140938

## 2019-05-22 23:03:56 findRegions: segmenting information

## 2019-05-22 23:03:56 findRegions: identifying candidate regions

## 2019-05-22 23:03:56 calculatePvalues: calculating F-statistics for permutation 16 and seed 20140939

## 2019-05-22 23:03:57 findRegions: segmenting information

## 2019-05-22 23:03:57 findRegions: identifying candidate regions

## 2019-05-22 23:03:57 calculatePvalues: calculating F-statistics for permutation 17 and seed 20140940

## 2019-05-22 23:03:57 findRegions: segmenting information

## 2019-05-22 23:03:57 findRegions: identifying candidate regions

## 2019-05-22 23:03:57 calculatePvalues: calculating F-statistics for permutation 18 and seed 20140941

## 2019-05-22 23:03:57 findRegions: segmenting information

## 2019-05-22 23:03:57 findRegions: identifying candidate regions

## 2019-05-22 23:03:57 calculatePvalues: calculating F-statistics for permutation 19 and seed 20140942

## 2019-05-22 23:03:58 findRegions: segmenting information

## 2019-05-22 23:03:58 findRegions: identifying candidate regions

## 2019-05-22 23:03:58 calculatePvalues: calculating F-statistics for permutation 20 and seed 20140943

## 2019-05-22 23:03:58 findRegions: segmenting information

## 2019-05-22 23:03:58 findRegions: identifying candidate regions

## 2019-05-22 23:03:58 calculatePvalues: calculating the p-values

## 2019-05-22 23:03:58 analyzeChr: Annotating regions

## No annotationPackage supplied. Trying org.Hs.eg.db.

## Loading required package: org.Hs.eg.db

## Loading required package: AnnotationDbi

## Loading required package: Biobase

## Welcome to Bioconductor
## 
##     Vignettes contain introductory material; view with
##     'browseVignettes()'. To cite Bioconductor, see
##     'citation("Biobase")', and for packages 'citation("pkgname")'.

## 

## Getting TSS and TSE.

## Getting CSS and CSE.

## Getting exons.

## Annotating genes.

## .............

##    user  system elapsed 
##   88.45    1.51   90.01

## Quick access to the results
regions <- results$regions$regions

## Annotation database to use
suppressPackageStartupMessages(library('TxDb.Hsapiens.UCSC.hg19.knownGene'))
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene

3.2 plotOverview()

Now that we have obtained the main results using derfinder, we can proceed to visualizing the results using derfinderPlot. The easiest to use of all the functions is plotOverview() which takes a set of regions and annotation information produced by bumphunter::matchGenes().

Figure 1 shows the candidate DERs colored by whether their q-value was less than 0.10 or not.

## Q-values overview
plotOverview(regions = regions, annotation = results$annotation, type = 'qval')

## 2019-05-22 23:05:20 plotOverview: assigning chromosome lengths from hg19!!!

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

Figure 2 shows the candidate DERs colored by the type of gene feature they are nearest too.

## Annotation overview
plotOverview(regions = regions, annotation = results$annotation,
    type = 'annotation')

## 2019-05-22 23:05:22 plotOverview: assigning chromosome lengths from hg19!!!

## Scale for 'x' is already present. Adding another scale for 'x', which will
## replace the existing scale.

In this particular example, because we are only using data from one chromosome the above plot is not as informative as in a real case scenario. However, with this plot we can quickly observe that nearly all of the candidate DERs are inside an exon.

3.3 plotRegionCoverage()

The complete opposite of visualizing the candidate DERs at the genome-level is to visualize them one region at a time. plotRegionCoverage() allows us to do this quickly for a large number of regions.

Before using this function, we need to process more detailed information using two derfinder functions: annotateRegions() and getRegionCoverage() as shown below.

## Get required information for the plots
annoRegs <- annotateRegions(regions, genomicState$fullGenome)

## 2019-05-22 23:05:23 annotateRegions: counting

## 2019-05-22 23:05:23 annotateRegions: annotating

regionCov <- getRegionCoverage(fullCov, regions)

## 2019-05-22 23:05:24 getRegionCoverage: processing chr21

## 2019-05-22 23:05:24 getRegionCoverage: done processing chr21

Once we have the relevant information we can proceed to plotting the first 10 regions. In this case, we will supply plotRegionCoverage() with the information it needs to plot transcripts overlapping these 10 regions (Figures 3, 4, 5, 6, 7, 8, 9, 10, 11, 12).

## Plot top 10 regions
plotRegionCoverage(regions = regions, regionCoverage = regionCov, 
    groupInfo = pheno$group, nearestAnnotation = results$annotation, 
    annotatedRegions = annoRegs, whichRegions=1:10, txdb = txdb, scalefac = 1,
    ask = FALSE, verbose = FALSE)

## Warning in matrix(c(-pad, -pad, pad, pad), nrow = nrow(d), ncol = 4, byrow =
## TRUE): data length exceeds size of matrix

The base-level coverage is shown in a log2 scale with any overlapping exons shown in dark blue and known introns in light blue.

3.4 plotCluster()

In this example, we noticed with the plotRegionCoverage() plots that most of the candidate DERs are contained in known exons. Sometimes, the signal might be low or we might have used very stringent cutoffs in the derfinder analysis. One way we can observe this is by plotting clusters of regions where a cluster is defined as regions within 300 bp (default option) of each other.

To visualize the clusters, we can use plotCluster() which takes similar input to plotOverview() with the notable addition of the coverage information as well as the idx argument. This argument specifies which region to focus on: it will be plotted with a red bar and will determine the cluster to display.

In Figure 13 we observe one large candidate DER with other nearby ones that do not have a q-value less than 0.10. In a real analysis, we would probably discard this region as the coverage is fairly low.

## First cluster
plotCluster(idx = 1, regions = regions, annotation = results$annotation,
    coverageInfo = fullCov$chr21, txdb = txdb, groupInfo = pheno$group,
    titleUse = 'pval')

## Parsing transcripts...

## Parsing exons...

## Parsing cds...

## Parsing utrs...

## ------exons...

## ------cdss...

## ------introns...

## ------utr...

## aggregating...

## Done

## Constructing graphics...

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

The second cluster (Figure 14) shows a larger number of potential DERs (again without q-values less than 0.10) in a segment of the genome where the coverage data is highly variable. This is a common occurrence with RNA-seq data.

## Second cluster
plotCluster(idx = 2, regions = regions, annotation = results$annotation,
    coverageInfo = fullCov$chr21, txdb = txdb, groupInfo = pheno$group,
    titleUse = 'pval')

## Warning in getIdeoGR(data): geom(ideogram) need valid seqlengths information for accurate mapping,
##                  now use reduced information as ideogram...

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Parsing transcripts...

## Parsing exons...

## Parsing cds...

## Parsing utrs...

## ------exons...

## ------cdss...

## ------introns...

## ------utr...

## aggregating...

## Done

## Constructing graphics...

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Warning: `panel.margin` is deprecated. Please use `panel.spacing` property
## instead

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

## Warning: Transformation introduced infinite values in continuous y-axis

## Warning: Transformation introduced infinite values in continuous y-axis

These plots show an ideogram which helps quickly identify which region of the genome we are focusing on. Then, the base-level coverage information for each sample is displayed in log2. Next, the coverage group means are shown in the log2 scale. The plot is completed with the potential and candidate DERs as well as any known transcripts.

3.5 vennRegions

derfinder has functions for annotating regions given their genomic state. A typical visualization is to then view how many regions overlap known exons, introns, intergenic regions, none of them or several of these groups in a venn diagram. The function vennRegions() makes this plot using the output from derfinder::annotateRegions() as shown in Figure 15.

## Make venn diagram
venn <- vennRegions(annoRegs)

## It returns the actual venn counts information
venn

##   exon intergenic intron Counts
## 1    0          0      0      0
## 2    0          0      1     11
## 3    0          1      0     13
## 4    0          1      1      0
## 5    1          0      0   1175
## 6    1          0      1    113
## 7    1          1      0      0
## 8    1          1      1      0
## attr(,"class")
## [1] "VennCounts"