A short introduction to the corto Package

Introduction

The corto (“correlation tool”) package [1] provides a pipeline to infer networks between “centroid” and “target” variables in a dataset, using a combination of Pearson correlation and Data Processing Inequality (DPI), first proposed in [2]. The main application of corto is in the field of Bioinformatics, Transcriptomics and Metabolomics, where co-occurrence between variables (transcripts, metabolites, or both) can be used as a means to infer regulatory mechanisms [3], gene functions [4], or metabolite abundances [5]. In this field, usually the tested features are genes (or rather, their expression profile across samples), metabolites, or any quantifiable biological molecule, whereas the centroids are Transcription Factors (TFs) and their targets are Target Genes (TGs). The TF-TG co-expression can hint at a causal regulatory relationship, as proven in many studies [6,7,8]. The corto tool follows the well-established pipeline of the ARACNe family of tools [9,10,11], using Pearson correlation instead of mutual information and, by default, a simplified DPI step (see below).

corto focuses on two aspects:

  1. Gene Network Inference, with DPI steps, Bootstrapping and optional CNV correction
  2. Gene Network Interrogation, via Master Regulator Analysis (MRA), with functions to visualize user-provided top deregulated networks on a specific signature.

The corto Gene Network Inference Algorithm

In brief, corto operates using the following steps:

  1. Load a gene expression matrix and a list of user-provided centroids (e.g. TFs). Data should be normalized via Variance Stabilizing Transformation (VST) for RNA-Seq, and Robust Multiarray Average (RMA) for microarrays, as described in [12].
  2. Removal of features with zero variance across the dataset.
  3. Calculate all pairwise Pearson correlation coefficients between centroid and target features, on the input matrix as provided. corto does not rank-transform the data: to work on ranks (i.e. Spearman correlation), a common procedure in co-expression based studies due to its benefits in reducing the effects of outliers [13], rank-transform each feature profile before running corto, e.g. inmat<-t(apply(inmat,1,rank)).
  4. Filter out all centroid-target features whose correlation coefficient absolute value is below the provided p-value threshold p.
  5. Apply DPI to the centroid-target edges. By default (dpi="target"), for each target only the edge from the centroid with the strongest absolute correlation is kept, i.e. the most likely regulator of the TG in the dataset (tied centroids are all kept). Alternatively, dpi="triangle" applies the ARACNe rule to all centroid-centroid-target triplets: an edge is removed only if it is the weakest of the three edges of a triplet, so that co-regulators of the same target that are not correlated with each other keep their edges. The triangle rule produces larger networks, closer to those of ARACNe, at a considerably higher computing time.
  6. All edges are tested for Data Processing Inequality, with the same rule, in a number of bootstrapped versions of the same input matrix (specified by the nbootstraps parameter, 100 by default, as in ARACNE-AP [10]). This step can be run in parallel by specifying the number of threads using the nthreads parameter. The number of occurrences of each edge (if they survived DPI) is calculated. This number can range from 0 (the edge is not significant in any bootstrap) to nbootstraps+1 (the edge is significant in all bootstraps, and also in the original matrix).
  7. A network is generated in the form of a list, where each element is a centroid-based list of targets. In order to follow the structure of the regulon class implemented by downstream analysis packages (e.g. VIPER [14]), each target is characterized by two parameters:
    • The tfmode, the mode of action, i.e. the Pearson correlation coefficient between the target and the centroid in the original input matrix
    • The likelihood of the interaction, calculated as the number of bootstraps in which the edge survives DPI, divided by the total number of bootstraps plus one (in this implementation, the presence of an edge in the original, non-bootstrapped matrix is considered as an additional evidence)

Running corto

Here is how to run corto. First, install the package:

install.packages("corto")

Then, load it:

library(corto)

Then, you can see how the input matrix looks like. For example, this dataset comes from the TCGA mesothelioma project [15] and measures the expression of 7000 genes across 87 samples:

load(system.file("extdata","inmat.rda",package="corto"))
inmat[1:5,1:5]
#         TCGA-3H-AB3K-01A TCGA-3H-AB3L-01A TCGA-3H-AB3M-01A TCGA-3H-AB3O-01A
# ZSCAN29         9.710627        10.045033         9.581725         9.966440
# MYOG            5.699747         5.699747         6.044523         5.699747
# PAX3            5.926491         5.699747         6.318039         6.667444
# BNC1           13.059897        13.289953        13.932597        11.860305
# SPIB            7.581019         7.818206         6.780734         7.629027
#         TCGA-3H-AB3S-01A
# ZSCAN29         9.317134
# MYOG            6.335294
# PAX3            6.303159
# BNC1           13.143449
# SPIB            5.985807
dim(inmat)
# [1] 7000   87

Another input needed by corto is a list of centroid features. In our case, we can specify a list of TFs generated from Gene Ontology with the term “Regulation of Transcription” [16].

load(system.file("extdata","centroids.rda",package="corto"))
centroids[15]
# [1] "BNC1"
length(centroids)
# [1] 294

Finally, we can run corto. In this example, we will run it with p-value threshold of 1e-30, 10 bootstraps and 2 threads

regulon<-corto(inmat,centroids=centroids,nbootstraps=10,p=1e-30,nthreads=2)
# Input Matrix has 87 samples and 7000 features
# Correlation Coefficient Threshold is: 0.889962633618839
# Removed 64 features with zero variance
# Calculating pairwise correlations
# Initial testing of triplets for DPI
# 115 edges passed the initial threshold
# Building DPI network from 30 centroids and 76 targets
# Running 10 bootstraps with 2 thread(s)
# Calculating edge likelihood
# Generating regulon object

The regulon object is a list:

regulon[1:2]
# $SCRT1
# $SCRT1$tfmode
#    ACTL6B    CACNG2     CDY1B 
# 0.9032878 0.8994689 0.9132632 
# 
# $SCRT1$likelihood
# [1] 0.2727273 0.3636364 0.3636364
# 
# 
# $ELF5
# $ELF5$tfmode
#      APOH  ARHGEF38       C8B   CEACAM5    CLDN18      CPB2    CXCL17       FGA 
# 0.9018898 0.9079885 0.9054840 0.9150803 0.9353695 0.9221696 0.8943170 0.9042021 
# 
# $ELF5$likelihood
# [1] 0.7272727 0.7272727 0.7272727 0.7272727 0.8181818 0.8181818 0.1818182
# [8] 0.7272727

The regulon in this dataset is composed of 28 final centroids with at least one target:

length(regulon)
# [1] 28
names(regulon)
#  [1] "SCRT1"   "ELF5"    "NEUROG1" "ZSCAN29" "BCL6B"   "MYOG"    "REST"   
#  [8] "MYBL2"   "FEZF2"   "PAX5"    "UHRF1"   "DMBX1"   "SOX7"    "PAX3"   
# [15] "BNC1"    "TBX21"   "IKZF3"   "EOMES"   "TFEC"    "ZNF831"  "FOXM1"  
# [22] "NKX2-1"  "SPIB"    "RFX4"    "SPI1"    "NOBOX"   "FOS"     "FOXO3"

Correcting with CNV data

As an additional, optional feature, corto gives the user the possibility to provide Copy Number Variation (CNV) data, if available, which can generate spurious correlations between TFs and targets [17]. The gene expression profiles for the target genes are corrected via linear regression and the residuals of the expression~cnv model substitute the original gene expression profile.

In this example, a CNV matrix is provided. The analysis will be run only for the features (rows) and samples (columns) present in both matrices

load(system.file("extdata","cnvmat.rda",package="corto",mustWork=TRUE))
regulon <- corto(inmat,centroids=centroids,nthreads=2,nbootstraps=10,verbose=TRUE,cnvmat=cnvmat,p=0.01)

Applying the network to another dataset

corto provides a function, mra(), to apply a previously calculated network to another dataset, in order to predict the relative levels of the centroids from their targets, measured in a different context. This can be done in three ways. Centroids can be any numerical variable, such as Transcription Factors, or Metabolites.

In all cases, mra() returns a list with the same four fields: nes (the Normalized Enrichment Scores), pvalue, sig (the signature used) and regulon (filtered for minsize).

Sample-by-sample centroid prediction

The network (regulon in the example below) is directly applied to a multi-sample dataset in the form of a matrix, where columns are samples and rows are targets (e.g. transcripts). For simplicity, all examples in this section use the same mesothelioma dataset (inmat) on which the network was generated. In this case, nes is a matrix with the same number of samples, including as rows the predicted centroids. The scores are intended as Normalized Enrichment Scores (NESs), based on a normalized T-test for one sample vs. the rest dataset. The NES is positive if the centroid network is higher in the sample vs the mean of the dataset, negative if lower.

predicted<-mra(inmat,regulon=regulon)
predicted$nes[1:5,1:5]

Contrast centroid prediction

As in the previous procedure, but the centroid score is provided as a differential NES between two sample groups (expmat1 and expmat2, the first two arguments of mra()). The resulting NES is positive if the centroid network is upregulated in expmat1 vs expmat2, negative if downregulated. In this example, the samples are split in two groups by the expression of BAP1, a tumor suppressor frequently lost in mesothelioma:

bap1<-inmat["BAP1",]
bap1low<-inmat[,bap1<median(bap1)]
bap1high<-inmat[,bap1>=median(bap1)]
predicted<-mra(bap1low,bap1high,regulon=regulon)

Signature centroid prediction

If a single named vector is provided (names are features, values are signature values such as T-test statistics or log2 fold changes), the centroid score is calculated on that signature alone, against a null model of shuffled signatures. Here, the signature is a Welch T-test statistic between the two groups above (features with zero variance give NaN values, which mra() ignores):

signature<-(rowMeans(bap1low)-rowMeans(bap1high))/
  sqrt(apply(bap1low,1,var)/ncol(bap1low)+apply(bap1high,1,var)/ncol(bap1high))
predicted<-mra(signature,regulon=regulon)

For the contrast and the signature versions of the mra() function, a mraplot() function can be applied to the output, in order to graphically show the most significant centroids (see Figure 1 in [1] for an example)

mraplot(predicted)

Funding

The development of this project (particularly the recent additions on metabolite-transcript network inference) was supported by the Italian Ministry of University and Research, under PRIN grant number 2022CEHEX8.

References

[1] Mercatelli, Daniele, et al. “corto: a lightweight R package for gene network inference and master regulator analysis.” Bioinformatics, Vol 36, Issue 12, June 2020. DOI: https://doi.org/10.1093/bioinformatics/btaa223

[2] Reverter, Antonio, and Eva KF Chan. “Combining partial correlation and an information theory approach to the reversed engineering of gene co-expression networks.” Bioinformatics 24.21 (2008): 2491-2497.

[3] D’Haeseleer, Patrik, Shoudan Liang, and Roland Somogyi. “Genetic network inference: from co-expression clustering to reverse engineering.” Bioinformatics 16.8 (2000): 707-726.

[4] Hansen, Bjoern O., et al. “Elucidating gene function and function evolution through comparison of co-expression networks of plants.” Frontiers in plant science 5 (2014): 394.

[5] Cavicchioli, Maria Vittoria, et al. “Prediction of metabolic profiles from transcriptomics data in human cancer cell lines.” International Journal of Molecular Sciences 23.7 (2022): 3867. DOI: https://doi.org/10.3390/ijms23073867

[6] Basso, Katia, et al. “Reverse engineering of regulatory networks in human B cells.” Nature genetics 37.4 (2005): 382-390.

[7] Amar, David, Hershel Safer, and Ron Shamir. “Dissection of regulatory networks that are altered in disease via differential co-expression.” PLoS computational biology 9.3 (2013): e1002955.

[8] Vandepoele, Klaas, et al. “Unraveling transcriptional control in Arabidopsis using cis-regulatory elements and coexpression networks.” Plant physiology 150.2 (2009): 535-546.

[9] Margolin, Adam A., et al. “ARACNE: an algorithm for the reconstruction of gene regulatory networks in a mammalian cellular context.” BMC bioinformatics 7.Suppl 1 (2006): S7.

[10] Lachmann, Alexander, et al. “ARACNe-AP: gene network reverse engineering through adaptive partitioning inference of mutual information.” Bioinformatics 32.14 (2016): 2233-2235.

[11] Khatamian, Alireza, et al. “SJARACNe: a scalable software tool for gene network reverse engineering from big data.” Bioinformatics 35.12 (2019): 2165-2166.

[12] Giorgi, Federico M., et al. “Comparative study of RNA-seq-and microarray-derived coexpression networks in Arabidopsis thaliana.” Bioinformatics 29.6 (2013): 717-724.

[13] Usadel, Björn, et al. “Co‐expression tools for plant biology: opportunities for hypothesis generation and caveats.” Plant, cell & environment 32.12 (2009): 1633-1651.

[14] Alvarez, Mariano J., Federico Giorgi, and Andrea Califano. “Using viper, a package for Virtual Inference of Protein-activity by Enriched Regulon analysis.” Bioconductor (2014): 1-14.

[15] Ladanyi, Marc, et al. “The TCGA malignant pleural mesothelioma (MPM) project: VISTA expression and delineation of a novel clinical-molecular subtype of MPM.” Journal of Clinical Oncology 36.15_suppl (2018): 8516.

[16] Gene Ontology Consortium. “The Gene Ontology (GO) database and informatics resource.” Nucleic acids research 32.suppl_1 (2004): D258-D261.

[17] Schubert, Michael, et al. “Gene networks in cancer are biased by aneuploidies and sample impurities.” Biochimica et Biophysica Acta - Gene Regulatory Mechanisms 1863.6 (2020): 194444. DOI: https://doi.org/10.1016/j.bbagrm.2019.194444