CEMiTool: Co-expression Modules Identification Tool

Basic usage

The CEMiTool R package provides users with an easy-to-use method to automatically run gene co-expression analyses. In addition, it performs gene set enrichment analysis and over representation analysis for the gene modules returned by the analysis.

For the most basic usage of CEMiTool only a data.frame containing expression data with gene symbols in the rows and sample names in the columns is needed, as following:

BiocManager::install("CEMiTool")
library("CEMiTool")
# load expression data
data(expr0)
head(expr0[,1:4])
##           X1913_d0  X1913_d3  X1913_d7  X1911_d0
## XIST     13.061894 13.290272 13.360468 13.178729
## DDX3Y     3.410819  3.164874  3.599792  3.400613
## RPS4Y1    6.326861  5.915121  6.341564  5.905167
## USP9Y     3.237749  3.362508  3.320674  3.365530
## CYorf15B  3.980988  4.201731  4.235020  4.046716
## EIF1AY    3.379857  3.229973  3.150274  3.196610

In this usage, the cemitool function receives the expression data, performs the co-expression modules analysis and returns a CEMiTool object:

cem <- cemitool(expr0)

To see a summary of the slots inside the CEMiTool, just call cem

cem
## CEMiTool Object
## - Number of modules: 4 
## - Modules:  (data.frame: 257x2): 
##   genes        modules
## 1  HBA1 Not.Correlated
## 2 RPS26 Not.Correlated
## 3   LYZ Not.Correlated
## - Expression file: data.frame with 4000 genes and 45 samples
## - Selected data: 257 genes selected
## - Gene Set Enrichment Analysis: null
## - Over Representation Analysis: null
## - Profile plot: ok
## - Enrichment plot: null
## - ORA barplot: null
## - Beta x R2 plot: null
## - Mean connectivity plot: null

The cemitool() function automatically executes some common analyses, depending on the input data. The following sections describes how to perform each of these analyses separately. Details on how to perform all analyses together are at the end of this vignette.

Gene filtering

As a default, the cemitool function first performs a filtering of the gene expression data before running the remaining analyses. This filtering is done in accordance to gene variance. In this example the filtering step has reduced the gene number to 257.

Module inspection

The module analysis has produced 4 modules and the allocation of genes to modules can be seen with the module_genes function:

# inspect modules
nmodules(cem)
## [1] 4
head(module_genes(cem))
##       genes        modules
## 1      HBA1 Not.Correlated
## 2     RPS26 Not.Correlated
## 3       LYZ Not.Correlated
## 4      PPBP             M3
## 5 abParts35 Not.Correlated
## 6     NAMPT             M2

Genes that are allocated to Not.Correlated are genes that are not clustered into any module.

If you wish to adjust the module definition parameters of your CEMiTool object, use find_modules(cem).

You can use the get_hubs function to identify the top n genes with the highest connectivity in each module: hubs <- get_hubs(cem,n). A summary statistic of the expression data within each module (either the module mean or eigengene) can be obtained using: summary <- mod_summary(cem)

Generating reports

The information generated by CEMiTool, including tables and images can be accessed by generating a report of the CEMiTool object:

generate_report(cem)

Also, you can create tables with the analyses results using:

write_files(cem)

Plots containing analysis results can be saved using:

save_plots(cem, "all")
## $M1
## 
## $M2
## 
## $M3
## 
## $Not.Correlated
## 
## $beta_r2_plot
## 
## $mean_k_plot
## [[1]]
## png 
##   2 
## 
## [[2]]
## png 
##   2 
## 
## [[3]]
## png 
##   2 
## 
## [[4]]
## png 
##   2 
## 
## [[5]]
## png 
##   2 
## 
## [[6]]
## png 
##   2 
## 
## [[7]]
## png 
##   2

Adding sample annotation

More information can be included in CEMiTool to build a more complete object and generate richer reports about the expression data. Sample annotation can be supplied in a data.frame that specifies a class for each sample. Classes can represent different conditions, phenotypes, cell lines, time points, etc. In this example, classes are defined as the time point at which the samples were collected.

# load your sample annotation data
data(sample_annot)
head(sample_annot)
##   SampleName Class
## 1   X1913_d0    g0
## 2   X1911_d0    g0
## 3   X1908_d0    g0
## 4   X1909_d0    g0
## 5   X1910_d0    g0
## 6   X1912_d0    g0

Now you can construct a CEMiTool object with both expression data and sample annotation:

# run cemitool with sample annotation
cem <- cemitool(expr0, sample_annot)

The sample annotation of your CEMiTool object can be retrieved and reassigned using the sample_annotation(cem) function. This function can also be used to define the columns with sample names and sample groupings (which are “SampleName” and “Class”, by default):

sample_annotation(cem, 
                  sample_name_column="SampleName", 
                  class_column="Class") <- sample_annot

Module enrichment

When sample annotation is provided, the cemitool function will automatically evaluate how the modules are up or down regulated between classes. This is performed using the gene set enrichment analysis function from the fgsea package.

You can generate a plot of how the enrichment of the modules varies across classes with the plot_gsea function. The size and intensity of the circles in the figure correspond to the Normalised Enrichment Score (NES), which is the enrichment score for a module in each class normalised by the number of genes in the module. This analysis is automatically run by the cemitool function, but it can be independently run with the function mod_gsea(cem).

# generate heatmap of gene set enrichment analysis
cem <- mod_gsea(cem)
cem <- plot_gsea(cem)
show_plot(cem, "gsea")
## $enrichment_plot

Expression patterns in modules

You can generate a plot that displays the expression of each gene within a module using the plot_profile function:

# plot gene expression within each module
cem <- plot_profile(cem)
plots <- show_plot(cem, "profile")
plots[1]
## $M1

Adding ORA analysis

CEMiTool can determine which biological functions are associated with the modules by performing an over representation analysis (ORA). To do this you must provide a pathway list in the form of GMT file. CEMiTool will then analyze how these pathways are represented in the modules.

You can read in a pathway list formatted as a GMT file using the read_gmt function. This example uses a GMT file that comes as part of the CEMiTool example data:

# read GMT file
gmt_fname <- system.file("extdata", "pathways.gmt", package = "CEMiTool")
gmt_in <- read_gmt(gmt_fname)

You can then perform ORA analysis on the modules in your CEMiTool object with the mod_ora function:

# perform over representation analysis
cem <- mod_ora(cem, gmt_in)

The numerical results of the analysis can be accessed with the ora_data function. In order to visualise this, use plot_ora to add ORA plots to your CEMiTool object. The plots can be accessed with the show_plot function.

# plot ora results
cem <- plot_ora(cem)
plots <- show_plot(cem, "ora")
plots[1]
## $M1
## $M1$pl

## 
## $M1$numsig
## [1] 9

Adding interactions

Interaction data, such as protein-protein interactions can be added to the CEMiTool object in order to generate annotated module graphs. Interaction files are formatted as a data.frame or matrix containing two columns for interacting pairs of genes.

# read interactions
int_fname <- system.file("extdata", "interactions.tsv", package = "CEMiTool")
int_df <- read.delim(int_fname)
head(int_df)
##   gene1symbol gene2symbol
## 1         DBH      REPIN1
## 2      RBFOX2       HERC5
## 3      ZNF460      CCDC22
## 4     SNRNP40        OAZ3
## 5       SRSF6        OAZ3
## 6      SPTAN1       ARL8A

You can add the interaction data to your CEMiTool object using the interactions_data function and generate the plots with plot_interactions. Once again, the plots can be seen with the show_plot function:

# plot interactions
library(ggplot2)
interactions_data(cem) <- int_df # add interactions
cem <- plot_interactions(cem) # generate plot
plots <- show_plot(cem, "interaction") # view the plot for the first module
plots[1]
## $M1

Putting it all together…

Finally, a CEMiTool object with all of the components mentioned above can also be constructed using just the cemitool function:

# run cemitool
library(ggplot2)
cem <- cemitool(expr0, sample_annot, gmt_in, interactions=int_df, 
                filter=TRUE, plot=TRUE, verbose=TRUE)
# create report as html document
generate_report(cem, directory="./Report")

# write analysis results into files
write_files(cem, directory="./Tables")

# save all plots
save_plots(cem, "all", directory="./Plots")

Troubleshooting

Sometimes, CEMiTool analyses can fail, usually due to problems in the input data. We provide a function, diagnostic_report which aims to try to assist in resolving these issues.

diagnostic_report(cem)

The function will return six different plots inside the report, all of which can be used to assess problems in the data. We will briefly discuss how each plot can be used to evaluate data problems.

Sample clustering tree

This plot aims to show if there are closely related groups within samples. If a sample annotation file is given, the plot will show different colors for each sample group, and any numerical data given in other columns as a heatmap. This information can be used in order to see how homogeneous/heterogeneous the input data are. Highly heterogeneous sample groups may be the cause of batch effects, which should be removed.

Mean x variance scatterplot

In this plot, the mean and the variance of each gene in the expression file is plotted as the x and y coordinates of the scatterplot, and a line is plotted in order to show the relationship between the two. Particularly for RNAseq data, if a strong R-squared value is found, one should set the apply_vst argument in the cemitool() function to TRUE in order to remove this correspondance.

Quantile-quantile plot and expression histogram

These plots are intended to highlight the distribution of expression values. A Q-Q plot is a mathematical approach to determine if data possibly arose from a theoretical distribution such as the normal distribution. This information can be used to guide the selection of an appropriate correlation coeffiecient for analyses, which can be changed via the cor_method argument of the cemitool() function. Currently accepted coefficients are “pearson” and “spearman”.

Beta x R2 plot and Mean connectivity plot

These plots can only be generated if the diagnostic_report() function is run after the cemitool() function. The Beta x R2 plot is used to visualize the selection of the soft-thresholding parameter Beta and its corresponding adherence to the scale-free topology model. Selected Beta values will be shown in red, unless NA. The Mean connectivity plot is intended to show the tradeoff between the network’s underlying connectivity and a higher adherence to the scale-free topology model (via higher values of the soft-threshold Beta parameter).