consensusDE aims to make differential expression (DE) analysis, with reporting of significance scores from multiple methods, with and without removal of unwanted variation (RUV) easy. It implements voom/limma, DESeq2 and edgeR and reports differential expression results seperately for each algorithm, as well as merging the results into a single table for determining consensus. The results of the merged table, are ordered by the summed ranks of p-values for each algorithm, the intersect at minimum p-value thresholds accross all methods is provided as the p_intersect, in addition to a number of statistics (see below).
consensusDE is simplified into two functions:
buildSummarized()
generate a summarized experiment that counts reads mapped (from bam files or htseq count files) against a transcriptomemulti_de_pairs()
perform DE analysis (all possible pairwise comparisons)Below the core functionality of consensusDE as well as how to plot results using the diag_plots
function.
Begin by first installing and loading the consensusDE
library. To illustrate functionality of consensusDE
, we will utilise RNA-seq data from the airway
and annotation libraries as follows. Begin by installing and attaching data from these libraries as follows:
A summarized experiment is an object that stores all relevant information for performing differential expression analysis. buildSummarized()
allows users to build a summarized experiment object by simply providing 1) a table of bam/htseq files (more below on format), 2) the directory of where to locate files and 3) a transcript database to map the reads to (either a gtf file or txdb). Below we will use bam files (from GenomicAlignments) as an example for creating a summarized experiment:
# build a design table that lists the files and their grouping
file_list <- list.files(system.file("extdata", package="GenomicAlignments"),
recursive = TRUE,
pattern = "*bam$",
full = TRUE)
# Prepare a sample table to be used with buildSummarized()
# must be comprised of a minimum of two columns, named "file" and "group",
# with one additional column: "pairs" if the data is paired
sample_table <- data.frame("file" = basename(file_list),
"group" = c("treat", "untreat"))
# extract the path to the bam directory - where to search for files listed in "sample_table"
bam_dir <- as.character(gsub(basename(file_list)[1], "", file_list[1]))
The minimum information is now ready to build a summarized experiment:
# NB. force_build = TRUE, is set to allow the Summarized Experiment to be built.
# This will report a Warning message that less than two replicates are present
# in the sample_table.
summarized_dm3 <- buildSummarized(sample_table = sample_table,
bam_dir = bam_dir,
tx_db = TxDb.Dmelanogaster.UCSC.dm3.ensGene,
read_format = "paired",
force_build = TRUE)
## strand_mode is defined as 0 (unstranded). This is appropriate for
## unstranded protocols, or if you wish to ignore strandedness when
## counting reads. See ?strandMode in GenomicAlignments for more
## information.
This will output a summarized object that has mapped the reads for the bam files that are listed in sample_table
, located in bam_dir
, against the transcript database provided: TxDb.Dmelanogaster.UCSC.dm3.ensGene
. Bam file format, whether “paired” or “single” end (the type of sequencing technology used) must be specified using the read_format
parameter.gtf formatted transcript databases can also be used instead of a txdb, by providing the full path to the gtf file using the gtf
parameter. To save a summarized experiment externally, for future use, specify a path to save the summarized experiment using output_log
.
strand_mode
is used to define how the stranded library prep protocol treated the strand. For paired data, this is used to indicate how the strand is inferred from the first and last fragments in the paired reads. If the protocol was unstranded or stranding should be ignored, strand_mode = 0
. ConsensusDE
calls strand_mode = 0
by default. If the protocol was stranded and the strand of the read is the strand of the first fragment (or read in single ended libraries), strand_mode = 1
. If the protocol was stranded and the strand of the read is the strand of the second fragment (or reverse of the read in single ended libraries) , strand_mode = 2
. For more information, see ?strandMode in the Genomic Alignments package.
To see details of all parameters see ?buildSummarized
.
Overview of the summarized experiment:
## class: RangedSummarizedExperiment
## dim: 15682 2
## metadata(2): gene_coords sample_table
## assays(1): counts
## rownames(15682): FBgn0000003 FBgn0000008 ... FBgn0264726 FBgn0264727
## rowData names(0):
## colnames(2): sm_treated1.bam sm_untreated1.bam
## colData names(2): file group
buildSummarized()
also allows users to filter out low read counts. This can be done when building the summarized experiment, or re-running with the summarized experiment output using buildSummarized()
. See “Performing Differential Expresssion” below with filter example.
Sometimes it will be convenient to first build a txdb
object and then pass this txdb
object to buildSummarized using the tx_db parameter. This can be done as follows:
txdb <- makeTxDbFromGFF("/path/to/my.gtf", format="gtf", circ_seqs=character())
For differential expression (DE) analysis we will use the airway
RNA-seq data for demonstration. See ?airway
for more details about this experiment. NOTE: the summarized meta-data must include the columns “group” and “file” to build the correct models. For illustration, we sample 1000 genes from this dataset.
# for compatability for DE analysis, add "group" and "file" columns
colData(airway)$group <- colData(airway)$dex
colData(airway)$file <- rownames(colData(airway))
# filter low count data
airway_filter <- buildSummarized(summarized = airway,
filter = TRUE)
# for illustration, we only use sa random sample of 1000 transcripts
set.seed(1234)
airway_filter <- sample(airway_filter, 1000)
# call multi_de_pairs()
all_pairs_airway <- multi_de_pairs(summarized = airway_filter,
paired = "unpaired",
ruv_correct = FALSE)
Running multi_de_pairs()
will perform DE analysis on all possible pairs of “groups” and save these results as a simple list of “merged” - being the merged results of “deseq”, “voom” and “edger” into one table, as well as the latter three as objects independently. The data frame is sorted by the rank_sum
. The following columns are included:
ID
- IdentifierAveExpr
- Average Expression (average of edgeR, DESeq2 and voom)LogFC
- Log2 Fold-Change, also known as a log-ratio (average of edgeR, DESeq2 and voom)LogFC_sd
- Log2 Fold-Change standard deviation of LogFC (average)edger_adj_p
- EdgeR p-value adjusted for multiple hypothesesdeseq_adj_p
- DESeq2 p-value adjusted for multiple hypothesesvoom_adj_p
- Limma/voom p-value adjusted for multiple hypothesesedger_rank
- rank of the p-value obtained by EdgeRdeseq_rank
- rank of the p-value obtained by DESeq2voom_rank
- rank of the p-value obtained by Limma/voomrank_sum
- sum of the ranks from edger_rank, voom_rank, rank_sump_intersect
- the largest p-value observed from all methods tested.
p_union
- the smallest p-value observed from all methods tested.
## [1] "untrt-trt"
# [1] "untrt-trt"
# to access data of a particular comparison
head(all_pairs_airway$merged[["untrt-trt"]])
## ID AveExpr LogFC LogFC_sd edger_adj_p deseq_adj_p
## 1 ENSG00000120129 11.38754 -2.811094 0.1118097 4.444781e-38 6.013628e-44
## 2 ENSG00000116584 11.02552 1.153092 0.1095189 1.508445e-16 9.643002e-42
## 3 ENSG00000139289 10.83881 1.131120 0.1098113 1.166207e-13 2.233808e-24
## 4 ENSG00000077684 10.46766 -1.066130 0.1135320 4.188089e-11 1.290408e-26
## 5 ENSG00000103196 11.26943 -2.540661 0.1367032 1.906080e-16 1.665175e-21
## 6 ENSG00000211445 13.09005 -3.598107 0.1163315 1.246261e-14 2.170616e-19
## voom_adj_p edger_rank deseq_rank voom_rank rank_sum p_intersect
## 1 4.189235e-05 1 1 1.0 3.0 4.189235e-05
## 2 2.107158e-04 2 2 2.0 6.0 2.107158e-04
## 3 5.341853e-04 5 5 3.5 13.5 5.341853e-04
## 4 5.341853e-04 8 3 3.5 14.5 5.341853e-04
## 5 9.963537e-04 3 6 6.5 15.5 9.963537e-04
## 6 9.963537e-04 4 9 6.5 19.5 9.963537e-04
## p_union
## 1 6.013628e-44
## 2 9.643002e-42
## 3 2.233808e-24
## 4 1.290408e-26
## 5 1.665175e-21
## 6 2.170616e-19
It is recommended to annotate with a GTF file byt providing the full path of a gtf file to the gtf_annotate parameter, in combination with a tx_db. If no tx_db is provided and the gtf path is provided, only gene symbol annotations will be performed.
Currently only ENSEMBL annotations are supported with the tx_db option.
It is often useful to add additional annotated information to the output tables. This can be achieved by providing a database for annotations via ensembl_annotate
. Annotations needs to be a Genome Wide Annotation object, e.g. org.Mm.eg.db
for mouse or org.Hs.eg.db
for human from BioConductor. For example, to install the database for the mouse annotation, go to http://bioconductor.org/packages/org.Mm.eg.db and follow the instructions. Ensure that after installing the database package that the library is loaded using library(org.Mm.eg.db)
. When running, “‘select()’ returned 1:many mapping between keys and columns” will appear on the command line. This is the result of multiple mapped transcript ID to Annotations. Only the first annotation is reported. See ?multi_de_pairs
for additional documentation.
An example of annotating the above filtered airway data is provided below:
# first ensure annotation database in installed
#library(org.Hs.eg.db)
#library(EnsDb.Hsapiens.v86)
# Preloaded summarized file did not contain meta-data of the tx_db. This is important if you want to extract chromosome coordinates. This can be easily updated by rerunning buildSummarized with the tx_db of choice.
airway_filter <- buildSummarized(summarized = airway_filter,
tx_db = EnsDb.Hsapiens.v86,
filter = FALSE)
## Warning in buildSummarized(summarized = airway_filter, tx_db = EnsDb.Hsapiens.v86, : No output directory provided. The se file and sample_table will not
## be saved
# call multi_de_pairs(),
# set ensembl_annotate argument to org.Hs.eg.db
all_pairs_airway <- multi_de_pairs(summarized = airway_filter,
paired = "unpaired",
ruv_correct = FALSE,
ensembl_annotate = org.Hs.eg.db)
# to access data of a particular comparison
head(all_pairs_airway$merged[["untrt-trt"]])
## ID AveExpr LogFC LogFC_sd edger_adj_p deseq_adj_p
## 275 ENSG00000120129 11.38754 -2.811094 0.1118097 4.444781e-38 6.013628e-44
## 245 ENSG00000116584 11.02552 1.153092 0.1095189 1.508445e-16 9.643002e-42
## 395 ENSG00000139289 10.83881 1.131120 0.1098113 1.166207e-13 2.233808e-24
## 90 ENSG00000077684 10.46766 -1.066130 0.1135320 4.188089e-11 1.290408e-26
## 149 ENSG00000103196 11.26943 -2.540661 0.1367032 1.906080e-16 1.665175e-21
## 807 ENSG00000211445 13.09005 -3.598107 0.1163315 1.246261e-14 2.170616e-19
## voom_adj_p edger_rank deseq_rank voom_rank rank_sum p_intersect
## 275 4.189235e-05 1 1 1.0 3.0 4.189235e-05
## 245 2.107158e-04 2 2 2.0 6.0 2.107158e-04
## 395 5.341853e-04 5 5 3.5 13.5 5.341853e-04
## 90 5.341853e-04 8 3 3.5 14.5 5.341853e-04
## 149 9.963537e-04 3 6 6.5 15.5 9.963537e-04
## 807 9.963537e-04 4 9 6.5 19.5 9.963537e-04
## p_union genename
## 275 6.013628e-44 dual specificity phosphatase 1
## 245 9.643002e-42 Rho/Rac guanine nucleotide exchange factor 2
## 395 2.233808e-24 pleckstrin homology like domain family A member 1
## 90 1.290408e-26 jade family PHD finger 1
## 149 1.665175e-21 cysteine rich secretory protein LCCL domain containing 2
## 807 2.170616e-19 glutathione peroxidase 3
## symbol kegg coords strand width
## 275 DUSP1 04010 chr5:172768090-172771195 - 3106
## 245 ARHGEF2 05130 chr1:155946851-156007070 - 60220
## 395 PHLDA1 <NA> chr12:76025447-76033932 - 8486
## 90 JADE1 <NA> chr4:128809623-128875224 + 65602
## 149 CRISPLD2 <NA> chr16:84819984-84920768 + 100785
## 807 GPX3 00480 chr5:151020438-151028993 + 8556
The following additional columns will now be present:
genename
- extend gene names (e.g. alpha-L-fucosidase 2)symbol
- gene symbol (e.g. FUCA2)kegg
- kegg pathway identifier (e.g. 00511)If metadata for the transcript database used to build the summarized experiment was included, the following annotations will also be included:
coords
- chromosomal coordinates (e.g. chr6:143494811-143511690)strand
- strand transcript is on (i.e. + or -)width
- transcript width in base pairs (bp) (transcript start to end) (e.g. 16880 bp)multi_de_pairs
provides options to automatically write all results to output directories when a full path is provided. Which results are output depends on which directories are provided. Full paths provided to the parameters of output_voom
, output_edger
, output_deseq
and output_combined
will output Voom, EdgeR, DEseq and the merged results to the directories provided, respectively.
consensusDE also provides the option to remove batch effects through RUVseq functionality. consensusDE currently implements RUVr which models a first pass generalised linear model (GLM) using EdgeR and obtaining residuals for incorporation into the SummarizedExperiment object for inclusion in the models for DE analysis. The following example, uses RUV to identify these residuals. To view the residuals in the model see the resisuals section below in the plotting functions. Note, that if ruv_correct = TRUE
and a path to a plot_dir
is provided, diagnostic plots before and after RUV correction will be produced. The residuals can also be accessed in the summarizedExperiment as below. These are present in the “W_1” column. At present only one factor of variation is determined.
# call multi_de_pairs()
all_pairs_airway_ruv <- multi_de_pairs(summarized = airway_filter,
paired = "unpaired",
ruv_correct = TRUE)
# access the summarized experiment (now including the residuals under the "W_1" column)
all_pairs_airway_ruv$summarized@phenoData@data
## SampleName cell dex albut Run avgLength Experiment
## SRR1039508 GSM1275862 N61311 untrt untrt SRR1039508 126 SRX384345
## SRR1039509 GSM1275863 N61311 trt untrt SRR1039509 126 SRX384346
## SRR1039512 GSM1275866 N052611 untrt untrt SRR1039512 126 SRX384349
## SRR1039513 GSM1275867 N052611 trt untrt SRR1039513 87 SRX384350
## SRR1039516 GSM1275870 N080611 untrt untrt SRR1039516 120 SRX384353
## SRR1039517 GSM1275871 N080611 trt untrt SRR1039517 126 SRX384354
## SRR1039520 GSM1275874 N061011 untrt untrt SRR1039520 101 SRX384357
## SRR1039521 GSM1275875 N061011 trt untrt SRR1039521 98 SRX384358
## Sample BioSample group file W_1
## SRR1039508 SRS508568 SAMN02422669 untrt SRR1039508 -0.08312388
## SRR1039509 SRS508567 SAMN02422675 trt SRR1039509 0.01853551
## SRR1039512 SRS508571 SAMN02422678 untrt SRR1039512 -0.16047843
## SRR1039513 SRS508572 SAMN02422670 trt SRR1039513 -0.26699805
## SRR1039516 SRS508575 SAMN02422682 untrt SRR1039516 0.55688461
## SRR1039517 SRS508576 SAMN02422673 trt SRR1039517 0.60106294
## SRR1039520 SRS508579 SAMN02422683 untrt SRR1039520 -0.29682660
## SRR1039521 SRS508580 SAMN02422677 trt SRR1039521 -0.36905611
# view the results, now with RUV correction applied
head(all_pairs_airway_ruv$merged[["untrt-trt"]])
## ID AveExpr LogFC LogFC_sd edger_adj_p deseq_adj_p
## 1 ENSG00000120129 11.38755 -2.795845 0.1043348 2.563603e-52 2.374465e-75
## 2 ENSG00000103196 11.26944 -2.498026 0.1089966 1.176625e-20 7.070363e-35
## 3 ENSG00000116584 11.02551 1.153351 0.1094772 7.233342e-19 1.343923e-30
## 4 ENSG00000211445 13.09006 -3.531512 0.0952666 1.176625e-20 1.113829e-31
## 5 ENSG00000077684 10.46767 -1.061481 0.1096889 2.006530e-12 7.017736e-35
## 6 ENSG00000139289 10.83879 1.130942 0.1093725 2.334717e-15 3.319270e-20
## voom_adj_p edger_rank deseq_rank voom_rank rank_sum p_intersect
## 1 2.077071e-05 1.0 1 1.0 3.0 2.077071e-05
## 2 3.970362e-04 2.5 3 4.5 10.0 3.970362e-04
## 3 2.121089e-04 4.0 5 2.0 11.0 2.121089e-04
## 4 4.617250e-04 2.5 4 8.0 14.5 4.617250e-04
## 5 3.970362e-04 11.0 2 4.5 17.5 3.970362e-04
## 6 3.970362e-04 6.0 8 4.5 18.5 3.970362e-04
## p_union
## 1 2.374465e-75
## 2 7.070363e-35
## 3 1.343923e-30
## 4 1.113829e-31
## 5 7.017736e-35
## 6 3.319270e-20
multi_de_pairs
supports DE with paired samples. Paired samples may include, for example, the same patient observed before and after a treatment. For demonstration purposes, we assume that each untreated and treated sample is a pair.
NB. paired analysis with more than two groups is not currently supported. If there are more than two groups, consider testing each of the groups and their pairs seperately, or see the edgeR, limma/voom or DESeq2 vignettes for establishing a multi-variate model with blocking factors.
First we will update the summarized experiment object to include a “pairs” column and set paired = "paired"
in multi_de_pairs
.
# add "pairs" column to airway_filter summarized object
colData(airway_filter)$pairs <- as.factor(c("pair1", "pair1", "pair2", "pair2", "pair3", "pair3", "pair4", "pair4"))
# run multi_de_pairs in "paired" mode
all_pairs_airway_paired <- multi_de_pairs(summarized = airway_filter,
paired = "paired",
ruv_correct = TRUE)
head(all_pairs_airway_paired$merged[["untrt-trt"]])
## ID AveExpr LogFC LogFC_sd edger_adj_p deseq_adj_p
## 1 ENSG00000211445 13.090146 -3.591562 0.11143193 2.013996e-191 2.107187e-261
## 2 ENSG00000120129 11.387754 -2.810740 0.11097119 9.093204e-136 7.801298e-140
## 3 ENSG00000103196 11.269573 -2.492442 0.11071604 2.411102e-113 6.650710e-116
## 4 ENSG00000253368 9.790702 -1.849633 0.11326685 3.576180e-67 9.160743e-57
## 5 ENSG00000137672 9.129518 -1.807741 0.11314394 2.082487e-56 4.560821e-44
## 6 ENSG00000180914 8.691937 -1.747506 0.09151028 2.133250e-41 2.155192e-29
## voom_adj_p edger_rank deseq_rank voom_rank rank_sum p_intersect
## 1 1.191186e-11 1 1.0 1 3.0 1.191186e-11
## 2 1.512098e-10 2 2.0 2 6.0 1.512098e-10
## 3 3.620473e-10 3 3.0 3 9.0 3.620473e-10
## 4 1.047378e-08 4 4.0 4 12.0 1.047378e-08
## 5 2.178982e-08 5 5.0 5 15.0 2.178982e-08
## 6 7.777989e-08 8 10.5 7 25.5 7.777989e-08
## p_union
## 1 2.107187e-261
## 2 7.801298e-140
## 3 6.650710e-116
## 4 3.576180e-67
## 5 2.082487e-56
## 6 2.133250e-41
The design matrix can be retrieved as follows (from e.g. the voom model fit)
## Intercept W_1 pair2 pair3 pair4 untrt
## SRR1039508 1 -0.25062331 0 0 0 1
## SRR1039509 1 0.24948639 0 0 0 0
## SRR1039512 1 -0.26981902 1 0 0 1
## SRR1039513 1 0.27385504 1 0 0 0
## SRR1039516 1 0.59944990 0 1 0 1
## SRR1039517 1 -0.59899109 0 1 0 0
## SRR1039520 1 -0.06879627 0 0 1 1
## SRR1039521 1 0.06543835 0 0 1 0
## attr(,"assign")
## [1] 0 1 2 2 2 3
## attr(,"contrasts")
## attr(,"contrasts")$pairs
## [1] "contr.treatment"
##
## attr(,"contrasts")$group
## [1] "contr.treatment"
consensusDE currently implements two main normalisation approaches in multi_de_pairs()
. These are specified with the norm_method
parameter, where options are: EDASeq
or all_defaults
. As per the parameter description, when all_defaults
is selected, this will use default normalisation methods for DE, EDASeq for QC (with control via EDASeq_method
), and edgeR “upperquantile” for determining RUV residuals (as per RUVSeq vignette). However, when EDASeq
is selected, this will use EDASeq normalisation and the specified EDASeq_method
throughout, for RUV, edgeR, DESeq2 and voom/limma. Using the EDASeq
allows for a standard normalisation approach to be used throughout, whereas all_defaults
, allows for variation of normalisation approach to also be modelled into the final merged results table.
When performing DE analysis, a series of plots (currently 10) can be generated and saved as .pdf files in a plot directory provided to multi_de_pairs()
with the parameter: plot_dir = "/path/to/save/pdfs/
. See ?multi_de_pairs
for description.
In addition, each of the 10 plots can be plotted individually using the diag_plots
function. See ?diag_plots
for description, which provides wrappers for 10 different plots. Next we will plot each of these using the example data.
Plot the number of reads that mapped to the transcriptome of each sample. The sample numbers on the x-axis correspond to the sample row number in the summarizedExperiment built, accessible using colData(airway)
. Samples are coloured by their “group”.
Residuals for the RUV model can be plotted as follows: