Gene Set Enrichment Analysis in RNA-Seq Data

Longitudinal gene expression profiles of subjects are collected in some clinical studies to monitor disease progression and understand disease etiology. The identification of gene sets that have coordinated changes with relevant clinical outcomes over time from these data could provide significant insights into the molecular basis of disease progression and lead to better treatments. In this article, we propose a Distance-Correlation based Gene Set Analysis (dcGSA) method for longitudinal gene expression data. dcGSA is a non-parametric approach, statistically robust, and can capture both linear and nonlinear relationships between gene sets and clinical outcomes. In addition, dcGSA is able to identify related gene sets in cases where the effects of gene sets on clinical outcomes differ across subjects due to the subject heterogeneity, remove the confounding effects of some unobserved time-invariant covariates, and allow the assessment of associations between gene sets and multiple related outcomes simultaneously. Through extensive simulation studies, we demonstrate that dcGSA is more powerful of detecting relevant genes than other commonly used gene set analysis methods. When dcGSA is applied to a real dataset on systemic lupus erythematosus, we are able to identify more disease related gene sets than other methods.

Cancer initiation and progression are likely caused by the dysregulation of biological pathways. Gene set analysis (GSA) could improve the signal-to-noise ratio and identify potential biological insights on the gene set level. However, platforms exploring cancer multi-omics data using GSA methods are lacking. In this study, we upgraded our GSCALite to GSCA (gene set cancer analysis, http://bioinfo.life.hust.edu.cn/GSCA) for cancer GSA at genomic, pharmacogenomic and immunogenomic levels. In this improved GSCA, we integrated expression, mutation, drug sensitivity and clinical data from four public data sources for 33 cancer types. We introduced useful features to GSCA, including associations between immune infiltration with gene expression and genomic variations, and associations between gene set expression/mutation and clinical outcomes. GSCA has four main functional modules for cancer GSA to explore, analyze and visualize expression, genomic variations, tumor immune infiltration, drug sensitivity and their associations with clinical outcomes. We used case studies of three gene sets: (i) seven cell cycle genes, (ii) tumor suppressor genes of PI3K pathway and (iii) oncogenes of PI3K pathway to prove the advantage of GSCA over single gene analysis. We found novel associations of gene set expression and mutation with clinical outcomes in different cancer types on gene set level, while on single gene analysis level, they are not significant associations. In conclusion, GSCA is a user-friendly web server and a useful resource for conducting hypothesis tests by using GSA methods at genomic, pharmacogenomic and immunogenomic levels.

The focus of analyzing data from microarray experiments has shifted from the identification of associated individual genes to that of associated biological pathways or gene sets. In bioinformatics, a feature selection algorithm is usually used to cope with the high dimensionality of microarray data. In addition to those algorithms that use the biological information contained within a gene set as a priori to facilitate the process of feature selection, various gene set analysis methods can be applied directly or modified readily for the purpose of feature selection. Significance analysis of microarray to gene-set reduction analysis (SAM-GSR) algorithm, a novel direction of gene set analysis, is one of such methods. Here, we explore the feature selection property of SAM-GSR and provide a modification to better achieve the goal of feature selection. In a multiple sclerosis (MS) microarray data application, both SAM-GSR and our modification of SAM-GSR perform well. Our results show that SAM-GSR can carry out feature selection indeed, and modified SAM-GSR outperforms SAM-GSR. Given pathway information is far from completeness, a statistical method capable of constructing biologically meaningful gene networks is of interest. Consequently, both SAM-GSR algorithms will be continuously revaluated in our future work, and thus better characterized.

Gene-set analysis (GSA) aims to identify sets of differentially expressed genes by a phenotype in DNA microarray studies. Challenges occur due to the salient characteristics of the data: (1) the number of genes is far larger than the number of observations; (2) gene expression measurements, especially within each gene set, can be highly correlated; and (3) the number of gene sets that can be examined is large and increasing rapidly. These challenges call for gene-set testing procedures that have both efficiency in computation for large GSAs and high power in the presence of the high correlation. We propose a new GSA approach called Linear Combination Test (LCT), incorporating the covariance matrix estimator of gene expression into the test statistic. The proposed LCT and two other GSA methods, a modification of Hotelling's T2 using a shrinkage covariance matrix and our SAM-GS (Dinu et. al. 2007), the two methods that have been reported by Tsai and Chen (2009) to perform best in terms of power, are evaluated in simulation studies and a real microarray study. The LCT method is more computationally efficient than the modified Hotelling's T2 and approximates the superb power of the modified Hotelling's T2. LCT is slightly faster than SAM-GS, but more powerful, due to incorporating the covariance matrix estimator. An extra step to enhance the interpretation of GSA results is also proposed in the form of a hierarchical LC (HLC) testing procedure, providing scientists useful hierarchical information on gene sets that LCT identified as differentially expressed.

Gene set analysis has been widely used to gain insight from high-throughput expression studies. Although various tools and methods have been developed for gene set analysis, there is no consensus among researchers regarding best practice(s). Most often, evaluation studies have reported contradictory recommendations of which methods are superior. Therefore, an unbiased quantitative framework for evaluations of gene set analysis methods will be valuable. Such a framework requires gene expression datasets where enrichment status of gene sets is known a priori. In the absence of such gold standard datasets, artificial datasets are commonly used for evaluations of gene set analysis methods; however, they often rely on oversimplifying assumptions that make them biased in favor of or against a given method. In this paper, we propose a quantitative framework for evaluation of gene set analysis methods by synthesizing expression datasets using real data, without relying on oversimplifying or unrealistic assumptions, while preserving complex gene–gene correlations and retaining the distribution of expression values. The utility of the quantitative approach is shown by evaluating ten widely used gene set analysis methods. An implementation of the proposed method is publicly available. We suggest using Silver to evaluate existing and new gene set analysis methods. Evaluation using Silver provides a better understanding of current methods and can aid in the development of gene set analysis methods to achieve higher specificity without sacrificing sensitivity.

Differential gene expression (DGE) studies often suffer from poor interpretability of their primary results, i.e., thousands of differentially expressed genes. This has led to the introduction of gene set analysis (GSA) methods that aim at identifying interpretable global effects by grouping genes into sets of common context, such as, molecular pathways, biological function or tissue localization. In practice, GSA often results in hundreds of differentially regulated gene sets. Similar to the genes they contain, gene sets are often regulated in a correlative fashion because they share many of their genes or they describe related processes. Using these kind of neighborhood information to construct networks of gene sets allows to identify highly connected sub-networks as well as poorly connected islands or singletons. We show here how topological information and other network features can be used to filter and prioritize gene sets in routine DGE studies. Community detection in combination with automatic labeling and the network representation of gene set clusters further constitute an appealing and intuitive visualization of GSA results. The RICHNET workflow described here does not require human intervention and can thus be conveniently incorporated in automated analysis pipelines.

Using Gene-Set Analysis To Gain Biological Knowledge Based On GWAS Results

In recent years, RNA-seq has become a very competitive alternative to microarrays. In RNA-seq experiments, the expected read count for a gene is proportional to its expression level multiplied by its transcript length. Even when two genes are expressed at the same level, differences in length will yield differing numbers of total reads. The characteristics of these RNA-seq experiments create a gene-level bias such that the proportion of significantly differentially expressed genes increases with the transcript length, whereas such bias is not present in microarray data. Gene-set analysis seeks to identify the gene sets that are enriched in the list of the identified significant genes. In the gene-set analysis of RNA-seq, the gene-level bias subsequently yields the gene-set-level bias that a gene set with genes of long length will be more likely to show up as enriched than will a gene set with genes of shorter length. Because gene expression is not related to its transcript length, any gene set containing long genes is not of biologically greater interest than gene sets with shorter genes. Accordingly the gene-set-level bias should be removed to accurately calculate the statistical significance of each gene-set enrichment in the RNA-seq.We present a new gene set analysis method of RNA-seq, called FDRseq, which can accurately calculate the statistical significance of a gene-set enrichment score by the grouped false-discovery rate. Numerical examples indicated that FDRseq is appropriate for controlling the transcript length bias in the gene-set analysis of RNA-seq data. To implement FDRseq, we developed the R program, which can be downloaded at no cost from http://home.mju.ac.kr/home/index.action?siteId=tyang.

Genome-wide association studies (GWAS) for epithelial ovarian cancer (EOC), the most lethal gynecologic malignancy, have identified novel susceptibility loci. GWAS for survival after EOC have had more limited success. The association of each single-nucleotide polymorphism (SNP) individually may not be well suited to detect small effects of multiple SNPs, such as those operating within the same biologic pathway. Gene set analysis (GSA) overcomes this limitation by assessing overall evidence for association of a phenotype with all measured variation in a set of genes. To determine gene sets associated with EOC overall survival, we conducted GSA using data from two large GWAS (N cases = 2,813, N deaths = 1,116), with a novel Principal Component-Gamma GSA method. Analysis was completed for all cases and then separately for high-grade serous histologic subtype. Analysis of the high-grade serous subjects resulted in 43 gene sets with P < 0.005 (1.7%); of these, 21 gene sets had P < 0.10 in both GWAS, including intracellular signaling pathway (P = 7.3 × 10(-5)) and macrolide binding (P = 6.2 × 10(-4)) gene sets. The top gene sets in analysis of all cases were meiotic mismatch repair (P = 6.3 × 10(-4)) and macrolide binding (P = 1.0 × 10(-3)). Of 18 gene sets with P < 0.005 (0.7%), eight had P < 0.10 in both GWAS. This research detected novel gene sets associated with EOC survival. Novel gene sets associated with EOC survival might lead to new insights and avenues for development of novel therapies for EOC and pharmacogenomic studies.

Microarray, developed at the end of the last century, is one of the most significant achievements in biotechnology history. One of the most important applications of microarrays is to uncover the signature genes characterizing different samples or biological states. This application was initially used in individual gene analysis, but recently has been applied to gene set analysis. We compare these two approaches and compare the performance using JProGO (Java Tool for the Functional Analysis of Prokaryotic Microarray Data using the Gene Ontology). JProGO is able to examine microarray data from more than 20 different samples of prokaryotic and thus practical analysis of gene behavior data with regards to Gene Ontology (GO) is made possible. Two Escherichia coli datasets are used to compare the performance of the two analyses. The results show that the gene set analysis is more powerful because of its ability to detect the GO terms according to biological expectation. This research will provide a new platform and opportunities for further research or studies in microarray based gene sets.

Gene set analysis methods have played a major role in generating biological interpretations from omics data such as gene expression datasets. However, most methods focus on detecting homogenous pattern changes in mean expression and methods detecting pattern changes in variance remain poorly explored. While a few studies attempted to use gene-level variance analysis, such approach remains under-utilized. When comparing two phenotypes, gene sets with distinct changes in subgroups under one phenotype are overlooked by available methods although they reflect meaningful biological differences between two phenotypes. Multivariate sample-level variance analysis methods are needed to detect such pattern changes. We use ranking schemes based on minimum spanning tree to generalize the Cramer-Von Mises and Anderson-Darling univariate statistics into multivariate gene set analysis methods to detect differential sample variance or mean. We characterize these methods in addition to two methods developed earlier using simulation results with different parameters. We apply the developed methods to microarray gene expression dataset of prednisolone-resistant and prednisolone-sensitive children diagnosed with B-lineage acute lymphoblastic leukemia and bulk RNA-sequencing gene expression dataset of benign hyperplastic polyps and potentially malignant sessile serrated adenoma/polyps. One or both of the two compared phenotypes in each of these datasets have distinct molecular subtypes that contribute to heterogeneous differences. Our results show that methods designed to detect differential sample variance are able to detect specific hallmark signaling pathways associated with the two compared phenotypes as documented in available literature. The results in this study demonstrate the usefulness of methods designed to detect differential sample variance in providing biological interpretations when biologically relevant but heterogeneous changes between two phenotypes are prevalent in specific signaling pathways. Software implementation of the developed methods is available with detailed documentation from Bioconductor package GSAR. The available methods are applicable to gene expression datasets in a normalized matrix form and could be used with other omics datasets in a normalized matrix form with available collection of feature sets.

IntroductionInfusion of glyceryl trinitrate (GTN), a donor of nitric oxide, induces immediate headache in humans that in migraineurs is followed by a delayed migraine attack. In order to achieve increased knowledge of mechanisms activated during GTN-infusion this present study aims to investigate transcriptional responses to GTN-infusion in the rat trigeminal ganglia.MethodsRats were infused with GTN or vehicle and trigeminal ganglia were isolated either 30 or 90 minutes post infusion. RNA sequencing was used to investigate transcriptomic changes in response to the treatment. Furthermore, we developed a novel method for Gene Set Analysis Of Variance (GSANOVA) to identify gene sets associated with transcriptional changes across time.Results15 genes displayed significant changes in transcription levels in response to GTN-infusion. Ten of these genes showed either sustained up- or down-regulation in the 90-minute period after infusion. The GSANOVA analysis demonstrate enrichment of pathways pointing towards an increase in immune response, signal transduction, and neuroplasticity in response to GTN-infusion. Future functional in-depth studies of these mechanisms are expected to increase our understanding of migraine pathogenesis.

Background: Bone metastasis-specific genes in breast cancer have been reported without considering the significant differences in ER status between bone and non-bone metastases. The aims of this study were to validate genes that had been reported as bone metastasis-specific genes using our data set and to identify bone metastasis-specific genes on the basis of biological differences between ERpositive and ER-negative breast cancers. Patients and Methods: We used Affymetrix GeneChip arrays to analyze tumor samples obtained from 365 primary invasive breast cancer patients who underwent surgery from 1999 to 2008. We excluded patients with HER2−positive breast cancer (normalized HER2 mRNA expression [probe set 216836_s_at] &amp;gt; 12.54). We classified the samples into 3 cohorts according to first metastatic site: bone, non-bone, or no metastasis. Differential expression of genes between bone and non-bone cohorts that were differentially expressed was identified using the Cox proportional hazards model, and gene sets was assessed using gene-set analysis. Results: Of the 365 patients, 34 (9.3%) were included in the bone cohort and 32 (8.8%) in the non-bone cohort. Two hundred fourteen (58.6%) had ER-positive and 151 (41.4%) had ER-negative breast cancer. First, we performed gene-set analysis using 5 gene sets that had been reported to be associated with bone metastasis. One gene set, which had been detected using an ER-negative breast cancer cell line, was validated as predicting bone metastasis in ER-positive breast cancer. None of the 5 gene sets predicted bone metastasis in ER-negative breast cancer. We then determined the levels of individual genes associated with bone metastasis by ER status using all 16,712 probe sets filtered by average gene expression level. When we analyzed all patients without any stratification by ER status, as in previous studies, 592 probe sets were significantly overexpressed in the bone cohort compared with the non-bone cohort, with a false discovery rate of ≤0.05. However, when we analyzed ER-positive and ER-negative breast cancers separately, no genes were found with significant differences between bone and non-bone cohorts. Finally, we used 2,246 functionally annotated gene sets assembled from Gene Ontology to examine possible biological differences between bone and non-bone cohorts. In the bone cohort, 151 and 125 gene sets were significantly overexpressed in ER-positive and ER-negative breast cancers, respectively (P ≤ 0.05). Ppathways related to Cellular growth and proliferation, and intracellular and second-messenger signaling were overexpressed in ER-positive breast cancer, whereas pathways related to nuclear receptor and cytokine signaling were overexpressed in ER-negative breast cancer. Most bone-metastasis-related pathways were different in ER-positive and ER-negative breast cancers (91.4% and 89.6% of the gene sets, respectively). Discussion: No genes were found that can predict bone metastasis. ER-positive and ER-negative breast cancers have different biological potentials for bone metastasis. Therefore, we need to assess the prediction model of bone metastasis based on the biological features for each ER status separately. Citation Information: Cancer Res 2011;71(24 Suppl):Abstract nr P4-16-02.

We present a novel multiple testing method for testing null hypotheses that are structured in a directed acyclic graph (DAG). The method is a top-down method that strongly controls the familywise error rate and can be seen as a generalization of Meinshausen's procedure for tree-structured hypotheses. Just as Meinshausen's procedure, our proposed method can be used to test for variable importance, only the corresponding variable clusters can be chosen more freely, because the method allows for multiple parent nodes and partially overlapping hypotheses. An important application of our method is in gene set analysis, in which one often wants to test multiple gene sets as well as individual genes for their association with a clinical outcome. By considering the genes and gene sets as nodes in a DAG, our method enables us to test both for significant gene sets as well as for significant individual genes within the same multiple testing procedure. The method will be illustrated by testing Gene Ontology terms for evidence of differential expression in a survival setting and is implemented in the R package cherry.

Protein kinase C (PKC) isozymes are key signal transducers involved in normal physiology and disease and have been widely implicated in cancer progression. Despite our extensive knowledge of the signaling pathways regulated by PKC isozymes and their effectors, there is essentially no information on how individual members of the PKC family regulate gene transcription. Here, we report the first PKC isozyme-specific analysis of global gene expression by microarray using RNAi depletion of diacylglycerol/phorbol ester-regulated PKCs. A thorough analysis of this microarray data revealed unique patterns of gene expression controlled by PKCα, PKCδ, and PKCε, which are remarkably different in cells growing in serum or in response to phorbol ester stimulation. PKCδ is the most relevant isoform in controlling the induction of genes by phorbol ester stimulation, whereas PKCε predominantly regulates gene expression in serum. We also established that two PKCδ-regulated genes, FOSL1 and BCL2A1, mediate the apoptotic effect of phorbol esters or the chemotherapeutic agent etoposide in prostate cancer cells. Our studies offer a unique opportunity for establishing novel transcriptional effectors for PKC isozymes and may have significant functional and therapeutic implications.