{'fraction_non_alphanumeric': 0.03824890601645928, 'fraction_numerical': 0.013902144747978606, 'mean_word_length': 5.097727023168991, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 3, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 12, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hypergraph Factorisation for Multi-Tissue Gene Expression Imputation

Integrating gene expression across scales and tissues is crucial for understanding the biological mechanisms that drive disease and characterise homeostasis. However, traditional multi-tissue integration methods cannot handle uncollected tissues or rely on genotype information, which is subject to privacy concerns and often unavailable. To address these challenges, we present HYFA (Hypergraph Factorisation), a novel method for joint imputation of multi-tissue and cell-type gene expression. HYFA imputes tissue-specific gene expression via a specialised graph neural network operating on a hypergraph of individuals, metagenes, and tissues. HYFA is genotypeagnostic, supports a variable number of collected tissues per individual, and imposes strong inductive biases to leverage the shared regulatory architecture of tissues. In performance comparison on data from the Genotype Tissue Expression project, HYFA achieves superior performance over existing transcriptome imputation methods, especially when multiple reference tissues are available. Through transfer learning on a paired single-nucleus RNA-seq (snRNA-seq) dataset, we further show that HYFA can accurately resolve cell-type signatures from bulk gene expression, highlighting the method’s ability to leverage gene expression programs underlying cell-type identity, even in tissues that were never observed in the training set. Using Gene Set Enrichment Analysis, we find that the metagenes learned by HYFA capture information about known biological pathways. Notably, the HYFA-imputed dataset can be used to identify regulatory genetic variations (eQTLs), with substantial gains over the original incomplete dataset. Our framework can accelerate effective and scalable integration of tissue and cell-type gene expression biorepositories.


Introduction
Sequencing technologies have enabled profiling the transcriptome at tissue and single-cell resolutions, with great potential to unveil intra-and multi-tissue molecular phenomena such as cell signalling and disease mechanisms. Due to the invasiveness of the sampling process, gene expression is usually measured independently in easy-to-acquire tissues, leading to an incomplete picture of an individual's physiological state and necessitating effective multi-tissue integration methodologies.
A question of fundamental biological significance is to what extent the transcriptomes of difficult-to-acquire tissues and cell types can be inferred from those of accessible ones [Basu et al., 2021, Consortium, 2020. Due to their ease of collection, accessible tissues such as whole blood could have great utility for diagnosis and monitoring of pathophysiological conditions through metabolites, signalling molecules and other biomarkers, including potential transcriptome-level associations [Yang et al., 2020]. Moreover, all human somatic cells share the same genetic information, which may regulate expression in a context-dependent and temporal manner, partially explaining tissueand cell-type-specific gene expression variation. Computational models that exploit these patterns could therefore be used to impute the transcriptomes of uncollected cell types and tissues, with potential to elucidate the biological mechanisms regulating a diverse range of developmental and physiological processes.
Multi-tissue imputation is a central problem in transcriptomics with broad implications for fundamental biological research and translational science. The methodological problem can powerfully influence downstream applications, including performing differential expression analysis, identifying regulatory (genetic) mechanisms, determining coexpression networks, and enabling drug target discovery. In practice, in experimental follow-up or clinical application, the task includes the special case of determining a good proxy or easily-assayed system for causal tissues and cell types. Multi-tissue integration methods can also be applied to harmonise large collections of RNA-seq datasets from diverse institutions, consortia, and studies [Xu et al., 2021] -each potentially affected by technical artifacts -and to characterise gene expression co-regulation across tissues. Reconstruction of unmeasured gene expression across a broad collection of tissues and cell types from available reference transcriptome panels may expand our understanding of the molecular origins of complex traits and of their context specificity.
Several methods have traditionally been employed to impute uncollected gene expression. Leveraging a surrogate tissue has been widely used in studies of biomarker discovery, diagnostics, and expression Quantitative Trait Loci (eQTLs), and in the development of model systems [Hoon et al., 2000, Cai et al., 2010, Istas et al., 2017, Gamazon et al., 2018, Kim et al., 2020. Nonetheless, gene expression is known to be tissue and cell-type specific, limiting the utility of a proxy tissue. Other related studies impute tissue-specific gene expression from genetic information [Zhou et al., 2020].  propose a mixed-effects model to infer uncollected data in multiple tissues from eQTLs. Sul et al. [2013] introduce a model termed Meta-Tissue that aggregates information from multiple tissues to increase statistical power of eQTL detection. However, these approaches do not model the complex, non-linear relationships between measured and unmeasured gene expression traits among tissues and cell types, and individual-level genetic information (such as at eQTLs) is subject to privacy concerns and often unavailable.
Computationally, multi-tissue transcriptome imputation is challenging because the data dimensionality scales rapidly with the number of genes and tissues, often leading to overparameterised models. TEEBoT [Basu et al., 2021] addresses this issue by employing principal component analysis (PCA) -a non-parametric dimensionality reduction methodto project the data into a low-dimensional manifold, followed by linear regression to predict target gene expression from the principal components. However, this technique does not account for non-linear effects and can handle only a single reference tissue, i.e. whole blood. Approaches such as standard multilayer perceptrons can exploit non-linear patterns, but are massively overparameterised and computationally infeasible.
To address these challenges, we present HYFA (Hypergraph Factorisation), a parameter-efficient graph representation learning approach for joint multi-tissue and cell-type gene expression imputation. HYFA represents multi-tissue gene expression in a hypergraph of individuals, metagenes, and tissues, and learns factorised representations via a specialised message passing neural network operating on the hypergraph. In contrast to existing methods, HYFA supports a variable number of reference tissues, increasing the statistical power over single-tissue approaches, and incorporates inductive biases to exploit the shared regulatory architecture of tissues. In performance comparison, HYFA attains improved performance over TEEBoT and standard imputation methods across a broad range of tissues from the Genotype-Tissue Expression (GTEx) project (v8) [Consortium, 2020]. Through transfer learning on a paired single-nucleus RNA-seq dataset (GTEx-v9) [Eraslan et al., 2021], we further demonstrate the ability of HYFA to resolve cell-type signatures from bulk gene expression. Thus, HYFA may provide a unifying transcriptomic methodology for multi-tissue imputation and cell-type deconvolution. In post-imputation analysis, application of eQTL mapping on the fully-imputed GTEx data yields a substantial increase in number of detected eQTLs. HYFA is publicly accessible at https://github.com/rvinas/HYFA.

Results
Hypergraph factorisation: a multi-tissue transcriptomic imputation approach We developed HYFA, a framework for inferring the transcriptomes of unmeasured tissues and cell-types from bulk expression collected in a variable number of reference tissues (Figure 1, Methods). HYFA receives as input gene expression measurements collected from a set of reference tissues, as well as demographic information, and outputs gene expression values in a tissue of interest (e.g. uncollected). The first step of the HYFA workflow is to project the input gene expression into low-dimensional metagene representations [Brunet et al., 2004, Raychaudhuri et al., 1999 for every collected tissue. Each metagene summarises abstract properties of groups of genes, e.g. sets of genes that tend to be expressed together [Svensson et al., 2020], that are relevant for the imputation task. In a second step, HYFA employs a custom message passing neural network [Gilmer et al., 2017] that operates on a 3-uniform hypergraph, yielding factorised individual, tissue, and metagene representations. Lastly, HYFA infers latent metagene values for the target tissue -a hyperedge-level prediction task -and maps these representations back to the original gene expression space. Through higher-order hyperedges, HYFA can also incorporate cell-type information and infer finer-grained cell-type-specific gene expression (Methods). Altogether, HYFA offers features to reuse knowledge across tissues, capture non-linear cross-tissue patterns of gene expression, learn rich representations of biological entities, and account for variable numbers of reference tissues.
Characterisation of cross-tissue relationships Characterising cross-tissue relationships at the transcriptome level can help elucidate coordinated gene regulation and expression, a fundamental phenomenon with direct implications on health homeostasis, disease mechanisms, and comorbidities [Roenneberg and Merrow, 2016, Davière and Achard, 2017, Bodine et al., 2021. We trained HYFA on bulk gene expression from the GTEx project (GTEx-v8; Methods) [Consortium, 2020] and assessed the cross-tissue gene expression predictability scores and quality of tissue embeddings ( Figure 2). The application of Uniform Manifold Approximation and Projection (UMAP) [McInnes et al., 2018] on the learnt tissue representations revealed strong clustering of biologically-related tissues (Figure 2a), including the gastrointestinal system (e.g. esophageal, stomach, colonic, and intestinal tissues), the female reproductive tissues (i.e. uterus, vagina, and ovary), and the central nervous system (i.e. the 13 brain tissues). The clustering properties were robust across UMAP runs and could be similarly appreciated using other dimensionality reduction algorithms such as t-distributed Stochastic Neighbor Embedding (t-SNE) [Van der Maaten and Hinton, 2008]. For every pair of reference and target tissues in GTEx, we then computed the Pearson correlation coefficient ρ between the predicted and actual gene expression, averaged the scores across individuals, and used a cutoff of ρ > 0.5 to depict the top pairwise associations (Figure 2b). We observed connections between most GTEx tissues and whole blood, which suggests that blood-derived gene expression is highly informative of (patho)physiological processes in other tissues [Ray et al., 2007]. Notably, brain tissues and the pituitary gland were strongly associated with several target tissues (ρ > 0.5), including gastrointestinal tissues (i.e. esophagus, stomach, and colon), the adrenal gland, and skeletal muscle, which may account for known disease comorbidities.
Predictability of tissue-specific gene expression from whole blood transcriptome Knowledge about tissue-specific patterns of gene expression can increase our understanding of disease biology, facilitate the development of diagnostic tools, and improve patient subtyping [Lage et al., 2008, Basu et al., 2021, but most tissues are inaccessible or difficult to acquire. To address this challenge, we studied to what extent HYFA can recover tissue-specific gene expression from whole-blood transcriptomic measurements ( Figure 3). For each individual in the test set with measured wholeblood gene expression, we employed the hypergraph neural network to predict tissue-specific gene expression in the individual's remaining collected tissues. We evaluated imputation performance using the Pearson correlation coefficient ρ between the high-dimensional inferred gene expression and the ground-truth target samples. We observed strong prediction scores for esophageal tissues (muscularis: ρ = 0.488, gastro: ρ = 0.461, mucosa: ρ = 0.356), heart tissues (left ventricle: ρ = 0.479, atrial: ρ = 0.455) and lung (ρ = 0.471), while Epstein Barr virus-transformed lymphocytes (ρ = 0.057), an accessible and renewable resource for functional genomics, was a notable outlier. We compared our method with TEEBoT [Basu et al., 2021] (without expression single-nucleotide polymorphism information), which first projects the high-dimensional blood expression data into a low-dimensional space through principal component analysis (30 components; 75-80% explained variance) and then performs linear regression to predict the gene expression of the target tissue. Overall, TEEBoT and HYFA attained comparable scores when a single tissue (i.e. whole blood) was used as reference ( Figure 3) and both methods outperformed standard imputation approaches (mean imputation, blood surrogate, and k nearest neighbours; Figure 4a).
Multiple reference tissues improve performance We hypothesised that using multiple tissues as reference would improve downstream imputation performance. To evaluate this, we selected individuals with measured gene expression both at the target tissue and 4 reference (accessible) tissues (whole blood, skin sun exposed, skin not sun exposed, and adipose subcutaneous) and employed HYFA to impute target expression values (Figures 5 and 6). We discarded underrepresented target tissues with less than 25 test individuals. Relative to using whole blood in isolation, using all accessible tissues as reference resulted in improved performance for 14 out of 18 target tissues (Figure 6), highlighting the benefit of leveraging multiple reference tissues. This particularly boosted imputation scores for esophageal tissues (muscularis: ∆ρ = 0.068, gastro: ∆ρ = 0.061, mucosa: ∆ρ = 0.048), colonic tissues (transverse: ∆ρ = 0.065, sigmoid: ∆ρ = 0.056), and artery tibial (∆ρ = 0.079) ( Figure 5). In contrast, performance for the pituitary gland (∆ρ = −0.011), lung (∆ρ = −0.003), and stomach (∆ρ = −0.002) remained stable or dropped slightly. Moreover, the performance gap between HYFA and TEEBoT (trained on the set of complete multi-tissue samples) further widened relative to the single-tissue scenario (Figures 4b and 5) -HYFA obtained better performance in all target tissues (49 out of 49). We attribute the improved scores to HYFA's ability to process a variable number of reference tissues, reuse knowledge across tissues, and capture non-linear patterns, resulting in increased statistical strength and better generalisation to unseen combinations of reference-target tissues.
Inference of cell-type signatures We next investigated the ability of HYFA to predict cell-type-specific signatures in a given tissue of interest. We first matched donors with collected bulk gene expression in GTEx-v8 Consortium [2020] and individuals with measured single-nucleus RNA-seq profiles from GTEx-v9 [Eraslan et al., 2021], resulting in 16 individuals. We selected the top 3000 highly variable genes using the Scanpy function scanpy.pp.highly_variable_genes with flavour setting seurat_v3 [Wolf et al., 2018, Stuart et al., 2019 (Methods). We then aggregated the snRNA-seq profiles for each cell-type, tissue, and individual in GTEx-v9, yielding 226 unique individual-specific signatures. We discarded under-represented cell-types with less than 10 signatures and split the signatures into train and test ( Figure  7). Next, we trained HYFA to infer the cell-type, tissue, and individual-specific signatures from the multi-tissue bulk expression profiles. To attenuate the small data size problem, we applied transfer learning on the model trained for the multi-tissue imputation task (Methods). We observed strong prediction scores for vascular endothelial cells (heart: ρ = 0.79; breast: ρ = 0.87, esophagus muscularis: ρ = 0.79) and fibroblasts (heart: ρ = 0.80; breast: ρ = 0.89, esophagus muscularis: ρ = 0.80). Strikingly, HYFA recovered the cell-type profiles of tissues that were never observed in the train set with high correlation (Figures 7 and 8), e.g. skeletal muscle (vascular endothelial cells: ρ = 0.80, fibroblasts: ρ = 0.76, pericytes/SMC: ρ = 0.70), demonstrating the benefits of the transfer learning approach and the factorised tissue representations. Our results highlight the potential of HYFA to impute unknown cell-type signatures even for tissues that were not considered in the original snRNA-seq dataset [Eraslan et al., 2021].
Multi-tissue imputation improves eQTL detection Gene expression acts as an intermediate molecular trait between DNA and phenotype and, therefore, genetic mapping of genome-wide gene expression can shed light on the genetic architecture and molecular basis of complex traits. The GTEx project has enabled the identification of numerous genetic associations with gene expression across a broad collection of tissues [Consortium, 2020], also known as expression Quantitative Trait Loci (eQTLs) [Nica and Dermitzakis, 2013]. However, eQTL datasets are characterised by small sample sizes, especially for difficult-to-acquire tissues and cell types, reducing the statistical power to detect eQTLs [Rockman and Kruglyak, 2006]. To address this problem, we employed HYFA to impute the transcript levels of every uncollected tissue for each individual in GTEx, yielding a complete gene expression dataset of 834 individuals and 49 tissues. We then performed eQTL mapping (Methods) on the original and imputed datasets and observed a substantial gain in the number of detected eQTLs ( Figure 9). Notably, this metric increased for tissues with low sample size (Spearman correlation coefficient ρ = −0.73) -which are most likely to benefit from borrowing information across tissues with shared regulatory architecture. Minor salivary displayed the largest gain with a 626% increase in number of eQTLs, while the increase for whole blood was only 5%. Our results generate a large catalog of new tissue-specific eQTLs (Data availability), with potential to enhance our understanding of how regulatory variation mediates variation in complex traits, including disease susceptibility.
Brain-gut axis The brain-gut axis is a bidirectional communication system of signalling pathways linking the central and enteric nervous systems. We investigated the extent to which the transcriptomes of tissues from the gastrointestinal system are predictive of gene expression in brain tissues. We selected all the unseen individuals with simultaneous measurements in gastrointestinal tissues (i.e. oesophago gastric junction) and brain tissues (i.e. frontal cortex, hippocampus, and anterior cingulate) and employed HYFA to predict the expression values of brain tissues ( Figure 10). We observed a small number of individuals with measurements in both brain and non-brain tissues (Supplementary Materials C). After ranking the genes according to their predicability scores and selecting the top 1000 genes for each brain tissue (Venn diagram; Figure 10), we found considerable overlap between the 3 brain tissues (153 common genes in the intersection). We then used Enrichr [Kuleshov et al., 2016] with the gene sets GO_Biological_Process_2021 and GO_Molecular_Function_2021 to identify the enriched Gene Ontology (GO) terms for the shared genes at the intersection. Overall, the best predicted genes were enriched in multiple signalling-related terms (e.g. cytokine receptor activity and interleukin-1 receptor activity). This aligns with studies that highlight that gut microbes communicate with the central nervous system through endocrine and immune signalling mechanisms [Martin et al., 2018]. Genes in the intersection were also notably enriched in the ciliary neurotrophic factor receptor activity (molecular function), which plays an important role in the survival of neurons [Davis et al., 1991], the development of the enteric nervous system [Liu, 2018], and the control of body weight [Xu and Xie, 2016]. Moreover, our results suggest an association with the Receptor for Advanced Glycation Endproducts (RAGE), which has been linked to inflammation-related pathological states such as vascular disease, diabetes, and neurodegeneration [Sparvero et al., 2009].
HYFA-learned metagenes capture known biological pathways A key feature of HYFA is that it reuses knowledge across tissues and metagenes, allowing to exploit shared regulatory patterns. We explored whether HYFA's inductive biases encourage learning biologically relevant metagenes. To determine the extent to which metagene-factors relate to known biological pathways, we applied Gene Set Enrichment Analysis (GSEA) [Subramanian et al., 2005] to the gene loadings of HYFA's encoder (Methods). Similar to [Zhao et al., 2021], for a given query gene set, we calculated the maximum running sum of enrichment scores by descending the sorted list of gene loadings for every metagene and factor. We then computed pathway enrichment p-values through a permutation test and employed the Benjamini-Hochberg method to correct for multiple testing. In total, we identified 18683 statistically significant enrichments (FDR < 0.05) of KEGG biological processes across all HYFA metagenes and factors ( Figure 11a). Among these, 2109 corresponded to signalling pathways and 1300 to pathways of neurodegeneration ( Figure 11b). We observed considerable overlap between several metagenes in terms of enriched biologically related pathways, e.g. factor 95 of metagene 11 had the lowest FDR for both Alzheimer's disease (FDR < 0.001) and Amyotrophic Lateral Sclerosis (FDR < 0.001) pathways ( Figure 11d and 11f). We then applied UMAP to project the latent values of metagene 11 and observed relatively strong clustering of individuals with Amyotrophic Lateral Sclerosis (spinal cord; Figure 11c) and Alzheimer's disease or Dementia (brain cortex, Figure 11e). Altogether, our analysis suggests that HYFA-learned metagenes capture information about known pathways and are therefore amenable to biological interpretation.

Discussion
Effective multi-tissue omics integration promises a system-wide view of human physiology, with potential to shed light on intra-and multi-tissue molecular phenomena. Such an approach challenges single-tissue and conventional integration approaches -often unable to model a variable number of tissues with sufficient statistical strength -necessitating the development of scalable, non-linear, and flexible methods. Here we developed HYFA (Hypergraph Factorisation), a parameter-efficient approach for joint multi-tissue and cell-type gene expression imputation that imposes strong inductive biases to learn entity-independent relational semantics and demonstrates excellent imputation capabilities. The hypergraph factorisation framework is flexible (it supports k-uniform hypergraphs) and may find application beyond computational genomics.
We performed extensive benchmarks on GTEx (v8) data [Consortium, 2020], the most comprehensive human transcriptome resource available, and evaluated imputation performance over a broad collection of tissues and cell types. In addition to standard transcriptome imputation approaches, we compared our method with TEEBoT [Basu et al., 2021], a linear method that predicts target gene expression from the principal components of the reference expression. In the single-tissue reference scenario, HYFA and TEEBoT attained comparable imputation scores, outperforming standard methods. In the multi-tissue reference scenario, HYFA consistently outperformed TEEBoT and standard approaches for all target tissues, demonstrating the effectiveness of HYFA's capabilities to borrow non-linear information across a variable number of tissues and exploit shared molecular patterns.
In addition to imputing tissue-level transcriptomics, we investigated the ability of HYFA to predict cell-type-level gene expression from multi-tissue bulk expression measurements. Through transfer learning, we trained HYFA to infer cell-type signatures from a cohort of single-nucleus RNA-seq [Eraslan et al., 2021] with matching GTEx-v8 donors. The inferred cell-type signatures exhibited a strong correlation with the ground truth despite the low sample size, indicating that HYFA's latent representations are rich and amenable to knowledge transfer. Strikingly, HYFA also recovered cell-type profiles from tissues that were never observed at transfer time, pointing to HYFA's ability to leverage gene expression programs underlying cell-type identity [Kotliar et al., 2019] even in tissues that were not considered in the original study of Eraslan et al. [2021]. For future work, we will investigate extensions of HYFA to resolve the cell-type composition of observed and imputed tissue-specific gene expression.
In post-imputation analysis, we studied whether the imputed data improves eQTL discovery. We employed HYFA to impute the gene expression levels of every uncollected tissue in GTEx-v8, yielding a complete dataset, and performed eQTL mapping. Compared to the original dataset, we observed a substantial gain in number of detected eQTLs, from 5% for whole blood to 626% for minor salivary. The increase peaked for tissues with low sample sizes, which are the ones expected to benefit the most from knowledge sharing across tissues. Our results uncover a large number of previously undetected tissue-specific eQTLs and highlight the ability of HYFA to exploit shared regulatory information across tissues.
In addition to imputing uncollected gene expression measurements, HYFA can provide insights on coordinated gene regulation and expression mechanisms across tissues. We analysed to what extent tissues from the gastrointestinal system are informative of gene expression in brain tissues -an important question that may shed light on the biology of the brain-gut axis -and identified enriched biological processes and molecular functions. Through Gene Set Enrichment Analysis [Subramanian et al., 2005], we further introspected the HYFA-learned metagenes and found a substantial amount of enriched pathways, opening the door to further biological interpretations. We believe that explainability methods could be integrated in the HYFA pipeline to generate interesting biological insights. Future work might also seek to impose stronger inductive bias to ensure that metagenes are identifiable and robust to batch effects.
We believe that HYFA, as a versatile graph representation learning framework, provides a novel methodology for effective integration of large-scale multi-tissue biorepositories. HYFA's implementation is available online at https: //github.com/rvinas/HYFA.

Problem formulation
Suppose we have a transcriptomics dataset of N individuals/donors, T tissues, and G genes. For each individual i ∈ {1, ..., N }, let X i ∈ R T ×G be the gene expression values in T tissues and define the donor's demographic information by u i ∈ R C , where C is the number of covariates. Denote by x (k) i the k-th entry of X i , corresponding to the expression values of donor i measured in tissue k. For a given donor i, let T (i) represent the collection of tissues with measured expression values. These sets might vary across individuals. LetX i ∈ (R ∪ { * }) T ×G be the measured gene expression values, where * denotes unobserved, so thatx

Multi-tissue model
An important challenge of modelling multi-tissue gene expression is that a different set of tissues might be collected for each individual. Moreover, the data dimensionality scales rapidly with the total number of tissues and genes. To address these problems, we represent the data in a hypergraph and develop a parameter-efficient neural network that operates on this hypergraph. Throughout, we make use of the concept of metagenes [Brunet et al., 2004, Raychaudhuri et al., 1999. Each metagene characterises certain gene expression patterns and is defined as a positive linear combination of multiple genes [Brunet et al., 2004, Raychaudhuri et al., 1999.

Hypergraph representation
We represent the data in a hypergraph consisting of three types of nodes: donor, tissue, and metagene nodes.
where V d is a set of donor nodes, V m is a set of metagene nodes, V t is a set of tissue nodes, and E is a set of multi-attributed hyperedges. Each hyperedge connects an individual i with a metagene j and a tissue k if k ∈ T (i), where T (i) are the collected tissues of individual i. The set of all hyperedges is defined as ij are hyperedge attributes that describe characteristics of the interacting nodes, i.e. features of metagene j in tissue k for individual i.
The hypergraph representation allows representing data in a flexible way, generalising the bipartite graph representation from You et al. [2020]. On the one hand, using a single metagene results in a bipartite graph where each edge connects an individual i with a tissue k. In this case, the edge attributes e On the other hand, using multiple metagenes leads to a hypergraph where each individual i is connected to tissue k through multiple hyperedges. For example, it is possible to construct a hypergraph where genes and metagenes are related by a one-to-one correspondence, with hyperedge attributes e (k) ij . The number of metagenes thus controls a spectrum of hypergraph representations and, as we shall see, can help alleviate the inherent over-squashing problem of graph neural networks.

Message passing neural network
Given the hypergraph representation of the multi-tissue transcriptomics dataset, we now present a parameter-efficient graph neural network to learn donor, metagene, and tissue embeddings, and infer the expression values of the unmeasured tissues. We start by computing hyperedge attributes from the multi-tissue expression data. Then, we initialise the embeddings of all nodes in the hypergraph, construct the message passing neural network, and define an inference model that builds on the latent node representations obtained via message passing.
Computing hyperedge attributes We first reduce the dimensionality of the measured transcriptomics values. For every individual i and measured tissue k, we project the corresponding gene expression values x where M , the number of metagenes, is a user-definable hyperparameter and W j ∀j ∈ 1..M are learnable parameters. In addition to characterising groups of functionally similar genes, employing metagenes reduces the number of messages being aggregated for each node, addressing the over-squashing problem of graph neural networks (Supplementary Materials G).
Initial node embeddings We initialise the node features of the individual V p , metagene V m , and tissue V t partitions with learnable parameters and available information. For metagene and tissue nodes, we use learnable embeddings as initial node values. The idea is that these weights, which will be approximated through gradient descent, should summarise relevant properties of each metagene and tissue. We initialise the node features of each individual with the available demographic information u i of each individual i (e.g. age and sex). Importantly, this formulation allows transfer learning between sets of distinct donors.
Message passing layer We develop a custom GNN layer to compute latent donor embeddings by passing messages along the hypergraph. At each layer of the GNN, we perform message passing to iteratively refine the individual node embeddings. We do not update the tissue and metagene embeddings during message passing, in a similar vein as knowledge graph embeddings [Bordes et al., 2013], because their node embeddings already consist of learnable weights that are updated through gradient descent. However, if we foresee that unseen entities will be present at test time (e.g. new tissue types), our approach can be extended by initialising their node features with constant values and introducing node-type-specific message passing equations.
where the functions ϕ e and ϕ h are edge and node operations that we model as multilayer perceptrons (MLP), and ϕ a is a function that determines the aggregation behavior. In its simplest form, choosing ϕ a (h m j , h t k , m ijk ) = 1 G|T (i)| m ijk results in average aggregation. Optionally, we can stack several message passing layers to increase the expressivity of the model.
The architecture is flexible to allow the following extensions: • Incorporation of information about the individual embeddings h d i into the aggregation mechanism ϕ a . • Incorporation of target tissue embeddings h t u , for a given target tissue u, into the aggregation mechanism ϕ a . • Update hyperedge attributes e (k) ij at every layer.
Aggregation mechanism In practice, the proposed hypergraph neural network suffers from a bottleneck. In the aggregation step, the number of messages being aggregated is M |T (i)| for each individual i. In the worst case, when all genes are used as metagenes (i.e. M = G; it is estimated that humans have around G ≈ 25, 000 protein-coding genes), this leads to serious over-squashing -large amounts of information are compressed into fixed-length vectors [Alon and Yahav, 2020]. Fortunately, choosing a small number of metagenes reduces the dimensionality of the original transcriptomics values which in turn alleviates the over-squashing and scalability problems (Supplementary Materials G). To further attenuate over-squashing, we propose an attention-based aggregation mechanism ϕ a that weighs metagenes according to their relevance in each tissue: where || is the concatenation operation and a and W are learnable parameters. The proposed attention mechanism, which closely follows the neighbour aggregation method of graph attention networks [Brody et al., 2022, Veličković et al., 2018, computes dynamic weighting coefficients that prioritise messages coming from important metagenes. Optionally, we can leverage multiple heads [Vaswani et al., 2017] to learn multiple modes of interaction and increase the expressivity of the model.
Hypergraph model The hypergraph model, which we define as f , computes latent individual embeddingsĥ d i from incomplete multi-tissue expression profiles asĥ

Downstream imputation tasks
The resulting donor representationsĥ d i summarise information about a variable number of tissue types collected for donor i, in addition to demographic information. We leverage these embeddings for two downstream tasks: inference of gene expression in uncollected tissues and prediction of cell-type signatures.

Inference of gene expression in uncollected tissues
Predicting the transcriptomic measurementsx (k) i of a tissue k (e.g. uncollected) is achieved by first recovering the latent metagene valuesê (k) ij for all metagenes j ∈ 1..M , a hyperedge-level prediction task, and then decoding the gene expression values from the predicted metagene representationsê (k) ij with an appropriate probabilistic model.

Prediction of hyperedge attributes
To predict the latent metagene attributesê (k) ij for all j ∈ 1..M , we employ a multilayer perceptron that operates on the factorised metagene h m j and tissue representations h t k as well as the latent , where the MLP is shared for all combinations of metagenes, individuals, and tissues. Negative binomial imputation model For raw count data, we use a negative binomial likelihood. To decode the gene expression values for a tissue k of individual i, we define the following probabilistic model p(x where NB is a negative binomial distribution. The mean µ (k) ij and dispersion θ (k) ij parameters of this distribution are computed as follows: are mean gene-wise proportions; W s , W θ , b s , and b θ are learnable parameters; and l (k) i is the library size, which is modelled with a log-normal distribution where W ν , W ω , b ν , and b ω are learnable parameters. Optionally, we can use the observed library size.
Gaussian imputation model For normalised gene expression data (i.e. inverse normal transformed data), we use the following Gaussian likelihood: where the mean µ ij and standard deviation σ ij are computed as follows: where W µ , W σ , b µ , and b σ are learnable parameters and softplus(x) = log 1 + exp(x) .
Optimisation We optimise the model to maximise the imputation performance on a dynamic subset of observed, pseudo-missing tissues, similar to Viñas et al. [2021]. For each individual i, we randomly select a subset C ⊂ T (i) of pseudo-observed tissues and treat the remaining tissues U = T (i) − C as unobserved. We then compute the individual embeddingsĥ d i using the gene expression of pseudo-observed tissues C and minimise the loss: which corresponds to the average negative log-likelihood across pseudo-missing tissues. Importantly, the pseudo-mask mechanism generates different sets of pseudo-missing tissues for each individual, effectively enlarging the number of training examples and regularising our model. We summarise the training algorithm in Algorithm 1 (Supplementary Materials E).
Inference of gene expression from uncollected tissues At test time, we infer the gene expression valuesx of an uncollected tissue v from a given donor i via the mean, i.e.x i . Alternatively, we can draw random samples from the conditional predictive distribution p(x

Prediction of cell-type signatures
We next consider the problem of imputing cell-type signatures in a tissue of interest. We define a cell-type signature as the sum of gene expression profiles across cells of a given cell-type in a certain tissue. Formally, let x (k,q) i be the gene expression signature of cell-type q in a tissue of interest k of individual i. Our goal is to infer x (k,q) i from the multi-tissue gene expression measurementsX i . To achieve this, we first compute the hyperedge features of a hypergraph consisting of 4-node hyperedges and then infer the corresponding signatures with a zero-inflated model.

Prediction of hyperedge attributes
We consider a hypergraph where each hyperedge groups an individual, a tissue, a metagene, and a cell-type node. For all metagenes j ∈ 1..M , we compute latent hyperedge attributesê (k,q) ij for a cell-type q in a tissue of interest k of individual i as follows: where h c q are parameters specific to each unique cell-type q and the MLP is shared for all combinations of metagenes, individuals, tissues, and cell-types.

Zero-inflated model
We employ the following probabilistic model: where ZINB is a zero-inflated negative binomial distribution. The mean µ , and dropout probability π (k,q) ij parameters are computed as: where W s , W θ , W π , b s , b θ , and b π are learnable parameters; n and b ω are learnable parameters. Optionally, we can use the observed library size.
Optimisation Single-cell transcriptomic studies typically measure single-cell gene expression for a limited number of individuals, tissues, and cell-types, so aggregating single-cell profiles per individual, tissue, and cell-type often results in small sample sizes. To address this challenge, we apply transfer learning by pre-training the hypergraph model f on the multi-tissue imputation task (see Section 4.3.1) and then fine-tuning the parameters of the signature inference module (Section 4.3.2) on the cell-type signature profiles. Concretely, we minimise the loss: , which corresponds to the negative log-likelihood of the observed cell-type signatures.
Inference of uncollected gene expression To infer the signature of a cell-type q in a certain tissue v of interest, we first compute the latent individual embeddingsĥ d i from the multi-tissue profilesX i and then compute the mean of the distribution p(x ). Alternatively, we can draw random samples from that distribution.

Data missingness assumption
By employing maximum likelihood inference on the observed data (Supplementary Materials), HYFA assumes that the data (i.e. tissues) are Missing At Random (MAR; Little and Rubin [2019]), that is, the missingness mechanism is independent of the unobserved data. Training HYFA via variational inference (Supplementary Materials) also necessitates the MAR assumption which, similar to Mattei and Frellsen [2019], arises from maximising the loglikelihood of the observed data through the Evidence Lower Bound (ELBO). The MAR assumption is less restrictive than the Missing At Completely at Random (MCAR) assumption -the missingness pattern is independent of the observed and unobserved data -of other methods such as mean imputation and GAIN [Yoon et al., 2018].
HYFA does not support data Missing Not At Random (MNAR), where the missingness mechanism depends on the unobserved data, i.e. the probability of being missing depends on unknown reasons [Van Buuren, 2018]. To handle this scenario, we would need to model the joint distribution p(X, R) of the observed dataX and the missingness mechanism R, that is, the missingness mechanism would be nonignorable and would need to be explicitly modelled. This could be achieved through selection modeling [Heckman, 1976], which factorises the joint distribution as p(X, R) = p(R|X)p(X), or pattern-mixture models [Glynn et al., 1986], which decompose the joint as p(X, R) = p(X)p(X|R). In general, it is impossible to test if MAR holds in a dataset [Schafer and Graham, 2002], but the impact of incorrectly assuming MAR is often minor [Collins et al., 2001].

Related work
Gene expression imputation methods A standard approach for imputation of uncollected transcriptomics values is to use a proxy tissue (e.g. whole blood) as a surrogate [Hoon et al., 2000]. However, gene expression is known to be tissue and cell-type specific, limiting the effectiveness of this technique. Other related studies infer tissue-specific gene expression from genetic information.  propose a mixed-effects model to infer uncollected data in multiple tissues from expression quantitative trait loci (eQTLs). Sul et al. [2013] introduce a model termed Meta-Tissue that aggregates information from multiple tissues to increase statistical power of eQTL detection. Nonetheless, these approaches do not model the complex relationships between measured and unmeasured gene expression traits among tissues and cell types, and individual-level genetic information (such as at eQTLs) is often unavailable and subject to privacy concerns. Instead, recent closely related work relies on linear factor analysis and dimensionality reduction techniques. TEEBoT (Tissue Expression Estimation using Blood Transcriptome) [Basu et al., 2021] projects gene expression from a single reference tissue (i.e. whole blood) into a low-dimensional space via principal component analysis (PCA), followed by linear regression to impute uncollected values. HYFA allows a departure from the linearity assumption of TEEBoT and also handles a variable number of reference tissues.
Knowledge graph embedding techniques Our framework leverages ideas from knowledge graph embedding techniques by using learnable embeddings for biological entities (i.e. tissues, cell-types, and metagenes). Since the advent of word embeddings [Mikolov et al., 2013], several approaches have emerged to learn vector representations of entities and relations in knowledge graphs [Bordes et al., 2013, Wang et al., 2014, Trouillon et al., 2016, Dettmers et al., 2018. TransE [Bordes et al., 2013] represents entities as low-dimensional vectors and relationships as translations in the embedding space, and optimises parameters through an energy-based objective. TransH [Wang et al., 2014] extends the TransE framework by projecting the entity embeddings into relation-specific hyperplanes. ComplEx [Trouillon et al., 2016] utilises complex vectors that capture antisymmetric entity relations. ConvE [Dettmers et al., 2018] models interactions between input entities and relations through convolutional and fully-connected layers. Despite all the recent advances, knowledge graph embeddings have been understudied for modelling higher-order structures (i.e. hyperedges). Moreover, while these methods are capable of link prediction, they are limited to a transductive setting, where the full set of entities (e.g. individuals) must be known at train time [Teru et al., 2020].
Graph representation learning Graph neural networks remedy this problem by leveraging the structure and properties of graphs to compute node features, allowing to handle unseen entities at inference time (e.g. individuals). The closest approach to our framework is GRAPE [You et al., 2020]. GRAPE represents tabular data as a bipartite graph, where observations and features are two types of nodes, and the observed feature values are attributed edges between the nodes [You et al., 2020]. Imputation of missing features then corresponds to an edge-level prediction task. HYFA subsumes GRAPE's bigraph in that our hypergraph becomes a bipartite graph when a single metagene is employed. This allows for a trade-off between feature granularity and over-squashing, which happens when information from a large receptive field is compressed into fixed-length node vectors [Alon and Yahav, 2020]. In terms of message passing, HYFA distinguishes between dynamic nodes (udpated during message passing) and static nodes (with learnable node features that are not updated during message passing), eliminating the dependency of tissue and metagene representations on donor features and, by transitivity, undesired dependencies between inviduals. HYFA is thus a hybrid and flexible approach that combines features from knowledge graph embedding and graph representation learning techniques.
Single-cell variational inference Our framework is related to single-cell variational inference (scVI) [Lopez et al., 2018] in that it can also be optimised via variational inference (Supplementary Materials), e.g. via a (zero-inflated) negative binomial likelihood, treating the individual representations as latent variables. In contrast to scVI, however, HYFA offers features to handle a variable number of reference tissues. It also incorporates inductive biases to reuse knowledge across tissues, allowing the model to scale to larger multi-tissue samples.

eQTL mapping
The breadth of tissues in the GTEx (v8) collection enabled us to comprehensively evaluate the extent to which eQTL discovery could be improved through the HYFA-imputed transcriptome data. We mapped eQTLs that act in cis to the target gene (cis-eQTLs), using all SNPs within ± 1 Mb of the transcription start site of each gene. For the imputed and the original (incomplete) datasets, we considered SNPs significantly associated with gene expression, at a false discovery rate ≤ 0.10. We applied the same GTEx eQTL mapping pipeline, as previously described [Consortium, 2015], to the imputed and original datasets to quantify the gain in eQTL discovery from the HYFA-imputed dataset.

Pathway enrichment analysis
Similar to [Zhao et al., 2021], we employed Gene Set Enrichment Analysis (GSEA) [Subramanian et al., 2005] to relate HYFA's metagene factors to known biological pathways. This is advantageous to over-representation analysis, which requires selecting an arbitrary cutoff to select enriched genes. GSEA, instead, computes a running sum of enrichment scores by descending a sorted gene list [Subramanian et al., 2005, Zhao et al., 2021.
We applied GSEA to the gene loadings in HYFA's encoder. Specifically, let W j ∈ R F ×G be the gene loadings for metagene j, where F is the number of factors (i.e. number of hyperedge attributes) and G is the number of genes (Equation 1). For every factor in W j , we employed blitzGSEA [Lachmann et al., 2022] to calculate the running sum of enrichment scores by descending the gene list sorted by the factor's gene loadings. The enrichment score for a query gene set is the maximum difference between p hit (S, i) and p miss (S, i) [Zhao et al., 2021], where p hit (S, i) is the proportion of genes in S weighted by their gene loadings up to gene index i in the sorted list [Zhao et al., 2021]. We then calculated pathway enrichment p-values through a permutation test (with N = 100 trials) by randomly shuffling the gene list. We employed the Benjamini-Hochberg method to correct for multiple testing.

GTEx bulk and single-nucleus RNA-seq data processing
The GTEx dataset is a public resource that has generated a broad collection of gene expression data collected from a diverse set of human tissues [Consortium, 2020]. We downloaded the data from the GTEx portal (Data availability). After the processing step, the GTEx-v8 dataset consisted of 15197 samples (49 tissues, 834 donors) and 12557 genes. The dataset was randomly split into 500 train, 167 validation, and 167 test donors. Each donor had an average of 18.22 collected tissues. The processing steps are described below.
2. Select set of overlapping genes across all tissues.

Normalise read counts across samples using the trimmed mean of M values (TMM) method [Robinson and
Oshlack, 2010].

Implementation and reproducibility
We report the selected hyperparameters in Supplementary Materials G. HYFA is implemented in Python [Van Rossum and Drake Jr, 1995]. Our framework and implementation are flexible (i.e. we support k-uniform hypergraphs), may be integrated in other bioinformatics pipelines, and may be useful for other applications in different domains. We used PyTorch [Paszke et al., 2019] to implement the model and scanpy [Wolf et al., 2018] to process the gene expression data. We performed hyperparameter optimisation with wandb [Biewald, 2020]. We employed blitzGSEA [Lachmann et al., 2022] for the pathway enrichment analysis. We also used NumPy [Harris et al., 2020], scikit-learn [Pedregosa et al., 2011], pandas [Wes McKinney, 2010], matplotlib [Hunter, 2007], and seaborn [Waskom, 2021]. Figure 1 was created with BioRender.com.

Data availability
The datasets analysed for this study, including bulk RNA-seq Consortium [2020]    Adipose Visceral Adr ena l Gla nd Ar te ry Ao rt a A rt e ry C o ro n a ry A m y g d a la A n t e r io r c in g u la t e C a u d a t e C e re b e ll a r C e re b e llu m Co rte x Fronta l Cortex Hippo camp us Hy po th ala m us N u cl e u s P u ta m e n S p in a l c o r d S u b s t a n t ia n ig r a B r e a s t M a m m a r y (a) UMAP representation of the tissue embeddings learnt by HYFA. Note that human body systems cluster in the embedding space (e.g. digestive system: stomach, small intestine, colon, esophagus; and central nervous system). (b) Network of tissues depicting the predictability of target tissues with HYFA. The dimension of each node is proportional to its degree. Edges from source to target tissues indicates an average Pearson correlation coefficient ρ > 0.5. Interestingly, central nervous system tissues strongly correlate with several non-brain tissues such as gastrointestinal tissues and skeletal muscles. Pearson correlation Prediction performance with accessible tissues as source tissues (whole blood, skin, and adipose subcutaneous) TEEBoT HYFA (accessible) Figure 5: Comparison between HYFA and TEEBoT when using all accessible tissues as reference. TEEBoT first projects the high-dimensional multi-tissue samples (each example corresponds to a single donor) into a low-dimensional space (30 components) and then predicts the gene expression of the target tissue via linear regression. Target tissues represented by less than 25 test individuals were discarded. Error bars indicate 95% confidence intervals. HYFA achieves better performance for all target tissue types. Pearson correlation Prediction performance with accessible tissues as source Whole Blood Skin Sun Epsd Skin Not Sun Epsd Adipose Subcutaneous Accessible All Figure 6: Prediction scores for different accessible tissues as reference. For each target tissue, we predicted the expression values based on accessible tissues (whole blood, skin sun exposed, skin not sun exposed, and adipose subcutaneous). We report the Pearson correlation coefficient between the predicted values and the actual gene expression values. For any given target tissue, we used the same set of individuals to evaluate performance, namely individuals in the test set with collected gene expression measurements in all the corresponding tissues. Target tissues represented by less than 25 test individuals were discarded. Error bars indicate 95% confidence intervals. The method attains the best performance in 14 out of 18 tissues when all accessible tissues are simultaneously used as reference. Enriched GO terms (intersection) Figure 10: Top predicted genes in multiple brain regions with the oesophago gastric junction as the reference tissue. (a) R 2 coefficient of top 20 predicted genes, ranked by average score. (b) Common genes in the top 1000 predicted genes for each brain tissue. (c) Enriched GO terms of the top shared genes at the interesection. The best predicted genes were enriched in signalling pathways (FDR < 0.05), consistent with studies reporting that gut microbes communicate to the central nervous system through endocrine and immune mechanisms. Moreover, these results highlight the potential association between the elements of the oesophago gastric junction and the ciliary neurotrophic factor, which has been linked to the survival of neurons [Davis et al., 1991] and the control of body weight [Xu and Xie, 2016].

D Connection with maximum likelihood
Let x obs be a random variable denoting the observed data (e.g. multi-tissue gene expression with missing values corresponding to uncollected tissues). Our optimisation procedure (Methods) splits x obs into pseudo-observedx obs and pseudo-missingx mis values, that is, x obs = (x obs ,x mis ). The log-likelihood of the observed data then corresponds to: log p(x obs ) = log p(x obs ,x mis ) = log p(x mis |x obs ) + log p(x mis ) and log p(x mis |x obs ) is precisely the quantity that our loss function is maximising through the pseudo-mask mechanism (Methods).

E Training algorithm
Optimisation. We optimise the model to minimise the mean squared error L between the normalised, ground-truth gene expression x (u) i and the imputed valuesx (u) i : where G is the number of genes. At train time, for any given individual, we dynamically mask out the expression values of a measured tissue type at random and treat them as uncollected, i.e. the ground truth. Algorithm 1 summarises the training algorithm. Parameters of inference model Given the updated donor representationsĥ p i , we compute the parameters of the inference model q(z i |X i , u i ) = N (z i ; µ i , diag(σ 2 i )) as follows: µ i = f µ (X i , u i ; ϕ) = MLP(ĥ p i ) log σ i = f σ (X i , u i ; ϕ) = MLP(ĥ p i ) where MLP denotes a multilayer perceptron. Parameters of generative model Assuming a Gaussian likelihood, for a given sample z i ∼ q(z i |X i , u i ), we compute the parameters of the generative model p(x   Table 2: Hyperparamter values considered for ablation studies. We used wandb [Biewald, 2020] to run Bayesian hyperparameter search over the variables and ranges considered. Note that with the Attentional GAT layer, the total dimension of the message m ijk is multiplied by the number of attention heads (here, 28 × 28 = 784).