RNA-seq is a revolutionary technology to measure genome-wide gene expression of biological samples at high resolution by sequencing messenger RNA (mRNA) molecules [1]. Common usages of RNA-Seq technology are computing transcription abundances [2], finding differentially expressed genes between two or more groups [3], de novo transcriptome assembly [4, 5] and finding novel genes and splicing patterns [6]. Among these usages, differential gene expression (DGE) analysis followed by functional enrichment is a common workflow in gene expression studies [2, 7, 8].
After RNA-seq reads are generated using a sequencing instrument, gene expression abundance is estimated by mapping the sequencing reads to a reference genome if there is an available reference genome or by building a transcriptome assembly de novo [9, 10]. DGE analysis is performed to compute statistically significant differentially expressed (DE) genes using tools such as DESeq2 [3], edgeR [11], limma [12] and Cufflinks [2]. DGE analysis could result thousands of genes, thus to better characterize the underlying biological functions of the DE genes, functional enrichment analysis is performed using tools such as GOSeq [8] and SeqGSEA [13].
However, particularly when biological samples are not well separated (e.g., mammalian tissue or human disease samples are often heterogeneous or heterocellular), a direct two-group DGE analysis can result in unmanageable lists of DE genes with uncertain significance [14]. Furthermore, batch effects may obscure the experimental signal or sample mishandling may generate outliers that perturb the experimental signal in ways unnoticed by the investigator.
In these scenarios, list of DE genes and even significantly enriched biological processes would be hard to interpret for biologists. Alternatively, after computing significantly enriched biological processes, samples could be visualized based on their activity for each of these biological processes. Per biological process visualization would enable biologists to have a deeper understanding of the samples’ activity with respect to each significant biological process.
To streamline the analysis of RNA-seq datasets to achieve the aforementioned goals, we developed a software toolkit GSEPD (gene set enrichment and projection display). GSEPD produces DE gene lists, significantly enriched gene ontology (GO) terms, and importantly their cross-product: a mapping of which genes are perturbed within each GO term, and how genes associated with those terms define the samples’ expression profiles in the context of the other RNA-Seq samples. GSEPD provides various plots and tables to summarize the results and give its users a comprehensive outlook of the underlying RNA-seq data.
We demonstrated the usage of GSEPD on a time series dataset of H1ESC cell lines [15]. GSEPD computed DE genes and significantly enriched GO terms between two time points, and clustered samples from all time points based on their activity in each significant GO term.
GSEPD is implemented as a Bioconductor package named rgsepd and freely available under GPL-3 license.
Implementation
We built GSEPD as a Bioconductor package named rgsepd to ensure that it is readily available, simple to install, and bundled with both test data and documentation. The system architecture of GSEPD is shown in Fig. 1. The interface to GSEPD is a short list of R commands and all the functions are fully automatic after providing the input data as a matrix. GSEPD generates all tables and figures for the input data within minutes.
GSEPD requires two types of input data to run: the multisample RNA-seq raw counts matrix and sample information matrix. Input should be loaded as a matrix in R with RefSeq ID numbers as row and sample identifiers as column names. The sample information matrix is used to link sample identifiers with test conditions and short labels (for plotting into figures). Given input data, GSEPD automatically computes DE genes between two groups with default parameters of DESeq2, adjusted if necessary for small sample counts [3]. GSEPD also utilizes GOSeq [8] for GO term enrichment analysis, once each for downregulated, upregulated and all genes in the DE gene list.
One of the novel features of GSEPD is to focus on each significantly enriched GO term and assess how samples are segregated with respect to the expression of genes in that GO term. In order to study if samples segregate in their original groups with respect to a particular GO term, GSEPD performs clustering of samples based on the expression of all genes in a significantly enriched GO term. GSEPD can also incorporate non-tested samples (i.e., samples that are not in the predefined groups) in clustering to enable investigators label unclassified or indeterminate samples by their expression profiles among GO terms relevant to the experiment.
GO term-based clustering of samples is performed by using k-means clustering where k = 2. Briefly, for a given GO term with N genes, each sample is represented as an N-dimensional vector of expression of all genes in the GO term. To avoid broad GO terms associated with thousands of genes, only GO terms with less than m (m = 31 by default) genes are evaluated by GSEPD for clustering.
To assess the quality of the clustering outcome, a validity score called V-measure [16] is computed. The V-measure computes the concordance between cluster assignments and actual class labels of the samples. The V-measure of a clustering is the harmonic mean of the cluster’s homogeneity and completeness. A cluster’s homogeneity is computed based on the entropy of class labels within the cluster, i.e., maximum homogeneity is achieved when all members of the cluster belongs to the same cluster. The completeness of a cluster is computed based on what percent of members of a class are assigned to the cluster. A cluster would have maximum completeness if it has all members of a class. In ideal cases, clusters should be homogenous and complete.
In order to assess the significance of V-measure score, GSEPD computes an empirical p-value for each GO term-based clustering by permuting sample group labels (i.e., class labels) and re-calculating the V-measure. The p-value is the proportion of random assignments that achieve a higher V-measure. By default, GSEPD performs adaptive permutation up to 400 times to resolve segregation by p < 0.01.
GSEPD visualizes significant GO terms in scatter plots and subspace principle component analysis (PCA) figures to allow further exploration of the results by the user. Vector projection of samples is performed based on gene set of the GO term to score each sample’s similarity to the centroid of each group and to highlight any outlier samples for the gene set.
In order to assess the concordance between group label of a sample and its localization in the clustering, GSEPD performs vector projection. First, we define the mean expression of the GO term gene set in samples of each group as the centroid of the group, and define an axis connecting both group centroids where one of the centroids is chosen as the origin in a N-dimensional Euclidean space (Fig. 2). Each sample is projected on this axis to compute two scores named alpha and beta. The alpha score is the distance between projected point on the axis to the origin and the beta score is the Euclidean distance between the sample and the projected point in the axis (Fig. 2). Beta score measures the goodness of fit and flag samples which do not fit the linear assumptions of the two-group comparison performed by DESeq2 whereas alpha score is used to measure the confidence of the cluster assignments. Alpha and beta scores are computed for samples from other groups and can help assess how samples from other groups “behave” for a given GO term.
GO term-based clustering and vector projection is performed for each significant GO term with gene sets ≤ m, creating an alpha and beta score for each sample and GO term pair. GSEPD produces heatmaps of gene expression for DE genes, heatmaps of alpha scores for significant GO terms, multi-panel scatterplots of genes in significant GO terms, PCA plots of samples and tables. All thresholds and parameters are configurable before runtime, and configurable output folders and formulaic file naming conventions ensure easy reproducibility or automated parameter sweeps. A tutorial and explanation of all outputs are available within the package vignette/manuals.
