Performing differential co-expression analysis using dcanr

Introduction

This package provides methods to perform differential co-expression analysis and to evaluate differential co-expression methods using simulated data. Differential co-expression analysis attempts to identify gene-gene associations that change across conditions. Currently, 10 methods that identify changes between binary conditions are included: 8 are novel Bioconductor implementations of previously published methods, and; 2 are accessed through interfaces to existing packages.

This vignette focuses on the application of differential co-expression inference methods to real data. Available methods, putative pipelines, and visualisations provided by the method are introduced.

Installation

Download the package from Bioconductor

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("dcanr")

Or install the development version of the package from Github.

BiocManager::install("DavisLaboratory/dcanr")

Load the installed package into an R session.

library(dcanr)

Available inference methods

The package implements 10 methods to infer differential co-expression networks across binary conditions. The list of available methods can be accessed by the dcMethods() function.

library(dcanr)
dcMethods()

A generic differential co-expression analysis pipeline

A differential co-expression pipeline generally consists of 4 steps:

  1. Generating scores/statistics for each pair of genes
  2. Assessing scores using statistical tests
  3. Correcting for multiple hypothesis testing
  4. Selecting differential associations

Not all methods follow this pipeline. EBcoexpress computes posterior probabilities therefore no statistical test needs to be performed and steps 2-3 can be skipped. Like-wise DiffCoEx does not perform any statistical tests and instead performs a soft-thresholding on the scores. FTGI performs a statistical test and $p$-values from this test are used as scores, therefore step 2 is skipped. A standard analysis with the z-score method using all 4 steps is shown here.

Load an example dataset (simulated)

We first load an example simulated dataset (included in the package) to extract the expression matrix and condition vector. Please note that multiple knock-down experiments are performed per simulation and we use one such knock-down as a condition here. The list of all knock-downs can be retrieved using getConditionNames().

#load data
data(sim102)
#get available conditions
getConditionNames(sim102)
#get expression data and conditions for 'UME6' knock-down
simdata <- getSimData(sim102, cond.name = 'UME6', full = FALSE)
emat <- simdata$emat
ume6_kd <- simdata$condition
print(emat[1:5, 1:5]) #149 genes and 406 samples
head(ume6_kd) #NOTE: binary conditions encoded with 1's and 2's

Step 1: Compute scores

All inference methods can be accessed using the same call therefore making it easier to change between methods. Method specific parameters can be passed to this function and will be managed accordingly. The default inference method is z-score therefore it does not need to be specified via dc.method. We recommend using the Spearman correlation as a measure of correlation as it is robust to outliers which may be present in RNA-seq data.

#apply the z-score method with Spearman correlations
z_scores <- dcScore(emat, ume6_kd, dc.method = 'zscore', cor.method = 'spearman')
print(z_scores[1:5, 1:5])

Step 2: Perform a statistical test

Appropriate statistical tests are automatically selected for the method applied. Tests are applied on the result of the dcScore() function (z-test for the z-score method and permutation tests for other methods). The testing function returns the score matrix (unmodified) if the method is either EBcoexpress, FTGI or DiffCoEx.

NOTE: Do NOT modify the result of the scoring method as this will result in failure of the testing function. This is intended as tests should be performed for all computed scores to prevent bias in the subsequent correction for multiple hypothesis testing. The same applies for the next step.

#perform a statistical test: the z-test is selected automatically
raw_p <- dcTest(z_scores, emat, ume6_kd)
print(raw_p[1:5, 1:5])

For methods such as MINDy that require a permutation test, the number of permutations can be specified by the B parameter. Permutation tests are computationally expensive therefore we also provide a parallelised implementation. See the help page of dcTest for examples.

Step 3: Correcting for multiple hypothesis testing

Since all pairwise combinations of genes are tested, $p$-values need to be adjusted. Given $n$ genes, the total number of hypothesis is $\frac{n(n-1)}{2}$ as the score matrices are symmetric. Adjustment is performed accordingly. The default adjustment function is stats::p.adjust with the 'fdr' method used, however, custom functions and their parameters can be specified instead. dcAdjust provides a wrapper to apply an adjustment method to the raw $p$-value matrix. Results from EBcoexpress and DiffCoEx remain unmodified.

#adjust p-values (raw p-values from dcTest should NOT be modified)
adj_p <- dcAdjust(raw_p, f = p.adjust, method = 'fdr')
print(adj_p[1:5, 1:5])

Step 4: Generating the differential co-expression network

The last step is thresholding the score/adjusted $p$-value matrix to select differential associations. Default adjusted $p$-value thresholds of 0.1 are applied where statistical tests are performed (to control for FDR at 0.1). Results are presented as an igraph object shown below where edges are coloured based on the score (negative to positive scores are represented using the purple to green gradient of colours).

library(igraph)

#get the differential network
dcnet <- dcNetwork(z_scores, adj_p)
plot(dcnet, vertex.label = '')
#convert to an adjacency matrix
adjmat <- as_adj(dcnet, sparse = FALSE)
print(adjmat[1:5, 1:5])
#convert to a data.frame
edgedf <- as_data_frame(dcnet, what = 'edges')
print(head(edgedf))

Session info {.unnumbered}

sessionInfo()


Try the dcanr package in your browser

Any scripts or data that you put into this service are public.

dcanr documentation built on Nov. 8, 2020, 5:48 p.m.