Skip to content

Models

Hierarchical Cell Type Deconvolution Model

hide_deconv.models.HIDE

===================================================== Model for hierarchical cell type deconvolution =====================================================

HIDE

Bases: Module

Hierarchical Cell Type Deconvolution Model.

  • Train the model with .train(...)
  • Predict with .predict(...)
Source code in src/hide_deconv/models/HIDE.py 
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
class HIDE(nn.Module):
    """
    Hierarchical Cell Type Deconvolution Model.

    - Train the model with *.train(...)*
    - Predict with *.predict(...)*
    """

    def __init__(
        self, X_l: list[pd.DataFrame], A_l: list[pd.DataFrame], lambdaNMSE: float = 0.0
    ):
        """
        Constructor of HIDE deconvolution model.

        Parameters
        ----------
        X_l : list[pd.DataFrame]
            List of reference profiles for each cell layer. First element should represent the finest coarsed resolution. (genes x celltypes_l)
        A_l : list[pd.DataFrame]
            List of projection matrices, projecting from the finest coarsed resolution to a higher one. The first element should always
            be an identity matrix. (celltypes_l x celltypes_fine)
        lambdaNMSE : float = 0.0
            Weighting of optional NMSE contribution in general loss.

        """
        super().__init__()

        self.L = len(A_l)
        self.lambdaNMSE = lambdaNMSE

        self.celltype_layer_labels = [X.columns for X in X_l]
        self.gene_labels = X_l[0].index

        self.p, _ = X_l[0].shape
        self.q_l = [len(X.columns) for X in X_l]

        self.A_l = [torch.tensor(A.to_numpy(), dtype=torch.float64) for A in A_l]
        self.X_l = [torch.tensor(X.to_numpy(), dtype=torch.float64) for X in X_l]

        self.g_l = nn.ParameterList(
            [
                nn.Parameter(
                    torch.empty(self.p, dtype=torch.float64).uniform_(0.001, 0.1),
                    requires_grad=True,
                )
                for _ in range(self.L)
            ]
        )

        self.progress_bar = Progress(
            TextColumn("[progress.percentage]{task.percentage:>3.0f}%"),
            BarColumn(),
            MofNCompleteColumn(),
            TextColumn("•"),
            TimeElapsedColumn(),
            TextColumn("•"),
            TimeRemainingColumn(),
        )

        self.epsilon = 1e-8  # small epsilon to prevent errors

    def get_loss(self, C: torch.float64, C_est: torch.float64):
        corr_terms = []
        nmse_terms = []

        # Ensure all datatypes are the same, as PyTorch tends to throw errors
        C = C.to(dtype=self.A_l[0].dtype)
        C_est = C_est.to(dtype=self.A_l[0].dtype)

        for layer in range(self.L):
            A = self.A_l[layer]

            C_l = A @ C
            C_est_l = A @ C_est

            mu_t = C_l.mean(dim=1, keepdim=True)
            mu_h = C_est_l.mean(dim=1, keepdim=True)

            vt = C_l - mu_t
            vh = C_est_l - mu_h

            num = (vt * vh).sum(dim=1)

            # Bugfix for the Intel MKI illegal variable problem: Calculate the denominator separately
            var_t = vt.pow(2).sum(dim=1)
            var_h = vh.pow(2).sum(dim=1)
            den = torch.sqrt(var_t * var_h + self.epsilon)

            r = num / den

            corr_terms.append(-r.mean())

            nmse_num = (C_l - C_est_l).pow(2).sum()
            nmse_den = (C_l).pow(2).sum() + self.epsilon
            nmse_terms.append(nmse_num / nmse_den)

        corr_loss = torch.stack(corr_terms).mean()
        nmse_loss = torch.stack(nmse_terms).mean()
        total_loss = corr_loss + self.lambdaNMSE * nmse_loss

        return total_loss, corr_loss, nmse_loss

    def train(self, Y: pd.DataFrame, C: pd.DataFrame, iter: int = 1000) -> list:
        """
        Trains the HIDE model.

        Parameters
        ----------
        Y : pd.DataFrame
            Bulk training data. (genes x samples)
        C : pd.DataFrame
            True cellular composition at finest coarsed cell resolution. (celltypes x samples)
        iter : int = 1000
            Number of training iterations.

        Returns
        -------
        list
            Loss per epoch.
        """
        Y = torch.tensor(Y.to_numpy(), dtype=torch.float64)
        C = torch.tensor(C.to_numpy(), dtype=torch.float64)

        optim = torch.optim.Adam(self.parameters(), lr=0.001)

        losses = []

        with self.progress_bar as p:
            for e in p.track(range(iter)):
                A_l = []
                B_l = []
                for layer in range(self.L):
                    g = torch.sqrt(self.g_l[layer] ** 2)
                    g = g * self.p / g.sum()
                    G_l = g.unsqueeze(1)
                    A_l.append(G_l * Y)
                    B_l.append(G_l * (self.X_l[layer] @ self.A_l[layer]))

                A_stack = torch.vstack(A_l)
                B_stack = torch.vstack(B_l)

                C_est = torch.linalg.lstsq(B_stack, A_stack).solution
                loss, _, _ = self.get_loss(C, C_est)

                loss.backward()
                optim.step()
                optim.zero_grad()

                losses.append(loss.item())

        # Norm weights and make non-negative
        # Option 1 => Might break relationship between layers
        for layer in range(self.L):
            gamma_l = self.g_l[layer] ** 2
            gamma_l = self.p * gamma_l / torch.sum(gamma_l)
            g_l = torch.sqrt(gamma_l)
            self.g_l[layer] = g_l

        return losses

    @torch.no_grad()
    def predict(
        self,
        Y: pd.DataFrame,
        norm: bool = False,
        library_sizes: pd.Series | None = None,
    ) -> dict[str, list[pd.DataFrame]]:
        """
        Predicts the cellular composition of a bulk.

        Parameters
        ----------
        Y : pd.DataFrame
            Bulk samples to be deconvoluted. (genes x samples)
        norm : bool = False
            Norm the results to one.
        library_sizes : pd.Series | None = None
            library_sizes at finest cell type layer. If provided, the results will be corrected for this.

        Returns
        -------
        dict[str,list[pd.DataFrame]]
            Dictionary where key "prediction" holds a list of compositions of each cell type layer.
            List elements are ordered in the same way as X_l and A_l. Element 0 corresponds to the finest
            coarsed cell type layer.
        """

        sample_names = list(Y.columns)
        Y = torch.tensor(Y.to_numpy(), dtype=torch.float64)

        A_l = []
        B_l = []
        for layer in range(self.L):
            G_l = self.g_l[layer].unsqueeze(1)
            A_l.append(G_l * Y)
            B_l.append(G_l * (self.X_l[layer] @ self.A_l[layer]))

        A_stack = torch.vstack(A_l)
        B_stack = torch.vstack(B_l)

        C_est = torch.linalg.lstsq(B_stack, A_stack).solution
        C_est[C_est < 0] = 0.0

        # Correct for library sizes if provided
        if library_sizes is not None:
            library_sizes = 1000000 / library_sizes
            C_est = torch.tensor(
                pd.DataFrame(
                    C_est.detach().cpu().numpy(),
                    index=self.celltype_layer_labels[0],
                    columns=sample_names,
                )
                .mul(library_sizes, axis=0)
                .to_numpy(),
                dtype=torch.float64,
            )

        if norm:
            C_est = C_est / C_est.sum(dim=0)

        prediction = []
        for layer in range(self.L):
            C_l = self.A_l[layer] @ C_est
            C_l = pd.DataFrame(
                C_l.detach().cpu().numpy(),
                index=self.celltype_layer_labels[layer],
                columns=sample_names,
            )
            prediction.append(C_l)

        return {"prediction": prediction}

__init__(X_l, A_l, lambdaNMSE=0.0)

Constructor of HIDE deconvolution model.

Parameters:

Name Type Description Default
X_l list[DataFrame]

List of reference profiles for each cell layer. First element should represent the finest coarsed resolution. (genes x celltypes_l)

 required
A_l list[DataFrame]

List of projection matrices, projecting from the finest coarsed resolution to a higher one. The first element should always be an identity matrix. (celltypes_l x celltypes_fine)

 required
lambdaNMSE float = 0.0

Weighting of optional NMSE contribution in general loss.

 0.0
Source code in src/hide_deconv/models/HIDE.py 
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
def __init__(
    self, X_l: list[pd.DataFrame], A_l: list[pd.DataFrame], lambdaNMSE: float = 0.0
):
    """
    Constructor of HIDE deconvolution model.

    Parameters
    ----------
    X_l : list[pd.DataFrame]
        List of reference profiles for each cell layer. First element should represent the finest coarsed resolution. (genes x celltypes_l)
    A_l : list[pd.DataFrame]
        List of projection matrices, projecting from the finest coarsed resolution to a higher one. The first element should always
        be an identity matrix. (celltypes_l x celltypes_fine)
    lambdaNMSE : float = 0.0
        Weighting of optional NMSE contribution in general loss.

    """
    super().__init__()

    self.L = len(A_l)
    self.lambdaNMSE = lambdaNMSE

    self.celltype_layer_labels = [X.columns for X in X_l]
    self.gene_labels = X_l[0].index

    self.p, _ = X_l[0].shape
    self.q_l = [len(X.columns) for X in X_l]

    self.A_l = [torch.tensor(A.to_numpy(), dtype=torch.float64) for A in A_l]
    self.X_l = [torch.tensor(X.to_numpy(), dtype=torch.float64) for X in X_l]

    self.g_l = nn.ParameterList(
        [
            nn.Parameter(
                torch.empty(self.p, dtype=torch.float64).uniform_(0.001, 0.1),
                requires_grad=True,
            )
            for _ in range(self.L)
        ]
    )

    self.progress_bar = Progress(
        TextColumn("[progress.percentage]{task.percentage:>3.0f}%"),
        BarColumn(),
        MofNCompleteColumn(),
        TextColumn("•"),
        TimeElapsedColumn(),
        TextColumn("•"),
        TimeRemainingColumn(),
    )

    self.epsilon = 1e-8  # small epsilon to prevent errors

train(Y, C, iter=1000)

Trains the HIDE model.

Parameters:

Name Type Description Default
Y DataFrame

Bulk training data. (genes x samples)

 required
C DataFrame

True cellular composition at finest coarsed cell resolution. (celltypes x samples)

 required
iter int = 1000

Number of training iterations.

 1000

Returns:

Type Description
list

Loss per epoch.
Source code in src/hide_deconv/models/HIDE.py 
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
def train(self, Y: pd.DataFrame, C: pd.DataFrame, iter: int = 1000) -> list:
    """
    Trains the HIDE model.

    Parameters
    ----------
    Y : pd.DataFrame
        Bulk training data. (genes x samples)
    C : pd.DataFrame
        True cellular composition at finest coarsed cell resolution. (celltypes x samples)
    iter : int = 1000
        Number of training iterations.

    Returns
    -------
    list
        Loss per epoch.
    """
    Y = torch.tensor(Y.to_numpy(), dtype=torch.float64)
    C = torch.tensor(C.to_numpy(), dtype=torch.float64)

    optim = torch.optim.Adam(self.parameters(), lr=0.001)

    losses = []

    with self.progress_bar as p:
        for e in p.track(range(iter)):
            A_l = []
            B_l = []
            for layer in range(self.L):
                g = torch.sqrt(self.g_l[layer] ** 2)
                g = g * self.p / g.sum()
                G_l = g.unsqueeze(1)
                A_l.append(G_l * Y)
                B_l.append(G_l * (self.X_l[layer] @ self.A_l[layer]))

            A_stack = torch.vstack(A_l)
            B_stack = torch.vstack(B_l)

            C_est = torch.linalg.lstsq(B_stack, A_stack).solution
            loss, _, _ = self.get_loss(C, C_est)

            loss.backward()
            optim.step()
            optim.zero_grad()

            losses.append(loss.item())

    # Norm weights and make non-negative
    # Option 1 => Might break relationship between layers
    for layer in range(self.L):
        gamma_l = self.g_l[layer] ** 2
        gamma_l = self.p * gamma_l / torch.sum(gamma_l)
        g_l = torch.sqrt(gamma_l)
        self.g_l[layer] = g_l

    return losses

predict(Y, norm=False, library_sizes=None)

Predicts the cellular composition of a bulk.

Parameters:

Name Type Description Default
Y DataFrame

Bulk samples to be deconvoluted. (genes x samples)

 required
norm bool = False

Norm the results to one.

 False
library_sizes pd.Series | None = None

library_sizes at finest cell type layer. If provided, the results will be corrected for this.

 None

Returns:

Type Description
dict[str, list[DataFrame]]

Dictionary where key "prediction" holds a list of compositions of each cell type layer. List elements are ordered in the same way as X_l and A_l. Element 0 corresponds to the finest coarsed cell type layer.
Source code in src/hide_deconv/models/HIDE.py 
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
@torch.no_grad()
def predict(
    self,
    Y: pd.DataFrame,
    norm: bool = False,
    library_sizes: pd.Series | None = None,
) -> dict[str, list[pd.DataFrame]]:
    """
    Predicts the cellular composition of a bulk.

    Parameters
    ----------
    Y : pd.DataFrame
        Bulk samples to be deconvoluted. (genes x samples)
    norm : bool = False
        Norm the results to one.
    library_sizes : pd.Series | None = None
        library_sizes at finest cell type layer. If provided, the results will be corrected for this.

    Returns
    -------
    dict[str,list[pd.DataFrame]]
        Dictionary where key "prediction" holds a list of compositions of each cell type layer.
        List elements are ordered in the same way as X_l and A_l. Element 0 corresponds to the finest
        coarsed cell type layer.
    """

    sample_names = list(Y.columns)
    Y = torch.tensor(Y.to_numpy(), dtype=torch.float64)

    A_l = []
    B_l = []
    for layer in range(self.L):
        G_l = self.g_l[layer].unsqueeze(1)
        A_l.append(G_l * Y)
        B_l.append(G_l * (self.X_l[layer] @ self.A_l[layer]))

    A_stack = torch.vstack(A_l)
    B_stack = torch.vstack(B_l)

    C_est = torch.linalg.lstsq(B_stack, A_stack).solution
    C_est[C_est < 0] = 0.0

    # Correct for library sizes if provided
    if library_sizes is not None:
        library_sizes = 1000000 / library_sizes
        C_est = torch.tensor(
            pd.DataFrame(
                C_est.detach().cpu().numpy(),
                index=self.celltype_layer_labels[0],
                columns=sample_names,
            )
            .mul(library_sizes, axis=0)
            .to_numpy(),
            dtype=torch.float64,
        )

    if norm:
        C_est = C_est / C_est.sum(dim=0)

    prediction = []
    for layer in range(self.L):
        C_l = self.A_l[layer] @ C_est
        C_l = pd.DataFrame(
            C_l.detach().cpu().numpy(),
            index=self.celltype_layer_labels[layer],
            columns=sample_names,
        )
        prediction.append(C_l)

    return {"prediction": prediction}