Differential expression analysis using a model-based gene clustering algorithm for RNA-seq data

Abstract

Background: RNA-seq is a tool for measuring gene expression and is commonly used to identify differentially expressed genes (DEGs). Gene clustering is used to classify DEGs with similar expression patterns for the subsequent analyses of data from experiments such as time-courses or multi-group comparisons. However, gene clustering has rarely been used for analyzing simple two-group data or differential expression (DE). In this study, we report that a model-based clustering algorithm implemented in an R package, MBCluster.Seq, can also be used for DE analysis.

Results: The input data originally used by MBCluster.Seq is DEGs, and the proposed method (called MBCdeg) uses all genes for the analysis. The method uses posterior probabilities of genes assigned to a cluster displaying non-DEG pattern for overall gene ranking. We compared the performance of MBCdeg with conventional R packages such as edgeR, DESeq2, and TCC that are specialized for DE analysis using simulated and real data. Our results showed that MBCdeg outperformed other methods when the proportion of DEG (P_DEG) was less than 50%. However, the DEG identification using MBCdeg was less consistent than with conventional methods. We compared the effects of different normalization algorithms using MBCdeg, and performed an analysis using MBCdeg in combination with a robust normalization algorithm (called DEGES) that was not implemented in MBCluster.Seq. The new analysis method showed greater stability than using the original MBCdeg with the default normalization algorithm.

Conclusions: MBCdeg with DEGES normalization can be used in the identification of DEGs when the P_DEG is relatively low. As the method is based on gene clustering, the DE result includes information on which expression pattern the gene belongs to. The new method may be useful for the analysis of time-course and multi-group data, where the classification of expression patterns is often required.

Keywords: Differential expression; Gene clustering; Posterior probability; RNA-seq.

Fig. 1

Results for two-group simulated data (P_DEG ≤ 0.25). Boxplots of AUC values (100 trials) for five methods for a total of eight conditions, P₁ = 0.5 (left) to 1.0 (right) with P_DEG = 0.05 (upper) and 0.25 (lower). The performance of MBCdeg (with K = 3) was high in most trials. The explanation for a representative trial (AUC = 0.9295) and the worst trial (AUC = 0.633) using MBCdeg2 with P_DEG = 0.25 and P₁ = 0.5 is given in Table 1

Fig. 2

Effect of the different cluster numbers for MBCdeg (K = 2–4). Boxplots of AUC values (100 trials) for MBCdeg1 (K = 2–4; colored in black) and MBCdeg2 (K = 2–4; colored in red) are shown. The remaining legends are the same as in Fig. 1. The AUC values for MBCdeg with K = 3 were similar to those in Fig. 1 (different trials were used). Note that two trials result in AUC < 0.4 (not shown): MBCdeg1 with K = 2 with P_DEG = 0.05 and P₁ = 0.9 (AUC = 0.3498), and MBCdeg1 with K = 4 with P_DEG = 0.25 and P₁ = 1.0 (AUC = 0.3893)

Fig. 3

Effect on larger P_DEG values (P_DEG ≥ 0.45). Boxplots of AUC values (50 trials) for five methods for a total of 16 conditions, P₁ = 0.5 (left) to 1.0 (right) with P_DEG = 0.45 (top) to 0.75 (bottom). MBCdeg1 and MBCdeg2 were used for analysis with K = 3. The performance of MBCdeg2 partly depended on the of DEGES normalization (see Tables 2 and 3)

Fig. 4

Results for three-group simulated data. Boxplots of AUC values (50 trials) for eleven methods (DESeq2, edgeR, TCC, and MBCdeg with K = 2–5) are shown. The simulation conditions used were P_DEG = 0.25, FC = 4, and n₁ = n₂ = n₃ = 3