Package {diffwrap}

diffwrap: Differential Expression Analysis of RNA-Seq Data

Description

Functions for differential expression analysis of read counts from messenger RNA (mRNA) sequencing (RNA-Seq) data or micro RNA (miRNA) expression values generated by the Comprehensive Analysis Pipeline for microRNA Sequencing (CAP-miRSeq) 'expression_reports.sh' script. The workflow follows the 'edgeR'-'limma' expression data analysis pipeline providing options for different approaches, such as "pure" 'edgeR', voom or paired samples. The functions in the package generate text files with differential expression lists, optionally annotated with information from 'biomart', expression summary plots as well as several quality control (QC) plots. The main function, diffExpr(), is a convenience wrapper performing all steps automatically based on sensible defaults. Methods are described in Robinson, McCarthy and Smyth (2010) doi:10.1093/bioinformatics/btp616, Ritchie et al. (2015) doi:10.1093/nar/gkv007, Law et al. (2014) doi:10.1186/gb-2014-15-2-r29 and Sun et al. (2014) doi:10.1186/1471-2164-15-423.

Author(s)

Maintainer: Vidal Fey vidal.fey@gmail.com

Authors:

• Vidal Fey vidal.fey@gmail.com

• Meeri Pekkarinen meeri.pekkarinen@gmail.com

• Reija Hieta reija.hieta@geneviatechnologies.com

• Bogdan Iancu bogdan.iancu@abo.fi

• Adrien Janssens Adrien.RFM.Janssens@proton.me

Pipe operator

Description

See magrittr::%>% for details.

Usage

lhs %>% rhs

Arguments

Value

The result of calling 'rhs(lhs)'.

Function to create a correlogram pheatmap, i.e., a plot to check randomness in the data set.

Description

Function to create a correlogram pheatmap, i.e., a plot to check randomness in the data set.

Usage

correlogram_pheatmap(
  expr.mat,
  clinical.mat,
  scale.fl = "none",
  legend.fl = TRUE,
  row.clust = TRUE,
  col.clust = TRUE,
  signif.stars.fl = FALSE,
  cell.size = 8,
  font.size = 11,
  color.blind.pal = "PuOr",
  color.extremes = c("#3182BD", "#E6550D"),
  anno.color = NULL,
  main.correl = "Correlogram",
  sample.correl = FALSE
)

Arguments

Details

Calculates Pearson's rank correlation coefficients for all possible pairs of the input matrix and plots the resulting correlation matrix, either with significance stars and coloured according to correlation or without stars and coloured according to clinical annotation.

Value

Returns a pheatmap plot object used in the diffr_pheatmap() function.

Author(s)

Bogdan Iancu - Genevia Technologies Oy

Main wrapper for executing the entire pipeline from reading in expression data such as count files to producing text files and graphs

Description

diffExpr is a convenience wrapper performing all steps automatically. Most sub-functions are exported and can be called by the user, as well, if desired. These functions may be applicable to different kinds of data/input, rely, however, on the conventions set for this package.

Usage

diffExpr(
  expr.dat,
  samp.info,
  control,
  design = NULL,
  samples = NULL,
  sample.plot.names = NULL,
  groups = NULL,
  pairs = NULL,
  block = FALSE,
  contrasts = NULL,
  out.dir = NULL,
  analysis.name = NULL,
  filter.strict = TRUE,
  filter.min.samp = NULL,
  fit.voom = FALSE,
  fit.voom.fun = "voomLmFit",
  fit.use.weights = FALSE,
  fit.norm.method = c("tmm", "quantile"),
  fit.disp = c("gene", "trend", "common"),
  fit.bayes.trend = FALSE,
  fit.bayes.robust = FALSE,
  fit.quasi.likelihood = NA,
  p.thr = 0.05,
  fdr.thr = 0.05,
  logfc.thr = 1,
  numlab = 25,
  point.lab = TRUE,
  de.plot.base.size = 16,
  biom.use = FALSE,
  biom.data.set = "hsapiens_gene_ensembl",
  biom.mart = "ensembl",
  biom.host = "https://www.ensembl.org",
  biom.timeout = 10,
  biom.filter = "ensembl_gene_id",
  biom.attributes = c("ensembl_gene_id", "hgnc_symbol", "description", "entrezgene_id"),
  biom.force.ensg = FALSE,
  biom.cache = NULL,
  biom.use.cache = FALSE,
  biom.sym.col = "hgnc_symbol",
  biom.rm.dups = FALSE,
  qc.top.n = 500,
  qc.gene.selection = "common",
  qc.pc = c(1, 2, 3),
  qc.type = c("both", "uncorrected", "pseudo-corrected"),
  qc.ellipse = TRUE,
  qc.ellipse.groups = NULL,
  qc.ellipse.legend = NA,
  qc.label.samples = TRUE,
  qc.point.size = 2,
  qc.label.size = 5,
  qc.circle = TRUE,
  qc.varname.size = 0,
  qc.var.axes = FALSE,
  hm.topn = 100,
  hm.p.thr = 0.05,
  hm.fdr.thr = 0.05,
  hm.logfc.thr = 1,
  hm.split.expr = FALSE,
  hm.pal.blind = "PuOr",
  hm.pal.n = 11,
  hm.pal.extremes = c("#3182BD", "#E6550D"),
  hm.pal.length = NULL,
  hm.anno.color = NULL,
  hm.anno.name = "Sample Class",
  enr.do = TRUE,
  enr.methods = c("clusterProfilerGO", "clusterProfilerKEGG", "gProfileR", "topGO"),
  enr.plot = FALSE,
  enr.plot.fdr.thr = fdr.thr,
  enr.plot.logfc.thr = logfc.thr,
  enr.plot.num.terms = 5,
  out.plots = TRUE,
  out.tables = TRUE,
  out.filtered.tables = TRUE,
  verbose = TRUE,
  log.file = NULL,
  dry.run = FALSE
)

Arguments

Details

For experimental designs involving comparisons within as well as between subjects inter-subject needs to be computed. In this case, the column specified in the 'pairs' argument must assign the subjects to the treatment/tissue/etc groups. For example, if we have two treatments the effects of which are to be observed in each two tissues, this design would apply. The 'pairs' factor is passed to the functions 'duplicateCorrelation()' and 'lmFit()'. The 'block' argument is used to specify whether the comparisons are to be made within AND between subjects or in the case of technical replicates, i.e., if the samples are not independent, in other words, correlated. That correlation is addressed by means of the 'duplicateCorrelation()' function in the limma package. If 'block' is set to TRUE, the (selected) 'voom' function is enforced. As of version 0.4, the 'edgeR::voomLmFit()' function is incorporated, which replaces 'voom()', 'lmFit()' and 'voomWithQualityWeights()'. voomLmFit()' ensures unbiased estimation of the residual variances and automates the estimation of sample weights and intrablock correlations. In edgeR, it is recommended to remove features without at least 1 read per million in n of the samples, where n is the size of the smallest group of replicates (determined from the 'groups' vector). The 'min.samp' argument is used to specify the number of samples in which a feature needs to be covered by at least one read per million. Quasi-likelihood pipeline: While the likelihood ratio test is a more obvious choice for inferences with GLMs, the QL F-test is preferred as it reflects the uncertainty in estimating the dispersion for each gene. It provides more robust and reliable error rate control when the number of replicates is small. The QL dispersion estimation and hypothesis testing is done by using the functions glmQLFit() and glmQLFTest().

Value

A list of all relevant objects generated in the course of the workflow:

• if 'voom=TRUE', voomed counts, otherwise DGEList object with TMM-normalisation factors

• if 'voom=TRUE', lmFit object, otherwise DGEList object with estimated dispersion

• if 'voom=TRUE', eBayes output of contrasted groups, otherwise glmFit object

• annotated topTable/topTags output

However, the function is first and foremost called for its side effects of generating results tables and plots.

Argument groups

The arguments are grouped by name prefix:

• no prefix: core experiment definition, plus the primary differential-expression thresholds (p.thr, fdr.thr, logfc.thr, numlab, point.lab);

• filter.*: count pre-filtering;

• fit.*: model fitting (voom / edgeR options);

• biom.*: 'biomart' gene annotation;

• qc.*: quality-control plots (MDS / PCA);

• hm.*: heatmaps, including their own significance thresholds and colour palette;

• enr.*: enrichment analysis;

• out.*, verbose, log.file, dry.run: output and run control.

Examples

out.dir <- file.path(tempdir(), "diffwrap_demo")
dir.create(out.dir, showWarnings = FALSE)
res <- diffExpr(expr.dat = diffwrap_counts,
                samp.info = diffwrap_samp_info,
                samples = "SampleName", groups = "Group",
                control = "control", analysis.name = "demo",
                out.dir = out.dir, enr.do = FALSE)
names(res$contrasts)

Function to generate a 3D scatterplot

Description

Function to generate a 3D scatterplot

Usage

diff_expr_3d_scatterplot(
  PCA,
  samp.name = NULL,
  groups,
  grp.nam = NULL,
  PC = c(1, 2, 3),
  main = NULL
)

Arguments

Value

No return value. Called for its side effect of drawing a three-dimensional PCA scatterplot on the current graphics device.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
pca <- diff_expr_PCA(edgeR::cpm(counts, log = TRUE), n = 100)
diff_expr_3d_scatterplot(pca, groups = groups)

Function to do PCA using 'stats::prcomp'

Description

Function to do PCA using 'stats::prcomp'

Usage

diff_expr_PCA(counts, n = 500, scale. = FALSE)

Arguments

Value

An object of class prcomp as returned by prcomp, holding the principal component decomposition of the (optionally variance-filtered) count matrix.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
pca <- diff_expr_PCA(edgeR::cpm(counts, log = TRUE), n = 100)
pca$sdev[1:3]

Function to generate a PCA biplot using 'ggbiplot.n', a version of 'ggbiplot' from https://github.com/vqv/ggbiplot.

Description

Function to generate a PCA biplot using 'ggbiplot.n', a version of 'ggbiplot' from https://github.com/vqv/ggbiplot.

Usage

diff_expr_PCA_ggbiplot(
  PCA,
  groups,
  grp.nam = NULL,
  ellipse = TRUE,
  circle = TRUE,
  varname.size = 0,
  var.axes = FALSE,
  main = NULL,
  fix.aspect = FALSE,
  tweak = FALSE,
  ...
)

Arguments

Value

A ggplot object containing the PCA biplot. The plot is returned rather than drawn, so it has to be printed to appear on a device.

See Also

[ggbiplot.n()]

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
pca <- diff_expr_PCA(edgeR::cpm(counts, log = TRUE), n = 100)
g <- diff_expr_PCA_ggbiplot(pca, groups = groups)
class(g)

Function to generate an ordinary two-dimensional PCA plot using 'ggplot2'

Description

Function to generate an ordinary two-dimensional PCA plot using 'ggplot2'

Usage

diff_expr_PCA_ggplot(
  PCA,
  samp.name = NULL,
  groups,
  grp.nam = NULL,
  PC = c(1, 2),
  main = NULL,
  ellipse = TRUE,
  ellipse.mapping.groups = NULL,
  ellipse.grp.nam = NULL,
  label.samples = TRUE,
  geom.point.size = 2,
  label.font.size = 5,
  plot.ellipse.legend = NA,
  do.plot = TRUE
)

Arguments

Value

A ggplot object containing the labelled PCA scatterplot.

See Also

[stat_ellipse()]

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
pca <- diff_expr_PCA(edgeR::cpm(counts, log = TRUE), n = 100)
g <- diff_expr_PCA_ggplot(pca, samp.name = NULL, groups = groups, do.plot = FALSE)
class(g)

Main wrapper function for QC plots

Description

Main wrapper function for QC plots

Usage

diff_expr_QC_plots(
  counts,
  samp.info,
  control,
  out.l,
  grp.nam = NULL,
  PC = c(1, 2, 3),
  sample.plot.names = NULL,
  ellipse = TRUE,
  ellipse.mapping.groups = NULL,
  ellipse.grp.nam = NULL,
  label.samples = TRUE,
  geom.point.size = 2,
  label.font.size = 5,
  plot.ellipse.legend = NA,
  circle = TRUE,
  varname.size = 0,
  var.axes = FALSE,
  pairs = NULL,
  pairs.name = NULL,
  gene.selection = "common",
  n = 500,
  type = NULL,
  analysis.name = NULL,
  out.dir
)

Arguments

Value

The input list out.l, with a QCplots element added (or extended) holding the generated quality control plots as ggplot objects, named after the plot type and the value of type.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
lcpm <- edgeR::cpm(counts, log = TRUE)
out.l <- diff_expr_QC_plots(counts = lcpm, samp.info = si, control = "control",
                            out.l = list(), grp.nam = "Group",
                            sample.plot.names = colnames(lcpm),
                            analysis.name = "demo", out.dir = tempdir())
names(out.l$QCplots)

Function to retrieve additional information from biomart

Description

Function to retrieve additional information from biomart

Usage

diff_expr_biomart(
  d3,
  biom.data.set = "hsapiens_gene_ensembl",
  biom.mart = "ensembl",
  host = "https://www.ensembl.org",
  biom.filter = "ensembl_gene_id",
  biom.attributes = c("ensembl_gene_id", "hgnc_symbol", "description"),
  biom.cache = NULL,
  use.cache = FALSE,
  sym.col = "hgnc_symbol",
  rm.dups = FALSE,
  force.ensg = FALSE,
  verbose = FALSE
)

Arguments

Value

A data.frame: the input table d3 merged with the annotation retrieved from biomart. The column holding gene symbols is renamed to gene_symbol.

Examples

# Needs network access to Ensembl BioMart. The query is wrapped in tryCatch() so that an
# unreachable or slow server reports the problem instead of failing the example.
if (requireNamespace("biomaRt", quietly = TRUE)) {
  d3 <- data.frame(ID = c("ENSG00000141510", "ENSG00000012048"))
  ann <- tryCatch(diff_expr_biomart(d3, biom.data.set = "hsapiens_gene_ensembl"),
                  error = function(e) {
                    message("BioMart not reachable: ", conditionMessage(e))
                    NULL
                  })
  if (!is.null(ann)) nrow(ann)
}

Function to generate dendrogram plots based on hierarchical clustering

Description

Function to generate dendrogram plots based on hierarchical clustering

Usage

diff_expr_dendro_plot(
  counts,
  groups,
  grp.nam = NULL,
  main = NULL,
  col.grps = FALSE
)

Arguments

Value

No return value. Called for its side effect of drawing a hierarchical clustering dendrogram on the current graphics device.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
diff_expr_dendro_plot(edgeR::cpm(counts, log = TRUE), groups = groups)

Function to extract contrasts and generate top tables and plots

Description

Function to extract contrasts and generate top tables and plots

Usage

diff_expr_extract_contrasts(
  contrasts = NULL,
  fit,
  fit2 = NULL,
  normcnt,
  out.l,
  do.voom = TRUE,
  quasi.likelihood = TRUE,
  out.dir,
  analysis.name = NULL,
  biomart = FALSE,
  biom.data.set = "hsapiens_gene_ensembl",
  biom.mart = "ensembl",
  host = "https://www.ensembl.org",
  biom.filter = "ensembl_gene_id",
  biom.attributes = c("ensembl_gene_id", "hgnc_symbol", "description"),
  biom.force.ensg = FALSE,
  biom.cache = NULL,
  use.cache = FALSE,
  sym.col = "hgnc_symbol",
  rm.dups = FALSE,
  p.thr = 0.05,
  fdr.thr = 0.05,
  logfc.thr = 1,
  numlab = 25,
  point.lab = TRUE,
  heatmap.topn = 100,
  hm.p.thr = 0.05,
  hm.fdr.thr = 0.05,
  hm.logfc.thr = 1,
  heatmap.split.expr = FALSE,
  color.blind.pal = "PuOr",
  n.pal.cols = 11,
  color.extremes = c("#3182BD", "#E6550D"),
  palette.length = NULL,
  anno.color = NULL,
  anno.name = "Sample Class",
  heatmap.main = NULL,
  font.size = 5,
  plots = TRUE,
  lists = TRUE,
  filtered.lists = TRUE,
  samp.info = NULL,
  samples = NULL,
  groups = NULL,
  sample.plot.names = NULL,
  de.plot.base.size = 16
)

Arguments

Value

The input list out.l, extended with a contrasts element holding one annotated result data.frame per contrast and, if plots were requested, with the MAplots, volcanoPlots and heatmapPlots elements. Called also for its side effects of writing result tables, plot files and Venn section spreadsheets below out.dir.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
d <- edgeR::calcNormFactors(edgeR::DGEList(counts, group = groups))
design <- diff_expr_make_design(si, groups)
contrasts <- diff_expr_make_contrasts(design, groups)
fit <- diff_expr_fit(counts, d, design, do.voom = FALSE, quasi.likelihood = TRUE)
normcnt <- edgeR::cpm(fit$d2, log = TRUE)
out <- diff_expr_extract_contrasts(contrasts = contrasts, fit = fit$fit,
                                   normcnt = normcnt, out.l = list(),
                                   do.voom = FALSE, out.dir = tempdir(),
                                   analysis.name = "demo", plots = FALSE,
                                   samp.info = si, samples = "SampleNames",
                                   groups = groups)
names(out$contrasts)

Function to filter counts

Description

Function to filter counts

Usage

diff_expr_filter_counts(counts, samp.info, strict = TRUE, min.samp = NULL)

Arguments

Details

In edgeR, it is recommended to remove features without at least 1 read per million in n of the samples, where n is the size of the smallest group of replicates (determined from the 'groups' vector).

Value

A matrix of counts with weakly expressed features and the non-informative htseq-count summary rows removed.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_read_counts(diffwrap_counts, si)
filtered <- diff_expr_filter_counts(counts, si, strict = TRUE)
c(before = nrow(counts), after = nrow(filtered))

Function to compute linear model fit and optionally apply 'voom' beforehand

Description

Function to compute linear model fit and optionally apply 'voom' beforehand

Usage

diff_expr_fit(
  counts,
  d,
  design,
  do.voom = TRUE,
  voom.fun = edgeR::voomLmFit,
  norm.method = c("quantile", "tmm"),
  use_weights = FALSE,
  quasi.likelihood = TRUE,
  bayes.trend = FALSE,
  bayes.robust = FALSE,
  pairs = NULL,
  pairs_col = NULL,
  block = FALSE,
  contrasts = NULL,
  disp = "tagwise.dispersion"
)

Arguments

Details

For experimental designs involving comparisons within as well as between subjects inter-subject needs to be computed. In this case, the column specified in the 'pairs' argument must assign the subjects to the treatment/tissue/etc groups. For example, if we have two treatments the effects of which are to be observed in each two tissues, this design would apply. The 'pairs' factor is passed to the functions duplicateCorrelation() and lmFit(). Quasi-likelihood pipeline: While the likelihood ratio test is a more obvious choice for inferences with GLMs, the QL F-test is preferred as it reflects the uncertainty in estimating the dispersion for each gene. It provides more robust and reliable error rate control when the number of replicates is small. The QL dispersion estimation and hypothesis testing is done by using the functions glmQLFit() and glmQLFTest().

Value

A named list. If do.voom=TRUE the elements are v (the EList of voom-transformed counts), fit (the MArrayLM linear model fit) and fit2 (the MArrayLM fit after eBayes(), with contrasts applied if these were supplied). Otherwise they are d (the input DGEList), d2 (the DGEList with estimated dispersions) and fit (the DGEGLM generalised linear model fit).

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
d <- edgeR::calcNormFactors(edgeR::DGEList(counts, group = groups))
design <- diff_expr_make_design(si, groups)
fit <- diff_expr_fit(counts, d, design, do.voom = FALSE, quasi.likelihood = TRUE)
class(fit$fit)

Function to standardize samp.info sample information data frame

Description

Function to standardize samp.info sample information data frame

Usage

diff_expr_get_samp_info(
  samp.info,
  samples,
  groups,
  ellipse.mapping.groups = NULL
)

Arguments

Value

A data.frame holding the sample sheet standardised to the conventions of this package: the sample and grouping columns renamed to SampleNames and Groups (and, where supplied, the ellipse grouping column to Ellipse), all coerced to factors with unused levels dropped, and the rows ordered by sample name.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info,
                              samples = "SampleName", groups = "Group")
head(si)

Function to generate a MDS plot using 'ggplot2'

Description

Function to generate a MDS plot using 'ggplot2'

Usage

diff_expr_ggplot_mds(
  counts,
  samp.name,
  groups,
  grp.nam = NULL,
  pairs = NULL,
  pairs.name = NULL,
  gene.selection = "common",
  dim.plot = c(1, 2),
  main = NULL
)

Arguments

Value

A ggplot object containing the multidimensional scaling plot.

See Also

[plotMDS()]

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
lcpm <- edgeR::cpm(counts, log = TRUE)
g <- diff_expr_ggplot_mds(lcpm, samp.name = colnames(lcpm), groups = groups)
class(g)

Function to generate a M-A plot using 'ggplot2'

Description

Function to generate a M-A plot using 'ggplot2'

Usage

diff_expr_ma_plot(
  dat,
  contr,
  id = NULL,
  sym.col = "gene_symbol",
  p.thr = 0.05,
  fdr.thr = 0.05,
  logfc.thr = 1,
  numlab = 25,
  out.dir,
  analysis.name = NULL,
  point.lab = TRUE,
  biom.attributes = c("ensembl_gene_id", "hgnc_symbol", "description"),
  font.size = 5,
  lists = TRUE,
  base.size = 16
)

Arguments

Value

A named list of ggplot objects with the elements FDR and Pval, holding the M-A plot with points highlighted by false discovery rate and by p-value, respectively.

Examples

set.seed(1)
d3 <- data.frame(gene_symbol = paste0("G", 1:100),
                 logFC = rnorm(100), PValue = runif(100), FDR = runif(100),
                 AveExpr = rnorm(100, 5))
rownames(d3) <- d3$gene_symbol
g <- diff_expr_ma_plot(d3, contr = "treated-control", sym.col = "gene_symbol",
                       out.dir = tempdir(), lists = FALSE)
names(g)

Function to make contrast matrix

Description

Function to make contrast matrix

Usage

diff_expr_make_contrasts(
  design,
  groups,
  pairs = NULL,
  block = FALSE,
  contrasts = NULL
)

Arguments

Details

The 'block' argument is used to specify whether the comparisons are to be made within AND between subjects or in the case of technical replicates, i.e., if the samples are not independent, in other words, correlated. @seealso [makeContrasts()]

Value

A contrast matrix as produced by makeContrasts, or NULL when no contrast matrix is needed because the comparisons are already inherent to the design matrix.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
groups <- stats::relevel(si$Groups, ref = "control")
design <- diff_expr_make_design(si, groups)
diff_expr_make_contrasts(design, groups)

Function to create design matrix

Description

Function to create design matrix

Usage

diff_expr_make_design(
  samp.info,
  groups,
  pairs = NULL,
  block = FALSE,
  use_weights = FALSE
)

Arguments

Details

The 'block' argument is used to specify whether the comparisons are to be made within AND between subjects or in the case of technical replicates, i.e., if the samples are not independent, in other words, correlated. If sample-specific quality weights are to be estimated by means of 'voomWithQualityWeights()' or 'voomLmFit()' and 'sample.weights' set to TRUE, 'use_weights' will be TRUE, enforcing a design matrix containing an 'intercept' column, i.e., where the columns reflect contrasts. The choice of the design matrix type impacts the estimated weights due to the effect of the intercept on the residual variance per sample in more complex designs, e.g., involving blocking factors, interactions or continuous covariates. The recommendation by the limma authors is to use the default design matrix, i.e., with intercept. With simple designs, the type of design matrix is not relevant. NOTE: This argument is not yet functional but a mere place-holder for future versions allowing for readily implemented more complex designs.

Value

A design matrix as produced by model.matrix: a means model without intercept (~0 + groups) for unpaired or blocked designs, or an additive model with intercept (~pairs + groups) when pairs enters the model as a fixed effect.

See Also

[model.matrix()]

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
groups <- stats::relevel(si$Groups, ref = "control")
diff_expr_make_design(si, groups)

Wrapper around 'limma::plotMDS' to generate a MDS plot

Description

Wrapper around 'limma::plotMDS' to generate a MDS plot

Usage

diff_expr_mds_plot(
  d,
  groups,
  n = 500,
  sample.plot.names = NULL,
  analysis.name = NULL,
  do.pdf = FALSE,
  out.dir
)

Arguments

Value

No return value. Called for its side effect of drawing a multidimensional scaling plot, optionally into a PDF file below out.dir.

See Also

[plotMDS()]

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
d <- edgeR::calcNormFactors(edgeR::DGEList(counts, group = groups))
diff_expr_mds_plot(d, groups = groups, do.pdf = FALSE, out.dir = tempdir())

Function to calculate pseudo counts representing batch-corrected normalised but untransformed values

Description

Function to calculate pseudo counts representing batch-corrected normalised but untransformed values

Usage

diff_expr_pseudo_counts(
  d,
  design,
  pairs = "pairs",
  disp = "tagwise.dispersion",
  do.cpm = TRUE
)

Arguments

Value

A matrix of pseudo counts in which the effect of the blocking variable has been removed; on the log2 counts per million scale if do.cpm=TRUE and on the count scale otherwise.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
si$Subject <- diffwrap_samp_info$Subject[order(diffwrap_samp_info$SampleName)]
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
d <- edgeR::calcNormFactors(edgeR::DGEList(counts, group = groups))
design <- diff_expr_make_design(si, groups, pairs = "Subject")
pc <- diff_expr_pseudo_counts(d = d, design = design, pairs = "pairs")
dim(pc)

Function to generate a histogram of the P-Value distribution

Description

Function to generate a histogram of the P-Value distribution

Usage

diff_expr_pval_hist_plot(d3)

Arguments

Value

Invisibly, an object of class histogram as returned by hist. Called for its side effect of drawing the p-value distribution.

Examples

d3 <- data.frame(ID = paste0("g", 1:200), PValue = runif(200))
diff_expr_pval_hist_plot(d3)

Function to read counts as produced by htseq-count

Description

Function to read counts as produced by htseq-count

Usage

diff_expr_read_counts(expr.dat, samp.info, miRSEQ = FALSE)

Arguments

Value

A matrix of raw counts with features in rows and samples in columns, restricted to the samples listed in samp.info and ordered as they are there.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_read_counts(diffwrap_counts, si)
dim(counts)

Function to generate a Volcano plot using 'ggplot2'

Description

'diff_expr_volcano_plot' generates two Volcano plots highlighting genes that are differentially expressed beyond custom thresholds for significance (set by parameters 'p.thr' and fdr.thr') and differential expression level (set by parameter 'logfc.thr').

Usage

diff_expr_volcano_plot(
  d3,
  id,
  sym.col = "gene_symbol",
  main = NULL,
  p.thr = 0.05,
  fdr.thr = 0.05,
  logfc.thr = 1,
  numlab = 25,
  point.lab = TRUE,
  base.size = 16
)

Arguments

Value

A named list of ggplot objects with the elements FDR and Pval, holding the volcano plot with points highlighted by false discovery rate and by p-value, respectively. Both plots are additionally printed, so that they are captured when the function is called with an open plotting device.

Examples

set.seed(1)
d3 <- data.frame(gene_symbol = paste0("G", 1:100),
                 logFC = rnorm(100), PValue = runif(100), FDR = runif(100))
g <- diff_expr_volcano_plot(d3, id = "gene_symbol", sym.col = "gene_symbol")
names(g)

Helper function to generate an output table with only most relevant columns

Description

Helper function to generate an output table with only most relevant columns

Usage

diffr_expr_generate_cleaned_de_table_output(
  contrast,
  annotated.normcnt,
  out.dir,
  samp.name.and.group.key,
  analysis.name = NULL,
  filtered.lists = TRUE,
  fdr.thr = 0.05,
  logfc.thr = 1
)

Arguments

Value

No return value. Called for its side effect of writing a reduced and optionally significance-filtered differential expression table to a tab-separated file below out.dir.

Examples

set.seed(1)
key <- data.frame(group = rep(c("control", "treated"), each = 4),
                  row.names = paste0("S", 1:8))
d3 <- data.frame(ID = paste0("g", 1:50), gene_symbol = paste0("G", 1:50),
                 logFC = rnorm(50), FDR = runif(50))
diffr_expr_generate_cleaned_de_table_output(
  contrast = "treated-control", annotated.normcnt = d3,
  samp.name.and.group.key = key, out.dir = tempdir(),
  analysis.name = "demo")

Function to create a heatmap from differential gene expression values

Description

Function to create a heatmap from differential gene expression values

Usage

diffr_pheatmap(
  expr.mat,
  clinical.mat,
  scale.fl = "none",
  legend.fl = TRUE,
  row.clust = TRUE,
  col.clust = TRUE,
  biserial.fl = FALSE,
  quantile.breaks.fl = FALSE,
  signif.stars.fl = FALSE,
  cell.size = 8,
  font.size = 10,
  color.blind.pal = "PuOr",
  n.pal.cols = 11,
  color.extremes = c("#3182BD", "#E6550D"),
  palette.length = NULL,
  anno.color = NULL,
  main = NULL,
  add.main = NULL,
  filt.info = NULL
)

Arguments

Details

The plot is produced as a side effect of function 'pheatmap_plots()' by printing the individual plot objects generated by 'diffr_pheatmap()'.

Value

Returns a list of pheatmap plot objects used in the pheatmap_plots() function.

Author(s)

Bogdan Iancu - Genevia Technologies Oy

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
expr <- edgeR::cpm(counts, log = TRUE)[1:30, ]
clin <- data.frame(Group = si$Groups, row.names = si$SampleNames)
hm <- diffr_pheatmap(expr, clin)

Function to produce a Venn diagram of differentially expressed gene tables

Description

Function to produce a Venn diagram of differentially expressed gene tables

Usage

diffr_venn(list.comp.tables, join.vec = "gene_symbol", .log = FALSE)

Arguments

Details

The actual plot is produced in the main plotting function by means of 'grid::grid.draw()' using the plot object as input.

Value

A list containing the Venn diagram grid object and the intersected (or joined) input tables.

Author(s)

Bogdan Iancu - Genevia Technologies Oy

Examples

mk <- function(sig) data.frame(gene_symbol = paste0("G", sig), fdr = 0.01)
tabs <- list(A = mk(1:20), B = mk(10:30))
v <- diffr_venn(tabs, join.vec = "gene_symbol")
names(v)

Differential expression analysis of RNA-Seq data

Description

This package provides functionality for differential expression analysis of read counts from mRNA sequencing data or miRNA expression values generated by the CAP-miRSEQ 'expression_reports.sh' script. The workflow follows the 'edgeR'-'limma' expression data analysis pipeline providing options for different approaches, such as "pure" 'edgeR', 'voom' or paired samples. The functions in the package generate text files with differential expression lists, optionally annotated with information from 'biomart', expression summary plots as well as several QC plots.

Details

The package contains utilities and methods for differential expression analysis of RNA-Seq data. It comes with a wrapper function, "diffExpr", performing all steps of the analysis. Most sub-functions are exported and can be called by the user, as well, if desired. These functions may be applicable to different kinds of data/input, rely, however, on the conventions set for this package.

Acknowledgements

This work was largely supported by university-level strategic and profiling area funding schemes at University of Turku and Tampere University. Downstream analysis functionality including all pathway analyses and visualisations as well as some DGE visualisations and helper functions were created or further developed at Genevia Technologies Oy.

Author(s)

Vidal Fey <vidal.fey@gmail.com>, Meeri Pekkarinen <meeri.pekkarinen@gmail.com>, Reija Hieta <reija.hieta@geneviatechnologies.com>, Bogdan Iancu <bogdan.iancu@abo.fi>, Adrien Janssens <Adrien.RFM.Janssens@proton.me> Maintainer: Vidal Fey <vidal.fey@gmail.com>

Simulated RNA-Seq read counts for the package examples

Description

A small simulated count matrix used throughout the examples, tests and the vignette. The data are not real measurements: they were generated from a negative binomial (gamma-Poisson) model so that the package can be demonstrated and tested without any external data or network access.

Usage

diffwrap_counts

Format

An integer matrix with 405 rows (400 genes plus 5 summary rows) and 8 columns (samples S01-S08). Row names are gene identifiers, column names sample names.

Details

The experiment consists of eight samples in two groups of four ("control" and "treated"). Each of the four subjects (P1-P4) contributes one control and one treated sample, so the data support paired and blocked analyses as well as the simple unpaired comparison. A subject-specific offset is built in, giving the paired designs a real effect to remove.

Sixty of the 400 genes are truly differentially expressed (30 up-, 30 down-regulated) with absolute log2 fold changes between 1.2 and 3. The biological dispersion is 0.15 (biological coefficient of variation of about 0.39) and library sizes are roughly 2-3 million reads.

Forty genes are given deliberately low expression so that they straddle the default strict filtering threshold of diff_expr_filter_counts() (more than 5 counts per million in at least half of the samples, which corresponds to about 13 counts here); 28 of them are removed by that filter.

The matrix additionally carries the five htseq-count summary rows (__no_feature, __ambiguous, __too_low_aQual, __not_aligned and __alignment_not_unique) that diff_expr_filter_counts() is expected to strip. They are retained on purpose so that the filtering step can be demonstrated.

Gene identifiers are well-formed but fictitious Ensembl gene IDs; they do not correspond to real genes and will not return annotation from 'biomart'.

Source

Simulated by ‘data-raw/make_example_data.R’ from ‘inst/extdata/example_counts.tsv’.

See Also

diffwrap_samp_info for the matching sample sheet.

Examples

data(diffwrap_counts)
dim(diffwrap_counts)
head(diffwrap_counts[, 1:4])
# the htseq-count summary rows that get filtered out:
diffwrap_counts[grep("^__", rownames(diffwrap_counts)), 1:4]

Sample sheet accompanying the simulated example counts

Description

The sample information table matching diffwrap_counts, in the layout expected by diffExpr() and diff_expr_get_samp_info().

Usage

diffwrap_samp_info

Format

A data.frame with 8 rows and 4 columns:

SampleName

character. Sample identifier, matching the column names of diffwrap_counts.

Group

character. Experimental group, "control" or "treated".

Subject

character. Subject identifier (P1-P4); pass to pairs for paired or blocked designs.

PlotName

character. Human-readable sample label for plots.

Details

The Subject column pairs each control sample with a treated sample from the same subject and is the column to pass to the pairs argument. PlotName supplies prettier labels for plotting and is the column to pass to sample.plot.names.

Source

Simulated by ‘data-raw/make_example_data.R’ from ‘inst/extdata/example_samp_info.tsv’.

See Also

diffwrap_counts for the matching count matrix.

Examples

data(diffwrap_samp_info)
diffwrap_samp_info
table(diffwrap_samp_info$Group, diffwrap_samp_info$Subject)

Close a graphics device only if it is actually open

Description

Defensive wrapper around dev.off(). Plotting helpers occasionally close the device they were handed, and a bare dev.off() with nothing but the null device left fails with “cannot shut down device 1 (the null device)”. This closes the requested device when it is still in dev.list() and is a silent no-op otherwise, so cleanup code and on.exit() handlers cannot error.

Usage

dw_dev_off(which = NULL)

Arguments

Value

Invisibly, TRUE if a device was closed and FALSE otherwise.

Convert a numeric value of seconds to human-readable duration

Description

Convert a numeric value of seconds to human-readable duration

Usage

fmt_dur(secs)

Arguments

Value

Character string with run time information.

Helper function for formatting the gene ID column of enrichment data frame.

Description

Helper function for formatting the gene ID column of enrichment data frame.

Usage

format_ensembl_ids_annotated_to_term(result, species, which.split = ",")

Arguments

Details

Enrichment tools report all genes annotated to a certain term in one string, separated by comma or similar. The function splits each row into individual IDs before converting and later re-collapses converted IDs.

Value

The input enrichment result with the Ensembl gene identifiers in the gene column replaced by the corresponding gene symbols.

Function to define breaks to be used for changing the palette of the heatmap.

Description

Function to define breaks to be used for changing the palette of the heatmap.

Usage

get_hm_breaks(
  expr.mat,
  scale.fl = "row",
  palette.length = 100,
  quantile.breaks.fl = FALSE,
  n = 11
)

Arguments

Value

a numeric vector containing the breaks to be used in the heatmap (see 'breaks' in [pheatmap()]).

See Also

[quantile()]

Function to compute colour palettes to be used in the heatmap.

Description

Function to compute colour palettes to be used in the heatmap.

Usage

get_hm_colors(
  palette.length = 100,
  color.extremes = c("#3182BD", "#E6550D"),
  breaks = NA,
  color.blind.pal = NULL,
  n.pal.cols = 11
)

Arguments

Value

a character vector containing the colour codes to be used in the heatmap.

See Also

[colorRampPalette()]

Make a biplot of PCA output data using ggplot2.

Description

Make a biplot of PCA output data using ggplot2.

Usage

ggbiplot.n(
  pcobj,
  choices = 1:2,
  scale = 1,
  pc.biplot = TRUE,
  obs.scale = 1 - scale,
  var.scale = scale,
  groups = NULL,
  grp.nam = NULL,
  ellipse = FALSE,
  ellipse.prob = 0.68,
  labels = NULL,
  labels.size = 3,
  alpha = 1,
  var.axes = TRUE,
  circle = FALSE,
  circle.prob = 0.69,
  varname.size = 3,
  varname.adjust = 1.5,
  varname.abbrev = FALSE,
  point.size = 1,
  axes.title.size = 1,
  legend.text.size = 1,
  legend.title.size = 1,
  ellipse.lwd = 1,
  main = NULL,
  fix.aspect = TRUE,
  tweak = TRUE,
  tidy = TRUE,
  ...
)

Arguments

Value

The final plot object returned by ggplot.

Examples

si <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
pca <- diff_expr_PCA(edgeR::cpm(counts, log = TRUE), n = 100)
g <- ggbiplot.n(pca, groups = groups)
class(g)

Function to create the annotation colour list used in the heatmap

Description

Function to create the annotation colour list used in the heatmap

Usage

make_pheatmap_anno_color(clinical.mat)

Arguments

Value

A named list of named character vectors of colours, one element per annotation variable, in the form expected by the annotation_colors argument of pheatmap.

Function to generate heatmap of gene expression values

Description

Function to generate heatmap of gene expression values

Usage

pheatmap_plots(
  d3,
  id,
  sym.col = "gene_symbol",
  samp.info,
  samples,
  groups,
  sample.plot.names = NULL,
  main = NULL,
  add.main = NULL,
  color.blind.pal = "PuOr",
  n.pal.cols = 11,
  color.extremes = c("#3182BD", "#E6550D"),
  palette.length = NULL,
  anno.color = NULL,
  anno.name = "Sample Class",
  p.thr = 0.05,
  fdr.thr = 0.05,
  logfc.thr = 1,
  topn = 100,
  split.expr = FALSE
)

Arguments

Details

The actual heatmap is produced as a side effect by printing the individual plot components.

This is the convenience entry point used by diffExpr() to build a contrast's heatmaps; it is exported so the same heatmaps can be regenerated from a stored result table (e.g. res$contrasts[[<contrast>]]) with different thresholds or colours, without re-running the pipeline. The lower-level work horse it calls is diffr_pheatmap.

Value

Returns a named list of heatmap objects: an element fdr (built from the FDR-filtered genes) and/or pval (used only when nothing passes the FDR threshold), each a list holding the regular heatmap and the gene/sample correlogram objects. Also called for its side effect of drawing those objects on the active graphics device.

Examples

si     <- diff_expr_get_samp_info(diffwrap_samp_info, "SampleName", "Group")
counts <- diff_expr_filter_counts(diff_expr_read_counts(diffwrap_counts, si), si)
groups <- stats::relevel(si$Groups, ref = "control")
d      <- edgeR::calcNormFactors(edgeR::DGEList(counts, group = groups))
design <- diff_expr_make_design(si, groups)
contr  <- diff_expr_make_contrasts(design, groups)
fit    <- diff_expr_fit(counts, d, design, do.voom = FALSE, quasi.likelihood = TRUE)
## a per-contrast table like the ones diffExpr() stores in res$contrasts
de <- edgeR::glmQLFTest(fit$fit, contrast = contr[, 1])
d3 <- merge(edgeR::cpm(fit$d2, log = TRUE),
            as.data.frame(edgeR::topTags(de, n = Inf)), by = "row.names")
names(d3)[1] <- "ID"; d3$gene_symbol <- d3$ID
## regenerate the heatmap at a stricter cut-off and a different palette, no refit
hm <- pheatmap_plots(d3, id = "ID", samp.info = si, samples = "SampleNames",
                     groups = groups, fdr.thr = 0.01, topn = 30, color.blind.pal = "RdBu")

Function for making network visualisation based on enrichment result, DE gene table and thresholding. Can take the input tables either as data frames or Excel files

Description

Function for making network visualisation based on enrichment result, DE gene table and thresholding. Can take the input tables either as data frames or Excel files

Usage

plot_enrichment_network(
  enrichment.result,
  DE.result,
  plot.filename,
  show.terms = 5,
  logfc.thr = NULL,
  fdr.thr = NULL,
  pdf.width = 11,
  pdf.height = 8,
  legend.cex.main = 0.8,
  legend.cex.text = 0.7
)

Arguments

Value

No return value. Called for its side effect of writing a network visualisation of the enrichment result to the file given by plot.filename.

Examples

# builds a network plot from an enrichment result and writes it to 'plot.filename'
if (requireNamespace("igraph", quietly = TRUE)) {
  enr <- data.frame(Description = c("Pathway A", "Pathway B"),
                    genes = c("G1,G2,G3", "G3,G4,G5"))
  de <- data.frame(gene_symbol = paste0("G", 1:5),
                   logFC = c(2.5, -1.8, 1.2, -2.1, 0.9),
                   FDR = c(0.001, 0.002, 0.01, 0.003, 0.02))
  plot_enrichment_network(enr, de,
                          plot.filename = file.path(tempdir(), "network.pdf"))
}

Helper function for enrichment visualisations: prepare plot legend y coordinates and labels

Description

Helper function for enrichment visualisations: prepare plot legend y coordinates and labels

Usage

prepare_scale_for_legend(
  scale.minimum,
  scale.maximum,
  int.values.for.ticks = NULL
)

Arguments

Details

The coordinates are used in a network graph and minimum and maximum scale values correspond to minimum and maximum fold-changes, by default.

Value

A named list with the elements scale.y.coordinates (the vertical positions of the legend tick labels) and scale.labels (the corresponding label values).

Examples

prepare_scale_for_legend(scale.minimum = -3, scale.maximum = 4)

Helper function that returns a volcano plot.

Description

Builds a single volcano plot for one significance measure. Points are coloured by an 'EnhancedVolcano'-style, colour-blind safe (Okabe-Ito) four-level scheme distinguishing genes that pass neither threshold ("not sign."), the fold-change threshold only, the significance threshold only, or both. Only the doubly-significant genes are eligible for text labels.

Usage

prepare_volcano_of_given_property(
  data.df,
  property.to.plot = c("fdr", "p"),
  property.column,
  property.thr = 0.05,
  logfc.thr = 1,
  main = NULL,
  numlab = 25,
  point.lab = TRUE,
  sym.col = "gene_symbol",
  pretty.breaks = FALSE,
  base.size = 16
)

Arguments

Value

A ggplot object containing the volcano plot for the requested significance measure.

Examples

set.seed(1)
d3 <- data.frame(gene_symbol = paste0("G", 1:100),
                 logFC = rnorm(100), PValue = runif(100))
g <- prepare_volcano_of_given_property(d3, property.to.plot = "p",
                                       property.column = "PValue",
                                       property.thr = 0.05, logfc.thr = 1)
class(g)

Function to define quantile breaks to be used for changing the palette of the heatmap.

Description

Function to define quantile breaks to be used for changing the palette of the heatmap.

Usage

quantile_breaks(xs, n = 20)

Arguments

Value

A numeric vector of unique quantile break points, of length at most n.

Note

Function by: Kamil Slowikowski, https://github.com/slowkow/slowkow.com/blob/master/_rmd/2017-02-16-heatmap-tutorial.R,

See Also

[quantile()]

Function which reorders the levels of a column of a data frame specified as a factor

Description

Function which reorders the levels of a column of a data frame specified as a factor

Usage

reorderFactors(df, column = "my_column_name", desired_level_order)

Arguments

Details

The desired order is a vector containing the levels of the factor in the desired order.

Value

The input data.frame with the levels of the specified factor column reordered.

Note

Written by https://stackoverflow.com/users/1701600/boern.

Function to run GO term enrichment analysis using the 'topGO' package.

Description

Function to run GO term enrichment analysis using the 'topGO' package.

Usage

run.topGO(
  background,
  foreground,
  ontologies = c("BP"),
  organism,
  ID_type = "ENSEMBL",
  pAdjustMethod = "BH"
)

Arguments

Details

The columns in the final table produced by topGO (their description is fetched from topGo documentation): - Annotated : number of genes in org.Hs.eg.db which are annotated with the GO-term. - Significant : number of genes belonging to your input which are annotated with the GO-term. - Expected : show an estimate of the number of genes a node of size Annotated would have if the significant genes were to be randomly selected from the gene universe. - pvalues : pvalue obtained after the test

Column 'p.adj.weight01' represents the adjusted p-value for the weight01 algorithm.

Value

A data frame with enriched terms and p-values from Fisher's exact test, using different algorithms: elim, classic, weight01.

Author(s)

Bogdan Iancu - Genevia Technologies Oy

Examples

# needs org.Hs.eg.db; builds a genome-wide GO universe, so it is slow
if (requireNamespace("topGO", quietly = TRUE) &&
    requireNamespace("org.Hs.eg.db", quietly = TRUE) &&
    requireNamespace("AnnotationDbi", quietly = TRUE)) {
  fg <- c("ENSG00000141510", "ENSG00000012048", "ENSG00000139618")
  # a small background keeps the example quick; a real analysis uses the full gene universe,
  # e.g. AnnotationDbi::keys(org.Hs.eg.db::org.Hs.eg.db, keytype = "ENSEMBL")
  bg <- unique(c(fg, utils::head(
    AnnotationDbi::keys(org.Hs.eg.db::org.Hs.eg.db, keytype = "ENSEMBL"), 2000)))
  go <- tryCatch(run.topGO(background = bg, foreground = fg, ontologies = "BP",
                           organism = "org.Hs.eg.db", ID_type = "ENSEMBL"),
                 error = function(e) {
                   message("topGO enrichment could not be run: ", conditionMessage(e))
                   NULL
                 })
}

Wrapper for executing various enrichment analyses

Description

runEnrichmentAnalyses enables the auto-run of some over-representation analysis (ORA) and gene set enrichment analysis (GSEA) functions for the output of diffExpr main wrapper. Functions in various R-packages (clusterProfiler, topGO, gProfileR) are integrated. Currently supports human or mouse!

Usage

runEnrichmentAnalyses(
  diffr.wrapper.output,
  analysis.name = "enrichment",
  use.background.from.diffr.output = TRUE,
  use.pval.in.DE.filtering.if.no.sign.fdrs = FALSE,
  out.dir,
  species = "human",
  p.thr = 0.05,
  fdr.thr = 0.05,
  logfc.thr = 1,
  do.plot = FALSE,
  plot.fdr.thr = fdr.thr,
  plot.logfc.thr = logfc.thr,
  plot.num.terms = 5,
  enrichment.methods = c("clusterProfilerGO", "clusterProfilerKEGG", "gProfileR",
    "topGO"),
  clusterProfilerGO.params = list(analysis.approach = "ORA", do.similarity.filtering =
    FALSE, min.gene.set.size = 10, max.gene.set.size = 1000, ontology = "BP", min.overlap
    = 2, p.adjust.method = "BH"),
  clusterProfilerKEGG.params = list(analysis.approach = "ORA", min.gene.set.size = 10,
    max.gene.set.size = 1000, ontology = "BP", min.overlap = 2, p.adjust.method = "BH"),
  gProfileR.params = list(data.sources = "GO:BP", show.only.significant = TRUE,
    measure_underrepresentation = FALSE, evidence_codes = TRUE, domain_scope =
    "annotated", highlight = TRUE),
  topGO.params = list(ontologies.used = c("BP"), org = "hsapiens")
)

Arguments

Value

A list of all relevant objects generated in the course of the enrichment analyses

Examples

# needs annotation packages; run on real data with mappable gene IDs
if (requireNamespace("clusterProfiler", quietly = TRUE) &&
    requireNamespace("org.Hs.eg.db", quietly = TRUE) &&
    requireNamespace("AnnotationDbi", quietly = TRUE)) {
  out.dir <- file.path(tempdir(), "diffwrap_demo")
  dir.create(out.dir, showWarnings = FALSE)
  res <- diffExpr(diffwrap_counts, diffwrap_samp_info, samples = "SampleName",
                  groups = "Group", control = "control", analysis.name = "demo",
                  out.dir = out.dir, enr.do = FALSE)
  # wrapped in tryCatch(): enrichment depends on annotation databases and, for some
  # methods, on remote services, neither of which may be available on a check machine
  enr <- tryCatch(
    runEnrichmentAnalyses(res, analysis.name = "demo", out.dir = out.dir,
                          species = "human", enrichment.methods = "clusterProfilerGO"),
    error = function(e) {
      message("Enrichment could not be run: ", conditionMessage(e))
      NULL
    })
}

Runs clusterProfiler GO enrichment function for a DEG list or for a ranked gene list.

Description

Runs clusterProfiler GO enrichment function for a DEG list or for a ranked gene list.

Usage

run_clusterProfiler_GO(
  input_genes,
  background_genes = "",
  file_name = NULL,
  ordered_query = FALSE,
  id_type = "ENSEMBL",
  ontology = "BP",
  OrgDb = "org.Hs.eg.db",
  pvalueCutoff = 0.05,
  min_set_size = 10,
  max_set_size = 1000,
  min_overlap = 2,
  pAdjustMethod = "BH",
  similarity_filtering = FALSE,
  rng.seed = NULL
)

Arguments

Value

A table listing statistically significant enrichment results according to threshold set in pvalueCutoff The table is also saved in xlsx format with user-specified name.

Examples

if (requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
  genes <- c("ENSG00000141510", "ENSG00000012048", "ENSG00000139618")
  # wrapped in tryCatch(): the result depends on the installed annotation database
  go <- tryCatch(run_clusterProfiler_GO(input_genes = genes, ontology = "BP",
                                        OrgDb = "org.Hs.eg.db", id_type = "ENSEMBL"),
                 error = function(e) {
                   message("GO enrichment could not be run: ", conditionMessage(e))
                   NULL
                 })
}

Runs clusterProfiler KEGG enrichment function for a DEG list or for a ranked gene list.

Description

Runs clusterProfiler KEGG enrichment function for a DEG list or for a ranked gene list.

Usage

run_clusterProfiler_KEGG(
  input_genes,
  background_genes = "",
  file_name = NULL,
  ordered_query = FALSE,
  id_type = "kegg",
  organism = "hsa",
  pvalueCutoff = 0.05,
  min_set_size = 10,
  max_set_size = 1000,
  min_overlap = 2,
  pAdjustMethod = "BH",
  rng.seed = NULL
)

Arguments

Value

A table listing statistically significant enrichment results according to threshold set in pvalueCutoff The table is also saved in xlsx format with user-specified name.

Examples

# enrichKEGG downloads pathway data from the KEGG web service, so this needs network access.
# Wrapped in tryCatch() so that an unreachable service reports the problem instead of
# failing the example.
if (requireNamespace("clusterProfiler", quietly = TRUE) &&
    requireNamespace("org.Hs.eg.db", quietly = TRUE)) {
  # KEGG works on Entrez identifiers
  kegg <- tryCatch(run_clusterProfiler_KEGG(input_genes = c("7157", "672", "675"),
                                            organism = "hsa", id_type = "kegg"),
                   error = function(e) {
                     message("KEGG service not reachable: ", conditionMessage(e))
                     NULL
                   })
}

Runs gprofiler function for a DEG list or for a ranked gene list.

Description

Runs gprofiler function for a DEG list or for a ranked gene list.

Usage

run_gprofiler(
  input_genes,
  organism = "hsapiens",
  background_genes = "",
  file_name = NULL,
  ordered_query = FALSE,
  multi_query = FALSE,
  show_only_significant = FALSE,
  evidence_codes = FALSE,
  exclude_iea = FALSE,
  measure_underrepresentation = FALSE,
  max_p_value = 0.5,
  correction_method = "fdr",
  domain_scope = "annotated",
  data_sources = NULL,
  highlight = FALSE
)

Arguments

Value

A table listing statistically significant enrichment results according to threshold set in padjCutoff. The table is also saved in xlsx format with user-specified name.

Examples

# Needs network access to the g:Profiler service. Wrapped in tryCatch() so that an
# unreachable server reports the problem instead of failing the example.
if (requireNamespace("gprofiler2", quietly = TRUE)) {
  gp <- tryCatch(run_gprofiler(input_genes = c("ENSG00000141510", "ENSG00000012048"),
                               organism = "hsapiens"),
                 error = function(e) {
                   message("g:Profiler not reachable: ", conditionMessage(e))
                   NULL
                 })
}