inst/doc/README.md
In sbcblab/geva: Gene Expression Variation Analysis (GEVA)
Gene Expression Variation Analysis (GEVA)
Itamar José Guimarães Nunes, Murilo Zanini David, Bruno César Feltes, and Marcio Dorn
Introduction
GEVA is a package for the analysis of differential gene expression in
multiple experimental comparisons. It takes into account the
fold-changes and p-values from previous differential expression (DE)
results that use large-scale data (e.g., microarray and RNA-seq) and
evaluates which genes would react in response to the distinct
experiments. This evaluation involves an unique pipeline of statistical
methods, including weighted summarization, quantile detection, cluster
analysis, and ANOVA tests, in order to classify a subset of relevant
genes whose DE is similar or dependent to certain biological factors.
This guide introduces the basic usage of geva package and focuses on
its main features to perform the entire analysis from the input to the
final classification. However, for more detailed specifications
regarding classes, functions, and arguments from geva, please check
the “Reference Guide” available in our GitHub repository. Alternatively,
the local documentation can be accessed by typing ?geva in the R
console.
Before proceeding to the current methodology, it is assumed that the
user already knows how to manipulate datasets and perform DE analyses
using Bioconductor packages or any external tool that is capable to
produce results from DE comparisons. For users with less familiarity
about this subject, please read the tutorials described to the available
R packages for DE analyses, such as limma for microarrays and
DESeq2 for RNA-seq. In addition, some standalone applications employ
the equivalent methods from R to achieve the same results, including
GEAP (for microarrays) and Chipster (for RNA-seq) , both of which
provide a graphical user interface and do not require programming
knowledge.
Installation
This package is available on GitHub and can be installed through the
following command:
BiocManager::install("geva")
Note that this command requires the BiocManager package (installed via
install.packages('BiocManager')). After downloading and installing the
sources, use the following command to load geva from the local package
library:
library(geva)
Data input
The input data is essentialy two or more tables produced by DE analyses
that include logFC and (adjusted) p-value columns in association to the
genes (row names). For microarrays, particularly, the probes may be used
as row names along with a Gene Symbol column, which can be attached to
the final results at the end of the analysis. Moreover, although only
two tables are required for GEVA’s minimal usage, the inclusion of
several columns is strongly recommended to achieve a resonable
statistical precision. Note that experiments can be grouped and analyzed
in multiple contexts at once in GEVA, and likewise in this case, each
group should include several experiments to attain better results from
the statistical tests.
GEVA gives some data input alternatives so that users can provide
objects from the local R environment or from external table files. These
alternatives are described in the sub-sections below, whereas only one
of them is required to accomplish the same desired output.
Alternative 1 – Tab-delimited Text Files
Programs that feature DE analysis usually output a table of DE results
which is exported as a plain text file. By convention, the saved file
should be formatted as one row per line and one tab-delimited value per
column, but other formats may be used as well. For the conventional
format, the geva.read.tables function can be called using default
parameters as demonstrated below:
# Replace the file names below with actual file paths
filenms <- c("cond_A_2h.txt", "cond_B_2h.txt", "cond_C_2h.txt",
"cond_A_4h.txt", "cond_B_4h.txt", "cond_C_4h.txt")
ginput <- geva.read.tables(filenms)
The code above will produce a GEVAInput object, which stores all the
relevant information regarding the input. It reads each file as a table
by calling read.table internally and extracting the columns containing
logFC and adj.P.Val columns, then merging all columns into two
tables (one for logFC values and one for weights).
In addition, the geva.read.tables function has some handful optional
parameters to be considered. For instance, if the dirname parameter is
used instead of filenames, all files inside the directory dirname
that match the pattern given by the files.pattern argument (default is
"\\.txt$" or TXT files) will be included. Other relevant arguments are
col.values (by default, "logFC") and col.pvals (by default,
"adj.P.Val"), used to indicate which columns names are used for
logFC and (adjusted) p-values. Vectors of multiple character
elements can be passed to these arguments if the column names differ
among the table files so that the first matching column is included.
Furhermore, if one wants to append additional columns in the analysis
(e.g., gene names or gene symbols) to associate them to the final
results, the column names can be specified at the col.other argument.
Alternative 2 – Multiple table objects
Table objects, particularly of matrix and data.frame types, can be
used as input to GEVA as long as they include the logFC and p-value
columns. The geva.merge.input function receives two or more table
arguments and extracts their corresponding columns to include in the
final merge. For example, given two data.frame objects defined as
dt1 and dt1 in the global environment, the command for this step
becomes:
# dt1 and dt2 are examples of input data.frames
# containing logFC and adj.P.Val columns
ginput <- geva.merge.input(dt1, dt2)
The code above will produce a GEVAInput object, which stores all the
relevant information regarding the input. Arguments are passed
individually and can also be named to define the columns in the final
merge (e.g., cond1=dt1, cond2=dt2 to append the extracted columns as
"cond1" and "cond2"). Note that some arguments from
geva.read.tables, including col.values, col.pvals, and
col.other, have the same principle as in geva.merge.input (see
Alternative 1).
Alternative 3 – Results from limma
If the DE analysis is being performed from a specific R package such as
limma, the results can be converted to a matrix or data.frame and
passed as arguments to geva.merge.input as demonstrated in the
previous section (see Alternative 2). For example, if limma was used
to produce two MArrayLM objects (i.e., DE results using linear model
fit), these can be converted to data.frame using limma::topTable,
then passed to geva.merge.input as demonstrated below:
# malm1 and malm2 are MArrayLM objects produced by
# limma (e.g., using eBayes)
dt1 <- topTable(malm1, number=999999, sort.by="none")
dt2 <- topTable(malm2, number=999999, sort.by="none")
ginput <- geva.merge.input(dt1, dt2)
The code above will produce a GEVAInput object, which stores all the
relevant information regarding the input. Since both dt1 and dt2
already include "logFC" and "adj.P.Value" columns,
geva.merge.input can be called using the defaults parameters.
Alternative 4 – Ideal input data (for tests only)
Be it due the abscence of experimental data or merely for didatical
reasons, there may be some situations where the features in this package
have to be immediately accessed and tested without needing to provide
any real data, since two or more DE analyses must be performed before
using GEVA. In this sense, the geva.ideal.example function can be used
to generate a random input that simulates real processed inputs by GEVA.
The function is called as follows:
# Generates a random GEVAInput with 10000 probes
# and 6 columns
ginput <- geva.ideal.example()
The code above will generate a GEVAInput object with random values
within a normal distribution and some random outliers to simulate the
relevant results. In addition, all columns are grouped into experimental
condition groups (factors) so that factor-dependent and
factor-specific results could be produced by the end of the analysis
test. Note that although the output is essentially “random”, the same
result can be reproduced by using set.seed before
geva.ideal.example.
Input data post-processing (optional)
Considering that the final results will strongly depend on the input
values in the concatenated tables, some tweaks in the obtained
GEVAInput can be done to improve them in terms of statistics and
presentation. Some features implemented in GEVA that allow this kind of
post-processing of the GEVAInput objects are presented over the
following sub-sections.
Numeric table correcting
First off, one may want to eliminate primary sources of errors from the
numeric tables before proceeding to the next steps. The calculations
become prone to bias when missing values (NA) or infinite numbers
(Inf) are present, so except in rare cases where their inclusion is
intentional, removing them is a reasonable choice.
In this sense, the geva.input.correct function will remove all missing
(NA), not-a-number (NaN), and infinite (Inf or -Inf) values from
GEVAInput upon calling the following command:
# Removes the rows containing missing and infinite values
ginput <- geva.input.correct(ginput)
The validation is only applied to the numeric tables in GEVAInput
(i.e., @values and @weights slots). As a result, if any invalid
values were found, their rows are removed. However, there is an
exceptional case where one column is entirely made of invalid values
which would cause all rows to be marked as invalid, so
geva.input.correct removes such columns in advance to prevent the
exclusion of the entire table.
Filtering values below statistical significance
The GEVAInput stores a table of transformed p-values as weights
(@weights slot, called by inputweights function) employed in some
calculations during the summarization step (discussed in the next
section). While the inclusion of weights is used to minimize statistical
errors, it also follows the assumption that all rows have at least one
significant p-value. In this sense, the geva.input.filter function can
be used to remove rows whose p-values are all above a certain threshold
(e.g., 0.05), as demonstrated below:
# Removes the rows that are entirely composed by
# insignificant values
ginput <- geva.input.filter(ginput, p.value.cutoff = 0.05)
The correction above is applied using a threshold of (\alpha < 0.05)
for (corrected) p-values. Just like any other statistical procedure, the
value of 0.05 given to the p.value.cutoff argument is arbitrary and it
is upon the user’s choice to define the best delimiter of significance.
Renaming the row names
Although large-scale experimental data is usually targeted to the
context of each gene, it is particularly common in microarrays to use
multiple probes that detect the expression levels for one or more genes.
If one desires to use gene names as primary row identifiers instead of
probes, these genes must replace the probes names accordingly. However,
multiple genes per probe become duplicates, so one of them must chosen
to provide unique identifiers for row names. In this sense, the
geva.input.rename.rows function is used to perform the renaming while
also solving such duplicates as demonstrated below:
# Replaces the row names with the "Symbol" column while
# selecting the most significant duplicates
ginput <- geva.input.rename.rows(ginput,
attr.column = "Symbol")
In the example above, the ginput has an additional column called
"Symbol" (accessed by featureTable(ginput)$Symbol) which is used to
replace the row names, but the attr.column argument could also be a
character vector with the same length of the number of rows. By default,
the above code will select the duplicates which have the least p-values
(i.e., lowest error probability), which is also specified by applying
the dupl.rm.method="least.p.vals" parameter. Alternatively, the
parameter dupl.rm.method="order" can be used to select the duplicated
value that appears first in the row order.
SV Analyses
By concluding the input step, a GEVAInput object that stores logFC
values and weights becomes available in the current session. The next
step will be to calculate the summarization and variation (SV) from
the concatenated input data to produce the SV points, which are used in
intermediate steps before the final classification.
Summarization
The geva.summarize function takes a GEVAInput object and performs
the summarization, as demonstrated below:
# Summarizes ginput to find the SV points
gsummary <- geva.summarize(ginput)
The code above uses the default parameters for summary.method and
variation.method ( and , respectively) but other methods are available
such as "median" and "mad" (median absolute deviation, or MAD). In
this context, they could be specified as follow:
# Summarizes ginput using median and MAD
gsummary <- geva.summarize(ginput,
summary.method = "median",
variation.method = "mad")
In addition, all the summarization methods specified in summary.method
and variation.method are implemented to take weights into account
(except if not available or when weights are equivalent).
As a result, geva.summarize returns a GEVASummary object storing the
table of S and V values. From this point, all objects defined by
intermediate steps can be plotted as a SV-plot, a type of scatter plot
where each point (called SV-point) represents a gene’s central logFC
value (S) and logFC variation (V). For instance, a plot can be
produced by calling the plot function on a GEVASummary object:
plot(gsummary)
Delimitation by quantiles
After obtaining the GEVASummary object, the next step will be
calculating the quantiles for every SV-point. That can be done by
calling the geva.quantiles function as shown below:
# Calculates the quantiles from a GEVASummary object
gquants <- geva.quantiles(gsummary)
The code above produces a GEVAQuantiles object which stores the
relevant partitions where the SV-points belong to. These partitions can
be viewed by calling plot on the produced object:
plot(gquants)
By default, the quantile detection is performed automatically using the
parameter (for more methods, call ?geva.quantiles). However, the
quantile delimiters can also be specified in the initial.thresholds
argument like the following example:
# Calculates the quantiles from a GEVASummary object
# using custom delimiters
gquants <- geva.quantiles(gsummary,
initial.thresholds = c(S=1, V=0.5))
In this second example, thresholds of 1 and 0.5 were defined for S
and V axes. As it can be noted from the SV-plot below, the results are
purposely exaggerated and may not represent a good separation between
relevant points, but this option is particularly useful to fine-tuning
the quantile delimiters in situations where the automatic methods did
not present a satisfatory outcome.
Note that the quantile detection does not define an absolute cutoff, but
partitionizes the SV space into estimated regions containing qualitative
classifications for the SV points. These classifications may change
after combining the GEVAQuantiles with the results from the next
steps.
Clustering
In this step, a cluster analysis is applied to separate relevant points
from the agglomeration of non-differentially expressed genes. Such
agglomeration is mostly proeminent at the bottom-center region of a SV
plot and essentially portraits the least relevant portion of the
results.
Basic usage of the wrapper function
The geva.cluster function is the top-level function for clusters
analysis and acts as a wrapper for more specific functions used to group
SV points. The inner function is specified by the cluster.method
argument with one of the following parameters: (i) "hierarchical",
calls the geva.hcluster function for hierarchical clustering; (ii)
"density", calls the geva.dcluster function for density-based
clustering; and (iii) "quantiles", calls the geva.quantiles function
shown in the previous section. Likewise, optional parameters from the
top funtion are passed to these calls.
In this section, only hierarchical and density-based clustering methods
are going to be discussed. Both methods use the resolution argument, a
single numeric between 0 and 1 that defines the ratio of output
clusters. If the resolution is 0.0 (zero), the least number of
clusters is assigned (i.e., usually one or two), while if 1.0 then
the maximum amount of clusters is assigned (i.e., aproximately one
cluster per point for hierarchical clustering). For example, to apply
geva.cluster using hierarchical clustering at 30% of the resolution,
the function is called as follows:
# Applies cluster analysis (30% resolution)
gcluster <- geva.cluster(gsummary,
cluster.method="hierarchical",
resolution=0.3)
The returned cluster data can be plotted using the generic plot
function:
plot(gcluster)
Combining clusters with summarized data (Optional)
Apart from its usage as a wrapper, the geva.cluster function can also
concatenate the summarized and grouped data into a single object by
setting grouped.return=TRUE in the arguments. With this setup, the
function will return a GEVAGroupedSummary object, which is a
GEVASummary that includes the list of group sets (GEVACluster or
GEVAQuantile objects). The code below illustrates this specific use
case:
# Applies cluster analysis with default parameters and
# returns a GEVAGroupedSummary
ggroupedsummary <- geva.cluster(gsummary,
grouped.return = TRUE)
Alternatively, multiple group sets (clusters and quantiles) can be
combined directly to the summarized data by appending each of them with
groupsets<-, which also converts the GEVASummary to a
GEVAGroupedSummary object. For example, assuming that gquants and
gcluster are output values from the previous quantiles
(geva.quantiles) and clustering (geva.cluster) steps, respectively,
the code would be:
# Makes a safe copy of the summary data
ggroupedsummary <- gsummary
# Appends the quantiles data
groupsets(ggroupedsummary) <- gquants
# Appends the clustered data
groupsets(ggroupedsummary) <- gcluster
# Draws a SV plot with grouped highlights (optional)
plot(ggroupedsummary)
Attaining and accessing the results
After obtaining the quantiles and clusters from the summarized data in
the previous step, now the entire data can be taken together to prospect
the final classifications for each gene. This section presents the final
steps to obtain the results table and some basic method to access it.
Final concatenation and factor analysis
In this final analysis step, the geva.finalize function takes a
GEVASummary object as argument in addition to the other values
returned from the intermediate steps, including the GEVAQuantiles and
GEVACluster objects. Alternatively, a GEVAGroupedSummary object
containing these intermediate results can be provided alone. The
function will correct the quantiles based on the clustered points and
return a classification that fits better both group assignments.
Furthermore, if factors (groups of experimental conditions) were defined
for the input columns, geva.finalize will also look for DE variations
according to these factors, thereby unlocking two additional possible
classifications ("factor-dependent" and "factor-specific"). The
possible use cases are discussed in the following sub-sections:
Alternative 1 – Without factors
If factors were not included, no additional steps are required. The
function call is done by passing the GEVASummary, GEVAQuantiles and
GEVACluster from previous steps:
# Calculates the final classifications based on the
# intermediate results from previous steps
gresults <- geva.finalize(gsummary, gquants, gcluster)
Or, if a GEVAGroupedSummary object is provided:
# Calculates the final classifications based on the
# intermediate results from previous steps
gresults <- geva.finalize(ggroupedsummary)
Note that, without factors, the only relevant classification is
"similar" (i.e., genes with similar logFC values among all
experiments).
Alternative 2 – With factors
Factors can be accessed and assigned to a GEVAInput object using
factors and factors<-, respectively, and both accessors are valid
for GEVASummary as well. The factors being set must be a factor or
character vector whose length is equivalent to the number of columns,
and it must contain at least two values per level to be considered since
the factor analysis is based on ANOVA.
For instance, considering a GEVASummary object that stores a
GEVAInput with 6 columns (experimental results), if one wants to
separate these columns into 3 factors (‘g1’, ‘g2’, and ‘g3’), the
following code could be applied:
# Assigning factors to an example input with 6 columns
# Example with GEVAInput
factors(ginput) <- c('g1', 'g1', 'g2', 'g2', 'g3', 'g3')
# Example with GEVAInput (using factor class)
factors(ginput) <- factor(c('g1', 'g1', 'g2',
'g2', 'g3', 'g3'))
# Example with GEVASummary
factors(gsummary) <- c('g1', 'g1', 'g2', 'g2', 'g3', 'g3')
By including factors in the current analysis, some optional arguments
related to the factor analysis become available in geva.finalize. The
p.value, for instance, determines the significance cutoff employed in
ANOVA tests (by default, this value is 0.05 for (\alpha < 0.05)). In
this case, the function call becomes:
# Calculates the final classifications based on the
# intermediate results from previous steps
gresults <- geva.finalize(gsummary, gquants, gcluster,
p.value=0.05)
Or, if a GEVAGroupedSummary object is provided:
# Calculates the final classifications based on the
# intermediate results from previous steps
gresults <- geva.finalize(ggroupedsummary, p.value=0.05)
The results can be plotted into a SV plot similarly as in the previous
steps, but now only the relevant points will be colored while the rest
are painted in black or gray:
plot(gresults)
Accessing and extracting the results
The returned GEVAResults object from geva.finalize represents the
concatenation of all previous steps in addition to the results table
and, if applicable, the intermediate steps from the factor analysis. The
results table stores the final gene classifications, including the
relevant ("similar", "factor-dependent", and "factor-specific")
and irrelevant ("sparse" and "basal") ones. Each classification can
be briefly described as follows:
basal: Genes with similar but mild logFC that approximates to
zero. Note that despite this name they not necessarily represent
basal levels of gene expression, especially if the control group
from DE analysis is not under normal conditions;sparse: Genes with high logFC variation but lacking any
relationship to the experimental conditions or the factors;similar: Genes with relevant logFC (far from zero) and low
logFC variance;factor-dependent: Genes with low logFC variance within the
specified factors, but high variance between diferent factors;factor-specific: Genes with low logFC variance within one
specific factor.
The function results.table can be used to return the table of final
gene classifications:
tail(results.table(gresults), 10)
| | classification | specific.factor |
| :----------- | :-------------- | :-------------- |
| probe_9991 | basal | NA |
| probe_9992 | basal | NA |
| probe_9993 | basal | NA |
| probe_9994 | basal | NA |
| probe_9995 | basal | NA |
| probe_9996 | basal | NA |
| probe_9997 | basal | NA |
| probe_9998 | factor-specific | Cond_2 |
| probe_9999 | basal | NA |
| probe_10000 | basal | NA |
On the other hand, the top.genes function may be a rather practical
way to return the most relevant results. It extracts by default the
"similar", "factor-dependent", and "factor-specific" results, and
can attach additional columns (e.g., gene symbols) specified by the
add.cols arguments. The code below shows an usage example of
top.genes:
# Extracts the top genes only
dtgens <- top.genes(gresults)
# Extracts the top genes and appends the "Symbol" column
dtgens <- top.genes(gresults, add.cols = "Symbol")
# Prints the last lines of the top genes table (optional)
print(tail(dtgens, 10))
| | Symbol | classification | specific.factor |
| :---------- | :---------- | :--------------- | :-------------- |
| probe_8487 | GENE_K8487 | factor-dependent | NA |
| probe_8740 | GENE_D8740 | factor-dependent | NA |
| probe_8823 | GENE_I8823 | factor-specific | Cond_1 |
| probe_9136 | GENE_J9136 | similar | NA |
| probe_9312 | GENE_D9312 | factor-dependent | NA |
| probe_9495 | GENE_E9495 | factor-dependent | NA |
| probe_9601 | GENE_G9601 | factor-specific | Cond_3 |
| probe_9758 | GENE_H9758 | factor-specific | Cond_3 |
| probe_9893 | GENE_M9893 | factor-dependent | NA |
| probe_9998 | GENE_N9998 | factor-specific | Cond_2 |
The resulting table can then be exported using functions such has
write.table from the R base package.
Shortcut function and reanalysis
The geva.quick function accepts a GEVAInput object and performs all
intermediate functions from the summarization to the final
concatenation. Optional (...) arguments are passed to the internal
calls to geva.summarize, geva.quantiles, geva.cluster and
geva.finalize, ultimately returning a GEVAResults object. The basic
usage is described as follows:
# Generates a random GEVAInput example
ginput <- geva.ideal.example()
# Performs all intermediate steps with geva.quick
# The resolution is used by the call to geva.cluster
gresults <- geva.quick(ginput, resolution=0.25)
## > Found 4 clusters and 31 significant genes
gresults <- geva.quick(ginput, resolution=0.4)
## > Found 16 clusters and 116 significant genes
This function can be applied to a GEVAResults object as well to
restore the parameters that produced this result, whereas optional
(...) arguments can overwrite them:
# Generates a random GEVAInput example
ginput <- geva.ideal.example()
# Performs all intermediate steps with geva.quick
# The summary.method is used by the call to geva.summarize
gresults <- geva.quick(ginput, summary.method='mean')
## > Found 60 significant genes
gresults <- geva.quick(gresults, summary.method='median')
## > Found 95 significant genes
In the example above, the entire analysis was redone using the
overwritten summary.method argument. Therefore, by following this
pattern, users can tweak different parameters depending on their
statistical choice regarding the current biological context.
sbcblab/geva documentation built on March 15, 2021, 10:08 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Gene Expression Variation Analysis (GEVA)
Itamar José Guimarães Nunes, Murilo Zanini David, Bruno César Feltes, and Marcio Dorn
Introduction
GEVA is a package for the analysis of differential gene expression in
multiple experimental comparisons. It takes into account the
fold-changes and p-values from previous differential expression (DE)
results that use large-scale data (e.g., microarray and RNA-seq) and
evaluates which genes would react in response to the distinct
experiments. This evaluation involves an unique pipeline of statistical
methods, including weighted summarization, quantile detection, cluster
analysis, and ANOVA tests, in order to classify a subset of relevant
genes whose DE is similar or dependent to certain biological factors.
This guide introduces the basic usage of geva package and focuses on
its main features to perform the entire analysis from the input to the
final classification. However, for more detailed specifications
regarding classes, functions, and arguments from geva, please check
the “Reference Guide” available in our GitHub repository. Alternatively,
the local documentation can be accessed by typing ?geva in the R
console.
Before proceeding to the current methodology, it is assumed that the user already knows how to manipulate datasets and perform DE analyses using Bioconductor packages or any external tool that is capable to produce results from DE comparisons. For users with less familiarity about this subject, please read the tutorials described to the available R packages for DE analyses, such as limma for microarrays and DESeq2 for RNA-seq. In addition, some standalone applications employ the equivalent methods from R to achieve the same results, including GEAP (for microarrays) and Chipster (for RNA-seq) , both of which provide a graphical user interface and do not require programming knowledge.
Installation
This package is available on GitHub and can be installed through the following command:
BiocManager::install("geva")
Note that this command requires the BiocManager package (installed via
install.packages('BiocManager')). After downloading and installing the
sources, use the following command to load geva from the local package
library:
library(geva)
Data input
The input data is essentialy two or more tables produced by DE analyses that include logFC and (adjusted) p-value columns in association to the genes (row names). For microarrays, particularly, the probes may be used as row names along with a Gene Symbol column, which can be attached to the final results at the end of the analysis. Moreover, although only two tables are required for GEVA’s minimal usage, the inclusion of several columns is strongly recommended to achieve a resonable statistical precision. Note that experiments can be grouped and analyzed in multiple contexts at once in GEVA, and likewise in this case, each group should include several experiments to attain better results from the statistical tests.
GEVA gives some data input alternatives so that users can provide objects from the local R environment or from external table files. These alternatives are described in the sub-sections below, whereas only one of them is required to accomplish the same desired output.
Alternative 1 – Tab-delimited Text Files
Programs that feature DE analysis usually output a table of DE results
which is exported as a plain text file. By convention, the saved file
should be formatted as one row per line and one tab-delimited value per
column, but other formats may be used as well. For the conventional
format, the geva.read.tables function can be called using default
parameters as demonstrated below:
# Replace the file names below with actual file paths
filenms <- c("cond_A_2h.txt", "cond_B_2h.txt", "cond_C_2h.txt",
"cond_A_4h.txt", "cond_B_4h.txt", "cond_C_4h.txt")
ginput <- geva.read.tables(filenms)
The code above will produce a GEVAInput object, which stores all the
relevant information regarding the input. It reads each file as a table
by calling read.table internally and extracting the columns containing
logFC and adj.P.Val columns, then merging all columns into two
tables (one for logFC values and one for weights).
In addition, the geva.read.tables function has some handful optional
parameters to be considered. For instance, if the dirname parameter is
used instead of filenames, all files inside the directory dirname
that match the pattern given by the files.pattern argument (default is
"\\.txt$" or TXT files) will be included. Other relevant arguments are
col.values (by default, "logFC") and col.pvals (by default,
"adj.P.Val"), used to indicate which columns names are used for
logFC and (adjusted) p-values. Vectors of multiple character
elements can be passed to these arguments if the column names differ
among the table files so that the first matching column is included.
Furhermore, if one wants to append additional columns in the analysis
(e.g., gene names or gene symbols) to associate them to the final
results, the column names can be specified at the col.other argument.
Alternative 2 – Multiple table objects
Table objects, particularly of matrix and data.frame types, can be
used as input to GEVA as long as they include the logFC and p-value
columns. The geva.merge.input function receives two or more table
arguments and extracts their corresponding columns to include in the
final merge. For example, given two data.frame objects defined as
dt1 and dt1 in the global environment, the command for this step
becomes:
# dt1 and dt2 are examples of input data.frames
# containing logFC and adj.P.Val columns
ginput <- geva.merge.input(dt1, dt2)
The code above will produce a GEVAInput object, which stores all the
relevant information regarding the input. Arguments are passed
individually and can also be named to define the columns in the final
merge (e.g., cond1=dt1, cond2=dt2 to append the extracted columns as
"cond1" and "cond2"). Note that some arguments from
geva.read.tables, including col.values, col.pvals, and
col.other, have the same principle as in geva.merge.input (see
Alternative 1).
Alternative 3 – Results from limma
If the DE analysis is being performed from a specific R package such as
limma, the results can be converted to a matrix or data.frame and
passed as arguments to geva.merge.input as demonstrated in the
previous section (see Alternative 2). For example, if limma was used
to produce two MArrayLM objects (i.e., DE results using linear model
fit), these can be converted to data.frame using limma::topTable,
then passed to geva.merge.input as demonstrated below:
# malm1 and malm2 are MArrayLM objects produced by
# limma (e.g., using eBayes)
dt1 <- topTable(malm1, number=999999, sort.by="none")
dt2 <- topTable(malm2, number=999999, sort.by="none")
ginput <- geva.merge.input(dt1, dt2)
The code above will produce a GEVAInput object, which stores all the
relevant information regarding the input. Since both dt1 and dt2
already include "logFC" and "adj.P.Value" columns,
geva.merge.input can be called using the defaults parameters.
Alternative 4 – Ideal input data (for tests only)
Be it due the abscence of experimental data or merely for didatical
reasons, there may be some situations where the features in this package
have to be immediately accessed and tested without needing to provide
any real data, since two or more DE analyses must be performed before
using GEVA. In this sense, the geva.ideal.example function can be used
to generate a random input that simulates real processed inputs by GEVA.
The function is called as follows:
# Generates a random GEVAInput with 10000 probes
# and 6 columns
ginput <- geva.ideal.example()
The code above will generate a GEVAInput object with random values
within a normal distribution and some random outliers to simulate the
relevant results. In addition, all columns are grouped into experimental
condition groups (factors) so that factor-dependent and
factor-specific results could be produced by the end of the analysis
test. Note that although the output is essentially “random”, the same
result can be reproduced by using set.seed before
geva.ideal.example.
Input data post-processing (optional)
Considering that the final results will strongly depend on the input
values in the concatenated tables, some tweaks in the obtained
GEVAInput can be done to improve them in terms of statistics and
presentation. Some features implemented in GEVA that allow this kind of
post-processing of the GEVAInput objects are presented over the
following sub-sections.
Numeric table correcting
First off, one may want to eliminate primary sources of errors from the
numeric tables before proceeding to the next steps. The calculations
become prone to bias when missing values (NA) or infinite numbers
(Inf) are present, so except in rare cases where their inclusion is
intentional, removing them is a reasonable choice.
In this sense, the geva.input.correct function will remove all missing
(NA), not-a-number (NaN), and infinite (Inf or -Inf) values from
GEVAInput upon calling the following command:
# Removes the rows containing missing and infinite values
ginput <- geva.input.correct(ginput)
The validation is only applied to the numeric tables in GEVAInput
(i.e., @values and @weights slots). As a result, if any invalid
values were found, their rows are removed. However, there is an
exceptional case where one column is entirely made of invalid values
which would cause all rows to be marked as invalid, so
geva.input.correct removes such columns in advance to prevent the
exclusion of the entire table.
Filtering values below statistical significance
The GEVAInput stores a table of transformed p-values as weights
(@weights slot, called by inputweights function) employed in some
calculations during the summarization step (discussed in the next
section). While the inclusion of weights is used to minimize statistical
errors, it also follows the assumption that all rows have at least one
significant p-value. In this sense, the geva.input.filter function can
be used to remove rows whose p-values are all above a certain threshold
(e.g., 0.05), as demonstrated below:
# Removes the rows that are entirely composed by
# insignificant values
ginput <- geva.input.filter(ginput, p.value.cutoff = 0.05)
The correction above is applied using a threshold of (\alpha < 0.05)
for (corrected) p-values. Just like any other statistical procedure, the
value of 0.05 given to the p.value.cutoff argument is arbitrary and it
is upon the user’s choice to define the best delimiter of significance.
Renaming the row names
Although large-scale experimental data is usually targeted to the
context of each gene, it is particularly common in microarrays to use
multiple probes that detect the expression levels for one or more genes.
If one desires to use gene names as primary row identifiers instead of
probes, these genes must replace the probes names accordingly. However,
multiple genes per probe become duplicates, so one of them must chosen
to provide unique identifiers for row names. In this sense, the
geva.input.rename.rows function is used to perform the renaming while
also solving such duplicates as demonstrated below:
# Replaces the row names with the "Symbol" column while
# selecting the most significant duplicates
ginput <- geva.input.rename.rows(ginput,
attr.column = "Symbol")
In the example above, the ginput has an additional column called
"Symbol" (accessed by featureTable(ginput)$Symbol) which is used to
replace the row names, but the attr.column argument could also be a
character vector with the same length of the number of rows. By default,
the above code will select the duplicates which have the least p-values
(i.e., lowest error probability), which is also specified by applying
the dupl.rm.method="least.p.vals" parameter. Alternatively, the
parameter dupl.rm.method="order" can be used to select the duplicated
value that appears first in the row order.
SV Analyses
By concluding the input step, a GEVAInput object that stores logFC
values and weights becomes available in the current session. The next
step will be to calculate the summarization and variation (SV) from
the concatenated input data to produce the SV points, which are used in
intermediate steps before the final classification.
Summarization
The geva.summarize function takes a GEVAInput object and performs
the summarization, as demonstrated below:
# Summarizes ginput to find the SV points
gsummary <- geva.summarize(ginput)
The code above uses the default parameters for summary.method and
variation.method ( and , respectively) but other methods are available
such as "median" and "mad" (median absolute deviation, or MAD). In
this context, they could be specified as follow:
# Summarizes ginput using median and MAD
gsummary <- geva.summarize(ginput,
summary.method = "median",
variation.method = "mad")
In addition, all the summarization methods specified in summary.method
and variation.method are implemented to take weights into account
(except if not available or when weights are equivalent).
As a result, geva.summarize returns a GEVASummary object storing the
table of S and V values. From this point, all objects defined by
intermediate steps can be plotted as a SV-plot, a type of scatter plot
where each point (called SV-point) represents a gene’s central logFC
value (S) and logFC variation (V). For instance, a plot can be
produced by calling the plot function on a GEVASummary object:
plot(gsummary)
Delimitation by quantiles
After obtaining the GEVASummary object, the next step will be
calculating the quantiles for every SV-point. That can be done by
calling the geva.quantiles function as shown below:
# Calculates the quantiles from a GEVASummary object
gquants <- geva.quantiles(gsummary)
The code above produces a GEVAQuantiles object which stores the
relevant partitions where the SV-points belong to. These partitions can
be viewed by calling plot on the produced object:
plot(gquants)
By default, the quantile detection is performed automatically using the
parameter (for more methods, call ?geva.quantiles). However, the
quantile delimiters can also be specified in the initial.thresholds
argument like the following example:
# Calculates the quantiles from a GEVASummary object
# using custom delimiters
gquants <- geva.quantiles(gsummary,
initial.thresholds = c(S=1, V=0.5))
In this second example, thresholds of 1 and 0.5 were defined for S
and V axes. As it can be noted from the SV-plot below, the results are
purposely exaggerated and may not represent a good separation between
relevant points, but this option is particularly useful to fine-tuning
the quantile delimiters in situations where the automatic methods did
not present a satisfatory outcome.
Note that the quantile detection does not define an absolute cutoff, but
partitionizes the SV space into estimated regions containing qualitative
classifications for the SV points. These classifications may change
after combining the GEVAQuantiles with the results from the next
steps.
Clustering
In this step, a cluster analysis is applied to separate relevant points from the agglomeration of non-differentially expressed genes. Such agglomeration is mostly proeminent at the bottom-center region of a SV plot and essentially portraits the least relevant portion of the results.
Basic usage of the wrapper function
The geva.cluster function is the top-level function for clusters
analysis and acts as a wrapper for more specific functions used to group
SV points. The inner function is specified by the cluster.method
argument with one of the following parameters: (i) "hierarchical",
calls the geva.hcluster function for hierarchical clustering; (ii)
"density", calls the geva.dcluster function for density-based
clustering; and (iii) "quantiles", calls the geva.quantiles function
shown in the previous section. Likewise, optional parameters from the
top funtion are passed to these calls.
In this section, only hierarchical and density-based clustering methods
are going to be discussed. Both methods use the resolution argument, a
single numeric between 0 and 1 that defines the ratio of output
clusters. If the resolution is 0.0 (zero), the least number of
clusters is assigned (i.e., usually one or two), while if 1.0 then
the maximum amount of clusters is assigned (i.e., aproximately one
cluster per point for hierarchical clustering). For example, to apply
geva.cluster using hierarchical clustering at 30% of the resolution,
the function is called as follows:
# Applies cluster analysis (30% resolution)
gcluster <- geva.cluster(gsummary,
cluster.method="hierarchical",
resolution=0.3)
The returned cluster data can be plotted using the generic plot
function:
plot(gcluster)
Combining clusters with summarized data (Optional)
Apart from its usage as a wrapper, the geva.cluster function can also
concatenate the summarized and grouped data into a single object by
setting grouped.return=TRUE in the arguments. With this setup, the
function will return a GEVAGroupedSummary object, which is a
GEVASummary that includes the list of group sets (GEVACluster or
GEVAQuantile objects). The code below illustrates this specific use
case:
# Applies cluster analysis with default parameters and
# returns a GEVAGroupedSummary
ggroupedsummary <- geva.cluster(gsummary,
grouped.return = TRUE)
Alternatively, multiple group sets (clusters and quantiles) can be
combined directly to the summarized data by appending each of them with
groupsets<-, which also converts the GEVASummary to a
GEVAGroupedSummary object. For example, assuming that gquants and
gcluster are output values from the previous quantiles
(geva.quantiles) and clustering (geva.cluster) steps, respectively,
the code would be:
# Makes a safe copy of the summary data
ggroupedsummary <- gsummary
# Appends the quantiles data
groupsets(ggroupedsummary) <- gquants
# Appends the clustered data
groupsets(ggroupedsummary) <- gcluster
# Draws a SV plot with grouped highlights (optional)
plot(ggroupedsummary)
Attaining and accessing the results
After obtaining the quantiles and clusters from the summarized data in the previous step, now the entire data can be taken together to prospect the final classifications for each gene. This section presents the final steps to obtain the results table and some basic method to access it.
Final concatenation and factor analysis
In this final analysis step, the geva.finalize function takes a
GEVASummary object as argument in addition to the other values
returned from the intermediate steps, including the GEVAQuantiles and
GEVACluster objects. Alternatively, a GEVAGroupedSummary object
containing these intermediate results can be provided alone. The
function will correct the quantiles based on the clustered points and
return a classification that fits better both group assignments.
Furthermore, if factors (groups of experimental conditions) were defined
for the input columns, geva.finalize will also look for DE variations
according to these factors, thereby unlocking two additional possible
classifications ("factor-dependent" and "factor-specific"). The
possible use cases are discussed in the following sub-sections:
Alternative 1 – Without factors
If factors were not included, no additional steps are required. The
function call is done by passing the GEVASummary, GEVAQuantiles and
GEVACluster from previous steps:
# Calculates the final classifications based on the
# intermediate results from previous steps
gresults <- geva.finalize(gsummary, gquants, gcluster)
Or, if a GEVAGroupedSummary object is provided:
# Calculates the final classifications based on the
# intermediate results from previous steps
gresults <- geva.finalize(ggroupedsummary)
Note that, without factors, the only relevant classification is
"similar" (i.e., genes with similar logFC values among all
experiments).
Alternative 2 – With factors
Factors can be accessed and assigned to a GEVAInput object using
factors and factors<-, respectively, and both accessors are valid
for GEVASummary as well. The factors being set must be a factor or
character vector whose length is equivalent to the number of columns,
and it must contain at least two values per level to be considered since
the factor analysis is based on ANOVA.
For instance, considering a GEVASummary object that stores a
GEVAInput with 6 columns (experimental results), if one wants to
separate these columns into 3 factors (‘g1’, ‘g2’, and ‘g3’), the
following code could be applied:
# Assigning factors to an example input with 6 columns
# Example with GEVAInput
factors(ginput) <- c('g1', 'g1', 'g2', 'g2', 'g3', 'g3')
# Example with GEVAInput (using factor class)
factors(ginput) <- factor(c('g1', 'g1', 'g2',
'g2', 'g3', 'g3'))
# Example with GEVASummary
factors(gsummary) <- c('g1', 'g1', 'g2', 'g2', 'g3', 'g3')
By including factors in the current analysis, some optional arguments
related to the factor analysis become available in geva.finalize. The
p.value, for instance, determines the significance cutoff employed in
ANOVA tests (by default, this value is 0.05 for (\alpha < 0.05)). In
this case, the function call becomes:
# Calculates the final classifications based on the
# intermediate results from previous steps
gresults <- geva.finalize(gsummary, gquants, gcluster,
p.value=0.05)
Or, if a GEVAGroupedSummary object is provided:
# Calculates the final classifications based on the
# intermediate results from previous steps
gresults <- geva.finalize(ggroupedsummary, p.value=0.05)
The results can be plotted into a SV plot similarly as in the previous steps, but now only the relevant points will be colored while the rest are painted in black or gray:
plot(gresults)
Accessing and extracting the results
The returned GEVAResults object from geva.finalize represents the
concatenation of all previous steps in addition to the results table
and, if applicable, the intermediate steps from the factor analysis. The
results table stores the final gene classifications, including the
relevant ("similar", "factor-dependent", and "factor-specific")
and irrelevant ("sparse" and "basal") ones. Each classification can
be briefly described as follows:
basal: Genes with similar but mild logFC that approximates to zero. Note that despite this name they not necessarily represent basal levels of gene expression, especially if the control group from DE analysis is not under normal conditions;sparse: Genes with high logFC variation but lacking any relationship to the experimental conditions or the factors;similar: Genes with relevant logFC (far from zero) and low logFC variance;factor-dependent: Genes with low logFC variance within the specified factors, but high variance between diferent factors;factor-specific: Genes with low logFC variance within one specific factor.
The function results.table can be used to return the table of final
gene classifications:
tail(results.table(gresults), 10)
| | classification | specific.factor | | :----------- | :-------------- | :-------------- | | probe_9991 | basal | NA | | probe_9992 | basal | NA | | probe_9993 | basal | NA | | probe_9994 | basal | NA | | probe_9995 | basal | NA | | probe_9996 | basal | NA | | probe_9997 | basal | NA | | probe_9998 | factor-specific | Cond_2 | | probe_9999 | basal | NA | | probe_10000 | basal | NA |
On the other hand, the top.genes function may be a rather practical
way to return the most relevant results. It extracts by default the
"similar", "factor-dependent", and "factor-specific" results, and
can attach additional columns (e.g., gene symbols) specified by the
add.cols arguments. The code below shows an usage example of
top.genes:
# Extracts the top genes only
dtgens <- top.genes(gresults)
# Extracts the top genes and appends the "Symbol" column
dtgens <- top.genes(gresults, add.cols = "Symbol")
# Prints the last lines of the top genes table (optional)
print(tail(dtgens, 10))
| | Symbol | classification | specific.factor | | :---------- | :---------- | :--------------- | :-------------- | | probe_8487 | GENE_K8487 | factor-dependent | NA | | probe_8740 | GENE_D8740 | factor-dependent | NA | | probe_8823 | GENE_I8823 | factor-specific | Cond_1 | | probe_9136 | GENE_J9136 | similar | NA | | probe_9312 | GENE_D9312 | factor-dependent | NA | | probe_9495 | GENE_E9495 | factor-dependent | NA | | probe_9601 | GENE_G9601 | factor-specific | Cond_3 | | probe_9758 | GENE_H9758 | factor-specific | Cond_3 | | probe_9893 | GENE_M9893 | factor-dependent | NA | | probe_9998 | GENE_N9998 | factor-specific | Cond_2 |
The resulting table can then be exported using functions such has
write.table from the R base package.
Shortcut function and reanalysis
The geva.quick function accepts a GEVAInput object and performs all
intermediate functions from the summarization to the final
concatenation. Optional (...) arguments are passed to the internal
calls to geva.summarize, geva.quantiles, geva.cluster and
geva.finalize, ultimately returning a GEVAResults object. The basic
usage is described as follows:
# Generates a random GEVAInput example
ginput <- geva.ideal.example()
# Performs all intermediate steps with geva.quick
# The resolution is used by the call to geva.cluster
gresults <- geva.quick(ginput, resolution=0.25)
## > Found 4 clusters and 31 significant genes
gresults <- geva.quick(ginput, resolution=0.4)
## > Found 16 clusters and 116 significant genes
This function can be applied to a GEVAResults object as well to
restore the parameters that produced this result, whereas optional
(...) arguments can overwrite them:
# Generates a random GEVAInput example
ginput <- geva.ideal.example()
# Performs all intermediate steps with geva.quick
# The summary.method is used by the call to geva.summarize
gresults <- geva.quick(ginput, summary.method='mean')
## > Found 60 significant genes
gresults <- geva.quick(gresults, summary.method='median')
## > Found 95 significant genes
In the example above, the entire analysis was redone using the
overwritten summary.method argument. Therefore, by following this
pattern, users can tweak different parameters depending on their
statistical choice regarding the current biological context.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.