- Research article
- Open Access
Bimodal distribution of RNA expression levels in human skeletal muscle tissue
BMC Genomics volume 12, Article number: 98 (2011)
Many human diseases and phenotypes are related to RNA expression, levels of which are influenced by a wide spectrum of genetic and exposure-related factors. In a large genome-wide study of muscle tissue expression, we found that some genes exhibited a bimodal distribution of RNA expression, in contrast to what is usually assumed in studies of a single healthy tissue. As bimodality has classically been considered a hallmark of genetic control, we assessed the genome-wide prevalence, cause, and association of this phenomenon with diabetes-related phenotypes in skeletal muscle tissue from 225 healthy Pima Indians using exon array expression chips.
Two independent batches of microarrays were used for bimodal assessment and comparison. Of the 17,881 genes analyzed, eight (GSTM1, HLA-DRB1, ERAP2, HLA-DRB5, MAOA, ACTN3, NR4A2, and THNSL2) were found to have bimodal expression replicated in the separate batch groups, while 24 other genes had evidence of bimodality in only one group. Some bimodally expressed genes had modest associations with pre-diabetic phenotypes, of note ACTN3 with insulin resistance. Most of the other bimodal genes have been reported to be involved with various other diseases and characteristics. Association of expression with cis genetic variation in a subset of 149 individuals found all but one of the confirmed bimodal genes and nearly half of all potential ones to be highly significant expression quantitative trait loci (eQTL). The rare prevalence of these bimodally expressed genes found after controlling for batch effects was much lower than the prevalence reported in other studies. Additional validation in data from separate muscle expression studies confirmed the low prevalence of bimodality we observed.
We conclude that the prevalence of bimodal gene expression is quite rare in healthy muscle tissue (<0.2%), and is much lower than limited reports from other studies. The major cause of these clearly bimodal expression patterns in homogeneous tissue appears to be cis-polymorphisms, indicating that such bimodal genes are, for the most part, eQTL. The high frequency of disease associations reported with these genes gives hope that this unique feature may identify or actually be an underlying factor responsible for disease development.
RNA expression levels indicate the overall effect of a combination of genetic and environmental factors [1, 2], which may result in distinct phenotypic differences, including disease susceptibility [3–5]. Attempts to unravel these genetic and environmental roles have met with only limited success for some complex diseases . The emergence of new profiling technologies, including genome-wide expression arrays, has increased hope that these factors can soon be disentangled.
Our current study measured RNA expression in the largest genome-wide analysis of human skeletal muscle tissue of which we are aware [7–17]. In process of analyzing this data, we noticed that the expression levels for some genes showed a distinct bimodal distribution. This pattern was quite intriguing (a nearly ubiquitous assumption with genome-wide expression arrays is that expression-level data are normally or log-normally distributed), and led us to ascertain how often such patterns occurred, and whether such were associated with diabetes-related phenotypes. We also sought to determine the relative role that genetic factors had in the influence of these bimodal patterns.
Limited reports of bimodal expression patterns in other studies have also been made [18–20], particularly when heterogeneous tissues are being contrasted [21–26] such as between healthy and diseased tissue (and for which bimodal expression levels would be expected for perhaps many genes). However, in homogeneous tissue of apparently healthy individuals, such bimodality may be surprising, and hints at unrecognized heterogeneity which may involve precipitating factors of future disease development. Previous reports of bimodality [18, 21–23, 27] have suggested that 10-30% of transcripts show bimodal expression, yet with an apparent lack of replication of such findings . As part of our analysis, we also observed that batch effects could potentially cause the artifactual detection of bimodality, inflating the prevalence estimate, but that after proper control, the prevalence of genuine bimodality in muscle tissue was in fact quite rare.
While this discrepancy in the estimated prevalence of bimodality is large and invites further investigation, the ability to distinguish genuine bimodals can lead to the more accurate evaluation of their biological role. Bimodality is an important biological phenomenon as it implies the existence of discrete populations (e.g., with and without disease) or discrete genetic or environmental influences on a trait. As such, bimodality has historically been sought in genetic research, yet few reproducible bimodal patterns have ever been observed in this field. One might, therefore, hope that genuine bimodally expressed genes will have significant genetic importance. Hence, it is most intriguing that many of the bimodal genes we identify herein have reported associations with a variety of human diseases and phenotypes. This elicits a multitude of hypotheses regarding the role of bimodal expression and human characteristics, including disease development. Our analysis also evaluates the role of bimodality in one such disease - type 2 diabetes, by assessing the association of the identified bimodal transcripts with pre-diabetic phenotypes in this population of non-diabetic Pima Indians who are at high risk for developing diabetes in their lifetime.
Prevalence of bimodal expression
The prevalence of bimodal expression was estimated from skeletal muscle biopsies taken from 225 non-diabetic, healthy Pima Indians. These tissue samples were scanned on Affymetrix Human Exon 1.0 ST microarray chips which provide on average 49 expression level measurements per gene. Subject characteristics are shown in Table 1. To assess bimodality of RNA expression, we fit unimodal and bimodal distributions to the gene-level expression data coming from each of the 17,881 core genes. We identified four clusters of chips by scan date which exhibited batch-type effects, two of which could be identified for lot numbers and fluidics station sets used in processing (Additional file 1). These batch effects were found to potentially cause spurious detection of bimodality, and hence the analyses were performed separately in two groups of identifiable chips of size 71 and 47. The full set of size 225 was also assessed for bimodality when such was observed in either of these independent groups (see Figure 1).
The total number of bimodally expressed genes estimated in each of the chip groups using different requirements for genuine bimodality is shown in Table 2. This analysis by group found only a small number of genes to be bimodally expressed with little overlap between the two components of the distribution, 19 in group A and 21 in group B, comprising a total of 32 unique genes (Table 3). Of these, 8 genes (GSTM1, HLA-DRB1, ERAP2, HLA-DRB5, MAOA, ACTN3, NR4A2, and THNSL2) were found to meet the stringent bimodal criterion in group A, group B, and the total study of 225 samples. Figure 2 shows the log2 expression level distributions for each of these confirmed bimodal genes. An additional 10 genes identified as bimodal had significant association with gender, of which 9 were located on the Y chromosome (Additional file 2).
As bimodal expression may reflect either two distinct levels of transcript abundance or a single transcript abundance level with the lower mode reflecting no expression, we compared the median levels of the lower mode of the confirmed bimodal genes to the median expression levels of all other genes (n = 17,783). The levels of the lower mode were indicative of background (no gene expression) for half of these confirmed bimodals - having median expression levels at the 25th or lower percentile of all genes, while four of them (THNSL2, MAOA, GSTM1, and ACTN3) appeared to have RNA expression at two distinct abundance levels, each with lower modes corresponding to levels near or above the median of the high mode of ERAP2 (55th percentile).
Many genes initially classified as bimodal in group A had very few data points (often 1 or 2) lying to the left or right of the intersection of the two component normal distributions. As discussed in Methods, genuine bimodality should not reflect the potential influence of a small number of outliers. For a robust reassessment of bimodality, we deleted the lowest 5% and highest 5% of expression values by chip group for each gene. Bimodality was again assessed in these trimmed data sets (n = 65 and n = 43), with similar bimodal prevalence estimates found to the estimate which required a minimum of 10% of data points to be in each component (Additional files 3 and 4). A combined Fisher analysis of these trimmed data sets found the same set of 8 confirmed bimodals plus an additional 6 genes (HLA-C, SLC44A5, NR4A3, HSPC157, ABCC6, and C19orf62) to have a false discovery rate (FDR) < 0.05. The p-values for bimodality of all 17,881 genes in both independent batch groups are plotted in Figure 3.
Diabetes related phenotypes and mode of expression
We assessed the association of the mode of expression with pre-diabetic traits (body mass index/percent body fat, insulin resistance, and insulin secretion), in the eight confirmed bimodal genes (Table 4). We found insulin sensitivity as measured by a hyperinsulinemic, euglycemic clamp to be higher in those in the low mode of ACTN3 expression (p = 0.0076) and percent body fat to be higher in those in the high mode of MAOA expression (p = 0.018). Body mass index and acute insulin response to an intravenous glucose tolerance test were not found to be associated with any of the bimodal genes in this analysis.
Association of bimodal expression and cis determinants
To investigate the potential cause of such bimodal expression, a subset of 149 individuals who were also genotyped (Affymetrix Genome-wide Human SNP Array 6.0) were analyzed for association with muscle expression levels. Cis- acting SNPs (defined as lying within a 200 kb region to either side of the gene location) were sought in this population. Nearly all of the confirmed bimodally expressed genes and half of the entire bimodal set had highly significant (unadjusted p < 10-8) associations between muscle tissue expression and SNP genotype (Table 5 and Additional file 5), and were hence expression quantitative trait loci (eQTL). Using the intersection of the two underlying normal distributions as the threshold for assigning individuals to either component of the bimodal distribution, the coefficient of agreement between expression component and SNP genotype was high (κ > 0.80) for the majority of the confirmed bimodal genes (Table 5). Still, some genes had little or less than expected agreement for the density of coverage of the SNP chip. For example, NR4A2 had no agreement at all (κ = 0.0). We further investigated this gene for any trans-determinant cause across the entire genome, but no SNP approached a genome wide significance level. We also noted that while GSTM1 had the clearest separation of expression modes, it had the lowest kappa coefficient of the seven cis-associated bimodals, possibly indicating influences other than nearby SNPs on the expression levels, such as environmental exposure or more complex genetic models.
Prevalence of copy number variation in bimodal genes
As polymorphic deletions may greatly affect expression levels, we sought to determine the frequency with which such copy number variation (CNV) may be responsible for producing the bimodal expression patterns. The genes listed in Table 3 were assessed from the UCSC Genome Browser for structural variation (gain, loss, or both) that involved segments of DNA larger than 1 kb. Of the 8 confirmed bimodally expressed genes, copy number variation is known to be present in 3, while 18 of the total 32 possible bimodal genes have currently identified copy number variation. As a comparison, we also assessed a random sample of 32 genes with normally distributed gene expression, and found the same number of genes - eighteen, with reported CNV in this database, indicating no difference in the prevalence of CNV in bimodal as compared to unimodally expressed genes. Hence while CNV may well cause bimodal expression for some genes (most notably GSTM1), for the majority of confirmed bimodals, known CNV is likely not the explanation for such a distribution.
Validation in other populations
We further assessed bimodal RNA expression in two publically available data sets containing expression data from healthy muscle biopsies (GSE13070 ; n = 59 and GSE5086 ; n = 62). We found the prevalence of bimodal RNA expression in these populations to be 0.72% and 1.4% respectively. These estimates were similar to the prevalence we observed, providing further validation that bimodal RNA expression in healthy muscle tissue is a rare event. We also re-assessed bimodal RNA expression in healthy lymphoblastoid cells (GSE1485 ; n = 193) which had previously been estimated to have high bimodal prevalence . We found the prevalence of bimodality in the restricted genes of that analysis to be only 2.9% by our criteria, the difference largely reflecting a previously used statistical threshold which is associated with a much higher false discovery rate than that used in other bimodal expression studies. Of great interest, we found that many of the bimodal genes of our analysis were also highly significant for bimodality in these other three populations (Table 6). These additional validations are most compelling as clear bimodality was also rare in these other populations, with p < 10-8 of chance concordance of the bimodal gene sets in each of these populations with the Pima bimodals (Additional file 6). We also note that while the most clearly bimodal gene replicating in our analysis (GSTM1) did not show strict bimodality in these other populations, related genes were highly bimodal in these other data sets (e.g. GSTM2 in GSE1485 and GSTT1 in GSE5086, data not shown).
Our detection of bimodal gene expression in skeletal muscle tissue reveals an apparently rare (<0.2%), yet potentially new predictor of human disease and phenotypes. Of the eight confirmed bimodal genes, we found that seven were expression quantitative trait loci, as they were strongly to completely associated with possibly causal SNPs occurring within a 200-kb cis region of the gene. This finding suggests that the majority of bimodal genes which may be found in healthy human tissue may be the result of polymorphisms within such gene regulatory regions. In the remaining potential bimodal genes, highly significant SNP associations were less common (possibly reflecting misclassification of bimodal genes), though still more frequent than occurs genome-wide. Hence in healthy individuals, eQTL are prime (yet not exclusive) locations for finding bimodally expressed genes. We also found that known copy number variation - while a potential explanation for some bimodal expression, was not found in higher abundance than across the genome, and is hence likely not the main source of such distributions.
The current study also assessed the relationship of RNA expression mode with known predictors of diabetes. We observed a likely association of insulin sensitivity with mode of ACTN3 expression, with 42.2% of the individuals falling in the higher expressed mode which was associated with insulin resistance. While not meeting strict correction for multiple testing, this gene has previously been associated with slow- and fast-twitch muscle fiber, with correlation to athletic performance . The bimodal expression of ACTN3 could reflect a bimodal distribution in the proportion of type 2 fibers present in the subjects, though we are not able to assess this directly. As insulin resistance is a known predictor of diabetes , the role of ACTN3 expression deserves further investigation.
Many of the other bimodal genes identified have previously been reported to be associated with other diseases and characteristics ranging from cancer to arthritis [31–39], some of these also reporting marked dichotomies in expression levels [38, 39]. Perhaps the most intriguing disease-related bimodal genes are the HLA loci, as the HLA system is known not only to convey risk for a variety of immune-related diseases, but contains the main genetic risk factors for type 1 diabetes . That three of the HLA genes (DRB1, DRB5, and C) were found to have bimodal distributions - as well as ERAP2 which plays a final role in producing HLA peptides , is of great interest, and we are currently further investigating the role of these genes in relation to diabetes and various metabolic phenotypes in the Pima population.
The prevalence of bimodal gene expression we have observed in a healthy single tissue is approximately 100-times lower than what has been reported in mixed-tissue studies [21, 22]. Since these studies assessed distinct expression differences across heterogeneous tissues, a larger number of bimodally expressed genes would be expected. Our consistently low estimate gives confidence that at least in healthy muscle tissue, unambiguous bimodal expression is not common - providing a new reference prevalence for bimodality in healthy tissue to which studies of bimodality in diseased and mixed tissues may be compared.
Our estimate also differs greatly from the only other prevalence reported for bimodal expression in a population comprised solely of healthy persons - 28% in lymphoblastoid cells [18, 28]. As the threshold for assessing bimodality in that study was much more liberal, we re-analyzed the data from this population and found only 2.9% of genes of the limited gene set previously analyzed to show bimodality consistent with the thresholds used in our analysis (see Methods). We also assessed two other publically available data sets containing muscle biopsies, finding estimates of clear bimodality in these populations to be 0.72% and 1.4%, and again with <0.2% prevalence with replication of bimodality across these two data sets (data not shown). These validation sets used chips which had on average 1 or 2 probes per gene, whereas our present study is the first to assess bimodality in expression data from exon arrays which contained on average 49 probes per gene. Hence, the bimodality we observed represents a strong replicated signal which could be missed by a single probe, yet could also miss uniquely bimodal probes when averaged at the gene level. The additional validation of our confirmed bimodal genes in these other populations with low bimodal prevalence (Table 6) provides even more convincing evidence of their authenticity.
We highlight two technical issues for future studies to be aware of: use of a model that allows for unequal variances in the bimodal distribution and the determining of batch effects. We observed that a model that allows for unequal variances was clearly more appropriate than a constant variance model for approximately 20% of the genuine bimodal genes identified in our data, and hence the failure to allow for such as well as differences in data scaling may contribute to some variability in prevalence estimates. Even more important is the potential of batch effects to cause artifactual bimodal expression. Such confounding can occur when any set of experimental factors is not uniform - causing the expression levels to be uniformly over- or under-estimated for a given transcript for a portion of the samples. We noted that these batch effects were associated more frequently with genes of higher transcript abundance levels. We also identified that chips which had distinct lot number differences and which had been processed on distinctly different fluidics station sets gave rise to such artifactual bimodality, though it was not possible to determine which of these factors had a more primary role in this study. However, because this grouped source of heterogeneity was identified for most arrays, bimodality could be analyzed in homogeneous groups, eliminating not only the unwanted bias, but also providing distinct replication sets in which tests for bimodality could be confirmed.
In conclusion, the prevalence of bimodal RNA expression in homogeneous tissue is quite low, and hence has been unrecognized in most expression analyses. The few prior genome-wide studies which have reported bimodality may have over-estimated its prevalence, possibly due to the presence of batch effects and/or overly lenient statistical thresholds. Differences in cell type, microarrays, and populations may also be responsible for some of the large difference in estimated prevalence. An additional novel insight from this analysis is that the majority of genes with bimodal RNA expression apparently exhibited this pattern due to cis-polymorphisms, being expression quantitative trait loci. Future investigations of bimodality in homogeneous tissues from healthy persons will be of great interest in assessing similarities and differences across other tissues and populations, as well as their association with other human characteristics and diseases. In the Pimas, the association of a bimodal gene with insulin resistance and the presence of bimodality in genes from the HLA locus which carries the greatest risk for type 1 diabetes underscore the need to further investigate and understand this interesting phenomenon.
All subjects were part of our study of the etiology of type 2 diabetes among Pima Indians living in Arizona. Volunteers providing muscle biopsies were non-diabetic, healthy, and were not taking any medications. Subjects were assessed for diabetes-free status via a fasting blood draw as well as a 75 gram 2 hour oral glucose tolerance test (OGTT) to confirm that the subject was indeed non-diabetic (fasting glucose <126 mg/dl and 2 hour glucose <200 mg/dl). All subjects provided informed consent (which places some restrictions on access to data and intended research purposes) and all studies were approved by both the institutional review board of the National Institute of Diabetes and Digestive and Kidney Diseases and the council of the Gila River Indian Community.
Most muscle biopsies were performed within one week's time of the diabetes assessment (all biopsies were performed within 30 days). After a 12 hour overnight fast and monitored physical inactivity, percutaneous needle biopsies were carried out on the vastus lateralis muscle under local anesthesia with 1% lidocaine. The biopsy specimens were snap frozen and stored at -20 to -70°C. Variation in storage time ranged up to 16 years, however no effect of the amount of time frozen on RNA quality was found.
Muscle tissue samples were homogenized using a PT 1300D (Brinkmann, Westbury, NY, USA) homogenizer with 12 mm rotor/stator head in 2.0 ml of cold TRIzol reagent (Invitrogen/Life Technologies, Carlsbad, CA, USA). The samples were then heated to 65°C and 0.4 ml chloroform (Sigma, St. Louis, MO, USA) was added. 5 μL of 4 ng/μL Carrier RNA was added to each sample to aid precipitation with ETOH and further purified using RNeasy Micro Kit (Qiagen, Valencia, CA, USA). 1 μg Total RNA Labeling Protocol starts with a ribosomal RNA (rRNA) reduction procedure where the 28S and 18S rRNA portion is significantly reduced from the total RNA sample minimizing the background and increasing the array detection sensitivity and specificity. All RNA was isolated by a single technician.
The rRNA was prepared and the cDNA synthesized for the Human Exon 1.0 ST Array microarray chips (Affymetrix, Santa Clara, CA, USA) using the GeneChip Whole Transcript Sense Target Labeling Assay kits (Affymetrix). Chip washing and staining was performed in two distinct sets of fluidics stations (GeneChip Fluidics Station 450, Affymetrix) and scanned on as many as five different scanners (GeneChip Scanner 3000 7G, Affymetrix). Two technicians performed the scanning measurements in discrete sets of dates.
Quality control metrics and technician insight were used to determine the need for re-scans. When re-scanning did not resolve the problem, such chips were discarded and cDNA hybridization mix placed on completely new chips (this occurred in 9.8% of the samples). A total of 225 baseline samples which passed such quality control metrics resulted in expression results available for analysis. Only transcript clusters identified as "core" (genes which have been identified with a high level of confidence) were analyzed, comprising a total of 17,881 genes.
To account for varying intensity levels and other variability between different scans, we normalized the 225 chips using the Robust Multichip Average method with quantile normalization of the log2 intensities with prior GC correction (GC-RMA [40, 41]), using the mean of all probes associated with the given gene. Visual inspection of chip intensities found some chips with limited smears or mars (a common occurrence in such experiments), however, the bimodal methods employed are likely to be robust to such occasional random technical events, and hence no chips or regions of chips were excluded.
In analyzing the 225 muscle files, we discovered a batch effect which was associated with distinctly different expression values by clusters of sample processing dates for some genes (Additional file 1). Investigating this batch effect, we determined that the final two clusters consisted of microarrays with distinctly different chip lot numbers and which had been processed on different fluidics stations. Initially (i.e. before discovering that a batch effect was present), use of all available chips led to the false detection of many bimodal genes - as combined expression values of normally distributed genes from different batches produced an apparent overall bimodal distribution. Hence, we conducted our analyses on the expression data from each batch separately - providing not only an independent confirmation set, but also an analysis free from this bias. While common linear adjustments may remove batch effects for normally distributed data, due to natural unequal proportions of data lying in each component of bimodal data, similar adjustments could result in bimodal data from two batches having 4 modes. Hence, we restrained our analysis to unadjusted data in separate batches.
The remaining 107 microarray chips could not be identified for lot number or for the fluidics station set they were processed on. Hence, we did not explore bimodality in this data separately, but constrained our initial analysis to the two identifiable clusters that had also been processed by a single technician (a cluster of 47 chips - group A, and another cluster of 71 chips - group B). While differences between clusters were observed for some genes, other genes appeared to be free from this batch effect. For these unaffected genes, power to detect associations with true bimodality would be greater in the total data set (225 chips) due to increased sample size. Hence, we first confirmed bimodality in the separate clusters before using the entire data set for further confirmation and association analyses.
Unimodal and Bimodal Distribution Fitting
For each gene, we had 47, 71, or 225 data points corresponding to the expression levels coming from each participant in the two batch groups or in the entire study. To avoid skewness with a unimodal distribution resulting in an artifactual inference of bimodality, the log2 expression data for each gene were transformed by the Box-Cox method to reduce skewness . The Box-Cox parameter was systematically varied in increments of 0.01 over the range of -4 to 4 to select the most appropriate value for each gene. Subsequently, maximum likelihood methods were used to estimate the parameters of a unimodal Gaussian distribution and of a bimodal Gaussian distribution (allowing for unequal variances of the two modes) that were most consistent with the transformed data. The unimodal density distribution is described by two parameters, the mean (μ) and standard deviation (σ), while the bimodal involves the means and standard deviations in each component (μ1, σ1, μ2, σ2) and the proportion of area in one of the components (e.g., p1). Such methods for model fitting have been used previously [43, 44]. A potential difficulty with maximum likelihood techniques is that it is difficult to be certain that the model has converged at the global, rather than a local, maximum of the likelihood. To ensure such convergence, parameters were estimated using a novel automatic parameter estimation technique. Essentially this method calculates the likelihood at various plausible values of the parameters chosen systematically (200,000 estimations per gene for the present study); the optimal values are then taken as initial estimates for a final maximum likelihood estimation using the Newton-Raphson algorithm. The result of this method was to greatly increase the likelihood of converging to the global maximum, with a coinciding increase in computational time that was not prohibitively long for the bimodal model (time was less than 1 minute per gene).
There is no unequivocal selection test for bimodal versus unimodal distributions, despite the fact that the models are nested [22, 24, 45]. Hence, a number of different criterion have been proposed, which we adopted. These include using a p-value threshold of 0.001 from a chi-square distribution with six degrees of freedom (as opposed to three in order to provide a more stringent cutoff) on the minus 2 log likelihood difference from the models, as well as a simultaneous requirement of misclassification area <0.1 based on the best estimated bimodal distribution [21, 22, 45]. This area is calculated as the minimum area under either curve to the right or left of the intersection of the two separate modes. The number of intersections can be 0, 1, or 2. When 0, the mean values were used as the dichotomous intersection estimates. When multiple, both intersections were used, with the intersection providing the minimal misclassification selected.
Additionally, we introduced a new requirement that the minimum number of data points to the left and right of the selected intersection point be at least 10%. This prevented bimodal model selection when a single or small number of outliers would otherwise have caused bimodal selection. Such a requirement is necessary to eliminate potentially large numbers of bimodal misclassifications for genes in which the best fit distribution's misclassification area is not significant, yet the bimodal fit is clearly the result of a single or small number of outlier points. As a consequence, this additional criterion will permit only the detection of bimodal genes that are dichotomously expressed in a non-trivial percentage of individuals, yet with greater confidence that these genes are not mis-assessed due simply to outlier points.
We confirmed the need for such a strategy in a separate analysis in which the lowest 5% and highest 5% of expression values were dropped. Such a strategy would prevent misclassification of bimodality due to a small number of outliers, though it would reduce the power to detect such in clean data. We noted this method to give very similar end-result prevalence estimates.
As some genes are expressed differently by gender, we fit a separate linear model on GC-RMA normalized expression levels to assess this association, using the resulting two-sided p-value as an additional tool to determine which genes exhibit bimodal expression due to different expression levels in men and women. If this p-value was ≥ 0.05, it was assumed that the bimodality was not due to gender.
A total of 149 individuals with expression data were also genotyped with Affymetrix Genome-wide Human SNP Array 6.0 chips. Potential cis-acting SNPs (expression quantitative trait loci, ) were identified as falling within 200 kb regions upstream or downstream of the start and stop locations of each bimodal transcript. A linear model was used to assess the association with adjustment for age, gender, and heritage with sibship included as a random effect to account for relations. A total of 1,221,921 SNP-expression associations were tested, from which a conservative Bonferroni estimate of significance would be 4.09 × 10-8. Expression levels were rank transformed with subsequent back-transformation to a normal distribution prior to fitting the linear models.
Assessment of pre-diabetic traits was made as part of the research protocol corresponding to the investigation nearest to the muscle biopsy, usually performed on the same day, and never greater than 30 days from the biopsy date. Insulin action in vivo was assessed using the hyperinsulinemic, euglycemic clamp as described in . Acute insulin secretion was assessed during an intravenous glucose tolerance test with analysis restricted to those with normal glucose tolerance (OGTT <140 mg/dl) and who were full heritage Pima . Percent body fat was performed via underwater weighing and total body dual energy X-ray absorptiometry .
We assessed bimodality in three additional data sets accessed from the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/). We re-analyzed bimodality in the expression of lymphoblastoid cells from 193 CEPH Utah residents (GSE1485, [18, 28]) and analyzed for the first time bimodality in 59 normal abdomen muscle biopsies obtained during gastrointestinal surgeries (GSE13070, ) and 62 skeletal muscle biopsies in a population with broad range of insulin sensitivity/resistance (GSE5086, ). The expression data files analyzed came from individuals prior to treatment when such occurred in the study designs and for which data were available. Criteria for bimodality were synonymous with those of the present study. When multiple probes for a gene were present in these data sets, the probe giving the best evidence of bimodality was used. Only 3,554 genes from GSE1485 which had less intra-individual variation (assessed from replicates) than inter-individual variation were analyzed, as was previously done [18, 28]. This would infer that our re-computed prevalence of 2.9% bimodality in this gene subset overestimates the total genome-wide prevalence of bimodal expression in that data. A total of 41,789 probes were analyzed in the GSE13070 and GSE5086 data sets. Probability of concordance of bimodal genes between the data sets was assessed via Fisher's exact test .
Fisher's method was used to combine p-values from independent chip groups. False Discovery Rates  were assessed in the 12,470 genes that were not associated with gender (p ≥ 0.05) and separately in the 5,411 genes that were associated with gender (p < 0.05). All other p-values listed are not adjusted for multiple comparisons. General linear models were used to determine association of traits with expression mode after adjustment for appropriate covariates with sibships treated as repeated measures. Coefficients of agreement between genotype and mode of expression were calculated as the maximum kappa-Cohen statistic of the four possible groupings of dominant models and mode pairings. Copy number variation was assessed from the UCSC Genome Browser (2006 Assembly) as structural variation (gain, loss, or both) that involved segments of DNA larger than 1 kb.
Gene scans were performed using Affymetrix GCOS software and Affymetrix Expression Console. GC-RMA normalization was performed using Partek Software Version 6.4 (6.09.0422). SAS Version 9.1 was used for Box-Cox transformations, Newton-Raphson maximum likelihood convergence, pre-diabetic phenotype associations, Fisher's exact test, as well as for all SNP association analyses. Code was also written in Microsoft Visual Studios C++ for data manipulation and maximum likelihood with automatic parameter estimation distribution fitting of unimodal and bimodal models.
Idaghdour Y, Storey JD, Jadallah SJ, Gibson GA: Genome-Wide Gene Expression Signature of Environmental Geography in Leukocytes of Moroccan Amazighs. PLoS Genetics. 2008, 4: e1000052-10.1371/journal.pgen.1000052.
Idaghdour Y, Czika W, Shianna KV, Lee SH, Visscher PM, Martin HC, Miclaus K, Jadallah SJ, Goldstein DB, Wolfinger RD, Gibson G: Geographical genomics of human leukocyte gene expression variation in southern Morocco. Nature Genetics. 2010, 42: 62-67. 10.1038/ng.495.
Emilsson V, Thorleifsson G, Zhang B, Leonardson AS, Zink F, Zhu J, Carlson S, Helgason A, Walters GB, Gunnarsdottir S, Mouy M, Steinthorsdottir V, Eiriksdottir GH, Bjornsdottir G, Reynisdottir I, Gudbjartsson D, Helgadottir A, Jonasdottir A, Jonasdottir A, Styrkarsdottir U, Gretarsdottir S, Magnusson KP, Stefansson H, Fossdal R, Kristjansson K, Gislason HG, Stefansson T, Leifsson BG, Thorsteinsdottir U, Lamb JR, et al: Genetics of gene expression and its effect on disease. Nature. 2008, 452: 423-428. 10.1038/nature06758.
Myers AJ, Gibbs JR, Webster JA, Rohrer K, Zhao A, Marlowe L, Kaleem M, Leung D, Bryden L, Nath P, Zismann VL, Joshipura K, Huentelman MJ, Hu-Lince D, Coon KD, Craig DW, Pearson JV, Holmans P, Heward CB, Reiman EM, Stephan D, Hardy J: A survey of genetic human cortical gene expression. Nature Genetics. 2007, 39: 1494-1499. 10.1038/ng.2007.16.
Dixon AL, Liang L, Moffatt MF, Chen W, Heath S, Wong KC, Taylor J, Burnett E, Gut I, Farrall M, Lathrop GM, Abecasis GR, Cookson WO: A genome-wide association study of global gene expression. Nature Genetics. 2007, 39: 1202-1207. 10.1038/ng2109.
Conrad DF, Pinto D, Redon R, Feuk L, Gokcumen O, Zhang Y, Aerts J, Andrews TD, Barnes C, Campbell P, Fitzgerald T, Hu M, Ihm CH, Kristiansson K, Macarthur DG, Macdonald JR, Onyiah I, Pang AW, Robson S, Stirrups K, Valsesia A, Walter K, Wei J, Wellcome Trust Case Control Consortium, Tyler-Smith C, Carter NP, Lee C, Scherer SW, Hurles ME: Origins and functional impact of copy number variation in the human genome. Nature. 2010, 464: 704-712. 10.1038/nature08516.
Palsgaard J, Brøns C, Friedrichsen M, Dominguez H, Jensen M, Storgaard H, Spohr C, Torp-Pedersen C, Borup R, De Meyts P, Vaag A: Gene expression in skeletal muscle biopsies from people with type 2 diabetes and relatives: differential regulation of insulin signaling pathways. PLoSOne. 2009 Aug 11, 4 (8): e6575-
Stepto NK, Coffey VG, Carey AL, Ponnampalam AP, Canny BJ, Powell D, Hawley JA: Global gene expression in skeletal muscle from well-trained strength and endurance athletes. Med Sci Sports Exerc. 2009, 41: 546-565.
Lee YH, Tokraks S, Pratley RE, Bogardus C, Permana PA: Identification of differentially expressed genes in skeletal muscle of non-diabetic insulin-resistant and insulin-sensitive Pima Indians by differential display PCR. Diabetologia. 2003, 46: 1567-1575. 10.1007/s00125-003-1226-1.
Yang X, Pratley RE, Tokraks S, Bogardus C, Permana PA: Microarray profiling of skeletal muscle tissues from equally obese, non-diabetic insulin-sensitive and insulin-resistant Pima Indians. Diabetologia. 2002, 45: 1584-1593. 10.1007/s00125-002-0905-7.
Sreekumar R, Panagiotis H, Schimke JC, Nair KS: Gene expression profile in skeletal muscle of type 2 diabetes and the effect of insulin treatment. Diabetes. 2002, 51: 1913-1920. 10.2337/diabetes.51.6.1913.
Frederiksen CM, Højlund K, Hansen L, Oakeley EJ, Hemmings B, Abdallah BM, Brusgaard K, Beck-Nielsen H, Gaster M: Transcriptional profiling of myotubes from patients with type 2 diabetes: no evidence for a primary defect in oxidative phosphorylation genes. Diabetologia. 2008, 51: 2068-2077. 10.1007/s00125-008-1122-9.
Yi Z, Bowen BP, Hwang H, Jenkinson CP, Coletta DK, Lefort N, Bajaj M, Kashyap S, Berria R, De Filippis EA, Mandarino LJ: Global relationship between the proteome and transcriptome of human skeletal muscle. J Proteome Res. 2008, 7: 3230-3241. 10.1021/pr800064s.
Zahn JM, Sonu R, Vogel H, Crane E, Mazan-Mamczarz K, Rabkin R, Davis RW, Becker KG, Owen AB, Kim SK: Transcriptional profiling of aging in human muscle reveals a common aging signature. PLoS Genetics. 2006, 2: e115-10.1371/journal.pgen.0020115.
Sears DD, Hsiao G, Hsiao A, Yu JG, Courtney CH, Ofrecio JM, Chapman J, Subramaniam S: Mechanisms of human insulin resistance and thiazolidinedione-mediated insulin sensitization. Proc Natl Acad Sci USA. 2009, 106: 18745-18750. 10.1073/pnas.0903032106.
Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, Puigserver P, Carlsson E, Ridderstråle M, Laurila E, Houstis N, Daly MJ, Patterson N, Mesirov JP, Golub TR, Tamayo P, Spiegelman B, Lander ES, Hirschhorn JN, Altshuler D, Groop LC: PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nature Genetics. 2003, 34: 267-273. 10.1038/ng1180.
Patti ME, Butte AJ, Crunkhorn S, Cusi K, Berria R, Kashyap S, Miyazaki Y, Kohane I, Costello M, Saccone R, Landaker EJ, Goldfine AB, Mun E, DeFronzo R, Finlayson J, Kahn CR, Mandarino LJ: Coordinated reduction of genes of oxidative metabolism in humans with insulin resistance and diabetes: Potential role of PGC1 and NRF1. Proc Natl Acad Sci USA. 2003, 100: 8466-8471. 10.1073/pnas.1032913100.
Yang Y, Tashman AP, Lee JY, Yoon S, Mao W, Ahn K, Kim W, Mendell NR, Gordon D, Finch SJ: Mixture modeling of microarray gene expression data. BMC Proceedings. 2007, 1 (Suppl 1): S50-10.1186/1753-6561-1-s1-s50.
Hellwig B, Hengstler JG, Schmidt M, Gehrmann MC, Schormann W, Rahnenführer J: Comparison of scores for bimodality of gene expression distributions and genome-wide evaluation of the prognostic relevance of high-scoring genes. BMC Bioinformatics. 2010, 11: 276-10.1186/1471-2105-11-276.
Bessarabova M, Kirillov E, Shi W, Bugrim A, Nikolsky Y, Nikolskaya T: Bimodal gene expression patterns in breast cancer. BMC Genomics. 2010, 11 (Suppl 1): S8-10.1186/1471-2164-11-S1-S8.
Ertel A, Tozeren A: Human and mouse switch-like genes share common transcriptional regulatory mechanisms for bimodality. BMC Genomics. 2008, 9: 628-10.1186/1471-2164-9-628.
Ertel A, Tozeren A: Switch-like genes populate cell communication pathways and are enriched for extracellular proteins. BMC Genomics. 2008, 9: 3-10.1186/1471-2164-9-3.
Gormley M, Tozeren A: Expression profiles of switch-like genes accurately classify tissue and infectious disease phenotypes in model-based classification. BMC Bioinformatics. 2008, 9: 486-10.1186/1471-2105-9-486.
Wang J, Wen S, Symmans WF, Pusztai L, Coombes KR: The Bimodality Index: A Criterion for Discovering and Ranking Bimodal Signatures from Cancer Gene Expression Profiling Data. Cancer Informatics. 2009, 7: 199-216.
McLachlan GJ, Bean RW, Peel D: A mixture model-based approach to the clustering of microarray expression data. Bioinformatics. 2002, 18: 413-422. 10.1093/bioinformatics/18.3.413.
Zhang J, Ji Y, Zhang L: Extracting three-way gene interactions from microarray data. Bioinformatics. 2007, 23: 2903-2909. 10.1093/bioinformatics/btm482.
Hsieh WP, Passador-Gurgel G, Stone EA, Gibson G: Mixture modeling of transcript abundance classes in natural populations. Genome Biol. 2007, 8: R98-10.1186/gb-2007-8-6-r98.
Morley M, Molony CM, Weber TM, Devlin JL, Ewens KG, Spielman RS, Cheung VG: Genetic analysis of genome-wide variation in human gene expression. Nature. 2004, 430: 743-747. 10.1038/nature02797.
Yang N, MacArthur DG, Gulbin JP, Hahn AG, Beggs AH, Easteal S, North K: ACTN3 genotype is associated with human elite athletic performance. Am J Hum Genet. 2003, 73: 627-631. 10.1086/377590.
Lillioja S, Mott DM, Spraul M, Ferraro R, Foley JE, Ravussin E, Knowler WC, Bennett PH, Bogardus C: Insulin resistance and insulin secretory dysfunction as precursors of non-insulin-dependent diabetes mellitus. Prospective studies of Pima Indians. N Engl J Med. 1993, 329: 1988-1992. 10.1056/NEJM199312303292703.
Engel LS, Taioli E, Pfeiffer R, Garcia-Closas M, Marcus PM, Lan Q, Boffetta P, Vineis P, Autrup H, Bell DA, Branch RA, Brockmöller J, Daly AK, Heckbert SR, Kalina I, Kang D, Katoh T, Lafuente A, Lin HJ, Romkes M, Taylor JA, Rothman N: Pooled analysis and meta-analysis of glutathione S-transferase M1 and bladder cancer: a HuGE review. Am J Epidemiol. 2002, 156: 95-109. 10.1093/aje/kwf018.
Todd JA, Bell JI, McDevitt HO: HLA-DQ beta gene contributes to susceptibility and resistance to insulin-dependent diabetes mellitus. Nature. 1987, 329: 599-604. 10.1038/329599a0.
Zhong H, Yang X, Kaplan LM, Molony C, Schadt EE: Integrating pathway analysis and genetics of gene expression for genome-wide association studies. Am J Hum Genet. 2010, 86: 581-591. 10.1016/j.ajhg.2010.02.020.
Saveanu L, Carroll O, Lindo V, Del Val M, Lopez D, Lepelletier Y, Greer F, Schomburg L, Fruci D, Niedermann G, van Endert PM: Concerted peptide trimming by human ERAP1 and ERAP2 aminopeptidase complexes in the endoplasmic reticulum. Nature Immunology. 2005, 6: 689-697. 10.1038/ni1208.
Schulze TG, Müller DJ, Krauss H, Scherk H, Ohlraun S, Syagailo YV, Windemuth C, Neidt H, Grässle M, Papassotiropoulos A, Heun R, Nöthen MM, Maier W, Lesch KP, Rietschel M: Association between a functional polymorphism in the monoamine oxidase A gene promoter and major depressive disorder. Am J Med Genet. 2000, 96: 801-803. 10.1002/1096-8628(20001204)96:6<801::AID-AJMG21>3.0.CO;2-4.
Meyer JH, Ginovart N, Boovariwala A, Sagrati S, Hussey D, Garcia A, Young T, Praschak-Rieder N, Wilson AA, Houle S: Elevated monoamine oxidase a levels in the brain: an explanation for the monoamine imbalance of major depression. Arch Gen Psychiatry. 2006, 63: 1209-1216. 10.1001/archpsyc.63.11.1209.
Aherne CM, McMorrow J, Kane D, FitzGerald O, Mix KS, Murphy EP: Identification of NR4A2 as a transcriptional activator of IL-8 expression in human inflammatory arthritis. Mol Immunol. 2009, 46: 3345-3357. 10.1016/j.molimm.2009.07.019.
Le WD, Xu P, Jankovic J, Jiang H, Appel SH, Smith RG, Vassilatis DK: Mutations in NR4A2 associated with familial Parkinson disease. Nature Genetics. 2002, 33: 85-89. 10.1038/ng1066.
Coles BF, Anderson KE, Doerge DR, Churchwell MI, Lang NP, Kadlubar FF: Quantitative analysis of interindividual variation of Glutathione S-Transferase expression in human pancreas and the ambiguity of correlating genotype with phenotype. Cancer Research. 2000, 60: 573-579.
Wu Z, Irizarry RA, Gentleman R, Murillo FM, Spencer F: A Model Based Background Adjustment for Oligonucleotide Expression Arrays. Journal of the American Statistical Association. 2004, 99: 909-917. 10.1198/016214504000000683.
Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, Speed TP: Exploration, Normalization, and Summaries of High Density Oligonucleotide Array Probe Level Data. Biostatistics. 2003, 4: 249-264. 10.1093/biostatistics/4.2.249.
MacLean CJ, Morton NE, Elston RC, Yee S: Skewness in commingled distributions. Biometrics. 1976, 32: 695-699. 10.2307/2529760.
Fan J, May SJ, Zhou Y, Barrett-Connor E: Bimodality of 2-h plasma glucose distributions in whites: the Rancho Bernardo study. Diabetes Care. 2005, 28: 1451-1456. 10.2337/diacare.28.6.1451.
Lim TO, Bakri R, Morad Z, Hamid MA: Bimodality in blood glucose distribution: is it universal?. Diabetes Care. 2002, 25: 2212-2217. 10.2337/diacare.25.12.2212.
McLachlan GJ: On bootstrapping the likelihood ratio test statistic for the number of components in a normal mixture. Applied Statistics. 1987, 36: 318-324. 10.2307/2347790.
Schadt EE, Monks SA, Drake TA, Lusis AJ, Che N, Colinayo V, Ruff TG, Milligan SB, Lamb JR, Cavet G, Linsley PS, Mao M, Stoughton RB, Friend SH: Genetics of gene expression surveyed in maize, mouse and man. Nature. 2003, 422: 297-302. 10.1038/nature01434.
Arnoff SL, Bennett PH, Gorden P, Rushforth N, Miller M: Unexplained hyperinsulinemia in normal and "prediabetic" Pima Indians compared with normal Caucasians. An example of racial differences in insulin secretion. Diabetes. 1977, 26: 827-840. 10.2337/diabetes.26.9.827.
Tataranni PA, Ravussin E: Use of dual-energy X-ray absorptiometry in obese individuals. Am J Clin Nutr. 1995, 62: 730-734.
Fury W, Batliwalla F, Gregersen PK, Li W: Overlapping probabilities of top ranking gene lists, hypergeometric distribution, and stringency of gene selection criterion. Conf Proc IEEE Eng Med Biol Soc. 2006, 1: 5531-5534.
Benjamini Y, Hochberg Y: Controlling the false discovery rate: a practical and powerful approach to multiple testing. Journal of the Royal Statistical Society, Series B (Methodological). 1995, 57: 289-300.
We gratefully thank the members and participants of the Gila River Indian Community, whose involvement made these findings possible. We acknowledge the additional technical help and data handling of Kristine Friesen and Sayuko Kobes as well as Drs. Paul Franks, Emilio Ortega, and Helen Looker. This research was supported in part by the Intramural Research Program of the National Institute of Diabetes and Kidney Diseases (NIDDK), and by a DCS award to C.B. from the American Diabetes Association (#7-04-DCS-02).
Data were acquired by CCM, CB, JK, LB, and VO; data were analyzed by CCM; interpretation of the data was made by CCM, RLH, and CB; CB and LJB provided funding; CCM conceived the bimodal study, designed the figures, and wrote the paper, with advice given from CB and RLH. All authors read and contributed to the manuscript.
Electronic supplementary material
Additional file 1:Figure illustrating how dichotomous gene level RNA expression can be artifact of batch effect. (DOC 228 KB)
Additional file 2:Summary of 16 genes (associated with gender) which met all criteria for bimodality in either group A or group B. (DOC 56 KB)
Additional file 3:Summary of bimodal genes found on trimmed bimodal data sets of non-gender associated genes. Listed are 28 genes which had a false discovery rate on the Fisher combined p-value from the two trimmed data sets <1.0. Fourteen of these have an FDR < 0.05. (DOC 94 KB)
Additional file 4:Summary of bimodal genes found on trimmed bimodal data sets of gender associated (p < 0.05) genes. Listed are 30 genes which had a false discovery rate on the Fisher combined p-value from the two trimmed data sets <1.0. Twelve of these have an FDR < 0.05. (DOC 94 KB)
Additional file 5:List of 211 SNPs associated with component of muscle expression with p < 7 × 10-8 in the 32 bimodally expressed genes. (DOC 278 KB)
About this article
Cite this article
Mason, C.C., Hanson, R.L., Ossowski, V. et al. Bimodal distribution of RNA expression levels in human skeletal muscle tissue. BMC Genomics 12, 98 (2011). https://doi.org/10.1186/1471-2164-12-98
- Copy Number Variation
- Batch Effect
- Expression Quantitative Trait Locus
- Bimodal Gene
- Transcript Abundance Level