\n<\/p>\n
TxPert predicts the transcriptional response \\({\\bf{y}}\\in {\\mathcal{Y}}\\subset {{\\mathbb{R}}}^{n}\\)<\/span> given a set of perturbation tokens \\(P\\subset {\\mathcal{P}}:=\\{1,\\ldots ,N\\}\\)<\/span> and a basal state representation derived from a control expression profile \\({\\bf{x}}\\in {\\mathcal{X}}\\subset {{\\mathbb{R}}}^{n}\\)<\/span>, which has been aligned with the predicted cell with respect to cell line and batch effect. Here, \\(n\\in {\\mathbb{N}}\\)<\/span> denotes the number of experimentally measured genes and \\(N\\in {\\mathbb{N}}\\)<\/span> denotes the total number of observed perturbations in the data. These perturbation tokens (or identifiers) are used to select node representations from a biological gene (product)\u2013gene (product) interaction KG whose embeddings are integrated with the basal state to produce the perturbed expression profile.<\/p>\n To combine the information of the basal state and the perturbations, we first learn latent representations of both, that is, \\({\\bf{x}}\\mapsto {\\bf{s}}\\in {{\\mathbb{R}}}^{d}\\)<\/span> and \\(P\\mapsto \\{{{\\bf{z}}}_{p}\\in {{\\mathbb{R}}}^{d}:p\\in P\\}\\)<\/span> for a chosen latent dimension \\(d\\in {\\mathbb{N}}\\)<\/span>. Then, we combine the information through latent shift, where a learned decoder g<\/i>\u03d5<\/i><\/sub> predicts the perturbation effect from the given context, that is, \\(\\widehat{{\\bf{y}}}={g}_{\\phi }({\\bf{s}}+{\\sum }_{p\\in P}{{\\bf{z}}}_{p})\\)<\/span>, using the mean squared error (MSE):<\/p>\n $${\\mathcal{L}}({\\bf{y}},\\widehat{{\\bf{y}}})=\\frac{1}{n}| | {\\bf{y}}-\\widehat{{\\bf{y}}}| {| }_{2}^{2}$$<\/span><\/p>\n<\/div>\n This setup naturally integrates with single, double and more general multiperturbation settings of order m<\/i>:\u2009=\u2009\u2223<\/span>P<\/i>\u2223<\/span> through the additive and compositional design. More sophisticated combination functions may be used to learn transition functions \\({{\\bf{s}}}^{{\\prime} }={T}_{\\psi }({\\bf{s}},{\\bf{z}})\\)<\/span> that allow sequential latent cell state modeling of subsequently applied perturbations, for example, \\(\\widehat{{\\bf{y}}}={g}_{\\phi }({T}_{\\psi }(\\ldots {T}_{\\psi }({\\bf{s}},{{\\bf{z}}}_{{{\\bf{p}}}_{{\\bf{1}}}})\\ldots ,{{\\bf{z}}}_{{{\\bf{p}}}_{{\\bf{m}}}}))\\)<\/span>. To obtain s<\/b> and z<\/b>p<\/i><\/sub>,\u2009p<\/i>\u2009\u2208<\/span>\u2009P<\/i>, we use learned encoders, namely the basal state model and the perturbation model, which are discussed below.<\/p>\n The basal state encoder is designed to capture intrinsic cellular attributes, such as cell-cycle stage, cell type and other baseline phenotypic features, by mapping a control cell\u2019s gene expression profile \\({\\bf{x}}\\in {{\\mathbb{R}}}^{n}\\)<\/span> to a compact, low-dimensional embedding, \\({\\bf{s}}={f}_{{\\mathtt{basal}}}({\\bf{x}})\\in {{\\mathbb{R}}}^{d}\\)<\/span>. The basal state model was subject to hyperparameter tuning, with an MLP performing best for unseen and double tasks (Supplementary Note 1 and Supplementary Fig. 7 for cross-cell type). The MLP learns a direct deterministic mapping from high-dimensional input gene expression data to a fixed-size embedding space. The MLP architecture offers a simple and computationally efficient framework for representation learning, while still retaining the capacity to model complex, nonlinear dependencies inherent in gene expression data.<\/p>\n An important aspect of modeling the basal state is the alignment of the control with the predicted perturbed cell. Note that experimental protocols can vary widely between different data sources. Therefore, we randomly sample this control according to the same cell line and dataset or experimental protocol. Furthermore, we explore basal state matching, where the control cell is selected to closely match the batch metadata of the perturbed sample. As this matching is not unique, we randomly sample one such appropriate control. Lastly, we experiment with basal state averaging, where instead of a single control, we compute the average expression profile across all matching controls for a given cell line and\/or batch. This produces a more stable estimate of the basal state. Both strategies consistently improved model performance in our experiments.<\/p>\n Beyond MLP encoders on raw gene expression profiles, we explored multiple transcriptomics foundation model embeddings to obtain a latent representation of the basal state. Specifically, we experiment with scGPT, and scVI pretrained on the Joung dataset39<\/sup>. We also include a variant with no basal state encoder, where the (latent) basal state space is represented directly by the raw expression profile (that is, s<\/b>: = x<\/b>). In this configuration, the perturbation decoder learns a \u0394<\/i> vector from the perturbation embedding, which is added to the control profile: \\({\\bf{x}}\\in {{\\mathbb{R}}}^{n}\\)<\/span>, that is, \\(\\widehat{{\\bf{y}}}={\\bf{x}}+{g}_{\\phi }({\\sum }_{p\\in P}{{\\bf{z}}}_{p})\\)<\/span>. This resembles the formulation of the general baseline, with a trained model predicting the \u0394<\/i> instead of a hand-crafted heuristic. Unsurprisingly, this variant shines most in settings with limited data availability for learning robust basal state representations, for example, perturbation effect prediction in cell lines without seen perturbations.<\/p>\n We rely on GNNs that use biological KGs capturing gene (product)\u2013gene (product) interactions to learn informative embeddings for gene perturbations. The GNN learns a matrix of node embeddings associated with the perturbation tokens \\(\\{1,\\ldots ,N\\}\\mapsto {\\bf{Z}}\\in {{\\mathbb{R}}}^{N\\times d}\\)<\/span>, where \\(N\\in {\\mathbb{N}}\\)<\/span> is the number of perturbations relevant to the investigated task. Each row (or node representation) of this matrix represents the latent encoding of a specific perturbation, that is, \\({{\\bf{z}}}_{p}\\in {{\\mathbb{R}}}^{d}\\)<\/span> required for the latent shift is the p<\/i>th row of Z<\/b>. More specifically, we first associate each perturbation p<\/i>\u2009\u2208<\/span>\u2009{1, \u2026, N<\/i>} with a randomly initialized input node embedding \\({{\\bf{h}}}_{p}^{0}\\in {{\\mathbb{R}}}^{{d}_{0}},{d}_{0}\\in {\\mathbb{N}}\\)<\/span> that are consolidated in the input node feature matrix \\({{\\bf{H}}}^{0}\\in {{\\mathbb{R}}}^{N\\times {d}_{0}}\\)<\/span>. During training, those node features are (1) treated as model parameters that are learned using backpropagation and (2) subsequently refined through the message passing of an L<\/i>-layer GNN, that is, \\({{\\bf{H}}}^{\\ell }={{\\mathtt{LAYER}}}_{{\\theta }_{\\ell }}({{\\bf{H}}}_{\\ell -1}),0\\le \\ell \\le L\\)<\/span> and Z<\/b>:\u2009=\u2009H<\/b>L<\/i><\/sup>. This allows the model to characterize perturbation effects on the basis of known relationships from KGs such as GO19<\/sup>, STRING18<\/sup> or proprietary data sources.<\/p>\n Real-world KGs present inherent imperfections; they often contain noisy or incorrect edges (false positives), suffer from missing connections (false negatives) and may originate from diverse sources offering multiple, sometimes conflicting, perspectives. The effective use of such complex data necessitates GNN architectures specifically chosen to address these challenges.<\/p>\n To this end, we select and adapt two fundamental GNN approaches. First, to handle noisy edges, we use attention-based models such as the GAT40,41<\/sup>. The ability of GAT to dynamically (re)weight neighbor importance provides much needed robustness by effectively downweighting less credible connections, which is a major difference relative to non-attention-based methods such as the simple graph convolution42<\/sup> used in GEARS. Second, to address graph incompleteness and capture long-range dependencies, we use GTs, specifically Exphormer16,17<\/sup>. Its capacity for attention beyond immediate graph neighbors allows it to potentially model implicit relationships and bridge structural gaps.<\/p>\n Furthermore, it proved crucial for the presented tasks to use multiple KGs that offer complementary and reinforcing perspectives on the task-related biology. For this, we explore architectures designed for synergistic learning from diverse sources. This includes extending GAT to GAT-Hybrid (allowing for node-level attention weighting of information from different KGs), introducing our provenance-aware Exphormer-MG variant and developing GAT-MLG, a multilayer extension of GAT that uses a supra-adjacency representation to effectively integrate information across multiple biological networks.<\/p>\n In Supplementary Note 3, we provide rigorous details on the relevant graph representation learning used for encoding perturbations, the proposed models and the techniques applied to take advantage of complementary information from multiple KGs. In our experimental setup, the learned embeddings for the perturbed gene(s) from the GNN were extracted and combined with a basal state representation to predict the resulting gene expression profile. This comparative analysis allowed us to investigate how different GNN strategies (attention, flexible connectivity and multigraph fusion) perform when learning from the complexities of biological KGs.<\/p>\n TxPert training algorithm<\/b><\/p>\n Require<\/b>: Pert. cells \\({\\mathcal{Y}}\\)<\/span>, control cells \\({\\mathcal{X}}\\)<\/span>, biological prior graph G<\/i>\u2009=\u2009(V<\/i>, E<\/i>) with perturbations \\({\\mathcal{P}}\\subset V\\)<\/span><\/p>\n Ensure<\/b>: Minimize MSE loss between predicted and true perturbed cell measurements<\/p>\n 1: Initialize input perturbation embeddings \\({\\{{{\\bf{h}}}_{v}^{0}\\}}_{v\\in V}\\)<\/span> randomly<\/p>\n 2: for<\/b> each training step do<\/b><\/p>\n 3: Sample a batch of perturbed cell profiles: \\({\\{{{\\bf{y}}}_{i}\\}}_{i=1}^{B}\\subset {\\mathcal{Y}}\\)<\/span><\/p>\n 4: Sample corresponding control cells from the same experimental batches: \\({\\{{{\\bf{x}}}_{i}\\}}_{i=1}^{B}\\subset {\\mathcal{X}}\\)<\/span><\/p>\n 5: Enrich perturbation embeddings using graph prior: \\({\\{{{\\bf{z}}}_{v}\\}}_{v\\in V}\\leftarrow {\\mathtt{GNN}}(G,{\\{{{\\bf{h}}}_{v}^{0}\\}}_{v\\in V})\\)<\/span><\/p>\n 6:\u2003for<\/b> each sample in batch do<\/b><\/p>\n 7: Encode control cells into basal latent space: b<\/b>i<\/i><\/sub>\u2009\u2190\u2009MLP<\/span>basal<\/sub>(x<\/b>i<\/i><\/sub>)<\/p>\n 8: Retrieve perturbations \\({P}_{i}\\subset {\\mathcal{P}}\\)<\/span> associated with target y<\/b>i<\/i><\/sub><\/p>\n 9:\u2003Combine control and perturbation embeddings: \\({\\widehat{{\\bf{z}}}}_{i}\\leftarrow {\\mathtt{COM}}({{\\bf{b}}}_{i},\\{{{\\bf{z}}}_{p}:p\\in {P}_{i}\\})\\)<\/span><\/p>\n 10: Decode to predicted perturbed profile: \\({\\widehat{{\\bf{y}}}}_{i}\\leftarrow {{\\mathtt{MLP}}}_{{\\rm{dec}}}({\\widehat{{\\bf{z}}}}_{i})\\)<\/span><\/p>\n 11: Compute loss for each sample: \\({{\\mathcal{L}}}_{i}\\leftarrow {\\mathtt{MSE}}({\\widehat{{\\bf{y}}}}_{i},{{\\bf{y}}}_{i})\\)<\/span><\/p>\n 12: end for<\/b><\/p>\n 13: Compute total loss over batch: \\({\\mathcal{L}}\\leftarrow {\\sum }_{i=1}^{B}{{\\mathcal{L}}}_{i}\\)<\/span><\/p>\n 14: Backpropagate and update model parameters<\/p>\n 15: end for<\/b><\/p>\n The data were split into training, validation and test sets through grouping by perturbation such that distinct sets of unseen perturbations were reserved both the validation and test sets with target ratios of 0.5625, 0.1875 and 0.25 for the training, validation and test sets, respectively. Moreover, for the cross-cell-type task, the test set was a reserved cell type with only control cells included during training with a breakdown into seen and unseen perturbations therein. As an exception, for the doubles task on the Norman dataset, predefined splits were loaded from the GEARS setup.<\/p>\n Optimal hyperparameters for each model were selected based on the validation Pearson \u0394<\/i> metric. Only metrics on the test set are reported.<\/p>\n All metrics are reported as weighted averages, that is, the mean of the mean across cells subjected to each unique perturbation, unless otherwise specified.<\/p>\n Where not otherwise specified, expression values \\({\\bf{x}}\\in {\\mathcal{X}}\\cup {\\mathcal{Y}}\\)<\/span>, are represented as log1p<\/i>-transformed and library size-normalized counts (with target library size of 4,000); that is, for a raw count \\({{\\bf{x}}}_{{\\rm{raw}}}\\in {{\\mathbb{R}}}^{n}\\)<\/span>, we define<\/p>\n $${\\bf{x}}:=\\log \\left(1+4000\\cdot \\frac{{{\\bf{x}}}_{\\mathrm{raw}}}{\\parallel {{\\bf{x}}}_{\\mathrm{raw}}{\\parallel }_{1}}\\right).$$<\/span><\/p>\n<\/div>\n The other representation used is a \u0394<\/i> representation, which is centered on batch-matched controls. Specifically, for each perturbed cell expression \\({{\\bf{y}}}_{i}\\in {\\mathcal{Y}}\\)<\/span> with cell line c<\/i> and batch b<\/i>, the expression is transformed to<\/p>\n $${\\delta }_{i}:={{\\bf{y}}}_{i}-{\\overline{{\\bf{x}}}}_{c,b},$$<\/span><\/p>\n<\/div>\n where \\({\\overline{x}}_{c,b}\\)<\/span> represents the mean expression of control cells \\({\\bf{x}}\\in {\\mathcal{X}}\\)<\/span> with batch b<\/i> and cell line c<\/i>.<\/p>\n Slightly modified from the metric \u2018Pearson correlation (\u0394<\/i> expression)\u2019 from the GEARS manuscript, Pearson \u0394<\/i> calculates the correlation between predicted and observed log fold change versus batch-matched control mean,<\/p>\n $${\\text{Pearson}}\\Delta (p):=\\mathrm{Pearson}({\\widehat{\\delta }}_{p},{\\delta }_{p}),$$<\/span><\/p>\n<\/div>\n where \\({\\widehat{\\delta }}_{p}\\)<\/span> and \u03b4<\/i>p<\/i><\/sub> are the batch-matched control centering of the prediction and ground truth, respectively, averaged across replicates of certain perturbation \\(p\\in {\\mathcal{P}}\\)<\/span>. For simplicity, we define this and following metrics for single perturbations \\(p\\in {\\mathcal{P}}\\)<\/span> but note that analogous formulations are appropriate for multiple perturbations \\(P\\subset {\\mathcal{P}}\\)<\/span>. The results across all predicted perturbation effects are then averaged to obtain an overall performance estimate.<\/p>\n Note that, for the GEARS model only, we report the exact \u2018Pearson correlation (\u0394<\/i> expression)\u2019 from the GEARS code base instead. We confirmed that any differences between \u2018Pearson correlation (\u0394<\/i> expression)\u2019 and our \u2018Pearson \u0394<\/i>\u2019 were much smaller in practice than the differences between models.<\/p>\n We use two variants of the retrieval rank metric that score a prediction\u2019s similarity to the ground truth not overall but relative to other perturbations. These metrics are the same as rank average from PerturBench10<\/sup>, except that they focus on similarity with a perfect score of 1, a random score of 0.5 and perfect anticorrelated prediction score of 0:<\/p>\n $$\\begin{array}{ll}\\mathrm{Retrieval}: & =\\displaystyle\\frac{1}{N}\\sum _{p\\in {\\mathscr{ \\mathcal P }}}\\,\\mathrm{rank}\\,({\\hat{\\delta }}_{p}),\\\\ \\mathrm{rank}\\,({\\hat{\\delta }}_{p}): & =\\displaystyle\\frac{1}{N-1}\\mathop{\\sum}\\limits_{{q\\in {\\mathscr{ \\mathcal P }}}\\atop {q\\ne p}}{{\\bf{1}}}_{\\{\\mathrm{Pearson}({\\hat{\\delta }}_{p},{\\delta }_{p})\\ge \\mathrm{Pearson}({\\hat{\\delta }}_{p},{\\delta }_{q})\\}}.\\end{array}$$<\/span><\/p>\n<\/div>\n For \u2018normalized\u2019 retrieval, the perturbation count \\(N:=| {\\mathcal{P}}|\\)<\/span> matches the original experiment, whereas, for \u2018fast retrieval\u2019, for computational efficiency, a seeded random reference set of only 100 perturbations is used, with the addition of the query perturbant p<\/i> when not in the reference set (thus, N<\/i>\u2009\u2208<\/span>\u2009{100, 101}. Similar to Person \u0394<\/i>, we report the averaged performance across all perturbations.<\/p>\n To establish a performance floor, we implement a nonlearned general baseline model that predicts expression profiles using mean values observed in the training data. This baseline uses an additive approach that combines the following:<\/p>\n The mean test cell type control expression profile<\/p>\n<\/li>\n Either the perturbation-specific mean changes (for seen perturbations) or the global perturbation mean (for unseen perturbations)<\/p>\n<\/li>\n When multiple cell lines are present in the training set, we either use a weighted average according to the number of samples per cell line or the perturbation-specific mean changes from the most similar cell line (nearest-cell-line baseline). Here, similarity is determined on the basis of mean correlation of shared perturbation \u0394<\/i> values between the test and candidate cell line.<\/p>\n<\/li>\n<\/ul>\n Consider a multiset of training samples \\(\\,\\mathrm{Train}\\,\\subset {{\\mathcal{P}}}_{\\mathrm{train}}\\times {{\\mathcal{C}}}_{\\mathrm{train}}\\times {{\\mathcal{B}}}_{\\mathrm{train}}\\)<\/span> consisting of combinations of perturbation(s), cell line(s) and batch effect(s) with a multiset test defined analogously. Consider a perturbation p<\/i> such that (p<\/i>,\u2009c<\/i>p<\/i><\/sub>,\u2009b<\/i>p<\/i><\/sub>)\u2009\u2208<\/span>\u2009test with cell line c<\/i>p<\/i><\/sub> and batch effect b<\/i>p<\/i><\/sub>. Implicitly, a (c<\/i>p<\/i><\/sub>,\u2009b<\/i>p<\/i><\/sub>) map is associated with a set of control cell profiles in that context.<\/p>\n If there exists (p<\/i>,\u2009c<\/i>,\u2009b<\/i>)\u2009\u2208<\/span>\u2009train, we have<\/p>\n $$\\begin{array}{ll}{\\widehat{{\\bf{y}}}}_{(p,{c}_{p},{b}_{p})} & =\\\\ & {\\bar{{\\bf{x}}}}_{({c}_{p},{b}_{p})}+\\displaystyle\\frac{1}{| \\{(q,c,b)\\in \\,\\mathrm{Train}:q=p\\}| }\\mathop{\\sum }\\limits_{{(q,c,b)\\in \\,\\mathrm{Train}}\\atop{q=p}}{{\\bf{y}}}_{(p,c,b)}-{\\bar{{\\bf{x}}}}_{(c,b)}.\\end{array}$$<\/span><\/p>\n<\/div>\n Otherwise, we use the global \u0394<\/i> across perturbations observed in the training set, that is,<\/p>\n $${\\widehat{{\\bf{y}}}}_{(p,{c}_{p},{b}_{p})}={\\bar{{\\bf{x}}}}_{({c}_{p},{b}_{p})}+\\frac{1}{| \\,\\mathrm{Train}| }\\mathop{\\sum }\\limits_{(q,c,b)\\in \\mathrm{Train}}{{\\bf{y}}}_{(p,c,b)}-{\\bar{{\\bf{x}}}}_{(c,b)}.$$<\/span><\/p>\n<\/div>\n For multiple perturbations, this baseline is implemented to initially attempt to use samples where the exact perturbation configuration is present. Otherwise, the perturbation is split into its components and each component is sequentially added to the test control mean according to the above method, adding a local \u0394<\/i> estimate if available and resorting to a global \u0394<\/i> otherwise.<\/p>\n As Perturb-seq is a destructive assay, we cannot observe the same cell in both perturbed and unperturbed states. This necessitates focusing on distribution means rather than individual cell accuracies. To approximate experimental reproducibility, we first use a split-half validation approach:<\/p>\n For each combination of perturbation(s), cell line context and batch, we apply three operations:<\/p>\n Divide test cells into two roughly equal halves<\/p>\n<\/li>\n Calculate mean expression profiles for each half<\/p>\n<\/li>\n Measure the agreement between these means using various metrics<\/p>\n<\/li>\n<\/ol>\n To account for the randomness in choosing the half-split, we repeat the experiment across multiple seeded runs and report average performance. This provides a performance benchmark analogous to human-level reproducibility, which is called accuracy in other machine learning domains.<\/p>\n Consider the set of expression profiles \\({\\mathcal{S}}\\subset {\\mathcal{Y}}\\)<\/span> for a fixed perturbation cell line context and batch (p<\/i>,\u2009c<\/i>,\u2009b<\/i>) in the test set:<\/p>\n $$\\begin{array}{rcl}{{\\mathcal{S}}}^{{\\prime} } & \\subseteq & {\\mathcal{S}}:\\,| {{\\mathcal{S}}}^{{\\prime} }| \\approx | {\\mathcal{S}}| \/2\\\\ {\\bar{{\\mathcal{S}}}}_{1} & = & \\frac{1}{| {{\\mathcal{S}}}^{{\\prime} }| }\\mathop{\\sum }\\limits_{{\\bf{y}}\\in {{\\mathcal{S}}}^{{\\prime} }}{\\bf{y}}\\\\ {\\bar{{\\mathcal{S}}}}_{2} & = & \\frac{1}{| {\\mathcal{S}}\\backslash {{\\mathcal{S}}}^{{\\prime} }| }\\mathop{\\sum }\\limits_{{\\bf{y}}\\in {\\mathcal{S}}\\backslash {{\\mathcal{S}}}^{{\\prime} }}{\\bf{y}}.\\end{array}$$<\/span><\/p>\n<\/div>\n We then report<\/p>\n $$\\mathrm{Reproduce}(p,c,b)=\\mathrm{Metric}({\\bar{{\\mathcal{S}}}}_{1},{\\bar{{\\mathcal{S}}}}_{2}),$$<\/span><\/p>\n<\/div>\n where metric<\/span> represents any of our evaluation metrics, for example, Pearson \u0394<\/i>, Retrieval<\/span> or MSE<\/span>. Theoretically, the split-half experimental reproducibility is not expected to establish an upper bound for performance of all models at test time because it operates on a different test set (only using half for prediction and testing, respectively). However, it empirically proves to be useful as a competitive (but still theoretically reachable) mark to beat.<\/p>\n As split-half reproducibility is likely an underestimate because of a reduced (halved) number of replicates, we also introduce a sample-based approach that gives an estimate for the reproducibility of the full-size dataset. Using the original count matrix, we calculate the per-batch probability distribution over genes (multinomial distribution maximum-likelihood estimator) for each perturbation. We then sample these distributions to generate two datasets that have the same number of observations (that is, cells) as the original dataset, but with stochastically resampled counts. These are then log1p<\/i>-normalized and subset to the HVGs in the same way as the original dataset, before calculating the experiment reproducibility as described above. A comparison of split-half and sampled reproducibility is shown in Supplementary Table 4.<\/p>\n We demonstrate the efficacy of our approach across a range of datasets, including CRISPRi (gene knockdown) of ~2,000 essential genes in K562 and RPE1 cell lines from a previous study13<\/sup> (also used in GEARS4<\/sup>) and similarly designed CRISPRi experiments in Jurkat and HEPG2 cell lines from another previous study24<\/sup>. Furthermore, we implement the Norman15<\/sup> dataset with 94 unique single and 110 unique double CRISPRa (gene overexpression) perturbations respectively in the K562 cell line.<\/p>\n The graphs used as inductive bias in this work can be classified into two main categories: (1) curated publicly available biological knowledge and (2) large-scale perturbation screens.<\/p>\n The curated graphs from category 1 include the GO graph, first used by GEARS, which is constructed by assigning edges between nodes that have a high Jaccard Index in their GO terms19<\/sup>, the STRING graph18<\/sup> and Reactome43<\/sup>.<\/p>\n Category 2 graphs are generated from large-scale perturbation screens including DepMap44<\/sup> and Perturb-seq45<\/sup>. These are extensive datasets linking genetic perturbation to either morphological or transcriptomic outcomes, which can offer particularly crucial insights into cellular responses to stimuli. To translate these experimental screens into graphs, we use derived embeddings to represent the genes and cell lines in a high-dimensional space, allowing for the analysis of relationships and identification of dependencies.<\/p>\n To curate these graphs, we first compute the pairwise similarity score between all combinations of genes. This means that, for each pair of genes (g<\/i>i<\/i><\/sub>,\u2009g<\/i>j<\/i><\/sub>), we compute the cosine similarity between their (aggregated) embeddings \\({{\\bf{x}}}_{{g}_{i}}\\)<\/span> and \\({{\\bf{x}}}_{{g}_{\\!j}}\\)<\/span>. Cosine similarity is computed as follows:<\/p>\n $$\\begin{array}{rcl}\\,{{\\rm{cosine}}\\; {\\rm {similarity}}}\\,({{\\bf{x}}}_{{g}_{i}},{{\\bf{x}}}_{{g}_{\\!j}}) & = & \\frac{{{\\bf{x}}}_{{g}_{i}}\\cdot {{\\bf{x}}}_{{g}_{\\!j}}}{\\parallel {{\\bf{x}}}_{{g}_{i}}\\parallel \\parallel {{\\bf{x}}}_{{g}_{\\!j}}\\parallel }\\\\ & = & \\frac{{\\sum }_{k=1}^{n}{{\\bf{x}}}_{{g}_{i}k}{{\\bf{x}}}_{{g}_{\\!j}k}}{\\sqrt{{\\sum }_{k=1}^{n}{{\\bf{x}}}_{{g}_{i}k}^{2}}\\cdot \\sqrt{{\\sum }_{k=1}^{n}{{\\bf{x}}}_{{g}_{\\!j}k}^{2}}}\\end{array}$$<\/span><\/p>\n<\/div>\n where<\/p>\n \\({{\\bf{x}}}_{{g}_{i}}\\cdot {{\\bf{x}}}_{{g}_{\\!j}}\\)<\/span> represents the dot product of the vectors<\/p>\n<\/li>\n \\(\\parallel {{\\bf{x}}}_{{g}_{i}}\\parallel\\)<\/span> and \\(\\parallel {{\\bf{x}}}_{{g}_{i}}\\parallel\\)<\/span> represent the Euclidean norms (magnitudes) of vectors \\({{\\bf{x}}}_{{g}_{i}}\\)<\/span> and \\({{\\bf{x}}}_{{g}_{\\!j}}\\)<\/span>, respectively<\/p>\n<\/li>\n \\({{\\bf{x}}}_{{g}_{i}k}\\)<\/span> and \\({{\\bf{x}}}_{{g}_{\\!j}k}\\)<\/span> are the individual components of vectors \\({{\\bf{x}}}_{{g}_{i}k}\\)<\/span> and \\({{\\bf{x}}}_{{g}_{\\!j}k}\\)<\/span>.<\/p>\n<\/li>\n<\/ul>\n These cosine similarities are converted to their absolute values because the difference between highly cosine negative and highly cosine positive does not translate literally to the signed weight of the edge in the graph.<\/p>\n We additionally use proprietary data from internal genome-wide perturbation screens, where we measure the similarity of perturbation effect using both microscopy imaging and transcriptomics in various cell types.<\/p>\n Filtering configurations were optimized empirically. We found that the most performant configuration involved selecting for the top 1% of edges by (absolute) weight for screen-based graphs. For all other graph types, we (additionally) filtered for no more than 20 incoming nodes by target. Edge direction was assigned arbitrarily for undirected edges.<\/p>\n Additional methods related to specific analyses are described below.<\/p>\n The Pharos initiative consolidates a variety of statistics relating to how researched and well known specific genes are26<\/sup>. Starting from this, we ranked knowledge levels as the mean of the rank of the Pharos Pubmed score and the rank of the Pharos negative log novelty score to create a single Pharos knowledge rank. We used this rank to break down and compare to the performance of models and understand potential bias. The \u2018knowledge levels\u2019 0, 1, 2 and 3 correspond to the following bins of the Pharos knowledge rank:<\/p>\n knowledge level 0 (least characterized): 0\u20130.2<\/p>\n<\/li>\n knowledge level 1: 0.2\u20130.4<\/p>\n<\/li>\n knowledge level 2: 0.4\u20130.6<\/p>\n<\/li>\n knowledge level 3: 0.6\u20130.8<\/p>\n<\/li>\n knowledge level 3 (most characterized): 0.8\u20131.<\/p>\n<\/li>\n<\/ul>\n In investigating the correlations between controls and mean baselines, we compared \u2018within\u2019-context correlation to \u2018across\u2019-context correlation. Generally, before calculating either, all examples were first split into two mutually exclusive halves, A and B, where within-context correlation is a comparison of A versus B, in each context, while across-context correlation is a comparison of an arbitrary half of one context to another context. The only exception is across-batch controls in Fig. 1a, for which, to make a conservative estimate of across-batch variance, full batches were aggregated without splitting. For batch comparison, individual control cells were split and aggregated; for the mean baselines, the \u03b4<\/i> of perturbant replicate cells was preaggregated and then split (such that the halves had nonoverlapping perturbations).<\/p>\n To achieve a descriptive biological summary of the actual gene expression changes in the mean baseline, we first calculated a meta mean baseline (mean of all cell types and datasets, using the intersect of provided expressed genes) and defined upregulated and downregulated genes as having a remaining \u03b4<\/i>\u2009>\u20090.05 or \u03b4<\/i>\u2009<\u2009\u22120.05, respectively. We then ran functional enrichment testing separately on each on these sets (versus the background of all genes in the dataset intersect) using the GOATOOLS package46<\/sup>.<\/p>\n To avoid confusion and distraction caused by reporting many similarly performing or perhaps slightly contradicting metrics, we first performed an empirical \u2018test of the test metrics\u2019. We adapted the evaluation method pioneered by a previous study12<\/sup> to work for our data and task. In short, this method uses cross-context retrieval of a perturbation as a way to judge whether a representation and metric together allow the retention and comparison of details necessary to distinguish perturbations. In our case, we modified the method to work on nonaggregated single-cell expression profiles (as this is the input and output of our model) and ran retrieval across the essential perturbation set on core cell types (K562, RPE1, HEPG2 and Jurkat). For each retrieval calculation, instead of aggregating, we first randomly sampled one cell of each perturbation. We ran three replicates on each of the cell\u2009\u00d7\u2009cell pairings for n<\/i>\u2009=\u20093\u2009\u00d7\u20096\u2009=\u200918 total estimates per perturbation. We report the 0.9 quantile, to focus on active perturbants; however, similar patterns were observed at other thresholds.<\/p>\n Note that the choice to focus on single cells excluded use of the representation selected by a previous study12<\/sup>, namely the signed Basal state model<\/h4>\n
Basal information matching and aggregation<\/h3>\n
Encoders<\/h3>\n
Perturbation model<\/h4>\n
Training and evaluation framework<\/h3>\n
Algorithm 1<\/h3>\n
Data splits<\/h4>\n
Metric definitions<\/h4>\n
Expression value representations and \u0394<\/i>
\n <\/h3>\nPearson \u0394<\/i>
\n <\/h3>\nRetrieval<\/h3>\n
Nonlearned general baseline<\/h4>\n
\n
Experimental reproducibility estimation: split-half validation and sample-based extension<\/h4>\n
\n
Data<\/h3>\n
Perturb-seq data sources<\/h4>\n
Graphs: sourcing and processing<\/h4>\n
\n
Data understanding<\/h4>\n
Pharos knowledge rank<\/h3>\n
\n
Within versus across<\/h3>\n
Functional enrichment<\/h3>\n
Metric selection through retrieval<\/h3>\n