dataesr/scientific-paragraphs-categorization huggingface dataset:
Gene expression is more constrained at the proteome level By combining these proteomic and transcriptomic datasets, we are in a position to simultaneously explore and compare the variation of these two gene expression layers at the population level. We therefore computed the absolute Log2(fold change) value (i. e. , |Log2(FC)|) for each gene in each pair of isolates and found that this value is 32% lower on average for the proteome (Figure  2A ), suggesting that protein abundance is less variable and more constrained than mRNA abundance. Furthermore, a higher correlation was observed between proteomes (ρ = 0. 92) compared to transcriptomes (ρ = 0. 83) (Figure  2B ). Finally, the variance observed for each gene was lower for the proteomic data (Figure  S3A ) and the Euclidean distances between each isolate were smaller when computed with the protein abundance dataset (Figure  S3B ).
The biological processes displaying a ratio lower than 0. 93 and a significant difference in terms of branch lengths (see Methods) were considered as strongly buffered. The biological processes displaying a ratio higher than 1 and a significant difference in term of branch lengths had an enhanced abundance variation at the proteome level. D. Biological processes detected as strongly buffered or with an enhanced variation using the procedure detailed in C. Using these datasets, we then sought to understand the main determinants shaping the proteome architecture at the population level. The S. cerevisiae yeast species exhibit a clear population structure, which potentially can impact the proteome landscape 43 (Table  S1 ). We performed a principal component analysis (PCA) with the protein abundance data and found that no clear grouping emerged from the subpopulations when plotting together the 6 first principal components (Figure  S4A-B-C ).
This observed bias is due to the presence of aneuploidies on these chromosomes in our population  43. These CNV-pQTL have also a higher impact on the protein abundance variation compared to the other CNV-pQTL, suggesting that aneuploidies represent a major source of proteome variation at the population level (Figure  S10 ). Only 24 local CNV-pQTL out of 4,528 were detected, and no significant effect size between local and distant CNV-pQTL was found (Figure  5C ). We then looked at the extent to which the genetic bases of protein abundance are common with those underlying the abundance of transcripts. We performed GWAS using the transcriptomic dataset and detected 596 SNP-eQTL and 4,877 CNV-eQTL (Figure  S11 , data file 2), which is of the same order of magnitude as the GWAS proteome results.
Discussion Quantifying transcripts and proteins expressed in a large natural population is fundamental for having a better understanding of the genotype-phenotype relationship. In this study, we have quantitatively analyzed the proteome of 942 natural isolates of S. cerevisiae, allowing in-depth exploration of protein abundance and precise characterization of the genetic origins of its variation at the species level. The S. cerevisiae species is characterized by a complex population structure, with domesticated and wild subpopulations  43. Structured populations are also observed in a large number of other species, such as humans, and their impact on the proteome remains unexplored. In our dataset, the population structure had no significant impact on the proteomic landscape. This observation is consistent with previous results obtained with the transcriptomes of S.
cerevisiae isolates  44, 64. In fact, most subpopulations are characterized by specific signatures related to a small set of genes but not to a general pattern. This dataset allowed us to have better insight into the architecture of the species-wide proteome variation. First, we found that the co-expression network captures main biological functions and is globally conserved across the species. Second, we detected differential protein expression signatures specific to subpopulations, reflecting an adaptation to specific ecological conditions, such as domesticated environments. Similar expression signatures can be also observed using transcriptomic data  44, 65  , highlighting that gene expression plasticity at both levels is a key mechanism of environmental adaptation. The species-wide proteomes and transcriptomes obtained in the same condition represent a unique opportunity to compare the gene regulation at both levels.
Although mass spectrometers are highly sensitive, it should be noted the limitation that proteomic methods are biased towards quantification of highly abundant proteins. Indeed, the fraction of the proteome quantified constitutes the vast majority of the total proteomic mass of a cell and is enriched for essential genes as well as in genes most connected in functional networks. Our dataset captures many of the fundamental processes. Yet, results related to low abundant proteins are missed by this approach. Overall, our study clearly highlights that the dependency between transcript and protein levels is complex, pointing to the importance of post-transcriptional regulation of protein abundance. Proteome and transcriptome are indeed two distinct layers of gene regulation, which need to be further explored to understand the genotype-phenotype relationship.
The CVs of QCs and wild isolates samples were calculated and had a median CV of 15. 15% and 34. 21%, respectively (Figure  S14 ; table  S22 ). Batch correction was carried out at the precursor level using median batch correction, which consists in bringing the median value of the precursors in the different batches to the same level. Proteins were then quantified from the peptide abundance using the maxLFQ 67 function implemented in the DIA-NN R package. The resulting dataset consists of 630 proteins for 942 strains. We imputed the missing value for further exploration using the KNN imputation method from the impute R package  68.
Unless specified, all the analyses described below use the quantile normalized transcriptomic and proteomic data. We recomputed the raw protein abundance coefficient of variation (CV) of each gene by dividing the standard deviation by the mean (using the non-normalized abundance) and transformed it to a percentage. Based on the CV, we performed a functional exploration by gene set enrichment analysis (GSEA)  69  using the fgsea R package  70  for the gene ontology annotation  71, 72  to detect cellular pathways with a conserved regulation across the population. The within-and across-gene mRNA-protein correlation was performed for each gene or each isolate using a Spearman correlation test. We selected the genes with a mRNA-protein correlation index higher than 0. 42 (> 95% percentile) and performed gene ontology (GO) enrichment analysis using the biological process (BP) database using the topGO R package  73.
We performed the GWAS using either the transcriptome log2 transformed TPM data or the protein abundance. For each dataset, we performed two separated GWAS, one based on SNP as genotype, and one based on the CNV as genotype. The SNP GWAS was run with total of 84,633 SNP displaying a minor allele frequency (MAF) > 5% and that were not located in the telomeric regions (< 20kb away from the chromosome ends). The CNV GWAS was run on a total of 1,019 CNV (MAF > 5%). We used the SNP matrix for both SNP and CNV GWAS, thus evaluating the kinship between the isolate to account for the population structure. We set a phenotype-specific p-value threshold using 100 permutation tests where the phenotypes were randomly permutated between the isolates. We use the 5% lowest p-value quantile from these permutation tests to define the significance threshold.
We finally scaled the significance thresholds of the CNV GWAS to account for the size difference between the SNP and CNV matrices. Regarding the SNP GWAS, the detected QTL were filtered to avoid false positives detection due to linkage disequilibrium among the SNP as described previously  44. This resulted in the filtration of 81 eQTL and 131 pQTL (out of respectively 677 and 729 QTL). The QTL were considered as "local" QTL when they were located 25 kb around their affected phenotype. We also sought to detect QTL hotspots in both transcriptome and proteome GWAS. We defined a hotspot as a concentration of at least 4 QTL in a 20 kb window. We compared the protein turnover rate  45  of obtained on 619 proteins encompassed in our dataset to see whether turnover rate had an impact on the overlap between SNP-eQTL and SNP-pQTL.
Acknowledgements This work was supported by a  National Institutes of Health (NIH)  grant  R01  ( GM147040-01 ) and a  European Research Council (ERC)   Consolidator grant  ( 772505 ) to J. S. And E. T. was supported by the  PhD Joint Programme CNRS & Weizmann Institute  and a fellowship from the medical association  la Fondation pour la Recherche Médicale  ( FDT202204014796 ). J. S. is a Fellow of the  University of Strasbourg Institute for Advanced Study (USIAS)  and a member of the  Institut Universitaire de France.
Deconvolution methods infer levels of immune and stromal infiltration from bulk expression of tumor samples. These methods allow projection of characteristics of the tumor microenvironment, known to affect patient outcome and therapeutic response, onto the millions of bulk transcriptional profiles in public databases, many focused on uniquely valuable and clinically-annotated cohorts. Despite the wide development of such methods, a standardized dataset with ground truth to evaluate their performance has been lacking. We generated and sequenced in vitro and in silico admixtures of tumor, immune, and stromal cells and used them as ground truth in a community-wide DREAM Challenge that provided an objective, unbiased assessment of six widely-used published deconvolution methods and of 22 new analytical approaches developed by international teams.
g. , memory, naïve, and regulatory CD4+ T cells (Fig. 1A ). To facilitate benchmarking and create a ground truth dataset, we generated in vitro and in silico expression profiles of cell populations admixed at predefined ratios. We isolated immune cells from healthy donors and obtained stromal, endothelial, and cancer cells from cell lines (Tables  S1  and S2 ; Methods). Cell type-specific marker expression was confirmed on the purified cells through RNA sequencing (Fig. S1 ). To robustly test algorithm performance across diverse conditions, we defined mixing proportions and generated admixtures from them, grouped into one of six datasets according to whether they: (1) included breast or colon cancer cells; (2) had proportions that were unconstrained or constrained by biologically-reasonable expectations ("biological" distribution); or (3) were created in silico or in vitro (Methods; Tables  S3-S8 ).
This resulted in a total of 96 in vitro admixtures and 140 in silico admixtures, with at least 18 admixtures in each dataset to ensure adequate sample size (Table  S9 ). We generated in silico admixtures as a linear combination of the mixing proportions and the purified expression profiles and in vitro admixtures by extracting RNA from the purified cells, mixing them at the specified proportions, and sequencing (Fig. 1B ). We provided participants with GEO microarray and RNA-seq sample identifiers curated as to inferred purified cell type for use in training (Tables  S10  and S11 ). Significantly, these were publicly available samples that did not include those used in generating our own admixtures. Methods were evaluated against the admixtures by correlating the predictions of cell type proportions with the predefined ("ground truth") proportions, independently for each cell type.
Since the aim of the Challenge concerns the microenvironment populations, we only assessed participants' predictions based on inference of immune and stromal cell content, with the admixed cancer cells effectively treated as contaminating noise. In order to rank methods, we defined an aggregate score that averaged correlations across cell types and validation datasets (Methods). Methods were first assessed using a "primary" Pearson correlation-based score. Statistical ties were resolved relative to the top-performing method (as determined by a Bayes factor; Methods) by a "secondary" Spearman correlation-based score. To account for sampling variability, we reported these Pearson-based (r) and Spearman-based (ρ) scores as their means across bootstraps (Methods). In silico admixtures were then defined as a linear combination of these signatures and specified ratios (unconstrained or constrained according to biologically reasonable expectation).
Several core algorithmic approaches were common across submissions, including those based on non-negative least squares (NNLS; 6 in the coarse-grained and 6 in the fine-grained sub-Challenge, respectively; Fig. 2A, B ) and summarization (SUM; 5 and 4). Nevertheless, there was wide methodological diversity amongst top performers. CIBERSORTx uses ν-SVR to simultaneously solve for all fractional abundances relating admixed and purified expression profiles using a signature matrix of ~525 differentially expressed genes spanning 22 immune cells types (LM22). 19, 32 Aginome-XMU, published subsequent to the Challenge 33 , utilizes a neural network composed of an input layer, five fully-connected hidden layers, and an output layer (Supplemental Methods; https://github. com/xmuyulab/DCTD_Team_Aginome-XMU; Table  S12 ).
The network effectively applies feature selection automatically and was trained here using synthetic admixtures. DA_505 applies a rank-based normalization, selects features by applying random forests to synthetic admixtures, and ultimately applies regression to predict abundance of each cell type independently (Supplemental Methods; https://github. com/martinguerrero89/Dream_Deconv_Challenge_Team_DA505; Fig. S2 ; Table  S13 ). mitten_TDC19 calculates a summarization score as the sum of the expression of selected markers, with the cell type-specific markers first nominated from expression profiles in purified bulk data or identified  34, 35  from single-cell data expression profiles and then prioritized according to their correlation with that cell type's proportion over synthetic admixtures (Supplemental Methods; https://github.
com/sdomanskyi/mitten_TDC19; Table  S14 ). Finally, Biogem, based on a previously published method,  36  uses robust linear modeling to perform deconvolution and differential expression-based feature selection to define the purified expression profiles (Supplemental Methods; https://github. com/giannimonaco/DREAMChallenge_Deconvolution; Fig. S3 ; Table  S15 ). Hence, despite their algorithmic differences, three of the top-performing methods were trained using synthetic admixtures, generated in silico from publicly available purified expression profiles (Table  S16 ). Importantly, the purified profiles that we created to generate Challenge admixtures were not made available to participants. Method performance improved for most teams over three allowed submissions (Fig. S4 ), as they were permitted to revise their method with each submission.
Deconvolution performance differs by cell type In addition to a ranked average performance across cell types, we assessed methods in their ability to predict individual cell type levels within an admixture (Fig. 3  and Figs S6-S15 ). All major lineage cell types could be predicted robustly by at least one method (max row of Fig. 3A  and Figs S8A, S12A , and S15A). For example, though CD4+ T cell levels were most difficult to predict on average, even these could be predicted with a Pearson correlation of 0. 86 by CIBERSORTx and xCell. Neutrophil levels were predicted best on average, with 18 of 28 methods having a Pearson correlation of at least 0. 90. Pearson correlation of method (left axis) prediction versus known proportion from admixture for each cell type (bottom axis). Pearson correlation is first averaged over validation dataset and then (n=1,000) over bootstraps (Methods) and subsequently averaged over coarse-and fine-grained sub-Challenges for cell types occurring in both.
Black entry indicates cell type not predicted by corresponding method. Bottom two rows ("mean" and "max") are the mean and maximum correlation, respectively, for corresponding cell type across methods. Rightmost column ("mean") is mean correlation for corresponding method across predicted cell types. Highest correlations for each cell type highlighted in bold italics. (B) Performance (Pearson correlation; x axis) of comparator baseline methods and participant methods ranking within the top three in either or both sub-Challenges (y axis) for each cell type (facet label). Distribution of Pearson correlations over bootstraps (n=1,000; Methods), computed as average over validation datasets and subsequently over coarse-and fine-grained sub-Challenges for cell types occurring in both. Blank row indicates cell type not reported by the corresponding method.
Our results provide the community cell type-level method performance to assist in making that decision, as well as novel datasets for benchmarking new approaches. Prediction performance was lowest for memory and naïve CD4+ T cells, macrophages, and monocytes, despite potential strides made by novel Challenge methods for these subpopulations. We can partially diagnose potential causes of this reduced performance by leveraging the sequenced purified populations to assess specificity, i. e. , "spillover" from the purified population to others, and sensitivity, assessed as the limit of detection of in silico spike ins. We observed that memory CD4+ T cells had poor specificity and sensitivity -ranking last in both metrics, with a median spillover across methods of 21%. Their best-case limit of detection was 22%, at which level they could be predicted above background by the most sensitive method.
Training data curation Participants were provided with a curated list of purified expression profiles in GEO. 40 GEO annotations were queried using regular expressions corresponding to cell populations of interest (e. g. , with patterns "T. reg", "regulatory", and "FOXP3" for regulatory T-cells). Specifically, GEO annotations for fields "source_name" or involving "characteristic" (e. g. , "characteristics_ch1") were accessed via GEOmetadb in R. 41  Cell type patterns are available at: https://github. com/Sage-Bionetworks/Tumor-Deconvolution-Challenge/blob/master/scripts/training-data-curation/phenotypes-to-query. tsv. Matches identified via grepl were manually curated, resulting in tables associating cell populations with GEO samples. These were further summarized according to dataset, listing the cell populations and the number of cell populations represented by each dataset.
This was intended to help participants prioritize datasets representing many or multiple cell types of interest. Per-sample and per-dataset tables were created separately for microarray (Table  S10   ), specifically the HiSeq (with pattern 'hiseq' in the 'title' field) or NextSeq (with pattern 'nextseq' in the 'title' field) platforms. Participants were provided only identifiers of GEO datasets (GSMxxx) and samples (GSExxx). In particular, cross dataset normalization to account for batch effects was not performed, but rather was left to participants.
Biological admixture generation Biologically-constrained admixtures were defined such that cell population proportions were within biologically-reasonable limits, in particular those detected by CyTOF in PBMCs and aggregated in the 10,000 Immunomes ("10KIP") database  42  (downloaded on July 9, 2018), observed in single-cell (sc)RNA-seq studies of breast cancer ("Azizi" 43 ) or melanoma ("Tirosh" 44 ), or inferred by CIBERSORT in the "Thorsson" TCGA pan-cancer study. 45  none of these resources included all (coarse-or fine-grained) cell types to be deconvolved, several were combined in a hierarchical fashion. Two such hierarchical models, one based on the Thorsson study and the second based on the Azizi study, were created. At each level of the hierarchy, the models defined a minimum and maximum proportion for each population relative to its parent population.
The minimum and maximum proportions for a particular cell type in a particular dataset were defined as two standard deviations above (~98th percentile) or below (~2nd percentile), respectively, the mean of proportions observed for that cell type in that dataset. The root of the model corresponds to the admixture of n cell populations. In both hierarchical models, the entire admixture was partitioned into cancer cell, leukocyte, and nonleukocyte stromal compartments, with minimum and maximum proportions for each compartment defined using the Thorsson study. Specifically, from the stromal fraction (SF), or total non-tumor cellular fraction, and the leukocyte fraction (LF) defined by Thorsson, we define the cancer cell proportion as 1 -SF, the leukocyte proportion as LF, and the non-leukocyte stromal proportion as 1 -SF -LF.
The biologically-constrained admixtures were generated using the "Hit and Run" Markov Chain Monte Carlo (MCMC) method for sampling uniformly from convex samples defined by linear (equality and inequality) constraints, as implemented in the hitandrun library in R (https://github. com/gertvv/hitandrun). The system of linear constraints included a variable for each of the n populations. As in the unconstrained admixtures, the corresponding n proportions sum to p <= 1, with p < 1 allowing the remaining 1-p proportion to be allocated to an (n+1) st cell type. The resulting equality constraint was passed to the solution. basis function, whose output was in turn passed to the createTransform function. 2n linear inequality constraints were defined from the minimum and maximum proportions of each of the n populations.
Code to generate the (unconstrained and biological) in silico validation admixtures is in analysis/validation_data/qc/generate-validation-in-silico-admixtures. R. The admixture expression data (i. e. , represented as TPM and as read counts) and (ground truth) admixtures are in the same Synapse folders as the corresponding in vitro data -i. e. , syn21821096 and syn21820011, respectively. They are the datasets designated 'AA', 'AB', 'AE', and 'AF. ' Participants were told the cancer cell type included in each dataset (BRCA or CRC), which was BRCA for AA and AE and CRC for AB and AF.
Deconvolution method scoring and comparison Pearson and Spearman correlation-based scores were calculated hierarchically for a given method : For each cell type and validation dataset (i. e. , DS1, DS2, DS3, DS4, AA, AB, AE, and AF), the correlation between the values predicted by and the ground truth was calculated. These correlations were then averaged over all cell types to define the score of method for dataset. These dataset-level scores were finally averaged over all datasets to define the aggregate score for method. To assess scoring differences in the primary metric between a top-performing method and another method , we computed a Bayes factor over 1000 bootstrap samples and considered as indicating a significant difference. More specifically, we bootstrap sampled (i. e. , sampled with replacement) prediction scores separately within each dataset (i.
e. , the mean of the column corresponding to cell type c in Fig. 7A  that excludes elements in which c is in the sample corresponding to the row). These spillovers were then averaged over sub-Challenges and the resulting distributions were plotted in Fig. 7B. Distributions of spillovers over cell types are plotted for each method in the coarse-(Fig. 7C ) and fine-grained (Fig. 7D ) sub-Challenges. Code to generate the format the purified expression profiles is in analysis/specificityanalysis/create-spillover-dataset. R. The processed expression data (i. e. , represented as TPM and as read counts) and the (ground truth) admixtures are in Synapse folder syn22392130.