Quantifying Variability in Gene Expression Profiles using Transcriptome and Isoform Diversity

1. Introduction

Following the central dogma of molecular biology, gene expression variability can aid in predicting and explaining the wide variety of protein products, functions, and, ultimately, variability in phenotypes. Over the past few decades, transcriptomic expression profiles have been assayed in many ways, with the two most common approaches being microarray-based and sequencing-based [1]. More recently, microarrays have been surpassed by Next-Generation Sequencing (NGS), also known as second-generation or short-read sequencing [2,3]. To assay gene expression with NGS technology (also known as RNA sequencing, RNA-Seq), RNA is first reverse-transcribed into cDNA, fragmented, and then constructed into an NGS library that is then read by the sequencer and then computationally mapped to the transcriptome for quantification [3]. The ability to barcode cells in combination with low-input protocols has also enabled single-cell RNA-Seq (scRNA-Seq) to measure the transcriptomic profiles of individual cells [4]. Most recently, third-generation or long-read sequencing popularized by Pacific Biosciences [5] (PacBio) and Oxford Nanopore Technologies [6] (ONT) has allowed the sequencing of contiguous reads of up to 2.3 million bases, considerably longer than the longest human messenger RNA (mRNA) transcripts [7]. These long reads are enabled by single-molecule real-time (SMRT) and ionic nanopore technology innovations by PacBio and ONT, respectively. In addition, long-read RNA-Seq (lrRNA-Seq) is capable of directly sequencing RNA and detecting its modifications, which was previously not possible with short-read technologies that require reverse transcription into cDNA [8]. Additionally, more technologies are being developed for the rapidly growing field of long-read sequencing at the single-cell level [9,10,11].

As both these and newer technologies continue to evolve with the goal of measuring transcriptomic profiles more precisely, the need to quantify and interpret those profiles continues to grow [9,10,11]. Dating back to microarray experiments, differential expression (DE) is the most common analysis of gene expression variability, and it is frequently used in both bulk RNA-Seq and scRNA-Seq. Generally, a basic DE analysis determines whether individual genes are up- or down-regulated (i.e., more highly or lowly expressed, respectively) between conditions, for example, across tissues [12,13] or disease states [14,15]. Popular R packages for differential gene expression include DESeq2 [16], EdgeR [17], and limma [18]. DE analysis is typically done at the gene level, collapsing all counts from a sequencing library that map to a single gene unit even though many genes undergo alternative splicing (AS) to produce several different mRNA molecules or isoforms. Alternatively, DE can also be examined at the isoform level by comparing reads mapped to each specific transcript as independent entities with the aforementioned [16,17,18] or other DE packages. Additional caveats of DE are that this analysis alone may miss biological complexity and information, that DE genes are often not causal (i.e., are disease-induced [19]), or that the function of causal genes changes but mRNA expression levels remain unchanged [20], and that interpretation must also account for biases associated with the expression level or abundance of reads for a particular gene [20]. Not all measured variability in gene expression profiles reflects changes in the conditions being studied because a portion may be due to technical variability (i.e., artifacts of sample processing and sequencing preparation). These caveats are also important considerations when performing DE analysis in scRNA-Seq [21].

An additional consideration for analyzing transcriptomic profiles is the need to quantify the complexity of biological systems because nothing in biology acts in isolation. To illustrate, a perturbed gene will likely have perturbed interactions with other genes and proteins in that biological system, which in turn may also contribute to phenotypic variability. There are different ways to assess these coordinated patterns of gene expression, such as grouping genes together by pathway or function in pathway analyses (e.g., KEGG [22] and GO [23]) to determine which pathways are up- or down-regulated between conditions. While this analytical approach considers genes in aggregate by function and interactions, it does not include all of their known or predicted interactions and genes can be part of multiple pathways, which confounds the analysis. Representing gene expression profiles as biological networks presents an alternative or complementary approach to differential gene expression analysis, and such networks have been shown to be dysregulated in disease [24]. For example, weighted gene co-expression network analysis (WGCNA) [25] examines which genes have coordinated co-expression. Combined with differential expression, WGCNA analysis can reveal differential co-expression associated with a disease state [26], like systemic lupus erythematosus (SLE) severity [27]. These network biology approaches, including gene regulatory networks (GRNs) and protein-protein interaction (PPI) networks, allow a researcher to incorporate additional information like the known targets of transcription factors or protein-protein interactions, respectively, to combine multiple information types and improve in silico predicted interactions within a condition [28].

However, the above approaches fail to describe how expression diversity changes between conditions. Transcriptional diversity quantifies transcriptomic changes as a measure of the variability in 1) the total expression of all genes or a gene across samples (transcriptome diversity) or 2) the isoform-specific expression of a given gene (isoform diversity) (Figure 1). For example, such approaches have identified unique disease-related genes across 16 human disease data sets compared to differential expression alone [29]. Here we review the causes and measurement of transcriptional diversity across samples, genes, and isoforms. We first overview modulators of gene expression variability and describe its quantification across genes. Then, we discuss the role AS plays in gene expression variability and how isoform diversity can be quantified. Additionally, we provide resources for calculating transcriptomic diversity for short- and long-read sequencing data. Finally, we discuss future applications for transcriptional diversity.

2. Transcriptome Diversity in Gene Expression Profiles

2.1. Biological Processes that Lead to Transcriptome Diversity

As gene expression analyses have become a critical tool for furthering phenotypic, mechanistic, and evolutionary interpretation, it is vital to understand the forces guiding gene expression variability [30]. Previous studies noted that genetic and environmental factors are the two main drivers of biological variability in gene expression [31]. However, additional intrinsic factors like cell cycle [32], circadian rhythm [33], and aging [34,35] (which are also influenced by genetic and environmental factors) also contribute to gene expression variability.

Promoters (DNA sequences where a protein binds and initiates transcription of a gene), enhancers (DNA sequences that bind proteins to alter transcription of a gene), and transcription factors (proteins that bind DNA sequences to control the rate of transcription of a gene) are key genetic features contributing to gene expression variability [36,37]. Though they each have a specific role, in combination they drive the overall gene expression variability observed across species, tissues, and cell types [36,37]. Known for their critical role in gene transcription, promoters have also been linked to gene regulation. The heavily studied RNA polymerase II core promoter directly regulates gene expression [38,39], and natural variations in promoter regions are linked directly to both gene expression and phenotype variability [40]. For example, a genome-wide association study (GWAS) recently identified the relationship between the promoter of the gene GSE5 to grain size variability of rice, where promoter deletions in specific rice varieties led to wider grains [40]. However, genes can also have multiple promoters, and coordinated usage of these various promoters, or alternative promoter usage, can drive transcriptional variability, for example, producing different transcripts of the brain-derived neurotrophic factor (BDNF) gene [41]. Promoters have also been used to control gene expression through technologies like CRISPR/Cas9 to further improve traits or understand the functionality of specific genes [42]. Likewise, the development and use of synthetic promoter libraries, which work by mutating regions (spacer sequences) near promoter consensus sequences to affect promoter strength, have exemplified the promoter’s role in driving gene expression variability [43,44].

Unlike promoters, enhancers are not necessary for transcription to occur but can increase the transcription of a gene, making it another key gene-expression regulatory element. By regulating transcription levels distally, enhancers also influence gene expression variability within specific cell types, tissues, and even species [45,46]. Similar to promoters, alteration in an enhancer region can lead to phenotype changes by impacting gene expression [47]. Naturally occurring variation among enhancers has been associated with diseases like Alzheimer’s disease and cardiomyopathy [37,48], further underscoring how enhancers can impact gene expression and disease phenotypes. Transcription factors are essential regulatory proteins that drive gene expression by interacting with DNA sequences like promoters and enhancers to control transcriptional processes [49,50]. For example, the transcription factors Sp1 and Sp3 control the expression of the IL-10 gene, which is vital in regulating immune responses [51]. Studies like the Encyclopedia of DNA Elements (ENCODE) project, which integrated over 450 experiments of 119 transcription factors, have demonstrated that transcription factors have dynamic regulatory networks that lead to measurable variability in homeostatic gene expression [52].

Additionally, epigenetic processes including DNA methylation, histone modifications, and other environmental or stress responses, can also drive gene expression variability. DNA methylation, notably mammalian m5C (methyl groups at the 5’ cytosine of a C-G dinucleotide) [53], regulates gene expression in multiple ways [54], including through transcription factor binding, the functionality of enhancers, insulator elements, and promoters, and by altering chromatin conformation [55]. Various studies have noted correlations between gene expression and DNA methylation, further supporting its role as a possible driver for gene expression variability [54,56]. Post-translational histone modifications (e.g., acetylation, methylation, phosphorylation, or ubiquitination) are also known to be correlated with gene expression [57,58,59] and can even be used to predict gene expression [58]. Environmental and stress-related effects, like hypoxia, can also impact the variability of gene expression. Many organismal studies have observed the impact of stress on producing a biological response and subsequent regulation of various genes to alleviate environmental damages (e.g., in oxidative stress) [60,61,62]. This environmental and stress response is an integral driver of the gene expression variability seen in agricultural, homeostatic, and disease-associated transcriptomic profiles. Finally, while not the focus of this review, additional biological processes contribute to the variability of gene products, such as post-transcriptional RNA modifications like m6A methylation [63] and other post-translational protein modifications like phosphorylation and acetylation.

2.2. Methods for Quantifying Transcriptome Diversity

Researchers have applied different approaches to empirically calculate gene expression variability for both bulk and single-cell transcriptome profiles, including standard deviation [64], coefficient of variation (CV) [29], variance [65], and others [66]. While the gene expression terms variation and diversity both describe changes in gene expression across samples, variation is more frequently associated with measures of dispersion (e.g., CV, variance), and diversity is more commonly associated with these probability-based measures, particularly Shannon or information entropy (Figure 2A). In fact, the application of CV and variance to gene expression profiling analysis is sometimes known as expression variance (EV) [66]. Originally described by Alemu et al., EV showed tissue-specific variation across gene expression profiles [67] and was later used to show expression variation associated with aging and methylation [68].

Standard deviation describes the dispersion of the data in relation to its mean. For example, standard deviation has been used to determine correlations between gene expression variation and tumor differentiation among individuals with prostate cancer [64]. Building on standard deviation, CV considers the standard deviation of the gene expression sample divided by its mean, and thus is a standardized measure [29]. Therefore CV can be used to compare across conditions or datasets to identify disease-associated genes that are not identified by DE alone [29]. Additionally, other studies have applied both technical CV and biological CV (BCV) to describe RNA-Seq gene expression variation associated with technical or biological variables, respectively, as well as [17] normalized CV to examine gene expression variation, for example, across neurological diseases [65]. Recent studies have also used CV to understand how gene expression variability among therapeutic targets determines drug effectiveness and safety, thus improving therapeutic development methodologies [69]. Another empirical measurement of gene expression is variance. In the Mar et al. study, variance measures the significance of the mean difference between groups by using a t-test or ANOVA [65], but the term has also been used synonymously with gene expression variability [65,68,70,71]. For example, Bachtiary et al. applied variance (here defined as standard deviation squared) to measure the variation of expression between and within cervical cancer patient samples [72].

Though CV, variance, and standard deviation are some of the most common methods for empirically calculating variation, there are a few less commonly used ways for describing variation across gene expression. The range of gene expression observed is one of the simplest measures of variability. Though generally not used in its simplest form (i.e., maximum value minus minimum value), a modified version of range has been used. For example, dynamic range, the log10 ratio between the maximum and minimum normalized gene expression counts, has been used to compare the expression of orthologous genes between humans and mice to determine genes constrained throughout early vertebrate evolution [73] as well as to describe gene expression variation patterns across organs and tissues [74]. Additionally, researchers have developed a metric based on a ratio of the percentage of reads covering a proportion of the genome to quantify gene expression variation [75]. When a large percentage of reads covers a smaller number of total genes in the genome, it indicates lower variability in that condition than when the percentage of reads spans over a larger set of genes in another condition. However, these metrics are biased towards longer genes if gene size is not properly accounted for during analysis.

In 1948, Shannon defined entropy as the probability of uncertainty of an outcome or the amount of choice in the outcome based on how much information [76]. The basis of Information Theory, Shannon entropy, is the log of the event probability so that an event with full certainty or a probability of one would have no surprise. Over the years, Shannon entropy has been applied to numerous biological processes, including gene expression [77]. When using Shannon entropy in this context, gene expression measurements for a specific gene are the information used to measure uncertainty, or as we describe it, diversity [78]. Previous studies have employed Shannon entropy to study diversity in drug targets [79], tissue-specificity [80], species-specificity [78], and even intraspecies genomic DNA information [81]. When used to compare gene expression in RNA-Seq data, differential Shannon entropy, compared to differential CV and DE, identified genes overlapping with CV-identified genes but also included unique disease-associated genes [29], underlining that Shannon entropy can identify biological signals that CV and DE do not. Shannon entropy has also been used in combination with WGCNA analyses by calculating entropy from the betweenness of networks [82]. Additionally, studies using adaptations of Shannon entropy, such as Tsallis entropy (also known as HCDT entropy), have divided transcriptome diversity into two categories: alpha and beta diversity [83], where alpha diversity represents diversity of a single profile, and beta diversity represents diversity between samples within a group. Example analytical packages that apply entropy in the context of gene expression diversity are described in Table 1.

Altogether, the aforementioned gene expression studies demonstrate not only the importance of further understanding the drivers of this gene expression diversity but also the importance of developing new and comprehensive ways to quantify this diversity through various methodologies. Quantifying transcriptome diversity is a salient part of ascertaining how biological processes lead to phenotypic manifestations, including in a disease context. Therefore, it is imperative to examine other sources of diversity, such as variability in mRNA transcripts due to AS.

3. Isoform Diversity in Gene Expression Profiles

3.1. Biological Processes that Lead to Isoform Diversity

Before the start of the human genome project, the human genome was expected to have approximately 100,000 genes [87] based on the approximated number of protein products. However, after completing the project, the human genome actually had between 20,000 and 25,000 genes, much less than projected [88]. While humans may not have more genes than all other organisms, their splicing patterns are more specific and complex [89]. Compared to other eukaryotic organisms, humans have the highest relative splicing abundance, and this abundance steadily decreases for species with a larger evolutionary divergence from humans [89]. Because 94% of human genes undergo AS [90], most genes have a variety of transcript isoforms that, in many cases, result in proteins with unique functions, therefore increasing protein diversity. There are between six and eight types of AS events depending on the classification used, and their abundance varies by species, with exon skipping being the most common in animals and intron retention more common in plants and fungi [91]. In addition to RNA splicing, differences in transcript usage in organisms can also be driven by promoter usage (as discussed above) or 3′ end usage [92]. In this section, we will focus on the variability driven by RNA splicing.

RNA splicing occurs as part of the process to produce mature mRNA from pre-RNA and was first described in 1977 [93]. RNA splicing happens in virtually all multi-exonic genes through either constitutive splicing (splicing out an intron that is always excluded from a final transcript) or AS (variably splicing out alternative exons and/or introns resulting in diverse mRNA sequences). Exons that are always incorporated in the final or mature transcript are described as constitutive exons, while alternatively-spliced exons are those that vary in usage from transcript to transcript [94]. This AS process is performed by the spliceosome, which contains approximately 170 proteins [95], including many RNA-protein complexes called small nuclear ribonucleoproteins (snRNPs), and recognizes splice sites to facilitate the transesterification reactions that lead to intron removal (further reviewed in [96,97,98]). A key property of AS is its high specialization to a given biological condition. For example, AS is species-specific [89], and as organisms gain more evolutionarily complexity, their AS patterns become more similar to humans. AS is also sex-specific and can lead to sex-specific traits. For example, in fruit flies (Drosophila), the sex splicing gene doublesex (dsx) controls sexual differentiation [99] and is regulated by AS. In addition to more extensive sex-specific splicing in fruit flies [100,101,102], sex-specific splicing has also been shown in fish [103], birds [104], non-human primates [105], and humans [106,107]. Moreover, AS is critical in developmental changes, particularly as coordinated AS changes help define tissue identity [98]. In fact, AS has been shown to be tissue-specific [108], driven at least in part by tissue-specific splicing factors [109]. These tissue-specific splicing factors govern complex splicing regulatory networks [110] that can influence protein interaction networks and thereby increase the functional diversity of proteins [111]. Further, AS is also highly cell-type specific [112], and the recent increase in single-cell studies has highlighted an increasing number of cases of AS that are cell-type and even cell-subtype specific [113]. This has been particularly well-documented in the brain during neuronal differentiation [114] (e.g., the cerebral cortex [112]) and in immune cells [115].

Changes in AS are also associated with many diseases [116,117]. Currently, an estimated 15% of human disease-causing point mutations result in an AS defect [118], and many diseases and disorders are associated with disrupted splicing patterns, like spinal muscular atrophy, cancer, and autism spectrum disorder (ASD) [119,120]. Because of all of these known changes in AS across numerous biological conditions, some pathogenic and others benign, it is critical for genomic researchers to quantify changes in AS using RNA sequencing.

3.2. Methods for Quantifying Isoform Diversity

Measuring alternatively-spliced transcript expression diversity requires first identifying and then quantifying transcripts from RNA-Seq data. One way to quantify alternatively spliced transcripts from gene expression data relies on identifying reads that cover splice junctions, genomic locus where two exons have been spliced together [121]. This process varies depending on the transcript quantification tool, as some tools only count if an exon is included at all (i.e., if any reads map to that exon), while others search for junctions to determine if an exon is actually spliced in, because reads mapping to a free exon (i.e., not including a junction) cannot resolve where in the transcript that exon has been spliced. A major limitation is that there must be sufficient read depth to detect all splice junctions from short-read data [122]. In some cases, junctions are specific to unique transcripts, so a read mapping to a unique junction could indicate that that transcript is being expressed without having any continuous reads capturing the entire transcript. One of these splice-junction-based methods is Splice Expression Variation Analysis (SEVA) [123], which compares the variability of the multivariate distribution of splice junction expression profiles between conditions. On the other hand, because short reads with few junctions usually do not match a unique transcript, probabilistic methods can be used to estimate exon inclusion [124] (also known as percent spliced in, PSI) [125], greatly reducing the precision of transcript quantification.

Differential transcript (or isoform) expression (DTE) is a data analysis method similar to differential gene expression as it identifies up- or down-expression of specific transcripts in one condition versus another. However, differential transcript usage (DTU), sometimes referred to as differential isoform usage, isoform switching, or differential splicing [126], can determine differential proportions of transcript expression within a gene across conditions [127]. Soneson et al. describe three methods for DTU: assembly-based, type of AS-based, and differential exon usage [127]. Due to limitations with short reads not covering entire transcripts that overlap, differential exon usage (DEU) can also be used in a similar way to measure shifts in functional unit expression (i.e., bins) across conditions, usually comparing PSI values [127]. Table 2 includes analysis packages comparing transcript expression across conditions and used for DTU, DEU (PSI), or other analyses.

However, short-read sequencing technology often fails to adequately resolve transcripts because the typical mRNA is over 1kb, whereas most short-read RNA-Seq data is only 100-200 bases in length [3], i.e., only long enough to cover an exon or less. The advent of long-read technologies has created an opportunity to capture more detailed gene expression profiles, especially for resolving transcript expression, but also for accurate sequencing and subsequent mapping of repetitive, hard-to-map, and/or duplicated gene regions [128]. Additionally, as lrRNA-Seq approaches can sequence full-length novel transcript isoforms, they are continuing to identify novel transcripts, including those that are lowly expressed [129]. While most of the short-read tools for transcript quantification can be used on long-read data, there are transcript quantification tools that are specialized for long-read sequencing data, like FLAIR [130] and BAMBU [131], which include steps for correcting misalignments that can result from less accurate reads. Applying variance and entropy-based diversity quantification approaches in combination with these lrRNA-Seq technologies, therefore captures transcriptomic changes across biological conditions and phenotypes.

Additionally, isoform diversity can be described by enumerating the total number of isoforms [132,133]. In contrast to only counting the number of transcripts, the distribution of isoform expression for a gene, such as in DTU, can also be considered. Similarly to the gene expression level, variance and Shannon entropy can be used to describe transcript isoform diversity (Figure 2) [134]. One way to measure isoform diversity is the Fano factor, or the squared variance over the mean [135], which describes the distribution of alternatively spliced transcripts while adjusting for the mean expression of that gene. Another method is by Shannon entropy, where a gene with many isoforms could be less diverse than a gene with few isoforms if the latter has a more even distribution of expression and perhaps equal usage of those gene products. Figure 2B provides an example of isoform diversity quantified with Shannon entropy.

The first instance of using Shannon entropy to describe diversity in three types of alternative transcription (AS, polyadenylation, and transcription initiation) using targeted microarray expression data was by Ritchie et al. in 2008 [134]. Ritchie et al.'s rationale was that Shannon entropy could capture aberrant transcription seen in cancer and was therefore used to compare patient cancerous tissue transcriptomic profiles with non-cancerous tissue transcriptomic profiles. The authors found that out of the three types of transcription studied, only AS had increased diversity in cancer tissues. They concluded that these changes in entropy are unlikely to reflect changes in gene function because they found it unlikely that, in a cancer context, shifts in isoform expression are functional or controlled. This general approach to measuring entropy has also been used to compare transcript diversity across conditions in the brain [128] and epithelial cells [136].

Several software approaches have been developed for quantifying isoform diversity with Shannon entropy, including Cuffdiff (from Cufflinks) [137,138,139], Whippet [140], and SplicingFactory [141]. Cuffdiff uses Jensen-Shannon divergence, which, like Shannon entropy, relies on probability to compare the distribution of transcript expression across conditions [139]. Whippet [140] applies Shannon entropy to define the entropy of individual AS events instead of at the gene level, meaning that each alternatively spliced exon is given a value based on PSI. SplicingFactory [141] is unique because it also includes multiple other methods for assaying diversity across isoforms of the same gene, like the Gini Index (originally developed for describing wealth inequalities) and the Simpson Index (originally developed to measure ecological diversity). The aforementioned Tsallis entropy could also potentially be used in an isoform context to describe biological variability since it has been shown to provide more information than Shannon entropy or Simpson Index alone by including a parameter (q) that can be manipulated to adjust the weight of highly-expressed genes [83], though it has not yet been applied to study alternatively-spliced isoform distributions.

To summarize, AS contributes to the gene expression diversity observed by increasing the number of products that can be produced by a single gene. Quantifying the isoform diversity of a given gene will not only identify which isoforms are highly expressed but how isoform expression shifts between conditions.

4. Conclusions

In conclusion, applying transcriptomic diversity when analyzing gene expression profiles is an integral analysis step to capturing biological information in tandem with the differential expression of individual genes. Critically, omitting isoform diversity from gene expression analysis can miss biological information since phenotypic diversity is partially driven by AS [91]. Quantifying transcriptional changes as a measure of the variability in 1) the total expression of all genes or a gene across samples (transcriptome diversity) or 2) the isoform-specific expression of a given gene (isoform diversity) are separate but complementary analyses. For example, transcriptome diversity can be quantified at the gene profile level across tissue samples in a population in order to identify the diversity and specificity of gene expression in tissue types [77]. Alternatively, if a research question relies on studying isoform usage and how specific transcript expression may be in a particular context, isoform diversity should be quantified. Combining both applications of diversity to better understand high-dimensional gene expression data can provide insight into the transcriptional differences between contexts and reveal gene expression differences that traditional DE analyses cannot like novel therapeutic targets [29].

This review provides an overview of the drivers of gene expression variability and diversity and describes how it has been quantified across all genes and at the isoform level. Additionally, we summarize resources for calculating diversity and variability. Applying additional variation and diversity measurements in transcriptomic analysis has the potential to capture additional gene expression profile changes between conditions and, in the future, could be adapted to additional gene expression profile analysis, such as spatial transcriptomic data. Clinically, these analyses may also elucidate new therapeutic targets for hard-to-treat disorders, further underscoring the importance of transcriptome diversity.