derfinder 1.24.2
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()
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.
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.
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)'.
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
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.
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"
)
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
)