dataesr/scientific-paragraphs-categorization huggingface dataset:
As any other constraint, the main reason to use it is providing the resolved profiles with chemical/physical meaning and decreasing the rotational ambiguity linked to the MCR solutions  [10]. As with the implementation of other constraints in the MCReALS framework, such as non-negativity, unimodality or closure, the sparseness constraint can be used on some or all the profiles of the different compounds in the concentration and/or spectral direction. The use of the sparseness constraint does not require that the original raw data are, by nature, sparse, but only some or all of the compound profiles in the underlying bilinear model. Moreover, the way the constraint is implemented does not assume the number of non-null coefficients, nor their position when reproducing the signal. It should be noted that the use of sparseness cannot and should not be generalized and it has to be applied only when the pure compound profiles in the underlying bilinear model fulfill this property.
Multivariate curve resolution e alternating least squares (MCReALS) The development and application of the MCReALS routine has covered the resolution of data sets of very diverse nature [23e26]. It assumes the expression of the Beer-Lambert law so that a bilinear model can be used to decompose a mixed data set into a model related to the pure individual contributions. Doing so allows us to write the system as in  (1) , in which the matrix D (m x n), containing m mixed spectra collected at n spectral channels can be reproduced by a bilinear model of k components, formed by concentration profiles, represented by matrix C (m x k), and their related spectral profiles, represented by the matrix S T (k x n). The matrix E (m x n) represents the error contribution of the model. D ¼ CS T þ E (1) MCReALS is an iterative alternating least-squares algorithm that resolves concentration and spectral profiles under the action of constraints and can be used on a single data set or multiset structures formed by several concatenated data matrices  [9, 26].
Constraints can be very diverse (e. g. , non-negativity, unimodality, closure, etc. ) and can be optionally applied to all or to some of the compounds in the system analyzed in the concentration and/or spectral direction. MCReALS has seen an increased use in the analysis of HSI [27e29]. The main characteristic of images is that they are defined by two spatial directions (x-coordinate and y-coordinate for each pixel) and a spectral direction, i. e. , as a data cube. The spectra are collected for each pixel and contain the chemical information of the constituents present at a given pixel position. Knowing this, it is indeed true that the bilinear model in (1) still holds for the data pixels. This implies that the three-dimensional data cube has to be transformed into a data matrix to be analyzed by MCReALS analysis, as shown in Fig.
This leads to a data size reduction from 56789 m/ z values per pixel to 8627. The data compression is then followed by a baseline correction by asymmetric least squares  [40] , spectral alignment by Correlation Optimized Warping  [41]  and balance of signal intensity variability by normalization (Euclidean norm) and scaling of the spectra. The scaling of the m/z spectral channels was performed by dividing each reading by the standard deviation of the channel readings plus 0. 1% of the maximum signal of the spectrum obtained by the sum of the selected purest mass spectra. After proper preprocessing, it can be seen that the ROI between 1150 and 2000 do not contain relevant information and, hence, this m/z range was removed from further analysis. The final dataset then considers 3750 m/z values per pixel.
Sparseness constraint in MCReALS The sparseness constraint was implemented in the MCReALS command line code (https://mcrals. wordpress. com/download/ mcr-als-command-line/) and all calculations were performed in MATLAB version R2014a (The MathWorks, Natick, MA, USA) on a computer with an Intel(R) Core(TM) i7-4770 CPU @ 3. 40 GHz CPU and 16 GB RAM memory (64-bit  Windows 7). The critical point of applying the constraint during the MCReALS analysis is the sparseness factor, l, which varies according to the degree of natural sparseness of each of the constituents analyzed, since the shape of the signal (e. g. intensity, number of peaks, etc. ) for each of them may be different. Different ways to select the optimal l have been proposed, as for example a cross- validation approach used in Ref. [12] , where the optimal l is selected as the one around the change in root mean square error (RMSE) slope.
Bearing in mind the nature of MS signals, we know that the resolved spectral profiles of the constituents must be sparse signals with only a few significant sharp peaks. Therefore, the use of the sparseness constraint should here improve the definition of m/z spectral signatures and, as a consequence, the overall description of the data set. After the preliminary MCReALS analysis of the bean MSI under non-negativity constraints, the analysis was carried out applying the sparseness constraint to the spectral profiles. Optimal sparseness factors are obtained as described in Section 3. 3. For a better understanding of the results, the sparseness factor chosen for each component (see caption Fig. 5 ) has been 'translated' to a sparseness ratio for the respective components, indicating the number of nonnull values in each spectral channel divided by the total number of spectral channels (in %).
Infrared spectroscopy imaging on a pharmaceutical sample For this data set, the number of components is known to be three and initial estimates were determined by using a SIMPLISMAbased method. An MCReALS analysis on the pre-processed data was first performed by using non-negativity as a constraint, on both the C and S T matrix and normalization of the S T profiles. The resolved distribution maps and spectra can be seen in Fig. 6a  (LOF: 7. 9%; explained variance: 99. 4%). The sparseness ratio for the component distribution maps is defined here as the ratio of the non-null pixels over the total number of pixels and is 55. 4%e 76. 3%e86. 2% for the respective components. When looking at the spectral profiles of the resolved components, the profiles of citric acid and acetyl salicylic acid (components 1 and 2) can be clearly distinguished, since the peak positions correspond with the IR spectra found in the database of the National Institute of Standards and Technology (NIST)  [44].
As a consequence of the presence of zones with lower rank or even being selective, an improvement of the spectral profiles recovered can be noticed. A closer look at the spectral profiles of citric acid and acetylsalicylic acid would indeed reveal an improvement in the intensity ratios of the different peaks with respect to the reference found in the NIST database  [44]. When focusing on the third component, an improvement can also be noticed. The resolution of the peaks in the spectrum has increased, clearly revealing two peaks (at ~1690 cm À1 and ~1725 cm À1 ) that can be assigned to caffeine. Additionally, it is worth noting that when the distribution map of this third component is not constrained by sparseness, the sparseness ratio increases to 93. 2% and the spectrum does not reveal the two caffeine peaks anymore (results not shown).
Concluding remarks A constraint that induces sparseness on the signals has been implemented in MCReALS for the resolution of complex multicomponent systems. It can be used on both spectral profiles and component distribution maps and is thus fully applicable to hyperspectral images. To induce sparseness in the signal fitted, a penalized least squares regression framework is used that constrains the number of non-null coefficients. The sparseness factor, l, is the only parameter to regulate within this framework, making it a flexible and powerful tool to adapt the degree of sparseness to the nature of the different signals and compounds within a data set. Sparseness can be used in combination with other appropriate constraints. The results that are shown correspond to cases where sparsity is desirable in either the resolved concentration profiles or the resolved spectra.
Acknowledgements The authors thank  Paul Eilers  for his involvement in the development of the algorithms and appreciate the fruitful discussions on the topic. C. R. thanks the  Agence National de la Recherche  ( ANR-15-CE09-0020-01 Blink ). C. B. and A. J have received funding from the  European Research Council  under the European Union's  Seventh Framework Program  (FP/2007-2013)/ ERC  Grant Agreement n. 32073 (CHEMAGEB project). S. P. acknowledges the support of the  Spanish government  through project  CTQ2015-66254-C2-2-P.
We used a single-cell RNA-seq dataset of ~11,000 cells from the tumor microenvironment to simulate bulk samples of known cell type proportions, and validated the results using independent, publicly available gold-standard estimates. This allowed us to analyze and condense the results of more than a hundred thousand predictions to provide an exhaustive evaluation across seven computational methods over nine cell types and ~1,800 samples from five simulated and real-world datasets. We demonstrate that computational deconvolution performs at high accuracy for well-defined cell-type signatures and propose how fuzzy cell-type signatures can be improved. We suggest that future efforts should be dedicated to refining cell population definitions and finding reliable signatures.
Moreover, cell type proportions might be biased in scRNA-seq data due to differences in single-cell dissociation efficiencies  (Lambrechts et al. (2018) ). At the same time, methods for gene expression profiling, RNA-seq, and microarrays, have been developed and optimized to be applicable in clinical settings, resulting in a plethora of transcriptomics datasets derived from patient tumor samples. However, while these methods provide transcriptomics data from heterogeneous samples considered as a whole, they do not provide information on their cellular composition, which therefore has to be inferred using computational techniques. Previously, available methods for estimating immune cell contents of tumor samples have been reviewed in terms of methodology  (Finotello and Trajanoski (2018) ;  Newman and Alizadeh (2016) ; Avila  Cobos et al.
quanTIseq comes as an entire pipeline that allows to directly analyze raw RNA-seq data (i. e. FASTQ files of sequencing reads), whereas the other methods only analyze pre-computed, normalized expression data. In this benchmark, we consider only its "deconvolution" module alone, acknowledging that the full quanTIseq pipeline can, in principle, provide higher performance. A limitation of TIMER and xCell is that the results of these methods depend on all samples that are analyzed in a single run, i. e. the same sample can have different scores when submitted together with different other samples. In particular, xCell uses the variability among the samples for a linear transformation of the output score. TIMER uses COM-BAT  (Johnson et al. (2007) ) to merge the input samples with reference profiles.
This is particularly problematic with small datasets and hampers comparability and interpretability of the score. Moreover, xCell does not detect any signal when ran with few non-heterogeneous samples  (Aran (2018) ). Nonetheless, all methods except CIBERSORT (relative mode) support between-sample comparisons and the performance of the methods can be readily compared using Pearson correlation with the gold-standard measure of cell fractions. Here we developed a systematic approach for benchmarking such computational methods and assessed the accuracy of estimating the abundance of nine immune and stromal cell types from bulk RNA-seq samples for seven state-of-the-art approaches specifically designed for immunooncology and provide immuno-oncologists guidance on which method they can apply in their studies.
We use a scRNA-seq dataset of more than 11,000 single cancer, stromal, and immune cells  (Schelker et al. (2017) ) to simulate bulk RNA-seq samples with known compositions and validate the benchmark results using three independent datasets that have been profiled with FACS  (Schelker et al. (2017) ;  Hoek et al. (2015) ;  Racle et al. (2017) ). We assess the performance of the methods by four metrics, namely: (i) predictive performance, (ii) minimal detection fraction, (iii) background predictions and (iv) spillover effects. We demonstrate, that by addressing spillover effects, 'fuzzy' signatures of still not well-understood cell-types can be significantly improved. Finally, we provide an R package, immunedeconv, integrating all of the investigated methods for userfriendly access. To facilitate benchmarking of newly developed methods and signatures, we provide a reproducible pipeline on GitHub, which can be straightforwardly extended to benchmark additional methods.
Results Benchmark on simulated and true bulk RNA-seq samples reveals differences in method performance between cell types We created 100 simulated bulk RNA-seq samples with known cell type proportions using the single cell dataset. For each sample, we individually retrieved and aggregated 500 random immune-and non-immune cells. This approach has been established and successfully applied by  Schelker et al. (2017)  for benchmarking CIBERSORT. Additional consistency checks support that simulated bulk RNA-seq data is not subject to systematic biases (supplementary figures 1-4). We applied the seven methods to these samples and compared the estimated to the known fractions. The results are shown in figure  1a. All methods obtained a high correlation on B cells (Pearson's r > 0. 71), cancer associated fibroblasts (r > 0.
72), endothelial cells (r > 0. 94) and CD8+ T cells (r > 0. 76). Most methods obtain a correlation of r > 0. 68 on macrophages/monocytes, NK cells and overall CD4+ T cells. We observed poor performance in distinguishing regulatory-from non-regulatory CD4+ T cells and in estimating dendritic cells (DCs). We validated the results using three independent datasets, which provide both bulk RNA-seq data and a gold-standard estimate of immune cell fractions using FACS (figure 1b-d and supplementary figure  5 ). In general, the results are highly consistent with the mixing benchmark, with the exception that the methods' performance on CD8+ T cells is worse on the validation data considered altogether (figure  1c ). This can be attributed to the fact that the variance of CD8+ T cell contents between the samples is low and the correlations are therefore not meaningful (see figure  1b  and supplementary figure  5 ).
(b) Performance of the methods on three independent datasets that provide immune cell quantification by FACS. Different cell types are indicated in different colors. Pearson's r has been computed as a single correlation on all cell types simultaneously. Note that only methods that allow both inter-and intra-sample comparisons (i. e. EPIC, quanTIseq, CIBERSORT absolute mode) can be expected to perform well here. (c-d) Performance on the three validation datasets per cell type. Spillover can be attributed to non-specific signature genes Motivated by the hypothesis that background predictions might be due to non-specific signature genes, we asked which cell types lead to methods erroneously predicting a higher abundance of another. We refer to this effect as 'spillover'. We assessed spillover effects using simulated bulk RNA-seq samples containing single immune cell types (figure  3 ) and validated the results using true bulk RNA-seq samples of FACS-purified immune cells (supplementary figure  6 ).
The former two spillover effects are conserved in the validation dataset (supplementary figure  6 ) supporting that the effects are not due to misannotated cells in the single cell dataset. The spillover between DCs and B cells could not be confirmed in the validation dataset, however, we demonstrate that in the single cell dataset, the B cell and DC clusters are both distinct and well-annotated (supplementary figure  7 ). Proceeding analogously to the CAF/macrophage spillover above, we identified six genes, TCL1A, TCF4, CD37, SPIB, IRF8, and BCL11A, that are expressed both in the B cell and DC subpopulations on the mRNA level. Indeed, in the LifeMap Discovery database  (Edgar et al. (2013) ), all six genes are annotated as being expressed in both B cells and plasmacytoid dendritic cells (pDCs).
When excluding these genes from the deconvolution matrices, the spillover effect is substantially reduced by 55 % for EPIC and 78 % for quanTIseq (figure  4b ). Given that the single cell dataset contains pDCs, quanTIseq has been trained on myeloid dendritic cells (mDCs) and EPIC does not consider any kind of DCs (supplementary table 2), it is not surprising that both methods view these genes as B cell markers and are not able to discriminate between pDCs and B cells.
Guidelines for method selection In table 2, we provide guidelines on which method to use for which cell type based on three criteria: interpretability of the score, overall performance, and possible limitations. It is very important to understand the implications of the different scoring strategies. EPIC and quanTIseq are the only methods providing an 'absolute score' that represents a cell fraction. All other methods provide scores in arbitrary units, that are only meaningful in relation to another sample of the same dataset. For this reason, and due to a robust overall performance, we recommend EPIC and quanTIseq for general purpose deconvolution. In practice, absolute scores are not always necessary. For instance, in a clinical trial, relative scoring methods can be used to infer fold changes between treatment and control group or to monitor changes in the immune composition in longitudinal samples.
Immunodeconv R package provides easy access to deconvolution methods We created an R package, immunedeconv, that provides a unified interface to the different deconvolution methods. The package is freely available from GitHub: https://github. com/grst/immunedeconv. We emphasize the reproducibility and extensibility of our benchmark. A snakemake pipeline  (Koster and Rahmann (2012) ), including a step-by-step instruction on how to benchmark an additional method is available from https://github. com/grst/immune_deconvolution_benchmark.
Discussion Quantification of the cellular composition of bulk tumor samples is a challenging computational problem that has been addressed by a variety of methods. For the utility of these methods, it is imperative to fully understand their capabilities and limitations. Here we provide a comprehensive overview and the first cell type-specific quantitative comparison of existing methods applied to RNA-seq data to guide researchers and clinicians in selecting the most suitable approach for the analysis of their samples. We assessed seven state-of-the-art computational methods by applying them to simulated bulk RNA-seq samples derived from a scRNA-seq dataset of ~11,000 single cells, integrated from different experiments  (Schelker et al. (2017) ). We validated the results using true bulk RNAseq samples profiled by FACS from three independent datasets  (Schelker et al.
(2017) ;  Racle et al. (2017); Hoek et al. (2015) ). FACS or IHC measurements are currently considered the gold standard for assessing cell type composition but in public databases only a few RNA-seq samples with matched FACS or IHC reference can be found. Further publiclyavailable matched datasets are needed to enhance our ability to reliably compare computational cell-type quantification methods to gold standard i i "output" -2019/1/27 -10:02 -page 7 -#7 i i i i i i.
Benchmark of cell-type quantification methods Table  2. Guidelines for method selection. For each cell type, we recommend which method to use, highlighting advantages and possible limitations. Overall performance: indicates how well predicted fractions correspond to known fractions in the benchmark. Absolute score: the method provides an absolute score that can be interpreted as a cell fraction. Methods that do not provide an absolute score only support inter-sample comparisons within the same experimental dataset, i. e. the score is only meaningful in relation to another sample. Background predictions: indicates, if a method tends to predict a cell type although it is absent.
Our analysis reveals substantial spillover between DCs and B cells in several methods caused by a subset of marker genes characterized by lower cell-specificity. We validated in a public knowledge base that these genes are indeed markers expressed in both pDCs and B cells on the mRNA level and could significantly reduce the spillover effect by removing the genes from the signature matrices. Similarly, we could substantially reduce the background prediction level of quanTIseq on macrophages/monocytes by excluding non-specific genes. By leveraging recently published large-scale single-cell RNA-seq datasets (Sade-Feldman et al. An open question is if profiles derived from one tumor type also apply to another. Schelker et al. to a certain extent, addressed this question, suggesting that cancer-type specific signatures are required.
However, their analysis is limited to CIBERSORT and two tumor types (melanoma and ovarian cancer). In our benchmark, other methods already achieve nearperfect correlations on some cell-types without requiring cancer-specific signatures, suggesting that universal cell-type signatures are possible. As a next step, recently published single-cell datasets of various cancer types  (Sade-Feldman et al. (2018) ;  Lambrechts et al. (2018) ;  Azizi et al. (2018) ;  Zheng et al. (2017) ;  Guo et al. (2018) ) could be combined with our pipeline to address this question comprehensively.
Single-cell data The dataset of 11,759 single cells from  Schelker et al. (2017)  was obtained through their figshare repository (https://figshare. com/ s/711d3fb2bd3288c8483a). We reproduced their analysis using their MATLAB source code and exported the final dataset to continue our analysis in R. We excluded all cells that were classified as "Unknown" from the downstream analysis. Moreover, we did not consider cells originating from PBMC as we were interested in the methods' performance on cancer samples. Gene expression values are expressed as transcripts per million (TPM).
Validation data We obtained preprocessed bulk RNA-seq and FACS data for 8 PBMC samples from  Hoek et al. (2015)  through personal communication from the quanTIseq authors  (Finotello et al. (2017) ). We obtained bulk RNAseq and FACS data for 4 metastatic melanoma samples from  Racle et al. (2017)  through GEO accession number GSE93722. We obtained bulk RNA-seq and FACS data for 3 ovarian cancer ascites samples  (Schelker et al. (2017) ) from the figshare repository (https://figshare. com/ s/711d3fb2bd3288c8483a). The 3 ovarian cancer ascites samples consist of two technical replicates each. The two replicates are highly consistent for each of the samples (Pearson's r ≥ 0. 98), we, therefore, took the straightforward approach of merging them by the arithmetic mean for each gene. All gene expression values are expressed as TPM.
Cell type mapping For comparing the methods, it is essential to map the cell types of the different methods to a controlled vocabulary (CV). It has to be taken into account, that the different methods resolve the cell types at different granularities. E. g. while CIBERSORT predicts naive-and memory CD8+ T cells, all other methods only predict CD8+ T cells. We address this issue by creating a hierarchy of cell types, and map each cell type from a dataset or method to a node in the cell type tree (see supplementary figure  8  and supplementary table  4 ). If a node (e. g. CD8+ T cell) is to be compared among the different methods, the estimated fraction is computed as follows: (1) if the method provides an estimate for the node, take that value. (2) If the method provides an estimate for all child nodes, sum up all children.