1 Basics

1.1 Install derfinder

R is an open-source statistical environment which can be easily modified to enhance its functionality via packages. derfinder 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 derfinder by using the following commands in your R session:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
      install.packages("BiocManager")
  }

BiocManager::install("derfinder")

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

1.2 Required knowledge

derfinder is based on many other packages and in particular in those that have implemented the infrastructure needed for dealing with RNA-seq data. That is, packages like Rsamtools, GenomicAlignments and rtracklayer that allow you to import the data. A derfinder user is not expected to deal with those packages directly but will need to be familiar with GenomicRanges to understand the results derfinder generates. It might also prove to be highly beneficial to check the BiocParallel package for performing parallel computations.

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 tag 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.

We would like to highlight the derfinder user Jessica Hekman. She has used derfinder with non-human data, and in the process of doing so discovered some small bugs or sections of the documentation that were not clear.

1.4 Citing derfinder

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

## Citation info
citation("derfinder")
## 
## 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>.
## 
## Frazee AC, Sabunciyan S, Hansen KD, Irizarry RA, Leek JT (2014).
## "Differential expression analysis of RNA-seq data at single-base
## resolution." _Biostatistics_, *15 (3)*, 413-426. doi:
## 10.1093/biostatistics/kxt053 (URL:
## https://doi.org/10.1093/biostatistics/kxt053), <URL:
## http://biostatistics.oxfordjournals.org/content/15/3/413.long>.
## 
## Collado-Torres L, Jaffe AE, Leek JT (2017). _derfinder:
## Annotation-agnostic differential expression analysis of RNA-seq data at
## base-pair resolution via the DER Finder approach_. doi:
## 10.18129/B9.bioc.derfinder (URL:
## https://doi.org/10.18129/B9.bioc.derfinder),
## https://github.com/lcolladotor/derfinder - R package version 1.24.2,
## <URL: http://www.bioconductor.org/packages/derfinder>.
## 
## To see these entries in BibTeX format, use 'print(<citation>,
## bibtex=TRUE)', 'toBibtex(.)', or set
## 'options(citation.bibtex.max=999)'.

2 Quick start to using to derfinder

Here is a very quick example of a DER Finder analysis. This analysis is explained in more detail later on in this document.

## Load libraries
library("derfinder")
library("derfinderData")
library("GenomicRanges")

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

## Load the data from disk -- On Windows you have to load data from Bam files
fullCov <- fullCoverage(files = files, chrs = "chr21", verbose = FALSE)

## Get the region matrix of Expressed Regions (ERs)
regionMat <- regionMatrix(fullCov, cutoff = 30, L = 76, verbose = FALSE)

## Get pheno table
pheno <- subset(brainspanPheno, structure_acronym == "AMY")

## Identify which ERs are differentially expressed, that is, find the DERs
library("DESeq2")

## Round matrix
counts <- round(regionMat$chr21$coverageMatrix)

## Round matrix and specify design
dse <- DESeqDataSetFromMatrix(counts, pheno, ~ group + gender)

## Perform DE analysis
dse <- DESeq(dse, test = "LRT", reduced = ~gender, fitType = "local")

## Extract results
mcols(regionMat$chr21$regions) <- c(mcols(regionMat$chr21$regions), results(dse))

## Save info in an object with a shorter name
ers <- regionMat$chr21$regions
ers

3 Introduction

derfinder is an R package that implements the DER Finder approach (Frazee, Sabunciyan, Hansen, Irizarry, and Leek, 2014) for RNA-seq data. Briefly, this approach is an alternative to feature-counting and transcript assembly. The basic idea is to identify contiguous base-pairs in the genome that present differential expression signal. These base-pairs are grouped into _d_ifferentially _e_xpressed _r_regions (DERs). This approach is annotation-agnostic which is a feature you might be interested in. In particular, derfinder contains functions that allow you to identify DERs via two alternative methods. You can find more details on the full derfinder users guide.

4 Sample DER Finder analysis

This is a brief overview of what a DER Finder analysis looks like. In particular, here we will be identifying expressed regions (ERs) without relying on annotation. Next, we’ll identify candidate differentially expressed regions (DERs). Finally, we’ll compare the DERs with known annotation features.

We first load the required packages.

## Load libraries
library("derfinder")
library("derfinderData")
library("GenomicRanges")

Next, we need to locate the chromosome 21 coverage files for a set of 12 samples. These samples are a small subset from the BrainSpan Atlas of the Human Brain (BrainSpan, 2011) publicly available data. The function rawFiles() helps us in locating these files.

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

Next, we can load the full coverage data into memory using the fullCoverage() function. Note that the BrainSpan data is already normalized by the total number of mapped reads in each sample. However, that won’t be the case with most data sets in which case you might want to use the totalMapped and targetSize arguments. The function getTotalMapped() will be helpful to get this information.

## Load the data from disk
fullCov <- fullCoverage(
    files = files, chrs = "chr21", verbose = FALSE,
    totalMapped = rep(1, length(files)), targetSize = 1
)

Now that we have the data, we can identify expressed regions (ERs) by using a cutoff of 30 on the base-level mean coverage from these 12 samples. Once the regions have been identified, we can calculate a coverage matrix with one row per ER and one column per sample (12 in this case). For doing this calculation we need to know the length of the sequence reads, which in this study were 76 bp long.

## Get the region matrix of Expressed Regions (ERs)
regionMat <- regionMatrix(fullCov, cutoff = 30, L = 76, verbose = FALSE)

regionMatrix() returns a list of elements, each with three pieces of output. The actual ERs are arranged in a GRanges object named regions.

## regions output
regionMat$chr21$regions
## GRanges object with 45 ranges and 6 metadata columns:
##      seqnames            ranges strand |     value      area indexStart
##         <Rle>         <IRanges>  <Rle> | <numeric> <numeric>  <integer>
##    1    chr21   9827018-9827582      * |  313.6717 177224.53          1
##    2    chr21 15457301-15457438      * |  215.0846  29681.68        566
##    3    chr21 20230140-20230192      * |   38.8325   2058.12        704
##    4    chr21 20230445-20230505      * |   41.3245   2520.80        757
##    5    chr21 27253318-27253543      * |   34.9131   7890.37        818
##   ..      ...               ...    ... .       ...       ...        ...
##   41    chr21 33039644-33039688      * |   34.4705 1551.1742       2180
##   42    chr21 33040784-33040798      * |   32.1342  482.0133       2225
##   43    chr21          33040890      * |   30.0925   30.0925       2240
##   44    chr21 33040900-33040901      * |   30.1208   60.2417       2241
##   45    chr21 48019401-48019414      * |   31.1489  436.0850       2243
##       indexEnd cluster clusterL
##      <integer>   <Rle>    <Rle>
##    1       565       1      565
##    2       703       2      138
##    3       756       3      366
##    4       817       3      366
##    5      1043       4      765
##   ..       ...     ...      ...
##   41      2224      17       45
##   42      2239      18      118
##   43      2240      18      118
##   44      2242      18      118
##   45      2256      19       14
##   -------
##   seqinfo: 1 sequence from an unspecified genome

bpCoverage is the base-level coverage list which can then be used for plotting.

## Base-level coverage matrices for each of the regions
## Useful for plotting
lapply(regionMat$chr21$bpCoverage[1:2], head, n = 2)
## $`1`
##   HSB113 HSB123 HSB126 HSB130 HSB135 HSB136 HSB145 HSB153 HSB159 HSB178 HSB92
## 1  93.20   3.32  28.22   5.62 185.17  98.34   5.88  16.71   3.52  15.71 47.40
## 2 124.76   7.25  63.68  11.32 374.85 199.28  10.39  30.53   5.83  29.35 65.04
##   HSB97
## 1 36.54
## 2 51.42
## 
## $`2`
##     HSB113 HSB123 HSB126 HSB130 HSB135 HSB136 HSB145 HSB153 HSB159 HSB178 HSB92
## 566  45.59   7.94  15.92  34.75 141.61 104.21  19.87  38.61   4.97   23.2 13.95
## 567  45.59   7.94  15.92  35.15 141.64 104.30  19.87  38.65   4.97   23.2 13.95
##     HSB97
## 566 22.21
## 567 22.21
## Check dimensions. First region is 565 long, second one is 138 bp long.
## The columns match the number of samples (12 in this case).
lapply(regionMat$chr21$bpCoverage[1:2], dim)
## $`1`
## [1] 565  12
## 
## $`2`
## [1] 138  12

The end result of the coverage matrix is shown below. Note that the coverage has been adjusted for read length. Because reads might not fully align inside a given region, the numbers are generally not integers but can be rounded if needed.

## Dimensions of the coverage matrix
dim(regionMat$chr21$coverageMatrix)
## [1] 45 12
## Coverage for each region. This matrix can then be used with limma or other pkgs
head(regionMat$chr21$coverageMatrix)
##         HSB113     HSB123      HSB126      HSB130      HSB135      HSB136
## 1 3653.1093346 277.072106 1397.068687 1106.722895 8987.460401 5570.221054
## 2  333.3740816  99.987237  463.909476  267.354342 1198.713552 1162.313418
## 3   35.3828948  20.153553   30.725394   23.483947   16.786842   17.168947
## 4   42.3398681  29.931579   41.094474   24.724736   32.634080   19.309606
## 5   77.7402631 168.939342  115.059342  171.861974  180.638684   93.503158
## 6    0.7988158   1.770263    1.473421    2.231053    1.697368    1.007895
##        HSB145       HSB153     HSB159      HSB178       HSB92        HSB97
## 1 1330.158818 1461.2986829 297.939342 1407.288552 1168.519079 1325.9622371
## 2  257.114210  313.8513139  67.940131  193.695657  127.543553  200.7834228
## 3   22.895921   52.8756585  28.145395   33.127368   23.758816   20.4623685
## 4   33.802632   51.6146040  31.244343   33.576974   29.546183   28.2011836
## 5   90.950526   36.3046051  78.069605   97.151316  100.085790   35.5428946
## 6    1.171316    0.4221053   1.000132    1.139079    1.136447    0.3956579

We can then use the coverage matrix and packages such as limma, DESeq2 or edgeR to identify which ERs are differentially expressed. Here we’ll use DESeq2 for which we need some phenotype data.

## Get pheno table
pheno <- subset(brainspanPheno, structure_acronym == "AMY")

Now we can identify the DERs using a rounded version of the coverage matrix.

## Identify which ERs are differentially expressed, that is, find the DERs
library("DESeq2")
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
## 
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
## 
##     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
##     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
##     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
##     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
##     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
##     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
##     colWeightedMeans, colWeightedMedians, colWeightedSds,
##     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
##     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
##     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
##     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
##     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
##     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
##     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
##     rowWeightedSds, rowWeightedVars
## 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")'.
## 
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
## 
##     rowMedians
## The following objects are masked from 'package:matrixStats':
## 
##     anyMissing, rowMedians
## Round matrix
counts <- round(regionMat$chr21$coverageMatrix)

## Round matrix and specify design
dse <- DESeqDataSetFromMatrix(counts, pheno, ~ group + gender)
## converting counts to integer mode
## Perform DE analysis
dse <- DESeq(dse, test = "LRT", reduced = ~gender, fitType = "local")
## estimating size factors
## estimating dispersions
## gene-wise dispersion estimates
## mean-dispersion relationship
## final dispersion estimates
## fitting model and testing
## Extract results
mcols(regionMat$chr21$regions) <- c(
    mcols(regionMat$chr21$regions),
    results(dse)
)

## Save info in an object with a shorter name
ers <- regionMat$chr21$regions
ers
## GRanges object with 45 ranges and 12 metadata columns:
##      seqnames            ranges strand |     value      area indexStart
##         <Rle>         <IRanges>  <Rle> | <numeric> <numeric>  <integer>
##    1    chr21   9827018-9827582      * |  313.6717 177224.53          1
##    2    chr21 15457301-15457438      * |  215.0846  29681.68        566
##    3    chr21 20230140-20230192      * |   38.8325   2058.12        704
##    4    chr21 20230445-20230505      * |   41.3245   2520.80        757
##    5    chr21 27253318-27253543      * |   34.9131   7890.37        818
##   ..      ...               ...    ... .       ...       ...        ...
##   41    chr21 33039644-33039688      * |   34.4705 1551.1742       2180
##   42    chr21 33040784-33040798      * |   32.1342  482.0133       2225
##   43    chr21          33040890      * |   30.0925   30.0925       2240
##   44    chr21 33040900-33040901      * |   30.1208   60.2417       2241
##   45    chr21 48019401-48019414      * |   31.1489  436.0850       2243
##       indexEnd cluster clusterL  baseMean log2FoldChange     lfcSE      stat
##      <integer>   <Rle>    <Rle> <numeric>      <numeric> <numeric> <numeric>
##    1       565       1      565 2846.2872     -1.6903182  0.831959  0.215262
##    2       703       2      138  451.5196     -1.1640426  0.757490  0.871126
##    3       756       3      366   29.5781      0.0461488  0.458097  3.132082
##    4       817       3      366   36.0603     -0.1866200  0.390920  2.225708
##    5      1043       4      765  101.6468     -0.1387377  0.320166  3.957987
##   ..       ...     ...      ...       ...            ...       ...       ...
##   41      2224      17       45 20.782035      -0.642056  0.427661 0.6047814
##   42      2239      18      118  6.410542      -0.634321  0.512262 0.5454039
##   43      2240      18      118  0.129717      -0.859549  3.116540 0.0206273
##   44      2242      18      118  0.702291      -0.628285  2.247378 0.5825105
##   45      2256      19       14  5.293293      -1.694563  1.252290 5.7895910
##         pvalue      padj
##      <numeric> <numeric>
##    1 0.6426743  0.997155
##    2 0.3506436  0.997155
##    3 0.0767657  0.863614
##    4 0.1357305  0.997155
##    5 0.0466495  0.862040
##   ..       ...       ...
##   41 0.4367595  0.997155
##   42 0.4602018  0.997155
##   43 0.8857989  0.997155
##   44 0.4453299  0.997155
##   45 0.0161213  0.725460
##   -------
##   seqinfo: 1 sequence from an unspecified genome

We can then compare the DERs against known annotation to see which DERs overlap known exons, introns, or intergenic regions. A way to visualize this information is via a Venn diagram which we can create using vennRegions() from the derfinderPlot package as shown in Figure 1.

## Find overlaps between regions and summarized genomic annotation
annoRegs <- annotateRegions(ers, genomicState$fullGenome, verbose = FALSE)

library("derfinderPlot")
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
venn <- vennRegions(annoRegs,
    counts.col = "blue",
    main = "Venn diagram using TxDb.Hsapiens.UCSC.hg19.knownGene annotation"
)
Venn diagram showing ERs by annotation class.

Figure 1: Venn diagram showing ERs by annotation class

We can also identify the nearest annotated feature. In this case, we’ll look for the nearest known gene from the UCSC hg19 annotation.

## Load database of interest
library("TxDb.Hsapiens.UCSC.hg19.knownGene")
## Loading required package: GenomicFeatures
## Loading required package: AnnotationDbi
txdb <- keepSeqlevels(TxDb.Hsapiens.UCSC.hg19.knownGene, "chr21")

## Find nearest feature
library("bumphunter")
## Loading required package: foreach
## Loading required package: iterators
## Loading required package: locfit
## locfit 1.5-9.4    2020-03-24
genes <- annotateTranscripts(txdb)
## No annotationPackage supplied. Trying org.Hs.eg.db.
## Loading required package: org.Hs.eg.db
## 
## Getting TSS and TSE.
## Getting CSS and CSE.
## Getting exons.
## Annotating genes.
annotation <- matchGenes(ers, subject = genes)

## Restore seqlevels
txdb <- restoreSeqlevels(txdb)

## View annotation results
head(annotation)
##       name
## 1     <NA>
## 2     LIPI
## 3 TMPRSS15
## 4 TMPRSS15
## 5      APP
## 6      APP
##                                                                                                                                                                                                                                                                                              annotation
## 1                                                                                                                                                                                                                                                                                                  <NA>
## 2                                                                                   NM_001302998 NM_001302999 NM_001303000 NM_001303001 NM_001379565 NM_001379566 NM_145317 NM_198996 NP_001289927 NP_001289928 NP_001289929 NP_001289930 NP_001366494 NP_001366495 NP_945347 XM_006723965 XP_006724028
## 3                                                                                                                       NM_002772 NP_002763 XM_011529654 XM_011529655 XM_011529656 XM_011529657 XM_011529658 XM_011529659 XP_011527956 XP_011527957 XP_011527958 XP_011527959 XP_011527960 XP_011527961
## 4                                                                                                                       NM_002772 NP_002763 XM_011529654 XM_011529655 XM_011529656 XM_011529657 XM_011529658 XM_011529659 XP_011527956 XP_011527957 XP_011527958 XP_011527959 XP_011527960 XP_011527961
## 5 NM_000484 NM_001136016 NM_001136129 NM_001136130 NM_001136131 NM_001204301 NM_001204302 NM_001204303 NM_001385253 NM_201413 NM_201414 NP_000475 NP_001129488 NP_001129601 NP_001129602 NP_001129603 NP_001191230 NP_001191231 NP_001191232 NP_001372182 NP_958816 NP_958817 XM_024452075 XP_024307843
## 6 NM_000484 NM_001136016 NM_001136129 NM_001136130 NM_001136131 NM_001204301 NM_001204302 NM_001204303 NM_001385253 NM_201413 NM_201414 NP_000475 NP_001129488 NP_001129601 NP_001129602 NP_001129603 NP_001191230 NP_001191231 NP_001191232 NP_001372182 NP_958816 NP_958817 XM_024452075 XP_024307843
##   description      region distance   subregion insideDistance exonnumber nexons
## 1 close to 3' close to 3'      815        <NA>             NA         NA      1
## 2  downstream  downstream   125774        <NA>             NA         NA      8
## 3    upstream    upstream   454170        <NA>             NA         NA     25
## 4    upstream    upstream   454475        <NA>             NA         NA     25
## 5 inside exon      inside   289903 inside exon              0         16     16
## 6 inside exon      inside   289899 inside exon              0         16     16
##     UTR strand  geneL codingL    Geneid subjectHits
## 1  <NA>      +     60      NA 100500815          15
## 2  <NA>      - 102077  101886    149998         149
## 3  <NA>      - 134537  133653      5651         579
## 4  <NA>      - 134537  133653      5651         579
## 5 3'UTR      - 290585  230434       351         354
## 6 3'UTR      - 290585  230434       351         354
## You can use derfinderPlot::plotOverview() to visualize this information

We can check the base-level coverage information for some of our DERs. In this example we do so for the first 5 ERs (Figures ??, ??, ??, ??, ??).

## Extract the region coverage
regionCov <- regionMat$chr21$bpCoverage
plotRegionCoverage(
    regions = ers, regionCoverage = regionCov,
    groupInfo = pheno$group, nearestAnnotation = annotation,
    annotatedRegions = annoRegs, whichRegions = seq_len(5), txdb = txdb,
    scalefac = 1, ask = FALSE, verbose = FALSE
)