ParTI Falsely Detects Pareto Fronts and Associated Tasks from the Transcriptomes of Yeast Gene Deletion Strains

To investigate the performance of ParTI on phenotypic data of organisms that are not adapted to their environments, we analyzed the transcriptomes of 1,484 yeast strains, each with a different nonessential gene deleted (Kemmeren et al. 2014). The gene deletion strains (and their transcriptomes) are apparently suboptimal because they are the products of artificial gene removals from a wild-type strain and would accumulate compensatory beneficial mutations when allowed to evolve in the lab (Szamecz et al. 2014). Furthermore, a number of past ParTI analyses reported Pareto optimality from transcriptomes of bacterial cells at different growth stages as well as transcriptomes of different eukaryotic cells, tissues, and cancers (Hart et al. 2015; Korem et al. 2015; Adler et al. 2019; Hausser et al. 2019). Hence, the transcriptomes of yeast gene deletion strains serve as a suitable negative control for PartTI. We preprocessed the yeast transcriptome data following a recent ParTI study of transcriptomes performed by the authors of ParTI (Adler et al. 2019).

Using ParTI, we found that the transcriptomes of 1,484 yeast gene deletion strains are well described by a triangle (P = 0.002; fig. 1D) with apparently meaningful archetypes (table 1 and supplementary table S1, Supplementary Material online). For example, archetype 1 is best described by mitochondrial functions, archetype 2 is best described by carbohydrate metabolism, and archetype 3 is best described by protein synthesis. We found the ParTI result to remain largely unchanged even after the removal of gene deletion strains that are fitter than the wild-type in the synthetic complete (SC) medium, where the transcriptome data were collected (fig. 1E). Finally, ParTI even identified a significant Pareto front from the 1,484 technical replicate transcriptomes of the same wild-type strain under the SC medium (fig. 1F). The archetypes inferred from the wild-type technical replicates (table 2 and supplementary table S2, Supplementary Material online) are similar to those inferred from the gene deletion strains. Furthermore, the PC coordinates of the wild-type technical replicates are highly correlated with those of the matched gene deletion strains in the microarray experiment (the correlation for PC1 coordinates is 0.96, whereas that for PC2 coordinates is 0.94). Therefore, the false-positive ParTI results found from the gene deletion strains are likely due, at least in part, to spurious signals from the microarray experiment. Clearly, ParTI can falsely detect Pareto fronts from phenotypes that are not Pareto optimal. This alarming finding prompts us to systematically study false-positive errors of ParTI.

Table 1.

Tasks Inferred from the Three Archetypes Identified by ParTI from the Transcriptomes of the Yeast Gene Deletion Strains.

Tasks	Corrected P Values
Archetype 1 (mitochondrial function)
Mitochondrial part	$$ 1.96\times {10}^{-23}$$
Envelope and organelle envelope	$$ 1.00\times {10}^{-23}$$
Mitochondrial envelope	$$ 5.57\times {10}^{-18}$$
Establishment of protein localization	$$ 2.53\times {10}^{-10}$$
Protein transport	$$ 2.53\times {10}^{-10}$$
Archetype 2 (carbohydrate metabolism)
Gluconeogenesis	$$ 4.73\times {10}^{-31}$$
Hexose biosynthetic process	$$ 4.73\times {10}^{-31}$$
Monosaccharide biosynthetic process	$$ 4.73\times {10}^{-31}$$
Glucose metabolic process	$$ 6.03\times {10}^{-29}$$
Monosaccharide metabolic process	$$ 6.03\times {10}^{-29}$$
Archetype 3 (protein synthesis)
Regulation of translational elongation	$$ 3.48\times {10}^{-42}$$
Regulation of cellular protein metabolic process	$$ 1.00\times {10}^{-40}$$
Regulation of protein metabolic process	$$ 1.00\times {10}^{-40}$$
Regulation of translation	$$ 6.24\times {10}^{-38}$$
Posttranscriptional regulation of gene expression	$$ 1.05\times {10}^{-36}$$

Note.—Some P values are identical to each other because of the highly redundant GO annotations of yeast genes. For each archetype, we present five GO categories with the largest median expression differences from the random expectation. The full table is provided in supplementary table S1, Supplementary Material online.

Table 2.

Tasks Inferred from the Three Archetypes Identified by ParTI from the Transcriptomes of the Yeast Wild-Type Technical Replicates.

Tasks	Corrected P Values
Archetype 1 (mitochondrial function)
Mitochondrial part	$$ 6.59\times {10}^{-32}$$
Chromatin organization	$$ 1.77\times {10}^{-21}$$
Organonitrogen compound biosynthetic process	$$ 6.03\times {10}^{-31}$$
Envelope and organelle envelope	$$ 3.00\times {10}^{-32}$$
Mitochondrial envelope	$$ 1.60\times {10}^{-24}$$
Archetype 2 (carbohydrate metabolism)
Gluconeogenesis	$$ 9.33\times {10}^{-31}$$
Hexose biosynthetic process	$$ 9.33\times {10}^{-31}$$
Monosaccharide biosynthetic process	$$ 9.33\times {10}^{-31}$$
Carbohydrate biosynthetic process	$$ 5.61\times {10}^{-32}$$
Glucose metabolic process	$$ 6.70\times {10}^{-27}$$
Archetype 3 (protein synthesis)
Regulation of translational elongation	$$ 1.86\times {10}^{-41}$$
Regulation of cellular protein metabolic process	$$ 9.87\times {10}^{-39}$$
Regulation of protein metabolic process	$$ 9.87\times {10}^{-39}$$
Regulation of biological quality	$$ 2.32\times {10}^{-37}$$
Coenzyme binding	$$ 5.47\times {10}^{-36}$$

Note.—Some P values are identical to each other because of the highly redundant GO annotations of yeast genes. For each archetype, we present five GO categories with the largest median expression differences from the random expectation. The full table is provided in supplementary table S2, Supplementary Material online.

Phylogenetic Relationships Cause False-Positive Errors

Because phylogenetic relationships may create spurious correlations among the traits of different species, it has been suggested that ParTI results could be affected by phylogenetic signals in the data (Edelaar 2013; Szekely et al. 2015; Tendler et al. 2015). However, this potential problem has not been explicitly studied, nor is the severity of the problem known. To assess the influence of phylogenetic relationships on PartTI’s performance, we simulated the evolution of 100 independent traits following a randomly generated phylogenetic tree of 100 species (see Materials and Methods). In the simulation, the traits evolved neutrally without adaptation. We then used the default setting of ParTI to infer Pareto fronts in the simulated data. This process of data generation (under different trees) and analysis was repeated 100 times.

To our surprise, ParTI inferred significant Pareto fronts (P < 0.05) from 75 of the 100 simulated data sets (fig. 2A). A typical example is shown in figure 2B. The supplementary materials of the ParTI method paper (Hart et al. 2015) noted that “If additional features are known about the data points, enrichment analysis can help determine whether the triangle is due to Pareto optimality or not. If no feature is enriched near the calculated archetypes, one may suspect that the Pareto approach does not apply.” To investigate whether the feature enrichment analysis is effective in guarding against false-positive errors in Pareto front identification, for each significant case, we also simulated the neutral evolution of 100 independent, binary traits along the tree according to a continuous-time Markov process (Harmon 2018) and used these traits to perform the feature enrichment analysis. Three to seven archetypes were inferred in the data with significant ParTI-identified Pareto fronts (fig. 2C). Importantly, in 72 of the 75 significant cases, there were uniquely enriched features for every archetype, and the frequency distribution of the minimum number of uniquely enriched features of an archetype across the 72 cases is shown in figure 2D. Clearly, the feature enrichment analysis of ParTI cannot curtail mistakenly called Pareto optimality.

Fig. 2.

Phylogenetic relationships among study organisms cause false-positive identifications of Pareto fronts by ParTI. Simulated data with neutrally evolving phenotypes are used. (A) Fraction of data sets with statistically significant Pareto fronts. The dotted line shows the expected fraction in an unbiased test. (B) An example of the inferred Pareto front projected to the first two principal components of the phenotypic space. (C) Frequency distribution of the number of inferred archetypes among data sets with significant Pareto fronts. (D) Frequency distribution of the minimum number of uniquely enriched features per archetype in the 72 data sets with significant Pareto fronts and uniquely enriched features for every archetype.

As a control, we generated data sets without phylogenetic structures using star trees, because all taxa are phylogenetically equally distant from one another in a star tree. We first simulated the neutral evolution of 100 independent traits following a star tree of 100 species with equal branch lengths. We repeated this simulation 100 times, but no significant Pareto front was discovered (fig. 2A). The same was found when a star tree with unequal branch lengths was considered (fig. 2A).

Population Structures Cause False-Positive Errors

A population may contain multiple subpopulations such that mating within subpopulations is more likely than that between subpopulations. Such population structures could mislead ParTI in identifying spurious Pareto fronts from organisms of the same species. To verify this hypothesis, we simulated the neutral evolution of the genotypes of 1,000 individuals according to a population structure (Gutenkunst et al. 2009) resembling a simplified version of the human population with three subpopulations (fig. 3A) and then computed the phenotype of each individual from its genotype using an additive model (see Materials and Methods). The simulation was repeated 100 times.

Fig. 3.

Population structure among study organisms cause false-positive identifications of Pareto fronts by ParTI. Simulated data with neutrally evolving phenotypes are used. (A) The out-of-Africa human demographic model is used in the simulation of neutral trait evolution. AF, African subpopulation; EU, European subpopulation; AS, Asian subpopulation. Double arrowheads indicate gene flows. (B) Fraction of data sets with statistically significant Pareto fronts. The dotted line shows the expected fraction in an unbiased test. The recombination and migrations rates indicated are relative to the base values provided in Materials and Methods. (C) An example of the inferred Pareto front projected to the first two principal components of the phenotypic space. The example is from the left-most bar of panel (B).

We found that ParTI inferred significant Pareto fronts with three archetypes more than 95% of the time despite that neutral evolution was simulated for the traits considered (left-most bar in fig. 3B). A typical example is shown in figure 3C. Furthermore, when using the subpopulation label (Africans, Europeans, and Asians) as discrete features in feature enrichment analysis, we found that in all cases, each archetype is enriched with one of the subpopulation labels, suggesting that the falsely identified archetypes reflect the subpopulations.

Because recombination between genotypes from different subpopulations reduces the population structure in the data, we assessed the impact of the recombination rate on our result. Interestingly, increasing the recombination rate to 10 times the human recombination rate does not substantially alter the ParTI result (second left bar in fig. 3B). Similarly, raising the effective recombination rate by breaking the single chromosome considered into 23 equal-size chromosomes does not qualitatively change the ParTI result (middle bar in fig. 3B). These observations are likely due to the low migration rates among subpopulations used in the simulation such that recombination primarily homogenizes genotypes within the same subpopulation. As expected, when we increased the migration rates among subpopulations by 10 fold, none of the 100 simulated data sets was found to be significantly Pareto optimal (second right bar in fig. 3B). Notwithstanding, in 39 of the 100 cases, each inferred archetype is significantly enriched with a different subpopulation, suggesting that the increased migration may not fully remove the population structure and that ParTI may still infer spurious Pareto fronts if the sample size is larger. On the contrary, when we reduced the migration rates to 10% of the original level, Pareto optimality was found to be significant at a frequency that is similar to that under the original migration rates (right-most bar in fig. 3B). Together, these findings reveal exceptionally frequent false-positive identifications of Pareto fronts by ParTI in the presence of population structure resembling that of humans.

The Flexibility of Data Analysis Inflates False Positives

In examining past studies that used ParTI and analyzing the real and simulated data, we noticed another potential source of false-positive errors—flexibility in data analysis, which provides an opportunity for p-hacking (Simmons et al. 2011). It has been well documented that the false-positive rate is substantially inflated when researchers are allowed to make many liberal decisions in data analysis (Simmons et al. 2011). The reason is simple: the probability that at least one of the many choices makes the result significant at the significance level of 0.05 is much higher than 0.05, the expected false-positive rate under a fixed procedure of data analysis.

Using ParTI to infer Pareto fronts and the associated tasks involves multiple steps, and the investigator has multiple decisions to make at each of these steps. As shown in figure 4A, to use ParTI, a researcher must first decide how to preprocess the data. Next, the researcher has to decide the number of potential archetypes by either using the default setting of ParTI (Hart et al. 2015) or manually specifying it (Hausser et al. 2019). In both cases, one has to first reduce the dimension of the data to a manually set number. Then, the researcher can choose one of five algorithms to compute Pareto fronts (Sisal, MVSA, MVES, SDVMM, and PCHA) (Hart et al. 2015). Among these five algorithms, Sisal, MVSA, and PCHA are most frequently used and have been claimed to show different advantages: Sisal is robust to outliers, MVSA is able to extract information from extreme points, and PCHA is stable with low iteration numbers (Hart et al. 2015). After the Pareto front is inferred, the researcher can perform feature enrichment analysis to infer the involved tasks. The features used can be any properties of the samples, discrete or continuous (Hart et al. 2015).

Fig. 4.

Flexibility in data analysis increases ParTI’s false-positive rate. (A) Pareto front inference involves multiple steps, each of which may be performed in multiple ways. (B) Fraction of data sets of neutral phenotypes that are detected by ParTI to have significant Pareto fronts. The number of ParTI analyses per raw data set equals four times the number of cutoffs considered (see Materials and Methods). The dotted line represents the expected fraction (5%) for an unbiased test.

Among the above steps, data preprocessing is probably the most flexible one. Let us take gene expression data preprocessing as an example. Starting from a raw data set, one may or may not normalize the data by the total expression level (Adler et al. 2019). Next, one can consider using log-transformation (Hart et al. 2015; Korem et al. 2015; Hausser et al. 2019), z-score transformation (Adler et al. 2019), or other transformations (Korem et al. 2015) to transform the data. Furthermore, one can filter the genes by their mean expression levels (Hart et al. 2015; Hausser et al. 2019) or the variability of gene expression levels typically measured by SD (Hart et al. 2015; Adler et al. 2019). Finally, one can decide to remove certain samples based on some flexible criteria. For example, the same breast cancer transcriptome data set used to detect Pareto fronts in cancer evolution was preprocessed differently in two studies; the normal tissue samples were retained in one study (Hart et al. 2015) but removed in the other (Hausser et al. 2019). To be fair to ParTI, most of the preprocessing techniques described above are also commonly used in other literature, especially in exploratory transcriptome analyses that are not associated with particular hypotheses. ParTI, however, is associated with a strong hypothesis that the phenotypes of the organisms examined have been optimized. As such, the flexibility in data analysis provides a huge space for p-hacking.

To illustrate this point, we simulated random phenotypes of 100 traits for each of 1,000 unrelated individuals and subjected the data to ParTI analysis (see Materials and Methods). Without data preprocessing, the false-positive rate of ParTI is 4% (fig. 4B). We then used various cutoffs to exclude genes from the analysis. For example, using the mean value cutoff of 20% excludes the 80% of traits with the lowest mean trait values, whereas using the SD cutoff of 20% excludes the 80% of traits with the lowest trait SDs. As the number of cutoffs explored in data preprocessing increases, the false-positive rate rises substantially (fig. 4B). When we explored 10 different cutoffs per data set (corresponding to 40 tests; see Materials and Methods), significant Pareto fronts were identified in 81% of the 100 replications (fig. 4B). Clearly, the flexibility in data preprocessing can cause a very high false-positive rate for ParTI. Note that the scenario considered here reflects a lower bound of the flexibility in data analysis because we only simulated part of the variations present in published ParTI studies (e.g., the variation in the number of archetypes specified was not simulated). With more researchers using ParTI, additional “creative” ways of data processing are expected.

Preprocessing of the Transcriptomes of Yeast Gene Deletion Strains

From the website http://deleteome.holstegelab.nl/ (last accessed December 24, 2020) (Kemmeren et al. 2014), we downloaded the two-channel microarray gene expression data of 1,484 yeast gene deletion strains grown in the SC medium. Let m_i be the expression level of gene i in a gene deletion strain, and w_i be the expression level of the same gene in the wild-type. The data from each gene deletion strain contained the information of M_i and A_i for gene i defined below:

Therefore, $$ m_{\mathrm{i}}=2^{A_{\mathrm{i}}+M_{\mathrm{i}}/2} \mathrm{and} w_{\mathrm{i}}=2^{A_{\mathrm{i}}-M_{\mathrm{i}}/2}$$, which provided the expression level estimates of gene i in the deletion strain as well as the wild-type. Because the expression levels measured here are meaningful only in a relative sense, we normalized the expression level of each gene in each strain by the sum of the expression levels of all genes in that strain. We retained the genes whose mean normalized expression levels are higher than 7.72 × 10⁻⁴, to strike a balance between the precision in the measurement of gene expression level and the number of genes that can be analyzed. Finally, we transformed the data into z-scores for each gene.

Some yeast gene deletion strains are fitter than the wild-type in the SC medium where the transcriptome data were collected, probably because of antagonistic fitness effects of genes across environments (Qian et al. 2012). To ensure that the transcriptomes used are nonoptimized, we also analyzed a subset of strains by excluding those gene deletion strains that are fitter than the wild-type in the SC medium, on the basis of published fitness of the gene deletion strains relative to the wild-type measured by barcode sequencing in SC (Qian et al. 2012).

To investigate whether Pareto fronts may be erroneously inferred from batch effects in quantifying the transcriptomes, we obtained the gene expression levels in all wild-type technical replicates by a similar procedure.

ParTI Analysis of the Transcriptomes of Yeast Gene Deletion Strains

Before analyzing the transcriptome data from the yeast gene deletion strains, we replicated several published results from ParTI analysis as well as examples in the ParTI package to ensure that we used ParTI as intended by its developers. Following the ParTI manual, we used the default “Sisal” algorithm to detect Pareto fronts from the preprocessed gene expression levels of the yeast gene deletion strains. The user-defined “dim” option, which specifies the number of PCs to retain in ParTI analysis, was set to 5 because the manual recommends it to be set between 5 and 10 (Hart et al. 2015). Similarly, the “binSize” option was set to 0.05 according to the manual (Hart et al. 2015). ParTI suggested three archetypes in the data, so the number of archetypes was set to 3. To perform the task inference, we used GO data from the GO2MSIG database (Powell 2014) (http://www.bioinformatics.org/go2msig/; last accessed December 24, 2020). We used all GO annotations, including IEA (GO terms inferred from electronic annotation). The genes were identified by the National Center for Biotechnology Information gene IDs. We similarly analyzed the gene deletion strains less fit than the wild-type and the wild-type replicates, respectively. We fixed the number of archetypes to 3 when analyzing these two data sets.

Simulation of Neutral Trait Evolution along Phylogenies and the Subsequent ParTI Analysis

We first generated phylogenetic trees with 100 tips using a birth–death process with the birth rate being 1 per unit time and the death rate being a uniform random variable sampled between 0 and 1 per unit time. We then simulated the neutral evolution of 100 independent traits according to the standard Brownian motion model (Harmon 2018). The rate of the Brownian motion was set at 1 per unit time for every trait, meaning that, starting from a mean trait value of μ, the mean trait value after t units of time of evolution is a normal random value with the mean equal to μ and variance equal to t. The simulation of trait evolution was performed using the “sim.traits” function in the “phylocurve” package (Goolsby 2015). For comparison, we simulated the neutral evolution of traits without phylogenetic relations (i.e., along a star phylogeny with equal branch lengths) (Harmon 2018). To achieve this goal, we used Pagel’s λ transformation, which stretches external branches of a tree relative to the internal branches, and set λ at 0. We also simulated neutral evolution along a star phylogeny with unequal branch lengths by randomly sampling the 100 branch lengths from a uniform distribution between 0 and 1 unit time and simulating trait evolution along these branches.

To simulate the evolution of 100 independent binary traits, we used the “sim.char” function in the “phylocurve” package. For each trait, the rate parameters were specified by the following Q matrix of the Markov process: $$ Q=\left[\begin{array}{cc}-0.05& 0.05\\ 0.05& -0.05\end{array}\right]$$.

We repeated the above process 100 times to generate 100 trees (the death rate varies among trees), each associated with a data set of 100 species, each having 100 continuous traits and, when the identified Pareto front is statistically significant, 100 additional binary traits.

We followed the ParTI manual to infer Pareto fronts and the associated tasks from the simulated data sets. The algorithm used was “Sisal.” The “dim” option was set to 8, which is in the middle of the manual-recommended range of [5, 10]; use of other dim values did not qualitatively alter the result. The “binSize” option was set to 0.2. We let ParTI automatically determine the number of archetypes in each data set.

Simulation of Neutral Trait Evolution According to a Human Demographic Model and the Subsequent ParTI Analysis

We used “msprime” (Kelleher et al. 2016) to generate the neutrally evolved genotypes of 100 diploid individuals. Limited by the amount of computational time required, we simulated only one of the 23 human chromosomes. The basic parameters were set according to a published population structure model describing the out-of-Africa origin of human populations with three main subpopulations—Africans, Europeans, and Asians—subject to intersubpopulation gene flow (Gutenkunst et al. 2009). We sampled 33 Africans (by randomly pairing 66 haplotypes sampled from the African subpopulation), 33 Europeans, and 34 Asians. We set the chromosome size to 10⁸ bases, roughly the size of an average human chromosome. The mutation rate was set to 10⁻⁸ per base per generation, and the recombination rate was set to 10⁻⁸ per base per generation, corresponding to the known parameters for humans (Milo et al. 2010). The above process was repeated 100 times to generate 100 independent data sets.

To assess the impact of recombination rate on our result, we generated 100 additional data sets with a recombination rate of 10⁻⁷ per base per generation. In addition, we generated 100 data sets by breaking the single chromosome into 23 chromosomes of 10⁸/23 bases to add random assortment. To assess the impact of the migration rate on our result, we generated 100 data sets with 10 times the original migration rates and another 100 data sets with 10% of the original migration rates.

After simulating the genotype of each individual, we generated the trait values for 100 traits from the genotype. In the simulation, every locus has two alleles in the population, specified by 0 and 1, respectively. The effect of allele 0 on any trait is 0, and the effect of allele 1 on a trait is a standard normal random variable. The trait value will be the sum of effects from all loci on the trait. This simulation is consistent with the universal pleiotropy model in which every locus affects every trait (Wagner and Zhang 2011). The simulated phenotypes were then analyzed by ParTI with default settings to detect three archetypes. The population labels of the output samples were supplied as features to perform the feature enrichment analysis.

Simulation of p-Hacking

To illustrate the point that flexibility in data analysis can lead to p-hacking, we simulated the phenotypes of 100 traits for 1,000 individuals. In each individual, the trait values were sampled from a multivariate normal distribution specified by a 100-dimension vector, whose entries are uniform random variables between 1 and 10 representing the mean phenotypic values of the 100 traits, and a 100 × 100 variance–covariance matrix. The variance–covariance matrix was based on a random trait–trait correlation matrix generated by the program “genPositiveDefMat” from the package “clusterGeneration” with default setting (Joe 2006), coupled with the variance of each trait generated by a random draw from the uniform distribution between 1 and 10. From this raw random phenotypic data set, we generated 40 data sets by commonly used preprocessing techniques. That is, we retained traits among the highest 10%, 20%, …, and 100% of traits in terms of the mean trait value or in terms of the SD, resulting in 20 data sets. From each of these 20 data sets, we converted the trait values to z-scores across individuals, generating 20 additional data sets. We then mimicked three different levels of p-hacking by considering two (50% and 100%), five (20%, 40%, 60%, 80%, and 100%), and all 10 cutoffs in terms of the mean trait value or SD, corresponding to 8, 20, and 40 ParTI analyses per raw data set, respectively. The algorithm used in the ParTI analysis was “PCHA,” another algorithm that was widely used in past ParTI studies. Because the raw phenotypic data are completely random, all Pareto fronts identified are false positives. To restrict the flexibility of ParTI in Pareto front inference, we fixed the number of archetypes to 3 by setting the variable “ForceNArchetypes” in ParTI to 3. The above data generation, preprocessing, and ParTI analysis were repeated 100 times.