DiCoExpress: a workspace to process multifactorial RNAseq experiments from quality controls to co-expression analysis through differential analysis based on contrasts inside GLM models.

Background RNAseq is nowadays the method of choice for transcriptome analysis. In the last decades, a high number of statistical methods, and associated bioinformatics tools, for RNAseq analysis were developed. More recently, statistical studies realized neutral comparison studies using benchmark datasets, shedding light on the most appropriate approaches for RNAseq data analysis. Nevertheless, performing an RNAseq analysis remains a challenge for the biologists. Results DiCoExpress is a workspace implemented in R that includes methods chosen based on their performance in neutral comparisons studies. DiCoExpress uses the pre-existing R packages as well as FactoMineR, edgeR and coseq, to perform quality control, differential, and co-expression analysis of RNAseq data. Users can perform the full analysis, providing a mapped read expression data file and a file containing the information on the experimental design. Following the quality control step, the user can move on to the differential expression analysis performed using generalized linear models with no effort thanks to the automated contrast writing function. DiCoExpress proposes a list of comparisons based on the experimental design, and the user needs only to choose the one(s) of interest for his research question. A co-expression analysis is implemented using the coseq package. Identified co-expression clusters are automatically analyzed for enrichment of annotations provided by the user, and several result outputs proposed. We used DiCoExpress to analyze a publicly available Bra ssica napus L. RNAseq dataset on the transcriptional response to silicon treatment in plant roots and mature leaves. This dataset, including two biological factors and three replicates for each condition, allowed us to demonstrate in a tutorial all the features of DiCoExpress.

Conclusions DiCoExpress is an R workspace to allow users without advanced statistical knowledge and programming skills to perform a full RNAseq analysis from quality controls to co-expression analysis through differential analysis based on contrasts inside generalized linear models . Hence, with DiCoExpress, the user can focus on the statistical modeling of gene expression according to the experimental design and on the interpretation of the results of such analysis in biological terms.

Background
During the last decades, Next-Generation Sequencing (NGS) technologies have developed at a fast pace with the improvement of data quality coupled with a reduction of experimental costs. Since the early years of NGS, the use of RNAseq to profile transcriptomes became the method of choice replacing in time microarraybased analyzes [1]. Plant biologists use RNAseq-based transcriptomic extensively, generating knowledge about transcriptional regulation in several biological processes [2][3][4][5]. Differential gene expression analysis across different experimental conditions is classically used to gain insight into gene regulation events and gene co-expression analysis to identify functional modules.
A classical analysis workflow will start with a data normalization step to account for technical biases that affect the number of reads mapped to a gene. Several methods are available, and among the most used, we can find RPKM (Reads Per Kilobase per Million mapped reads) [6], Upper quartile normalization [7], RLE (Relative Log Expression) [8] and TMM (Trimmed Mean of the M-values) [9]. Multiple methods, based on different statistical modeling of data, are available to perform differential expression analysis. Negative binomial-based models with robust meanvariance modeling, have been used extensively at the beginning, and they are available in the R-packages edgeR [10] and DESeq [8]. More recently, the linear models and their generalized extensions for negative binomial distributions (GLM) have been proposed to account for the versatility of multifactorial experiments.
They are available in the R-package limma [11] for the linear models and in the Rpackages edgeR [12] and DESeq2 [13] for the generalized linear models. Following differential gene expression analysis, several approaches to identify and group coexpressed genes have been in use over the years. Pearson's or Spearman's correlations, WGCNA (Weighted correlation network analysis) method [14], hierarchical clustering and K-means are the most conventional approaches found in the literature [15,16]. With these approaches, the number of clusters is chosen, either a priori or a posteriori,by the user. Mixture models offer a different approach by identifying an underlying structure which corresponds to clusters of co-expressed genes. Moreover a model selection criterion allows determinining the most appropriate cluster number [17,18].
To perform such analysis, tools associated with the methods are available in R to quite quickly get from data to results [19]. In parallel bioinformatics tools offering a Graphical User Interface (GUI), and interactive visualization tools were developed.
First-generation tools included the RNAseq read-mapping step [20][21][22], more often realized independently at present, depending on data and genome availability.
Several of these GUI tools [23][24][25][26][27][28][29] ease the use of the main RNAseq analysis Rpackages for normalization and differential expression analysis such as limma [12], DESeq [9], DESeq2 [14] and/or edgeR [13]. The role of these GUI is to realize Rbased RNAseq data analysis with little or no experience in the command line. More recent tools take advantage of the R-shiny framework that eases the creation of a GUI for R-packages and pipelines [30]. The majority of these GUI tools include a high number of data visualization options and the possibility to generate figures for publications.
Even with all these tools, a biologist is often in front of a real dilemma on how to analyze his dataset correctly. Indeed a characteristic shared by the majority of GUI tools developed up to date is to offer the user the possibility to choose among multiple statistical methods for each step of the analysis with no specific propositions. However, the RNAseq data specificities, such as heterogeneity of counts or overdispersion among biological replicates, represent a methodological challenge that has to be addressed by proper statistical modeling of the gene expression. It is worth noting that in the case of multifactorial experiments, if interaction terms are included in the modeling, the writing of the contrasts might become tricky, requiring a good understanding of some statistical concepts, not always mastered by a biologist. As a result, the large-scale data analysis of RNAseq data is not straightforward for a biologist.
DiCoExpress aims to offer a tool usable without advanced statistical knowledge and/or programming skills to analyze RNAseq projects with complete and balanced designs with at most two biological factors and one technical factor. To offer a validated set of tools, we based our choice on three neutral comparison studies [18,31,32]. The idea of such studies is to design and implement a framework to generate realistic benchmark datasets with known truth to make an objective and reproducible performance assessment. Comparing normalization methods, Dillies et al. [31] showed that the RLE method implemented in the package DESeq2 [14] and the TMM method implemented in the package edgeR [11] demonstrate satisfactory behavior in the presence of highly expressed genes. Both these methods maintain a reasonable false-positive rate without loss of power. The choice of both methods was confirmed by Reddy et al. and Evans et al. [33,34] even in experiments with slightly asymmetric differential expression or different amounts of mRNA/cell per condition. Based on these detailed evaluations, both RLE and TMM are suitable, but we decided to choose the TMM normalization as the default method and proposed RLE as an alternative for normalization due to the choice made for the differential analysis described below.
Rigaill et al. [32] made a neutral comparison study among differential gene expression methods, including negative binomial-based, generalized linear models, and linear models on transformed data. Performance analyzes based on the p-value distributions, ROC curves, and proportion of true and false-positive rates show a clear difference of behavior between negative binomial-based methods and the others. Linear models on transformed data or generalized linear models are consequently the most adapted for the differential analysis. Among these models, as also observed in Schurch et al. [35], when the proportion of differentially expressed genes is low, the results obtained with the method implemented in the edgeR package are more satisfying. We thus chose the statistical model implemented in the edgeR package as a method of choice for differential expression data analysis. Moreover, we propose automatic writing of a large number of contrasts in order to facilitate the comparisons between the biological conditions considered in the experimental design. This automatic writing is a real advantage because, in the available R-packages, most contrasts in GLM with interactions between two factors must be handwritten and require thus an excellent understanding of the statistical modeling.
For the co-expression analysis, we preferred mixture models to correlation-based approaches. Mixture models aim at identifying an underlying structure in modeling the unknown distribution by a weighted sum of parametric distributions, each one representing a group of co-expressed genes. Gaussian mixture models were relevant for microarray data and were applied with success on several datasets [36,37]. For RNAseq data, which are discrete, Rau et al. [18] first concluded that normalized expression profiles modeled with a Poisson mixture are relevant for coexpression analysis. However, in the Poisson mixture, the dependence structure between samples is not considered and can mislead the results. To tackle this problem, they proposed then a Gaussian mixture after a transformation of the normalized expression profiles [19]. This model seems to be more suitable for RNAseq co-expression analysis by providing a proper identification of the groups of co-expressed genes because it accounts for per-cluster correlation structures among samples. For these reasons, we chose this Gaussian mixture implemented in the coseq R-package.
In conclusion, using these neutral comparison studies, we combined the most adapted tools for each step of a standard RNAseq analysis. DiCoExpress is a workspace, illustrated in Figure 1, to be installed on a computer to create a userfriendly workspace for analyzing RNAseq datasets. The directory Data will store all the projects, and the directory Results will contain a subdirectory per project with all the results of the different steps. The directory Sources contains the R functions used by DiCoExpress. Finally, the directory Template_scripts will contain an R script file for each project, allowing a semi-automated data analysis where the user is guided through all the steps from normalization to co-expression analysis. Hence, with DiCoExpress, our objective is to focus on the statistical modeling of gene expression according to the experimental design and on the interpretation of the results of such analysis in biological terms.

Implementation
To create DiCoExpress, we use the R programming language and several R-packages from CRAN and Bioconductor [39]. Each step of the analysis has a dedicated function available in the directory Sources. Seven functions compose the core of DiCoExpress (Fig. 2), and they are combined in a script, stored in the directory Template_scripts for each project to specify the steps of the analysis and the parameters to use. A full description of these seven functions is available in Additional File 1.

Input files and data quality controls
To run DiCoExpress on a project, the user has to provide only two input files: one containing a count table summarizing the mapped reads for each gene, named Project_Name_COUNTS.txt, and a second one with a description of the project design according to the experimental factors, named Project_Name_TARGET.txt. If functional gene annotations are available in a file, the user has the option to upload it. This information will be integrated into the result tables and can also be used to perform enrichment tests.
The (1)  DiCoExpress performs analysis for a complete and balanced experimental design. If this condition is not verified, then an error message appears, and the script stops running. A filter option in Load_Data_Files is proposed to extract a subset of samples leading to complete and balanced design, thus avoiding manual modifications of the expression file. The filtering rules are described according to the Project_Name_TARGET.txt (see section Results for an example). The (2) Quality_Control function produces several representations of the dataset before and after filtering low expressed genes and correction of the library size effect. This step is optional, but we advise the users to perform it to evaluate the quality of the RNAseq data before further analyses.

Differential expression analysis
The differential analysis is based on a negative binomial GLM, where the log of the gene expression is modeled by all the factors describing the experiment. When the number of observations is twice greater than the number of parameters of the model, we advise to include interaction terms between the biological factors. Such terms in the gene expression definition might reveal meaningful interactions such as genotype-environment interaction and answer in a direct way to some biological questions [40][41][42]. The (3) GLM_Contrasts function will automatically write a list of contrasts based on the model specified by the user. We focused on contrasts involving the biological factors, and their names are sufficiently explicit to understand the associated biological question addressed. For example, we proposed automatic writing of the difference between two modalities of a biological factor averaged on the second factor or for a given modality of the second factor. Hence thanks to this list of proposed contrasts, the user will be able to choose the ones that are relevant for the project without worrying about the complexity of their statistical formulation. Running this function is a prerequisite to run the differential expression analysis. The (4) DiffAnalysis_edgeR function uses edgeR R-package to estimate the parameters of the GLM and then test every contrast chosen by the user. As proposed by Rigaill et al. [32], the distribution of raw p-values of each contrast is inspected to assess the quality of the statistical modeling of the gene expression. Since the distribution of raw p-values is theoretically dominated by a uniform distribution, the fit between the statistical model and the data can be observed on these raw p-value histograms. If the raw p-value distribution is not satisfactory (see example in the Results here below and Fig. 4), we advise repeating the analysis using a more stringent cut-off for the filtering step or another rule of filtering. If the raw p-value distribution remains unsatisfactory, the problem might come from a large number of parameters compared to the number of observations available to estimate them. In this case, we advise modifying the modeling of the gene expression removing, for example, the interaction term. For each contrast, a subdirectory is created to store the differentially expressed genes (DEGs) lists and other useful results files.

Co-expression analysis
The (5) Venn_Intersection_Union function helps the biologist in the interpretation of the results by comparing different DEG lists. This function also generates the union and/or the intersections of these DEG lists in order to perform a co-expression analysis with the (6) Coexpression_coseq function. This latter function uses coseq Rpackage [18] to transform the raw data into normalized expression profiles. We kept the filter function of coseq removing the genes with low mean normalized counts.
Those discarded genes are assigned in Cluster 0. A co-expression analysis is

Enrichment analysis
Once an RNAseq analysis is complete, the next step is to evaluate the coherence of the results by comparing them with biological knowledge. To this end, the (7) Enrichment function performs hypergeometric tests in order to find annotation terms that are specifically enriched or depleted in a given list of genes with respect to a reference specified by an annotation file. The enrichment analysis following the co-expression analysis is automatically performed on all the co-expressed gene clusters. This function can also be applied to any list, e.g., lists of differentially expressed genes from the GLM analysis.

Results and Discussion
We illustrate the use of DiCoExpress by analyzing a dataset associated with the publication of Haddad et al. [43]. This RNAseq dataset describes gene expression in roots and mature leaves of Brassica napus with or without silicon (Si) treatment.
Three biological replicates are available. The experimental design can be described by two biological factors Tissue and Treatment and a technical factor Replicate with three modalities (Fig. 3). To illustrate the outputs of enrichment tests, we used the annotation of B. napus v.5 from the Brassica genome database [44] to perform enrichment analyzes.
We tested DiCoExpress on the full dataset available in contrast to Haddad et al., who only focused on the root samples. The procedure is described in Additional File  Fig. 1A and 1B). A hierarchical clustering heatmap and principal component analysis graphs are generated to look at the sample similarities. In our analysis, we observe, as expected, a clear difference between the two tissues as well as an apparent clustering of mature leaf samples according to the treatment (Supplementary Fig.   1C and 1D).
We performed a differential expression analysis using a GLM with both biological factors and the technical replicate factor. We included an interaction between the two biological factors in the model. We checked the quality by looking at the raw p- observe an increase of the frequency around 1: this usually suggests that data are not properly filtered (Fig. 4A). Following this observation, we went back to the beginning of the analysis, setting a more stringent CPM_Cutoff = 5, and we obtained satisfying raw p-value histograms for the seven contrasts (Fig. 4B). We observed, as expected, that the highest number of differentially expressed genes is found for the comparison of both tissues with 28 Fig. 2A). An advantage of using a GLM with an interaction term is to identify straightforward genes that respond differently to the Silicon treatment in the two tissues using the [MatureLeaf_NoSi-MatureLeaf_Si]- [Root_NoSi-Root_Si] contrast. In this interaction analysis, we found 106 genes differentially expressed, and an example of a gene in this list is shown in Supplementary Fig. 2B. The hierarchical clustering on the top 50 DEGs ranking on their p-values for this contrast is also proposed by DiCoExpress (Supplementary Fig. 2C). On the bottom of this plot, we observe groups of genes with a clear opposite behavior between the two tissues. For the others, the behavior is more variable, but all these genes are declared to be the most impacted genes by the treatment and in different ways in the two tissues.
As users often need to compare DEG lists, in DiCoExpress, we propose the  3A). Within the Venn diagram, we can distinguish genes whose expression varies in response to treatment in a specific tissue or in both treatments with a similar or different behavior depending on the tissue and examples from each class is shown in Supplementary Fig. 3B. This grouping of genes using a Venn diagram is based only on the results of the single contrast differential analysis. However, by performing a co-expression analysis, we can go further in the interpretation by clustering these genes according to their average expression profile in all samples.
We applied the Coexpression_coseq function with default parameters to group the 945 DEGs from the union of the three contrasts. The convexity of the ICL curve has a clear minimum (a marker of a good quality clustering analysis), and we found seven clusters of co-expressed genes (Fig. 5). Three genes with low mean normalized counts were assigned in Cluster 0, i.e, coseq could not assign them to a cluster. Clusters 3 and 6 (71 and 106 genes, respectively) contain genes with low expression and no change in the roots, but their expression varies in response to the Si treatment in the mature leaves (over-expression in Cluster 3 and underexpression in Cluster 6). Conversely, Cluster 7 (201 genes) show low expression with no significant change in the mature leaves, but they are strongly expressed in the roots with a slight reduction following the treatment. Cluster 1 and 4 (203 and 146 genes, respectively) include genes more expressed in one tissue compared to the other one (higher expression in mature leaves for cluster 4 and higher expression in roots for cluster 1) but without significant Si treatment response. In Following the co-expression analysis, that finishes the statistical analysis of this dataset used to illustrate the use of DiCoExpress, we performed enrichment analyses on these 7 clusters, and also of the 106 DEGs for the interaction contrast, and they are available in the tutorial (Additional File 2). Further interpretation and discussion of the biology behind these enrichments are beyond the scope of our presentation of the DiCoExpress usage.

Conclusions
DiCoExpress is a user-friendly workspace for analyzing efficiently multifactorial RNAseq transcriptome experiments from quality controls to co-expression analysis through differential expression analysis. We based the development of DiCoExpress on neutral comparison studies combining the most performant statistical approaches for each step of a standard RNAseq analysis. In DiCoExpress, we used generalized linear models (GLM) implemented in the R-package edgeR for differential gene expression analysis and Gaussian mixture models implemented in the R-package coseq to perform the co-expression analysis. DiCoExpress simplifies the GLM analysis proposing automated writing of all possible contrasts and optimizes the co-expression analysis with the re-estimation of the collection of  Overview of DiCoExpress workflow DiCoExpress is composed of seven functions written in R p Histograms of raw p-values Brassica napus analysis Histograms of raw p-value for the contra