Loss-function learning for digital tissue deconvolution
Franziska Görtler Stefan Solbrig Tilo Wettig Peter J. Oefner Rainer Spang Michael Altenbuchinger
q-bio.QM

Symmetry breaking in planar and maximal outerplanar graphs 
    
Saeid Alikhani ^[Corresponding author]
Samaneh Soltani 

    December 30, 2023
=============================================================


Background: The gene expression profile of a tissue averages the expression profiles of all cells in this tissue. Digital tissue deconvolution (DTD) addresses the following inverse problem: Given the expression profile y of a tissue, what is the cellular composition c of that tissue?
If X is a matrix whose columns are reference profiles of individual cell types, the composition c can be computed by minimizing ℒ(y-Xc) for a given loss function ℒ. Current methods use predefined all-purpose loss functions. They successfully quantify the dominating cells of a tissue, while often falling short in detecting small cell populations. 


Results: Here we learn the loss function ℒ along with the composition c. This allows us to adapt to application-specific requirements such as focusing on small cell populations or distinguishing phenotypically similar cell populations. Our method quantifies large cell fractions as
accurately as existing methods and significantly improves the detection of small cell populations and the distinction of similar cell types.




§ INTRODUCTION

Different tissues of the body have different cellular compositions. The composition of tumor tissue is different from that of normal tissue. 
Also, when comparing two tumor tissues, their cellular composition can differ greatly. 
The relatively small populations of tumor-infiltrating immune cells are of particular importance. 
They affect progression of disease  <cit.> and success of treatment <cit.>. Immune therapies block communication lines between tumor cells and 
infiltrating immune cells. Whether they are successful or not depends on the presence, quantity, and molecular sub-type of the infiltrating 
immune cells <cit.>. Immune-cell populations are typically small, and their molecular phenotype can be difficult to observe under the microscope.
Single-cell technologies such as fluorescence-activated cell sorting (FACS; <cit.>), cytometry by time-of-flight 
(CyTOF; e.g. <cit.>), and single-cell RNA sequencing <cit.> assess molecular features on the single-cell level and can thus 
be used to determine the cellular tissue composition experimentally. 

A more cost- and work-efficient alternative to single-cell assays is a combination of bulk-tissue gene expression profiling with 
digital tissue deconvolution (DTD) <cit.>.
DTD addresses the following inverse problem: given the bulk gene expression profile y of a tissue, what is the cellular composition c of 
that tissue? 
Supervised DTD assumes that there is a matrix X whose columns are reference profiles of individual cell types. 
The composition c of y can be computed by minimizing ℒ(y-Xc) for a given loss function ℒ. 
Competing DTD methods use different predefined all-purpose loss functions ℒ and different estimation algorithms 
to distil c from y and X. 

The practical objective of DTD is to estimate c correctly, while the formal objective of common DTD algorithms is to estimate y correctly. If tissue expression profiles were exact mixtures of reference profiles, existing methods should work perfectly. They are not and this causes problems:




  (1)Collections of references profiles can be incomplete.
There might be cells in the tissue that are not represented by the reference profiles. In that case the global DTD problem is not solvable, and DTD-algorithms will compensate for the contributions of these cells by increasing the frequencies of other cell types.



  (2)Small cell fractions are hard to quantify.
From a practical point of view this is probably the most important point, and improvements are needed badly. Immunological cell populations in a tumor are small, 
but they may determine the reaction of a tumor to immunotherapy. Therefore, DTD algorithms must use faint signals from small cell populations more effectively.


  (3)Some cell types can hardly be distinguished by their expression profiles.
The profile of an epithelial cell differs greatly from that of a lymphoid cell. For two immunological sub-entities of CD8+ T cells the differences are more subtle. The more similar two cell types are, the more similar are their expression profiles, and the more difficult is their distinction.


In summary, different applications need different approaches. One way to adapt the estimation of c is to adapt the loss function  ℒ.
If the focus of an application is on a predefined set of cell types, genes that are informative to distinguish exactly these cells should dominate ℒ. This
is even more important if the focus is on small cell populations, the faint signals of which must not be suppressed. 
Unfortunately, it is not clear a priori which genes to ignore and which to focus on.  




§ METHODS




 §.§ Notations

Let X∈ℝ^p × q be a matrix with cellular reference profiles X_·, j in its columns, where the dot stands for all row indices. 
X_ij is the reference expression value of gene i in cells of type j, 
p the number of genes, and q the number of cell types in X, respectively. 
We further introduce a matrix Y∈ℝ^p× n with bulk profiles of n cell mixtures 
Y_·, k in its columns and a matrix C∈ℝ^q× n with the cellular compositions 
of the mixtures C_·, k as columns. 




 §.§ Loss-function learning

Following the established linear DTD algorithms, we approximate the mixture Y_·,k by a linear 
combination of reference profiles (the columns of X) with C_.,k as weights and estimate the 
composition of the k-th mixture C_·,k by minimizing 

    ℒ_g(Y_·,k-XC_·,k) ,

where 

    ℒ_g= || diag(g) (Y_·,k-XC_·,k) ||^2_2  .

In contrast to standard DTD algorithms, which determine g by prior knowledge or separate statistical analysis, 
we will learn g directly from data.
To this end we assume that we have a training set of mixtures Y_·,k 
from a specific application context with known cellular compositions 
C_·,k.  
The entries of g are the gene weights that define the 
loss function.
We want to learn g from the training data such that minimizing 
ℒ_g(y-Xc) with respect to c
yields accurate quantifications of cell populations for future samples with 
similar characteristics as those used for training.  

Our method has two nested objective functions: An outer function L(g) and 
an inner function ℒ_g, which is here given by equation (<ref>). 
L evaluates discrepancies 
between the estimated and the true cellular frequencies of cell types across samples:

    L(g) =-∑_j=1^qcor(C_j,·,Ĉ_j,·(g))
       g_i≥0  ||g||_2=1 ,

where the Ĉ_j,·(g) are the estimates of C_j,· given g.
To evaluate L( g)  we need to calculate all Ĉ_j,·(g), which requires 
optimizing ℒ_g with respect to all C_·, k.
Note that if  ĝ is a minimum of L, so is αĝ. The
constraint || g||_2=1 is thus needed to ensure unique solutions.

The minimum of ℒ_g can be calculated analytically, yielding

    Ĉ(g) = (X^T Γ X)^-1 X^T Γ Y

with Γ=diag(g). Inserting this term into L leaves us with a single optimization problem in g.
We minimize L  by a gradient-descent algorithm.
Let μ_j and σ_j be the mean and standard deviation of C_j,·, respectively. 
We obtain the gradient 

    ∂  L(g)/∂ g_i   =   ∑_j=1^q ∑_k=1^n 1/σ_jσ̂_j(cov(C_j,·,Ĉ_j,·)/nσ̂_j^2(Ĉ_jk - μ̂_j)
    	   - 1/n(C_jk- μ_j) )
    	   ∂Ĉ_jk(g)/∂ g_i
			
with 

    ∂Ĉ(g)/∂ g_i = 
     (X^TΓ X)^-1 X^T δ(i) (1-X  (X^T Γ X)^-1 X^T Γ)Y ,

where δ(i)∈ℝ^p× p is defined as

    δ(i)_jk= {[      1 i=j=k,;      0        ]

The constraints ||g||_2=1 and g_i≥0 were incorporated  by normalizing g by its length and by restricting the search space to g_i≥0. 





§ RESULTS



 §.§ DTD of melanomas

 
For both training and validation we need expression profiles of cellular mixtures of known composition. We used expression data
of melanomas whose composition has been experimentally resolved using single-cell RNAseq profiling <cit.>.
These data included 4,645 single-cell profiles from 19 melanomas. The cells were annotated as   
T cells (2,068), B cells (515), macrophages (126), endothelial cells (65), cancer-associated fibroblasts (CAFs) (61), 
natural killer (NK) cells (52), and tumor/unclassified (1,758). 
The first 9 melanomas  defined our validation cohort and the remaining 10 our training data. 

First, data were transformed into transcripts per million. Then, for each cell  cluster we sampled 20% of single-cell profiles in the training data, summed them up,  normalized them
to a common number of counts, and removed them from the training data. This yielded reference profiles X_·, j. The 1,000 genes with the highest variance across all reference profiles were used to train models. 

The sum of all single-cell profiles of a melanoma gave us bulk profiles. In addition, we generated a large number of artificial 
bulk profiles by randomly sampling single-cell profiles and summing them up. All bulk profiles 
were normalized to the same number of reads as those in X_·, j. 









 §.§ Loss-function learning improves DTD accuracy in the case of incomplete reference data


We generated 2,000 artificial cellular mixtures from our training cohort. For each of these mixtures, we randomly drew 100 
single-cell profiles, summed up their raw counts, and normalized them to a fixed number of total counts. Analogously, we 
generated 1,000 artificial cellular validation mixtures. 

Then, we restricted X to three cell types (T cells, B cells, and macrophages). Hence endothelial cells, CAFs, NK cells and tumor/unclassified cells in the mixtures are not represented in X. For standard DTD with g=(1,…,1), we observed correlation coefficients
of 0.70 (T cells), 0.39 (B cells), and 0.52 (macrophages) between true and estimated cell population sizes for the validation mixtures. These improved to 0.86 (T cells), 0.89 (B cells), and 0.83 (macrophages) for loss-function learning, after we ran 1000 iterations of the gradient descent algorithm on the training data.

To test the limits of the approach, we excluded all but the macrophages, which account for less than 3% of all cells, from 
the reference data X. We observed, that standard DTD broke down, while loss-function learning yielded a model that predicted macrophage abundances 
that still correlated well (r=0.84) with the true abundances (Figure  <ref>).  









 §.§ Loss-function learning improves the quantification of small cell populations

We generated data as above but this time controlling the abundance of B cells in the simulated mixtures at
0 to 5 cells, 5 to 15, 15 to 30, 30 to 50, and 50 to 75 out of 100 cells. Not surprisingly, small fractions of B cells were harder to quantify than large ones. Loss-function learning
improved the accuracy for all amounts of B cells, but the improvements were greatest for small amounts (Figure <ref>a). 
With only 0 to 5 cells in a mixture the accuracy improved from r=0.22 to r=0.79. Furthermore, we observed that loss-function learning on small B-cell proportions yielded a model that was highly predictive of B-cell contributions over the whole spectrum (Figure <ref>a green stars). 

If we compare the top-ranked genes of the model learned for the small B-cell population (Figure <ref>b) to that of the macrophage-focussed simulation (Figure <ref>c), we observe that the former still comprises marker genes to distinguish all cell types, while the latter focusses on genes that characterize macrophages. 







 §.§ Loss-function learning improves the distinction of closely related cell types


The cell types that were annotated by <cit.> displayed very different expression profiles.
If we are interested in T-cell subtypes such as CD8+ T cells, CD4+ T-helper (Th) cells, and regulatory T cells (Tregs),
reference profiles are more similar and DTD is more challenging.
We subdivided the fraction of annotated T-cell profiles as follows: 
all T cells with positive CD8 (sum of CD8A and CD8B) and zero CD4 count were labelled CD8+ T cells (1,130).
Vice versa, T cells with zero CD8 and positive CD4 count were labelled CD4+ T cells (527).
These were further split into Tregs if both their FOXP3 and CD25 (IL2RA) count was positive (64),
and CD4+ Th cells otherwise (463).  
T cells that fulfilled neither the CD4+ nor the CD8+ criteria (411) contributed to the mixtures, but were not assessed by DTD.
We augmented the reference matrix X, here consisting of T cells, B cells, macrophages, endothelial cells, CAFs and NK cells, by these cell types, replacing the original all T-cell profile with the more specific profiles for CD8+, CD4+ Th and Tregs. Then we simulated 2,000 training and 1,000 test mixtures as described above. 


For standard DTD with g=1 we observed correlation coefficients
of 0.19 (CD4+ Th), 0.53 (CD8+), and 0.08 (Tregs) between true and estimated cell population sizes.
These improved to 0.58 (CD4+ Th), 0.78 (CD8+), and 0.57 (Tregs) for our method 
(Figure  <ref>).







 §.§ Loss-function learning is beneficial even for small training sets, and the performance improves as the training dataset grows


 We repeated the simulation in subsection <ref>, but varied the size of the training dataset. We observed that loss-function learning improved accuracy for 
training datasets as small as 15 samples. Moreover, with more training data added the boost in performance grew and saturated only for training sets with more than 1,000 samples 
(Figure <ref>).






 §.§ HPC-empowered loss-function learning rediscovers established cell markers and complements them by new discriminatory genes for improved performance


Here, we introduce a final model, optimized on the 5,000 most variable genes. For this purpose, we generated 25,000 training mixtures from the melanomas of the training data. 
With standard desktop workstations the solution of this problem was computationally not feasible. A single computation of the gradient took 16 hours (2x Intel Xeon CPU [X5650; Nehalem Six Core, 2.67 GHz], 148 Gb RAM), and this needs to be computed several hundred times until convergence.
Therefore, we developed a High-Performance-Computing (HPC) implementation
of our code by parallelizing equations (<ref>) and (<ref>) with MPI, using the
pbdMPI library (<cit.>, <cit.>) as an interface. Furthermore, we linked R with
the Intel Math Kernel Library for threaded and vectorized matrix
operations. We ran the algorithm on 25 nodes of our QPACE 3 machine <cit.> with 8 MPI tasks per node and 32 hardware threads per task, where each thread can use two AVX512 vector units.
In 16 hours 5,086 iterations were finished, after which the loss (<ref>) was stable to within 1%.




The high-performance model includes several genes, whose expression is characteristic for the cells distinguished in the present study. These include, among others, 
the CD8A gene, which encodes an integral membrane glycoprotein essential for the activation of cytotoxic T-lymphocytes <cit.> and the protection of 
a subset of NK cells against lysis, thus enabling them in contrast to CD8- NK cells to lyse multiple target cells <cit.>. As evident from Figure <ref>, 
NK cells are clearly set apart from all the other cell types studied by the expression of the killer cell lectin like receptor genes KLRB1, KLRC1, and KLRF1 <cit.>. B cells, on the other hand, are clearly characterized by the expression of (i) CD19, which assembles with the antigen receptor of B lymphocytes and influences B-cell selection and differentiation <cit.>, (ii) CD20 (MS4A1), which is coexpressed with CD19 and functions 
 as a store-operated calcium channel <cit.>, (iii) B Lymphocyte Kinase (BLK), a src-family protein tyrosine kinase that plays an important role in B-cell receptor signaling and phosphorylates specifically (iv) CD79A at Tyr-188 and Tyr-199 as well as CD79B (not among the top 150 genes) at Tyr-196 and Tyr-207, which are required for the surface expression and function of the B-cell antigen receptor complex (<cit.>), and (v) BLNK,  which bridges BLK activation with downstream signaling pathways <cit.>. 
The expression of FOXP3 is also highly cell specific. FOXP3 distinguishes regulatory T cells from other CD4+ cells  and functions as a master regulator of their 
 development and function (<cit.>). Finally, CD4+ T-helper (Th) cells are distinguished indirectly from all the other aforementioned lymphocytes by the lack of expression of cell 
 type-specific genes. In contrast to lymphocytes, macrophages, cancer-associated fibroblasts (CAFs), and 
 endothelial cells, which line the interior surface of blood vessels and lymphatic vessels, are characterized each by a much larger number of
 genes. Exemplary genes include CD14, CD163, MSR1, STAB1, and CSF1R for macrophages. The monocyte differentiation antigen CD14, for instance, mediates the innate immune response 
 to bacterial lipopolysaccharide (LPS) by activating the NF-κB pathway and cytokine secretion (<cit.>), while the colony stimulating factor 1 receptor (CSF1R) acts as a receptor for the
 hematopoietic growth factor CSF1, which controls the proliferation and function of macrophages (<cit.>). CAFs, on the other hand, are distinguished by the expression of genes encoding 
 extracellular matrix proteins such as fibulin-2 (FBLN2) and fibulin-3 (EFEMP1), various collagens (COL1A1, COL3A1, COL6A1, COL6A3), versican (VCAN), a well known mediator of cell-to-cell and cell-to-matrix interactions (<cit.>) that plays critical roles in 
 cancer biology (<cit.>), as well as the matrix metalloproteinases MMP1 and MMP2, two collagen degrading enzymes that allow cancer cells to migrate out of the primary
 tumor to form metastases (<cit.>). 
 Noteworthy is also GREM1, an antagonist of the bone morphogenetic protein pathway. Its expression and secretion by stromal cells in tumor tissues promotes the survival and proliferation of cancer cells <cit.>. 
 Genes characteristic for endothelial cells include among others CDH5, a member of the cadherin superfamily essential for endothelial adherens junction 
 assembly and maintenance (<cit.>), the endothelial cell-specific chemotaxis receptor (ECSCR) gene, which encodes a cell-surface single-transmembrane domain glycoprotein that plays a role in endothelial cell migration, apoptosis and proliferation (<cit.>), claudin-5 (CLDN5), which forms the backbone of tight junction strands between endothelial cells <cit.>, and the von Willebrand factor (VWF), which mediates the adhesion of platelets to sites of vascular damage by binding to specific platelet membrane glycoproteins and to constituents of exposed connective tissue (<cit.>). 


We discussed 29 genes of the top 150 shown in Figure <ref>. These genes have a total weight of 28% of all 5,000 gene weights (calculated as ĝ_i ×var(X_i,·)). Our algorithm complements this gene set with additional genes, including many that were, to our knowledge, not yet used to characterize cell types. An interesting example is CXorf36 (DIA1R), which has been described as being expressed at low levels in many tissues and deletion and/or mutations of which have been associated with autism spectrum disorders <cit.>. However, nothing is known about its function to date. Therefore, its observed overexpression in endothelial cells may provide an important clue for future study on its function.



 §.§ Loss-function learning shows similar performance as CIBERSORT for the dominating cell populations
and improves accuracy for small populations and in the distinction of closely related cell types
 

Next we compared our model trained in subsection <ref> to a competing method. For this, we generated 1,000 test mixtures from our validation melanomas. We chose CIBERSORT <cit.> for comparison, because it was consistently among the 
best DTD algorithm in a broad comparison of five different algorithms on several benchmark datasets <cit.>. We ran CIBERSORT on the test mixtures, using two distinct approaches: first, we uploaded 
our validation data to CIBERSORT using their reference profiles. The performance is summarized in Figure <ref> as CIBERSORT^a (yellow). We observed that the large population of B cells was estimated accurately, while smaller populations 
were inaccurate (NK cells, Tregs). Next, we uploaded our reference profiles and used the CIBERSORT gene selection (CIBERSORT^b green). We found that highly abundant cell types (B cells and CD8+ T cells) were predicted with high accuracy. However, 
the distinction of similar cell types such as CD4+ T helper cells and Tregs was compromised, r=0.42 and r=0.42, respectively. 
Similarly, predictions for the small populations of CAFs were compromised. That might be explained 
by the fact that CIBERSORT does not take into account their distinction and thus appropriate marker genes might be missing. 
In a direct comparison to CIBERSORT our method showed similar or better performance.  














 §.§ Loss-function learning improves the decomposition of bulk melanoma profiles

All mixtures discussed so far were artificial because only 100 single-cell profiles were chosen randomly. They might differ significantly from mixtures in real tissue. Therefore, we generated 19 full bulk melanoma profiles by summing up the respective single-cell profiles. These should reflect bulk melanomas <cit.>. Our predictions are contrasted with the true proportions in Figure <ref>. Only the predictions for Tregs were compromised with r=0.48, while the predictions for all other cell types were reliable with correlations ranging from r=0.70 (CD4+ Th) to r=0.99 (CAFs) on the validation melanomas.







§ DISCUSSION


We suggest loss-function learning for digital tissue deconvolution to adapt the deconvolution
algorithm to the requirements of specific application domains. The concept is similar to an
embedded feature-selection approach in regression or classification problems. In both contexts
feature selection is directly linked to a prediction algorithm and not treated as an independent
preprocessing step.    

We described and tested a specific instance of loss-function learning using squared residuals
for ℒ_g. The concept is not limited to this type of inner loss function
and can also be used in combination with other loss functions such as those from
penalized least-squares regression <cit.>, l_1 regression, 
or support vector regression <cit.>.
However, the least-squares loss function allowed us to state the outer optimization problem in a closed analytical form,
reducing computational burden.  

The outer loss function L evaluates the fit of estimated and true cellular proportions
in the training samples. We chose the correlation of estimated versus true quantities across samples, and no 
absolute measure of deviation such as ||c-ĉ||_2^2. Moreover, we did not require the estimated 
proportions Ĉ_·,k for tissue k to sum up to one.
Consequently, when testing our method we did not look at absolute deviations of true versus estimated cell proportions
but only at their correlation. We cannot reliably infer how many cells of a specific type (e.g., T cells)
are in a tissue (Figure <ref>), nor can we infer whether they constituted 10% or 20% of the cells in this tissue.
However, if we had two tissues and estimated that there were more cells of that type in the first tissue compared to the second, this relation was also found in the true cell populations.
Calibrating the method to absolute cell numbers or percentages within a tissue appears to be difficult.
This might be due to the normalization of the libraries from single cells to a common size, which might not reflect
true biology if, for example, some cell types contain more RNA than others. Post-hoc normalization of the estimated quantities might be an option.

In summary, we introduced loss-function learning as a new machine-learning approach to the digital tissue deconvolution
problem. It allows us to adapt to application-specific requirements such as focusing on small cell populations or delineating similar cell types.
In simulations and in an application to melanoma tissues our method quantified large cell fractions as
accurately as existing methods and significantly improved the detection of small cell populations and the distinction of similar cell types.
 


§ ACKNOWLEDGEMENT

This work was supported by BMBF (eMed Grant 031A428A)  and DFG (FOR-2127  and SFB/TRR-55).



apalike
