MetaVolcanoR: Differential expression meta-analysis tool
June 5, 2019
Introduction
The measurement of how gene expression changes under different biological conditions is necessary to reveal the gene regulatory programs that determine the cellular phenotype.
Comparing the expresion of genes under a given condition against a reference biological state is usually applied to identify sets of differentially expressed genes (DEG). These DEG point out the genomic regions functionally relevant under the biological condition of interest.
Athough individual genome-wide expression studies have small signal/noise ratio, today’s genomic data availability usually allows to combine differential gene expression results from independent studies to overcome this limitation.
Databases such as GEO (https://www.ncbi.nlm.nih.gov/geo/), SRA (https://www.ncbi.nlm.nih.gov/sra), ArrayExpress, (https://www.ebi.ac.uk/arrayexpress/), and ENA (https://www.ebi.ac.uk/ena) offer systematic access to vast amounts of transcriptome data. There exists more than one gene expression study for many biological conditions. This redundancy could be exploit by meta-analysis approaches to reveal genes that are consistently and differentially expressed under given conditions.
MetaVolcanoR was designed to identify the genes whose expression is consistently perturbed across several studies.
Usage
Overview
The MetaVolcanoR R package combines differential gene expression results. It implements three strategies to summarize gene expression activities from different studies. i) Random Effects Model (REM) approach. ii) a vote-counting approach, and iii) a p-value combining-approach. MetaVolcano exploits the Volcano plot reasoning to visualize the gene expression meta-analysis results.
Installation
BiocManager::install("MetaVolcanoR", eval = FALSE)Load library
library(MetaVolcanoR)Input Data
Users should provide a named list of objects containing differential gene expression results. Each object of the list
must contain a gene name, a fold change, and a p-value variable. It is highly recommended to also include the variance or the confidence intervals of the fold change variable.
Take a look at the demo data. It includes differential gene expression results from five studies.
data(diffexplist)
class(diffexplist)## [1] "list"head(diffexplist[[1]])## Symbol Log2FC pvalue CI.L CI.R
## 1 A1BG -0.70126879 0.000140100 -1.0087857 -0.39375189
## 2 A1BG-AS1 -0.25106351 0.008694757 -0.4304790 -0.07164803
## 3 A1CF 0.03332573 0.615989488 -0.1036882 0.17033968
## 4 A2M 0.83504214 0.018550388 0.1568214 1.51326289
## 5 A2ML1 0.03942552 0.843222358 -0.3728473 0.45169836
## 6 A4GALT -0.20815882 0.282488068 -0.6025247 0.18620708length(diffexplist)## [1] 5Implemented meta-analysis approaches
Random Effect Model MetaVolcano
The REM MetaVolcano summarizes the gene fold change of several studies taking into account the variance. The REM estimates a summary p-value which stands for the probability of the summary fold-change is not different than zero. Users can set the metathr parameter to highlight the top percentage of the most consistently perturbed genes. This perturbation ranking is defined following the topconfects approach.
meta_degs_rem <- rem_mv(diffexp=diffexplist,
pcriteria="pvalue",
foldchangecol='Log2FC',
genenamecol='Symbol',
geneidcol=NULL,
collaps=FALSE,
llcol='CI.L',
rlcol='CI.R',
vcol=NULL,
cvar=TRUE,
metathr=0.01,
jobname="MetaVolcano",
outputfolder=".",
draw='HTML',
ncores=1)## index Symbol Log2FC_1 CI.L_1 CI.R_1 vi_1 Log2FC_2
## 1 4795 MXRA5 0.8150851 0.3109324 1.3192377 0.06616251 1.3001104
## 2 2166 COL6A6 -1.7480348 -2.5780749 -0.9179947 0.17934364 -0.8388366
## 3 2053 CIDEA NA NA NA NA NA
## 4 7115 SULT1A4 0.9689025 0.5103475 1.4274575 0.05473571 0.7513323
## 5 130 ACACB -0.8431142 -1.4708480 -0.2153804 0.10257437 -1.1119841
## 6 6528 SLC27A2 -0.6782948 -0.9931027 -0.3634869 0.02579759 -1.8916655
## CI.L_2 CI.R_2 vi_2 Log2FC_3 CI.L_3 CI.R_3 vi_3
## 1 0.6603306 1.9398901 0.10654886 1.1895480 0.8401301 1.5389659 0.031781777
## 2 -1.3578456 -0.3198277 0.07011930 -1.0300519 -1.4730328 -0.5870710 0.051080819
## 3 NA NA NA -1.0111528 -1.3226326 -0.6996729 0.025255027
## 4 0.4707021 1.0319624 0.02050012 NA NA NA NA
## 5 -1.7417389 -0.4822293 0.10323592 -0.5305046 -0.6957455 -0.3652637 0.007107599
## 6 -2.6822584 -1.1010726 0.16270229 -1.2126830 -1.6702908 -0.7550753 0.054509799
## Log2FC_4 CI.L_4 CI.R_4 vi_4 Log2FC_5 CI.L_5 CI.R_5
## 1 0.2188594 -1.052230 1.4899492 0.4205720 0.8051543 0.1367255 1.4735830
## 2 -1.3755263 -2.162453 -0.5885999 0.1611967 -0.7213490 -1.5714484 0.1287505
## 3 -1.7991026 -2.918939 -0.6792665 0.3264351 -0.8738120 -1.6373061 -0.1103179
## 4 NA NA NA NA NA NA NA
## 5 -0.7991042 -1.457868 -0.1403403 0.1129659 -0.5155929 -0.8606782 -0.1705076
## 6 -1.3554403 -2.288444 -0.4224370 0.2265970 -1.4905464 -2.5565023 -0.4245905
## vi_5 signcon ntimes randomSummary randomCi.lb randomCi.ub randomP
## 1 0.11630493 5 5 1.0333001 0.7882044 1.2783958 1.420312e-16
## 2 0.18811668 -5 5 -1.0649749 -1.3396138 -0.7903361 2.956522e-14
## 3 0.15173972 -3 3 -1.0417876 -1.3210774 -0.7624977 2.653168e-13
## 4 NA 2 2 0.8106154 0.5712566 1.0499741 3.187477e-11
## 5 0.03099851 -5 5 -0.5830624 -0.7212245 -0.4449003 1.324963e-16
## 6 0.29577833 -5 5 -1.2207058 -1.6760435 -0.7653680 1.484852e-07
## het_QE het_QEp het_QM het_QMp error se rank
## 1 4.179945 0.38220032 68.27752 1.420312e-16 FALSE 0.12504883 1
## 2 4.580708 0.33308457 57.76318 2.956522e-14 FALSE 0.14012185 2
## 3 1.980047 0.37156797 53.44957 2.653168e-13 FALSE 0.14249482 3
## 4 0.629179 0.42765661 44.05825 3.187477e-11 FALSE 0.12212181 4
## 5 4.317851 0.36469506 68.41455 1.324963e-16 FALSE 0.07049086 5
## 6 11.099093 0.02547263 27.60901 1.484852e-07 FALSE 0.23231516 6head(meta_degs_rem@metaresult, 3)## Symbol signcon randomSummary randomCi.lb randomCi.ub randomP het_QE
## 1 MXRA5 5 1.033300 0.7882044 1.2783958 1.420312e-16 4.179945
## 2 COL6A6 -5 -1.064975 -1.3396138 -0.7903361 2.956522e-14 4.580708
## 3 CIDEA -3 -1.041788 -1.3210774 -0.7624977 2.653168e-13 1.980047
## het_QEp het_QM het_QMp error rank
## 1 0.3822003 68.27752 1.420312e-16 FALSE 1
## 2 0.3330846 57.76318 2.956522e-14 FALSE 2
## 3 0.3715680 53.44957 2.653168e-13 FALSE 3meta_degs_rem@MetaVolcano
draw_forest(remres=meta_degs_rem,
gene="MMP9",
genecol="Symbol",
foldchangecol="Log2FC",
llcol="CI.L",
rlcol="CI.R",
jobname="MetaVolcano",
outputfolder=".",
draw="HTML")
draw_forest(remres=meta_degs_rem,
gene="COL6A6",
genecol="Symbol",
foldchangecol="Log2FC",
llcol="CI.L",
rlcol="CI.R",
jobname="MetaVolcano",
outputfolder=".",
draw="HTML")
 The REM MetaVolcano also allows users to explore the forest plot of a given gene based on the REM results.
Vote-counting approach
The vote-counting MetaVolcano identifies differential expressed genes (DEG) for each study based on the user-defined p-value and fold change thresholds. It displays the number of differentially expressed and unperturbed genes per study. In addition, it plots the inverse cumulative distribution of the consistently DEG, so the user can identify the number of genes whose expression is perturbed in at least 1 or n studies.
meta_degs_vote <- votecount_mv(diffexp=diffexplist,
pcriteria='pvalue',
foldchangecol='Log2FC',
genenamecol='Symbol',
geneidcol=NULL,
pvalue=0.05,
foldchange=0,
metathr=0.01,
collaps=FALSE,
jobname="MetaVolcano",
outputfolder=".",
draw='HTML')
head(meta_degs_vote@metaresult, 3)## Symbol deg_1 deg_2 deg_3 deg_4 deg_5 ndeg ddeg idx degvcount
## 1 ABCC3 1 1 1 1 1 5 5 25 2.Up-regulated
## 2 ABHD5 -1 -1 -1 -1 -1 5 -5 -25 0.Down-regulated
## 3 ACACB -1 -1 -1 -1 -1 5 -5 -25 0.Down-regulatedmeta_degs_vote@degfreq
The vote-counting MetaVolcano visualizes genes based on the number of studies where genes were identified as differentially expressed and the gene fold change sign consistency. It means that a gene that was differentially expressed in five studies, from which three of them it was downregulated, will get a sign consistency score of 2 + (-3) = -1. Based on user preference, MetaVolcano can highlight the top metathr percentage of consistently perturbed genes.
meta_degs_vote@MetaVolcano
Combining-approach
The combinig MetaVolcano summarizes the fold change of a gene in different studies by the mean or median depending on the user preference. In addition, the combinig MetaVolcano summarizes the gene differential expression p-values using the Fisher method. The combining MetaVolcano can highlight the top metathr percentage of consistently perturbed genes.
meta_degs_comb <- combining_mv(diffexp=diffexplist,
pcriteria='pvalue',
foldchangecol='Log2FC',
genenamecol='Symbol',
geneidcol=NULL,
metafc='Mean',
metathr=0.01,
collaps=TRUE,
jobname="MetaVolcano",
outputfolder=".",
draw='HTML')
head(meta_degs_comb@metaresult, 3)## Symbol metap metafc idx
## 1 MMP9 9.002947e-15 1.9693517 27.66076
## 2 ACVR1C 3.548802e-20 -1.2544105 -24.39818
## 3 ANG 5.674270e-26 -0.9364936 -23.64280meta_degs_comb@MetaVolcano