VCNet: vector-based gene co-expression network construction and its application to RNA-seq data

Abstract

1 Introduction

Large-scale quantification of gene expression has been made possible with advance in microarray and high throughput sequencing (HTS) technologies, and the importance of interactions among genes has been widely recognized. A network is a natural tool to represent functional interactions, and it also provides clues to the causes and mechanisms of biological phenomena (Wu et al., 2010; Yalamanchili et al., 2014a). At first, gene networks were built under a prior or given a framework of already recognized biological pathways. However, knowledge of biological pathway are limited and thus, such knowledge-driven network building approach limits the information gathered from the data obtained in an experiment. This drawback led to development of methods of network construction directly from experimental data without using prior biological knowledge. Before long, weighted network construction methods were also developed (Horvath, 2011; Zhang et al., 2005). Gene co-expression networks (GCN) constructed by these algorithms have provided new perspective to understand biological processes. In essence, GCN can detect the gene pairs who have similar and/or correlated expression patterns. These genes are often controlled by the same transcriptional regulatory mechanism, functionally related, or members of the same pathway or protein complex (Weirauch, 2011). Thus, GCN may reveal potentially new functional relationships of genes and biologic pathways, and help us to identify new targets for treatment (Cheng et al., 2012; Yang et al., 2014). Also, Hansen et al. (2014) extends the application to plants.

With the wide-spread application of RNA-seq to obtain expression data, building a weighted network based on GCN becomes more tricky. During the microarray era in the last decade, a genome-wide expression typically of 30 000 genes were obtained (Kerr et al., 2000). After filtering out genes with no/very low expression, experimental datasets would only have a reduced number of attributes (genes) in the range of 10 000–15 000. Nowadays, RNA-seq has become a routine method to analyse transcriptomes (Marguerat and Bähler, 2010). In RNA-seq, mRNA is converted to a library of cDNA fragments and sequenced in a high-throughput manner to produce short reads. Typically, 30–100 million reads are obtained in a sample. Then these reads are aligned to a reference genome to obtain counts of reads that are mapped onto transcribed regions of the genome. These regions could be genes, exons and transcripts (Wang et al., 2009). By this way, the expressions of transcriptional regions (typically over 100 000 exons or transcripts) are estimated, i.e. number of attributes increased by >10- to 20-folds from the previous microarray era. Several GCN algorithms that had been designed to handle microarray data were modified to accommodate RNA-seq data. However, limitations of these previous algorithms are apparent. There are two key hurdles which prevent a simple migration of old GCN algorithms to an application using RNA-seq data. First, the number of attributes increased by >10 times from 10 000 genes (in a microarray) to over 100 000 exons or transcripts (in RNA-seq). Second, expression of a gene quantified as a vector of exonic expression read counts for each gene in RNA-seq data instead of a scalar in microarray (one expression value per gene). Such exonic expression values may have correlation among themselves. For example, various exons of the same gene would have correlated expression levels. Although many microarray algorithms claimed that they could work for RNA-seq, they could only handle different entities of the same gene as separate entities and thus, could not produce a summative measure of these exons/transcripts. Other shortcomings include a requirement of large sample size, a lack of statistical basis and thus results are not robust, not applicable to a large dataset with such a large number of attributes. Therefore, new tailor made GCN methods are needed.

Only few tools of GCN reconstruction tailored for RNA-seq data exist (López-Kleine and González-Prieto, 2016). Two algorithms have been proposed to perform GCN analysis using a gene-vector basis to analyse RNA-seq data. Firstly, Hong et al. (2013) proposed a method called RNASeqNet to construct the GCN based on the read count vectors of genes. For a gene, each component of this vector is the expression value (read counts) of one exon of the gene, and the dimension of a vector is the number of exons. Then, exhaustive pairwise Canonical Correlation Analysis (CCA) are performed using vectors of all pairs of a given gene set. However, CCA is well known for its dimensional limitation, and it cannot work when the total number of exons of two genes is greater than the sample size. Otherwise, the results of canonical correlations are constrained to unity irrespective of the real population canonical correlations (Pezeshki et al., 2004). Yalamanchili et al. (2014b) tried to develop another algorithm termed SpliceNet to ease the requirement for sample size. They performed pairwise hypothesis tests on covariance matrices of all pairs of two genes using a large dimensional trace. SpliceNet procedure involves the exact computation of the inverse of two covariance matrices for gene pairs. However, this procedure is also constrained by dimensional restriction when the maximal number of exons of two genes is equal to or larger than the sample size. The sample covariance matrix is singular in this situation. In fact, according to TCGA annotation, ∼10% genes are with 25 or more exons, which is illustrated in Supplementary Figure S1 and Supplementary Table S7. However, few samples (e.g. <25 samples, N<25) are usually analysed in RNA-seq experiments by a typical researcher. TCGA represents the largest scale of RNA-seq analysis nowadays and even so, for some tumour no more than 25 samples are analysed. Therefore, this dimensional issue is a real problem, and a novel method which can handle vectors with a large number of exons in datasets with only a moderate sample size is needed.

In this paper, we propose a new algorithm called VCNet to construct GCN from RNA-seq data. VCNet performs a new hypothesis test based on transformed correlation matrix of a gene–gene pair using the Frobenius norm. The Frobenius norm is a commonly used metric of matrices, with the advantage of having an asymptotic distribution of test statistics. Similar to RNASeqNet, we use the read counts vector to represent the expression of a gene. VCNet alleviates the dimensional bottleneck of sample size because it does not require computation of inverse of any matrix. Moreover, VCNet used the asymptotic distribution of its statistics to determine the P-value. In the simulation study, the true positive rate of VCNet is higher than those of SpliceNet and RNASeqNet under the same false positive rate. VCNet is very efficient to detect gene interaction, even in situations where the other two methods failed. Besides, real RNA-seq datasets were also analysed. In the analysis of the normal breast tissue, we find the network inferred by VCNet contained more biologically meaningful edges which are in agreement with literature than those inferred by SpliceNet and RNASeqNet. In differential network analysis of kidney chromophobe cancer dataset, VCNet can identify the interactions associated with EP300 and CREBBP together with many components in the Hypoxia-Inducible Factor pathway, while SpliceNet and RNASeqNet fail to find them because the exon numbers of these two genes are larger than the sample size.

2 Materials and Methods

There are various approaches to represent gene expression data in RNA-seq datasets. Traditionally, the RPKM of gene is used to represent gene expression, which is a scalar. However, this approach does not take full advantage of the high resolution in RNA-seq. Therefore, vector-based expression is proposed. So far, there are three different ways to express RNA-seq data as vectors. First, one vector corresponds to one gene. Each entry represents counts of sequence reads mapped to corresponding exon of a gene. Therefore, the dimension of a vector is equal to the number of exons of that gene. Second, one vector corresponds to one isoform and it is similar to a exonic vector. However, each component is the read counts mapped to a corresponding exon that is specific to that isoform. Third, expression levels of various isoforms of a gene form the component of a gene-vector, in which the dimension is the number of isoforms belonging to the gene. Each component is the estimated expression of the corresponding isoform. PennSeq (Hu et al., 2014), RSEM (Li and Dewey, 2011) and WemIQ (Zhang et al., 2014) are commonly used software to estimate isoform expression. The illustration of input data is shown in Supplementary Figure S6. The detailed discussion about the four different formats of data input is in the Supplementary Material. After extensive evaluation, the performance of the exon expression of a gene is better than those of other inputs. We select the exonic expression of a gene as our input data for the main text.

2.1 Methods to determine pair-wise co-expression genes using gene vectors

In practice, we use the eigenvalues of $R^11$ and $R^22$ to estimate the $(λ1,…,λp)$ and $(μ1,…,μq)$⁠, where $R^11$ and $R^22$ are the sample correlation matrices of X and Y, respectively. A co-expression edge is kept between two genes if we reject the null hypothesis. We reject it when the observed value of $(n−2.5)TVCNet$ is larger than the critical value $Tα$⁠. The critical value $Tα$ is the αth upper quantile of $∑i=1p∑j=1qλiμjχij2$⁠. The asymptotic distribution of $(n−2.5)TVCNet$ under null hypothesis is a mixture of χ² distribution, which was approximated by Davies method (Davies, 1980). False Discovery Rate (FDR) is used to account for multiple testing.

Compared with SpliceNet and RNASeqNet, VCNet completely avoided the need to compute the inverse matrices of the covariance matrices. Therefore, it is tailored to work under large dimension and small sample size condition. In contrast, SpliceNet performs hypothesis testing on the covariance matrix $Σ12$ between X and Y, and is vulnerable to the heterogeneity of variances which is a common encounter in gene expression data. VCNet, as an general GCN algorithm, can be used with different inputs. The specific results are in the Supplementary Material.

3 Result

3.1 Simulation study

To make the simulation data more close to real data, we conduct the following simulation procedure. We first generated the isoform expression of genes. And then exonic expressions corresponding to this gene were obtained according to its isoform expressions. Specifically, $U=(U1,U2)$ and $V=(V1,V2)$ follow two-dimensional normal distributions $N2(0,Σ1)$ and $N2(0,Σ2)$⁠, respectively. $Σ1=(1c0c01)$ and $Σ2=(1−c0−c01)$⁠. Here, we assume that both of genes have two isoforms. $F(1)=(U1,V1)$ is the isoform expressions for the first gene and $F(2)=(U2,V2)$ is the isoform expressions for the second gene. The two genes are independent when $c0=0$ and are correlated when $c0≠0$⁠. We use equation $X=Ap×2(1)F(1)+ϵ$ and $Y=Aq×2(2)F(2)+ϵ$ to obtain the exonic expression corresponding to the two genes. F is the isoform expression of a gene. A is an indicator matrix, its row corresponding to exon and column corresponding to isoform. A_ij takes one if the j^th isoform contains the i^th exon and zero otherwise. ϵ follows an independent normal distribution with zero mean and one variance. These two equations are an approximation to the relationship between exonic and isoform expression. Many methods to estimate isoform expression are based on them (Jiang and Wong, 2009). Note that we use X and Y as input. For gene pair X and Y, we set the exon numbers (p and q) of two genes to 5-5, 15-20, 25-25, and more extreme values of 10-100, 90-100 and the sample size to 20, 30, and 40. The specific A under different exon number setting are shown in Supplementary Table S15. For c₀, we set it to 0.3, 0.4, and 0.5. In total, there are 45 combinations of settings in the simulation. For each combination, 2000 samples were generated. Half of them are from the null model (⁠$c0=0$⁠) and the remaining samples are from the alternative model (⁠$c0≠0$⁠). We use the Receiver Operating Characteristic (ROC) curves and area under curve (AUC) to evaluate the performances of VCNet, SpliceNet and RNASeqNet. ROC curves are displayed for the three methods in Figure 1.

Fig. 1

ROC curves comparing various algorithms using the simulation study. Six different sets of simulation experiments are shown with different number of exons in two genes under evaluation. p and q are the number of exons in the two genes, respectively. n represents the sample size and c₀ represents the dependency of two vectors. In some simulation settings, SpliceNet and/or RNASeqNet failed to work and their ROC are not shown. A vertical line is drawn at false positive rate = 0.1 to compare the corresponding true positive rates among the algorithms

Open in new tabDownload slide

On the top row of Figure 1, the simulation results of the setting of sample size = 40 and two genes with 5 exons. Under this condition, all three algorithms showed similar performance though the association strength c₀ had some effect on the performances among them. It is observed that AUCs of all the methods were significantly increasing with the increasing association strength. When c₀ is reduced to its lowest end, VCNet has higher true positive rate than the other two methods at the same false positive rate of 0.1. SpliceNet is better than RNASeqNet.

Next, we evaluate the influence of the sample size and the exon number on the performances of three methods. The second row of Figure 1, showed the simulation results of the setting of two genes with 15 exons and 20 exons, respectively. Effect of different sample size on the performance was shown when c₀ was set to 0.5. When the sample size is 20, the total number (15 + 20) and the maximal number (20) of the exons for two genes are larger than or equal to the sample size. Therefore, RNASeqNet and SpliceNet failed to analyse such datasets, and there was no ROC curves for them. In contrast, VCNet shows pretty good ROC curve. An increasing trend of performance for VCNet was also observed with sample size from 20 to 40. When the sample size increased to 30 or 40, VCNet showed the largest AUC among the three algorithms. It shows the ability of VCNet to handle small samples, which is a bottleneck of many existing algorithms. We also display AUCs for all methods under different conditions for a comprehensive comparison. The detailed result is shown in Table 1.

In summary, we can see that VCNet is optimal among the three methods in handling simulation data. Overall, VCNet needs fewer samples than SpliceNet and RNASeqNet to achieve the same true positive rate.

3.2 Analysis of real data

3.2.1 Network construction

In this section, we evaluate whether VCNet can reconstruct networks that is closed to the ground truth (known biological pathways). The exon expressions of genes (level 3, exon-level RPKM) of RNA-seq data from healthy breast tissue are downloaded from TCGA data portal. In total, there are 113 normal samples. We analysed two pathways: ‘resolution of D loop structure’, ‘homologous DNA pairing and strand exchange’ from Reactome Database (Croft et al., 2014; Milacic et al., 2012) for illustration. We use Reactome and Physical interaction from GeneMANIA as gold standard to evaluate the predictions of three methods, seperately. To display the effects of sample size on the performance of network construction, we randomly re-sampled N (⁠$N=20,30,40,100$⁠) samples to produce validation datasets. VCNet, SpliceNet and RNASeqNet were then analysed using these re-sampled datasets, separately. This procedure was replicated 20 times. Twenty networks were constructed by each algorithm under each sample size setting. For VCNet and SpliceNet, we obtain the edges under the cut-off 0.05 of FDR. For RNASeqNet, we select the 5% of edges with top weights as suggested by original paper. We compare the networks to the two ground truths. We use precision, recall to evaluate the performances. F-score is also employed and $F−score=2×Precision×RecallPrecison+Recall$⁠. F-score is a balance between precision and recall. A larger F-score indicates a good behaviour. The results for the two pathways are shown in Table 2. The pathway ‘resolution of D loop structure’ using Reactome Pathway as a gold standard is shown in Table 2. For all algorithms, the precision, recall and F-score increased with increasing sample sizes. For low sample size, VCNet had the highest precision, recall and F-score as well.

One reason for low F-score of SpliceNet is that the networks constructed by SpliceNet are almost empty when sample size is under 20. The maximum exon number of genes is 63 and there are 5 genes whose exon numbers are >20 in the ‘resolution of D loop structure’ pathway. Therefore, SpliceNet and RNASeqNet fail to infer some edges at sample size of this range. SpliceNet is getting better when the sample size is increasing. Next, when we take physical interaction as ground truth, the performances for all methods are worse than those of Reactome Pathway since co-expression network reflects the transcriptional regulation. Similar results are obtained for the ‘homologous DNA pairing and strand exchange’ pathway. There are 8 genes whose exon numbers are >20 in this pathway. Therefore, SpliceNet and RNASeqNet performed poorly when sample size was 20.

As some genes in above pathways have large numbers of exons, so SpliceNet and RNASeqNet failed to work. Next, we evaluated three methods on pathways whose genes contain small number of exons. Therefore, all methods are valid on these pathways. The sample size is 65. The pathway information is displayed in Supplementary Table S9. The precision, recall and F-score are shown in Supplementary Table S11. We can see that VCNet is better than the other two methods. In Table 2, when the sample size is 100, SpliceNet is slightly better than VCNet under pathway1. Covariance matrix contains more information than correlation matrix. Therefore, SpliceNet may perform over VCNet under some cases. Above analysis indicated that VCNet is comparable to SpliceNet in the pathways on which all methods are valid.

3.2.2 Differential network analysis between tumour and normal tissues

Comparison between a normal network and a tumour networks is an important means to identify key genes related to pathogenesis. We applied VCNet to Chromophobe Renal Cell Carcinoma (KICH) which was downloaded from TCGA data portal. There are 25 normal samples and 65 KICH tumour samples. We analyse the ‘regulation of hypoxia inducible HIF by oxygen’ pathway and ‘Tp53 Regulates Metabolic Genes’ pathway, which are related to pathogenesis of KICH (Haase, 2012). For each method, two co-expression networks are constructed using the normal and case samples, respectively. The FDR cut-off was set at 0.05.

Like the analysis in normal breast tissue, we compare the constructed networks from VCNet, SpliceNet and RNASeqNet to the two gold standards. Precision, recall and F-score are shown in Table 3 for different sample size and conditions. We first inspect the networks constructed from the 25 normal samples. For ‘regulation of hypoxia inducible HIF by oxygen’ pathway, under the Reactome Pathway as a gold standard, SpliceNet has the largest precision and VCNet has the biggest recall. VCNet is the highest when we turn to F-score. When the physical interactions are as golden standard, the same trend is observed. As we do not know the truth network under tumour condition, we still compare the networks from tumour samples to the two golden standards. To investigate the impact of sample size on the performances of all methods, we select 25 samples from the tumour data and apply three methods to them. We also apply three methods to 65 tumour samples. The results are shown in Table 3. Under 25 tumour samples, RNASeqNet has the highest precision under two criteria. VCNet has the largest recall under two criteria. In terms of F-score, VCNet is the best. Under the 65 tumour samples, SpliceNet is better than the other two methods in term of F-score. For ‘Tp53 Regulates Metabolic Genes’ pathway, we can see that VCNet is the best in terms of F-score under both normal and tumour samples.

We compare the networks of ‘regulation of hypoxia inducible HIF by oxygen’ pathway inferred by VCNet using normal data and tumour data, which is displayed in Figure 2. For SpliceNet and RNASeqNet, the differential networks are displayed in Supplementary Figures S2 and S3. In Figure 2, the dotted line represents the edges which were only detected in tumour tissue, and the solid line represents the edges which were only detected in normal tissue. The edges presenting in both networks are removed. We find that two genes EP300 and CREBBP have interactions with HIF1A in the tumour networks of three methods, which is in agreement with literature. And String Database also supports these interactions, which is shown in Supplementary Figure S4. It has been reported that EP300 protein can bind to CREBBP protein and activate the expression of HIF1A (Ema et al., 1999). HIF1A encodes the alpha subunit of transcription factor hypoxia-inducible factor (HIF-1) (Iyer et al., 1998). HIF-1 functions as a master regulator of cellular and systemic homeostatic response to hypoxia by activating transcription of many genes (Ietta et al., 2006). And the over-expression of HIF-1 is heavily implicated in promoting tumour growth and metastasis of kidney cancer (Bos et al., 2003). Above analysis indicates that VCNet, SpliceNet and RNASeqNet can detect the biologically meaningful correlations using the large case samples. We also showed the inferred networks of VCNet, SpliceNet and RNASeqNet under normal and cancer tissues. The six networks are displayed in Supplementary Figures S7–S9. The differential networks of ‘Tp53 Regulates Metabolic Genes’ pathway are shown in Supplementary Tables S12–S14.

Fig. 2

The differential network between normal and KICH cancer as constructed from results of VCNet. The co-expression network of normal and cancer were inferred separately. Solid lines represent those differentially co-expressed edges which were only detected in normal tissue and the dotted lines were those detected only cancer. Number of exons in two genes (EP300 and CREBBP) were larger than sample size and VCNet was able to handle their RNA-seq data. (c.f. Figures S2 and S3 in Supplementary document show differential networks of results produced by other algorithms)

Open in new tabDownload slide

EP300 and CREBBP have interaction with HIF3A instead of HIF1A in a normal network constructed by VCNet. These two edges were not present in the normal networks inferred by SpliceNet and RNASeqNet. As both EP300 and CREBBP have over 30 exons which exceed the sample size in normal tissue datasets. Both SpliceNet and RNASeqNet fail to detect them. HIF3A has a high similarity with HIF1A in the bHLH and PAS domains but lacks structures for transactivation found in the C-terminus of HIF1A (Gu et al., 1998). Since the C-terminal structure of HIF3A is different from that of HIF1A, the transfected HIF3A protein shows different characteristics from those of the transfected HIF1A. HIF3A suppresses hypoxia-inducible HIF-mediated gene expression and might be a negative regulator of hypoxia-inducible gene expression in the human kidney (Hara et al., 2001). Network constructed by VCNet suggested a differential co-expression relationship between HIF3A and other genes in cancers versus controls tissues, which underpins the ability of VCNet to identify meaningful interactions even in small samples.

As for the prediction for interplays between genes, CUL2, HIF1A have interactions in the tumour network constructed by VCNet. This interaction is also supported by String Database (Supplementary Figure S4). CUL2 is also required for HIF1A to mediate the expression of other genes in tumour samples. Some study shows that CUL2 is required for normal vasculogenesis in the zebrafish embryos, at least in part mediated by its regulation of HIF-mediated transcription (Maeda et al., 2008). Besides, both EP300 and CREBBP have interactions with CUL2. It indicates that EP300 and CREBBP may participate the expression of CUL2. Both EP300 and CREBBP are transcriptional co-activators (Grunstein, 1997; Gusterson et al., 2003) and have functions in cell differentiation (Giles et al., 1998), cell proliferation (Goodman and Smolik, 2000; Yao et al., 1998). Above analysis demonstrates that VCNet can provide rich molecular information to biological process, especially for genes whose exon numbers are large.

3.2.3 Comparison of computation time

Since SpliceNet is designed to run on Linux system and RNASeqNet is designed to run on Windows system, the running time of them cannot be compared directly. We run VCNet on Linux system and Windows system separately and then compared the running times to SpliceNet and RNASeqNet. CPU used in the Windows was Intel(R) Core(TM) i7-2600 3.40 GHz and the R version is 2.15.2. The CPU used in Linux system was Intel Xeon Processor E5-2670 2.60 GHz and the R version is 3.0.1. The dataset were same as that used in differential network analysis. And we analyse the 18 genes related to ‘regulation of hypoxia inducible HIF by oxygen’ pathway. The results are shown in Table 4. All the methods completed in a short time. Constructing network using all expressed genes (whole genome analysis) will be more demanding and take more time. Theoretically we may analyse the whole well annotated sets of 20 497 human genes, such complete pairwise network will be consists of 210 053 256 pairs. We attempted this challenge by modifying the code of VCNet for parallel computing with help of Rcpp package. The runtime of VCNet is ∼20–24 h using 43 CPUs on a Linux system. SpliceNet and RNASeqNet are not designed for parallel computing. Therefore, we did not apply them to the whole genome dataset. There are no technological constraints for VCNet. Its algorithm is simple. More importantly, VCNet has a test statistics of asymptotic distribution, therefore, we can obtain P-values from the distribution directly rather than the permutation strategy.

4 Discussion

With the advance of the next generation sequencing technology, we can obtain the most precise and accurate measurement of gene expressions possible so far, and it will enable more insights into biological processes in health and disease conditions. However, the increased data amount also created a new problem in term of computer storage, processing algorithm and not the least, challenges in the mathematical and statistical analysis of such massive data. In order to gain new insight into biological pathways, analysis of gene network represents a fundamental analysis of RNA-seq data, however, development of accurate co-expression analysis applicable to the high dimensional RNA-seq data is still in its infancy. Most methods developed previously during the microarray era suffer from the ‘curse of dimension’. As a consequence, determination of covariance matrix is singular, and its inverse matrix is undefined when elements of vectors exceed the number of samples. Considering this problem, we proposed a new approach using Frobenius norm to construct a new test statistics and provided its asymptotic distribution. This method can provide an adequate analysis of gene vectors co-expression network of RNA-seq data. This new test statistic computes squared Frobenius norm of the transformed correlation coefficient matrix between two random vectors. Therefore, it does not need the inverse of any matrix and provide the unique advantage for VCNet in the analysis of vectors of gene expression. We use Frobenius norm to evaluate whether the transformed R₁₂ is equal to zero matrix. There are several advantage in the application of Frobenius norm. One is that the test statistics based on the norm are relatively easier to determine than other norms, in particular, in term of the asymptotic distribution of the test statistics. Furthermore, Frobenius norm is robust against the noise. VCNet shows high AUC in the simulation study and high F-scores in real data analysis especially when the sample size is moderate. Only VCNet can detect the interactions between genes with large vector attributes (e.g. large exon numbers) while SpliceNet and RNASeqNet miss these edges. Moreover, we provide support for an asymptotic distribution of the new test statistic under null model, which reduces the computation complexity of P-values a lot.

When we prepare for this paper, another group proposed another method to test the independence of two random vectors (Li et al., 2015). Their method is similar to ours. However, they performed the hypothesis test using the covariance matrix. As we state in Section 2, the hypothesis test on the covariance matrix $Σ12$ is easily influenced by the heterogeneity of variances. The assumption of equal variances of exons vectors in the RNA-seq experiments would be frequently violated. We compared VCNet to their method and the results are shown in the Supplementary Material. Among most simulation conditions, VCNet has larger AUC than theirs. Specht et al. (2015) proposed a generalization of the Pearson’s correlation coefficient to construct GCN. However, their method is confined to the total read counts mapped onto a gene, instead of using vector of read counts mapped to individual exons. We think that the read count vectors of genes contain more information, in particular it would encompass the features of alternative splicing. We also investigate the effects of imbalance between p and q on the performances of three methods. The results are shown in Supplementary Table S10, which does not give a clear trends with the ratio of p and q.

VCNet also has another advantage of robustness even when the data do not follow a multivariate normal distribution. The simulation procedure can be found in the Supplementary Material, and the results are shown in Supplementary Figure S5. It showed that the ROC curves are similar to those in Figure 1, which demonstrate the test statistic is robust even under non-Gaussian distribution.

Like the case in all co-expression analysis, VCNet cannot distinguish if the interaction between two genes is direct interaction or indirect interaction. The inference of direct interaction is another field of causative analysis. Traditionally, direct interaction network may be inferred by the inverse of the covariance matrix in graphic model theory. In the case of microarray data, the (i, j) element in the inverse matrix corresponds to the gene pair of (i, j). If some element of the inverse matrix is non-zero, we say that the corresponding gene pair is directly associated. Feizi et al. (2013) proposed a network deconvolution method to obtain the direct interactions, while their method is still based on the gene expression of a scalar. On the contrary, in the case of RNA-seq data, the expression of a gene is denoted by a vector. How to define the direct interaction between two vectors is still unsettled and is a topic of further research.

In conclusion, VCNet based on an asymptotic distribution of Frobenius norm provides a good algorithm to compare two vectors of gene expressions. In this report, the elements of the vectors are read counts of each exons. It can be extended to analyse the isoform–isoform co-expression network, where various exons of an isoform form the vector. Moreover, 95% of human genes undergo alternative splicing (Pan et al., 2008), it will be a suitable data for analysis. Further application of VCNet may include lncRNA–mRNA interaction or microRNA–mRNA interaction. Prediction of interactions between lnRNA and its target mRNAs will help us uncover the mechanisms of biological processes.

Acknowledgement

Part of the analysis was performed on the Computing Platform of the Center for Life Science.

Funding

This work was supported by the National Key Basic Research Project of China [grant numbers 2015CB910303, 2016YFA0502303] and the National Natural Science Foundation of China [grant numbers 31171262, 31428012, 31471246].

Conflict of Interest: none declared.

References

Author notes