# A robust kernel machine regression towards biomarker selection in multi-
omics datasets of osteoporosis for drug discovery

Md Ashad Alam, Hui Shen, and Hong-Wen Deng
Tulane Center for Biomedical Informatics and Genomics
Division of Biomedical Informatics and Genomics
Deming Department of Medicine
Tulane University, New Orleans, LA 70112, USA

###### Abstract

Many statistical machine approaches could ultimately highlight novel features
of the etiology of complex diseases by analyzing multi-omics data. However,
they are sensitive to some deviations in distribution when the observed
samples are potentially contaminated with adversarial corrupted outliers
(e.g., a fictional data distribution). Likewise, statistical advances lag in
supporting comprehensive data-driven analyses of complex multi-omics data
integration. We propose a novel non-linear M-estimator-based approach, “robust
kernel machine regression (RobKMR)," to improve the robustness of statistical
machine regression and the diversity of fictional data to examine the higher-
order composite effect of multi-omics datasets. We address a robust kernel-
centered Gram matrix to estimate the model parameters accurately. We also
propose a robust score test to assess the marginal and joint Hadamard product
of features from multi-omics data. We apply our proposed approach to a multi-
omics dataset of osteoporosis (OP) from Caucasian females. Experiments
demonstrate that the proposed approach effectively identifies the inter-
related risk factors of OP. With solid evidence (p-value $\leq 0.00001$),
biological validations, network-based analysis, causal inference, and drug
repurposing, the selected three triplets ((DKK1, SMTN, DRGX), (MTND5, FASTKD2,
CSMD3), (MTND5, COG3, CSMD3)) are significant biomarkers and directly relate
to BMD. Overall, the top three selected genes (DKK1, MTND5, FASTKD2) and one
gene (SIDT1 at p-values $\leq 0.001$) significantly bond with four drugs-
Tacrolimus, Ibandronate, Alendronate, and Bazedoxifene out of $30$ candidates
for drug repurposing in OP. Further, the proposed approach can be applied to
any disease model where multi-omics datasets are available.

##  1 Introduction

Biomedical technology has accelerated the cycle of multi-omics data discovery
for complex disease control and prevention. Multi-omics data can facilitate
our comprehensive understanding of the biological insight into the fundamental
mechanism of complex traits and diseases (e.g., osteoporosis (OP)), which
represent a significant burden in the global healthcare system [47, 16, 44,
53]. These are often a result of the composite interplay between multiple
layers of multi-omics data [28, 30, 48]. An individual-omics technique may
detect a subset of biomarkers of a complex disease and thus can only capture
changes in a small portion of the biological cascade [30, 59]. However, a
robust integrated risk factors analysis and understanding of comprehensive
interactions between various omics data are still rare and challenging.

In the last decade, several statistical methods have been used to detect gene-
gene interactions [21, 41, 13]. Logistic regression, multifactor
dimensionality reduction, linkage disequilibrium, entropy-based statistics,
and the sequence kernel association test are the example of such methods [55,
20]. While most of these methods are based on the unit association of single
nucleotide polymorphisms (SNPs), testing the associations between the
phenotype and SNPs has limitations. It is not sufficient for interpreting
gene-gene interactions [40, 20]. However, many researchers have shown that
alterations at other levels (i.e., transcriptome, epigenome, proteome, etc.)
also play a significant role in complex traits [51]. Thus, only focusing on
genomic data is not sufficient to identify the related risk factors for
complex traits. To address this concern, researchers have extended the model-
based kernel machine method proposed initially for detecting gene-gene
interactions to analyze the interaction of genomic and multiple extra-genomic
data to select discriminatory biomarkers [7, 24, 37, 1].

Statistical machine learning approaches (e.g., kernel-based methods) offer
beneficial ways to study an extensive collection of genetic variants
associated with complex traits. It helps to explore the relationship between
genetic markers and a disease state [55, 3, 45, 57]. A kernel machine method
has been proposed to jointly model the genetic and non-genetic features and
their interactions [24]. While these methods could ultimately highlight novel
features of the etiology of complex diseases, they cannot be reliably used for
genomic data with multiple extra-genomic data. Recently, researchers have
extended these methods to examine the higher-order interaction effect of
multimodal ($\leq$ 2 data types) datasets [8]. However, the major limitation
of these methods is their ineffectiveness in presenting outliers or
adversarial attacks that often occur in genomics and functional genomics
datasets [19, 56].

To date, almost all genomic analysis methods apply to genome data sets with
the assumption that the data sets come from a non-fictional data distribution
(a normal distribution only) [43]. This assumption can lead to an inaccurate
inference in real-world genomic data analysis. For example, Figure 1 presents
a density (a) and a volcano (b) plot of the genomic (SNP) and extra genomic
(RNA-seq and reduced representative bisulfite sequencing (RRBS)) data. This
figure clearly shows a fictional data distribution and outliers in the genomic
data. However, standard genomic data science approaches are sensitive to such
deviations in distribution when the observed samples are potentially
contaminated with adversarial corrupted outliers [6, 32, 23, 17]. Robust
learning approaches thus are the critical aspect in achieving true biomedical
genomic data integration and are necessary parameters that should not be
ignored if our approach to modeling complex diseases is to grow and evolve
[5]. The novelty is in the use of a sophisticated methodology to examine the
interaction and composite effects for network analysis of genomic systems in
complex diseases.

Figure 1:  The density (a) and volcano plot (b) of raw multi-omics data
clearly shows that a fictional data distribution and outliers in the dataset.

In biomedical data analysis, robustness is a critical and challenging issue
where outliers or adversarial attacks (targeted and untargeted attacks) often
occur. Robustness of efficiency means that even if the observations do not
match the distributional assumptions properly, the standard error of the
statistic will barely be affected. On the other hand, in classical statistical
methods, outliers can often cause havoc. To overcome this problem, since the
1960s, many robust methods have been developed, which are less sensitive to
the outliers [29, 26]. On the other hand, big data-powered machine learning
and deep learning approaches are sensitive to fictional data distribution. The
advanced methods also need a strong assumption that train and test data came
from the same distribution [43, 4]. Many robust/adversarial learning
techniques have been studied in biomedical data analysis, which is less
sensitive to contaminated data and distribution assumptions [52, 31, 2]. As a
consequence, robust non-linear integrated approaches (e.g., robust kernel-
based machine) are an essential feature in the comprehensive analysis of
multi-omics biomedical datasets [39, 11, 36]. The robust positive definite
kernel-based machine approach can overcome the non-linearity problem and
inconsistent results of multi-omics biomedical datasets.

Researchers have investigated the issue of robustness for the support vector
machine, [15, 18], kernel density estimation (Kim and Scott, 2012), kernel
canonical correlation, and kernel principal analysis, yet no well-founded
robust kernel machine regression method has been proposed. Motivated by these
considerations, we introduce a robust kernel machine approach to identify
composite effects in multi-omics datasets. Herein, we introduce a novel non-
linear M-estimator-based approach, “robust kernel machine regression", and
apply it to identify composite effects in multi-omics data of OP-related
traits (e.g., Bone mineral density (BMD)). To achieve robustness, we
formulated an empirical optimization problem by combining empirical
optimization with the idea of Huber’s and Hampel’s M-estimation model. The
proposed robust kernel machine regression can be computed efficiently via
kernelized iteratively re-weighted least squares (KIRWLS) [8].

To examine the robustness of our model, we compare the performance of the
proposed method with existing methods using synthesized and real OP datasets.
OP is a bone disorder that increases bone resorption and (or) decreases bone
formation by osteoclasts and osteoblasts. We used three-omics datasets of OP.
For more information on the data, we refer the reader to the data construction
section 2.1. This study explores novel genomes, epigenomes, and transcriptomes
that robustly and efficiently identify corresponding factors. To validate the
results, we performed network-based analysis, causal inference, GeneMANIA
prediction analysis 111http://genemania.org, GeneHancer analysis, and drug
repurposing. RobKMR shows that the selected three triplets ((DKK1, SMTN,
DRGX), (MTND5, FASTKD2, CSMD3), (MTND5, COG3, CSMD3)) are significant
biomarkers directly related to BMD. In general, the top three selected genes
(DKK1, MTND5, FASTKD2) and one gene (SIDT1, at p-values $\leq$ $0.001$) are
significantly bonded with four drugs- Tacrolimus, Ibandronate, Alendronate,
and Bazedoxifene out of $30$ candidates for drug repurposing in OP. The
following aspects make this paper highly novel.

  1. i.

We proposed a novel non-linear M-estimator-based robust kernel machine
regression for improving the robustness of statistical machine regression and
the diversity of fictional multi-omics data.

  2. ii.

We examined the robust test statistic of the different effects of multi-omics
data, including overall interactions and composites.

  3. iii.

We analyzed a simulated and a real multi-omics OP dataset to reveal that the
proposed method is highly efficient.

  4. iv.

To validate the results, we performed a comprehensive GeneMANIA, GeneHancer,
pathway, gene-gene-gene network, and causal analysis. The selected biomarkers
are significant and directly related to BMD.

  5. v.

For drug repurposing, we used molecular docking. Using molecular docking, we
observe that the top selected genes significantly bond with candidates lead
four drugs- Tacrolimus, Ibandronate, and Bazedoxifene for drug repurposing in
OP.

##  2 Materials and methods

Statistical machine-based approaches utilized in the integrated analysis of
multi-omics data that provide a general view of the biological insights of
complex diseases and traits [12, 35]. These integrated approaches facilitated
practical ways to identify interrelated effects, including interactions and
composites effects of multi-omics data. Moreover, the critical concern for
these approaches is the fictional data distribution and outliers in the
dataset. Such outliers or adversarial attacks often occur in genomics and
functional genomics and present a critical challenge for robust data science
methods for interactive analyses [43]. Consequently, this paper proposes a
robust kernel machine regression via a robust kernel Gram matrix. To that end,
different loss functions (e.g., Huber’s, Hampel’s, and Tukey’s bi-weight
Cauchy’s, Welsch’s, and Andrews’) is used in M-estimation(16, 59, 60). Unlike
the mean square error loss function, the derivative of these robust loss
functions is bounded. We propose an algorithm to estimate a centered robust
kernel Gram matrix using the weight of robust kernel mean element. The
proposed robust kernel machine regression can be computed efficiently via
kernelized iteratively re-weighted least squares (KIRWLS) [7, 38]. This method
examines the robust test statistic of the different effects, including joint
and higher-order interaction and composite effects for identifying risk
factors leading to recurrence of the desired complex disease and enhancing
disease prediction.

###  2.1 Dataset

We apply the proposed method to our generated multi-omics dataset from an
osteoporosis study. In this paper, we conducted integrate study of genome
(3,997,535 SNPs which annotated to 25, 442 genes), epigenome (46,690 CpG
methylations which annotated to 4,676 genes.), and transcriptome (22,682 genes
expression profiles) data from 57 Caucasian females with high BMD and 51 with
low BMD. This dataset (our Louisiana osteoporosis study) is available on
request at our lab and some of data has been already deposited in dbGaP
(phs001960.v1.p1) [8].

###  2.2 Robust model setting and estimation

Robustness to outliers, noisy samples, and heavy-tailed distributions is an
essential issue for statistical machine learning approaches, including kernel
machine regression. In kernel machine regression, we can reduce the effect of
outliers, noisy samples, or heavy-tailed distribution using robust
$M$-estimation. Consequently, we propose a robust kernel machine regression
via robust kernel Gram matrix.

###  2.3 Robust loss function

Different loss functions with bounded derivative are used in M (maximum
likelihood-type)-estimation. Huber’s, Hampel’s, and Tukey’s biweight Cauchy’s,
Welsch’s , and Andrews’ loss function are common loss functions for the
M-estiamtion [29, 27, 50, 54]. The basic notions of classical and robust loss
function which can be used for standard and robust kernel machine approach are
as follows:

The least-squares loss function is defined as

$\displaystyle\rho(t)=t^{2}/2.$ (1)

This mean square error loss function is a standard and well-known loss
function which is used the sum of all the squared differences between the
actual value and the estimated value in statistics, statistical machine
learning as well as in data science. The influence function of this function
is unbounded. Hence, least-square estimators are not robust.

The least-absolute loss function is defined as

$\displaystyle\rho(t)=|t|$ (2)

The sum of all the absolute deviations between the estimated value and the
true value is used to minimize the error. The least absolute value is an
unstable estimator because this loss function $|x|$ is not strictly convex in
x. While estimators reduce the influence of large errors than the least-
squares estimators, they still have an influence because the influence
function has no cut-off point.

The Huber’s loss function is defined as

$\displaystyle\rho(t)=\begin{cases}t^{2}/2,\qquad\qquad 0\leq t\leq c\\\
ct-c^{2}/2,\qquad c\leq t,\end{cases}$

where c ($c>0$) is a tuning parameter. This is a hybrid approach between
squared and absolute error losses functions. It is a parabola around the area
of zero and increases linearly at a given level (e.g., $|x|>c$). For almost
all situations, the Huber estimator is acceptable; very rarely, it has been
found to be inferior to some other loss -functions. However, due to the lack
of stability in the gradient values of the function ( e.g., its discontinuous
second derivative), this estimator is possibly insufficient.

The Hampel’s loss function is defined as:

$\displaystyle\rho(t)=\begin{cases}t^{2}/2,\qquad\qquad\qquad\qquad 0\leq
t\leq c_{1}\\\ c_{1}t-c_{1}^{2}/2,\qquad\qquad c_{1}\leq t<c_{2}\\\
-\frac{c_{1}}{2(c_{3}-c_{2})}(t-c_{3})^{2}+\frac{c_{1}(c_{2}+c_{3}-c_{1})}{2},\qquad
c_{2}\leq t<c_{3}\\\ \frac{c_{1}(c_{2}+c_{3}-c_{1})}{2},\qquad\qquad c_{3}\leq
t,\end{cases}$

where $c_{1}<c_{2}<c_{3}$ the non-negative free parameters that allow us to
control the degree of suppression.

The Tukey’s biweight loss functions is defined as:

$\displaystyle\rho(t)=\begin{cases}\frac{c^{2}}{6}(1-(1-(t/c)^{2})^{3}),\qquad\qquad
0\leq t\leq c\\\ \frac{c^{2}}{6},\qquad c\leq t,\end{cases}$

where $c>0$. Since the loss incurred by large residuals is constant, it is
even more insensitive to outliers. But, like Hubar loss, it reveals quadratic
behavior near the origin.

The Cauchy’s loss function is defined as:

$\displaystyle\rho(t)=\frac{c^{2}}{2}\rm{log}(1+(t/c)^{2}),$

where $c$ the non-negative free parameter. Unlike the mean square error and
least-absolute loss functions, the Cauchy’s loss function can alleviate the
influence of a considerable noise with a sample for estimating the residuals.
Thus, this loss function has less dependence on the noise distribution and is
more robust to the noise.

The Welsch’s function is defined as:

$\displaystyle\rho(t)=\frac{c^{2}}{2}[1-\rm{exp}(-(t/c)^{2})],$ (3)

where $c$ the non-negative free parameter. The Welsch’s loss functions try to
reduce the effect of significant errors further even suppress the outlier. The
windows of this loss function are like a bell. It emphasizes the impact of
data with $c$ close to zero and gradually reduces the weights of those further
from zero.
The Geman-MeClure is defined as:

$\displaystyle\rho(t)=\frac{t^{2}/2}{1+t^{2}}.$

The Geman-McClure loss function, which behaves almost quadratically for small
values and saturates for large ones, is similar to the truncated least squares
loss. It produces accurate results and should be preferred over the other
ones.

These loss functions hold basic assumptions of the loss functions (i) non-
decreasing, (ii) weight function exists and is finite, (iii) influence and
weight functions are continuous and bounded, and (iv) $\varphi(t)$ is
Lipschitz continuous [31].

Figure 2: Family of (a) Least-squares, (b) Least-absolute, (c) Huber’s, (d)
Hampel’s, (e) Tukey’s, (f) Cauchy, (g) Welsch, and (h) Geman-MeClure loss
functions (OF: original function ($\rho(\cdot))$; IF: influence function,
$\psi(\cdot)$; WF: weight function; and DIF: the derivative of influence
function).

###  2.4 Robust kernel Gram matrix

We propose a robust kernel Gram matrix using a robust kernel mean element. The
classical kernel mean element is the solution to the empirical risk
optimization problem, which is a least-square class of estimators. The method
of least squares is sensitive to the presence of outliers in the data. To
reduce the effect of outliers, we can use $M$-estimation for estimating kernel
mean element. Kim and Scot (2012) has been proposed the robust kernel ME for
density estimation [31]. The robust kernel ME, based on a robust loss function
$\zeta(t)$ on $t\geq 0$, is defined as

$\displaystyle\widehat{\mathcal{M}}_{R}=\operatornamewithlimits{argmin}_{f\in\mathcal{H}_{X}}\sum_{i=1}^{n}\zeta(\|\Phi(X_{i})-f\|_{\mathcal{H}_{X}}).$
(4)

Essentially Eq. (4) does not have a closed form solution, but using KIRWLS,
the solution of robust kernel mean is,

$\widehat{\mathcal{M}}_{R}^{(h)}=\sum_{i=1}^{n}w_{i}^{(h-1)}k_{X}(\cdot,X_{i}),$

where
$w_{i}^{(h)}=\frac{\varphi(\|\tilde{\Phi}(X_{i})-f^{(h)}\|_{\mathcal{H}_{X}})}{\sum_{b=1}^{n}\varphi(\|\Phi(X_{b})-f^{(h)}\|_{\mathcal{H}_{X}})}\,,\rm{and}\,\varphi(x)=\frac{\zeta^{\prime}(x)}{x}.$
the centered robust Gram matrix is

$\tilde{K}_{ij}=\langle\tilde{\Phi}(X_{i}),\tilde{\Phi}(X_{j})\rangle=(\mathbf{C}\mathbf{K}_{X}\mathbf{C}^{T})_{ij},$

where $\mathbf{K}_{X}=(k_{X}(X_{i},X_{j}))_{i=1}^{n}$ is a Gram matrix,
$\mathbf{1}_{n}=[1_{1},1_{2},\cdots,1_{n}]^{T}$ and
$\mathbf{C}=\mathbf{I}-\mathbf{1}_{n}\mathbf{w}^{T}$.

Given the weights of the robust kernel ME,
$\mathbf{w}=[w_{1},w_{2},\cdots,w_{n}]^{T}$, of a set of observations
$X_{i},\cdots,X_{n}$, the points

$\tilde{\Phi}(X_{i}):=\Phi(X_{i})-\sum_{a=1}^{n}w_{a}\Phi(X_{a})$

are centered. The algorithm of estimating robust Gram matrix is given in
Figure 1.

Input: $D=\\{\mathbf{X}_{1},\mathbf{X}_{2},\ldots\mathbf{X}_{n}\\}$ in
$\mathbb{R}^{m}$. The kernel matrix $\mathbf{K}$ with kernel $k$ and
$\mathbf{K}_{\mathbf{X}_{i}}=k(\cdot,\mathbf{X}_{i})$. Threshold $TH$, (e.g.,
$10^{-8}$). The objective function of robust mean element is

$M_{R}=\arg\min_{f\in\mathcal{H}}J(f),\qquad\rm{where}\,J(f)=\frac{1}{n}\sum_{i=1}^{n}\rho(\|K_{X_{i}}-f\|_{\mathcal{H}})$

  1. Do the following steps until:

$\frac{|J(M_{R}^{(h+1)})-J(M_{R}^{(h)})|}{J(M_{R}^{(h)})}<TH,$

where$M_{R}^{(h)}=\sum_{i=1}^{n}w_{i}\mathbf{K}_{\mathbf{X}_{i}},\,w_{i}^{(h)}=\frac{\varphi(\|\mathbf{K}_{\mathbf{X}_{i}}-M_{R}^{(h)}\|_{\mathcal{H}})}{\sum_{i=1}^{n}\varphi(\|\mathbf{K}_{\mathbf{X}_{i}}-M_{R}^{(h)}\|_{\mathcal{H}})}\,,\rm{and}\,\varphi(x)=\frac{\xi^{\prime}(x)}{x}$

     * (1)

Set $h=1$ and $w_{i}^{(0)}=\frac{1}{n}$.

     * (2)

Solve
$w_{i}^{(h)}=\frac{\varphi(\epsilon_{i}^{[h]})}{\sum_{i}^{n}\varphi(\epsilon_{i}^{[h]})}$
and make a vector $\mathbf{w}$ for $i=1,2,\cdots n$.

     * (3)

Update the mean element, $M_{R}^{(h+1)}=[\mathbf{w}^{(h)}]^{T}\mathbf{K}$.

     * (4)

Update error,
$\epsilon^{[h+1]}=(\rm{diag}(\mathbf{K})-2[\mathbf{w}^{(h)}]^{T}\mathbf{K}+[\mathbf{w}^{(h)}]^{T}\mathbf{K}[\mathbf{w}^{(h)}]^{T}\mathbf{1}_{n})^{\frac{1}{2}}$.

     * (5)

Update $h$ as $h+1$.

Output: the centered robust kernel matrix,
$\tilde{\mathbf{K}}_{R}=\mathbf{H}\mathbf{K}\mathbf{H}^{T}$ where
$\mathbf{H}=\mathbf{I}_{n}-\mathbf{1}_{n}\mathbf{w}^{T}$

Algorithm 1 The algorithm of estimating centered kernel matrix using robust
kernel mean element.

###  2.5 Robust kernel machine regression

Robustness theory helps us to understand the behavior of statistical genomic
procedures in real-life situations of genomic data without imposing
assumptions on the data. While many researchers have been studying the
robustness issue in a machine learning setting (e.g., support vector machine
for classification and regression, kernel PCA, kernel CCA, etc.), a well-
founded robust machine learning method has yet to be proposed for multi-omics
[14, 15, 18]. The standard kernel machine regression associate the output
$y_{i}$ $(i=1,2,\cdots,n)$ with $(q-1)$ covariates
$X_{i}=[X_{i1},X_{i2},\cdots,X_{i(q-1)}]^{T}$ and $m$-modal datasets,
$\mathbf{M}_{i}^{(1)},\mathbf{M}_{i}^{(2)},\cdots,\mathbf{M}_{i}^{(m)}$ of $n$
independent identical distributed (IID) subjects that obey:

$\displaystyle
y_{i}=X_{i}^{T}\mathbf{\beta}+f(\mathbf{M}_{i}^{(1)},\cdots,\mathbf{M}_{i}^{(m)})+\epsilon_{i},$
(5)

where $X_{i}$ is a vector of $q$ covariates including intercept for the $i-$th
subject, $\beta$ is a vector of $q$ fixed effects, $f$ is an unknown function
on the product domain,
$\mathcal{M}=\mathcal{M}^{(1)}\otimes\mathcal{M}^{(2)}\otimes,\cdots,\otimes\mathcal{M}^{(m)}$
with $\mathbf{M}_{i}^{(\ell)}\in\mathcal{M}_{\ell},\ell=1,2,\cdots m$ and the
error, $\epsilon_{i}\sim\rm{NIID}(0,\sigma^{2})$. By the ANOVA decomposition,
$f$ can be decomposed into the main effects, pairwise interactions effects,
the interactions effects of the respective dataset, and so on. Similarly, we
can also decompose the functional space, RKHS [7, 8].

The kernel machine regression is not robust to contaminated data and fictional
data distribution. Several researchers have investigated the issue of
robustness for the support vector machine, [15, 18], kernel density estimation
kernel canonical correlation, kernel principal analysis [5]. There are no
well-founded robust kernel machine regression method has been proposed yet.
Herein, we will introduce a novel non-linear M-estimator-based approach,
“robust kernel machine regression," and apply it to identify composite effects
in Genomic and Multiple extra-Genomic Data of disease-related traits. To
achieve robustness, we will formulate an empirical optimization problem by
combining empirical optimization with the idea of Huber’s and Hampel’s
M-estimation model.

The key goal of this paper is to propose robust KMR and apply this to three
views (m=3): genome, epigenome, transcriptome, along with BMD information of
subjects. To that end, assume that we have $n$ IID subjects under
investigation; $y_{i}\,(i=1,2,\cdots n)$ is a quantitative phenotype for the
$i$-th subject. We also associate the clinical covariates (e.g., age, weight,
and height) with these five views. Under this setting, Eq. (5) becomes:

$\displaystyle
y_{i}=X_{i}^{T}\beta+f(\mathbf{M}_{i}^{(1)},\mathbf{M}_{i}^{(2)},\mathbf{M}_{i}^{(3)})+\epsilon_{i}$
(6)

and

$f(\mathbf{M}^{(1)}_{i},\mathbf{M}^{(2)}_{i},\mathbf{M}^{(3)}_{i})=h_{\mathbf{M}^{(1)}}(\mathbf{M}_{i}^{(1)})+\cdots+\\\
h_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}}(\mathbf{M}_{i}^{(1)},\mathbf{M}_{i}^{(2)})+\cdots+\\\
h_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}}(\mathbf{M}_{i}^{(1)},\mathbf{M}_{i}^{(2)},\mathbf{M}_{i}^{(3)}),$
(7)

Using the representer theorem [34, 46] and the fact that the reproducing
kernel of a product of an RKHS is the product of the reproducing kernels [9],
the expanded functions of $f$ for arbitrary
$\tilde{\mathbf{M}}^{(\ell)}\in\mathcal{M}^{(ell)}$,
$\ell,\zeta\in\\{1,2,3\\}$ can be written as:

$h_{\mathbf{M}^{(\ell)}}=\sum_{i=1}^{n}\alpha^{(1)}_{i}k^{(1)}(\tilde{\mathbf{M}}^{(\ell)},\mathbf{M}_{i}^{(\ell)}),$

$h_{\mathbf{M}^{(\ell)}\times\mathbf{M}^{(\zeta)}}=\sum_{i=1}^{n}\alpha^{(1\times
2)}_{i}k^{(\ell)}(\tilde{\mathbf{M}}^{(\ell)},\mathbf{M}_{i}^{(\ell)})k^{(\zeta)}(\tilde{\mathbf{M}}^{(\zeta)},\mathbf{M}_{i}^{(\zeta)}),$

and

$h_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}}=$
$\sum_{i=1}^{n}\alpha^{(1\times 2\times
3)}_{i}k^{(1)}(\tilde{\mathbf{M}}^{(1)},\mathbf{M}_{i}^{(1)})k^{(2)}(\tilde{\mathbf{M}}^{(2)},\mathbf{M}_{i}^{(2)})k^{(3)}(\tilde{\mathbf{M}}^{(3)},\mathbf{M}_{i}^{(3)}).$

We can define the robust kernel Gram matrix for all data set:
$\mathbf{K}^{(1)}=(k^{(1)}(M_{i}^{1},\mathbf{M}_{j}^{1}))_{ij}$,
$\mathbf{K}^{(2)}=(k^{(2)}(M_{i}^{2},\mathbf{M}_{j}^{2}))_{ij}$,
$\mathbf{K}^{(3)}=(k^{(3)}(M_{i}^{3},\mathbf{M}_{j}^{3}))_{ij}$,
$\mathbf{K}^{(1\times 2)}=\mathbf{K}^{(1)}\odot\mathbf{K}^{(2)}$,
$\mathbf{K}^{(1\times 3)}=\mathbf{K}^{(1)}\odot\mathbf{K}^{(3)}$,
$\mathbf{K}^{(2\times 3)}=\mathbf{K}^{(2)}\odot\mathbf{K}^{(3)}$ and
$\mathbf{K}^{(1\times 2\times
3)}=\mathbf{K}^{(1)}\odot\mathbf{K}^{(2)}\odot\mathbf{K}^{(3)}$, where $\odot$
denotes the element-wise product of two matrices. Now we have

$\mathbf{h}_{\mathbf{M}^{(\ell)}}=\mathbf{K}^{(\ell)}\alpha^{(\ell)}$
$\mathbf{h}_{\mathbf{M}^{(\ell)}\times\mathbf{M}^{(\zeta)}}=\mathbf{K}^{(\ell\times\zeta)}\alpha^{(\ell\times\zeta)},$
$\mathbf{h}_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}}=\mathbf{K}^{(1\times
2\times 3)}\alpha^{(1\times 2\times\times 3)},$

where
$\alpha^{(\ell)}=[\alpha^{(\ell)}_{1},\alpha^{(\ell)}_{2},\cdots,\alpha^{(\ell)}_{n}]^{T},$
$\alpha^{(\ell\times\zeta)}=[\alpha^{(\ell\times\zeta)}_{1},\alpha^{(\ell\times\zeta)}_{2},\cdots,\alpha^{(\ell\times\zeta)}_{n}]^{T},$
and $\alpha^{(1\times 2\times 3)}=[\alpha^{(1\times 2\times
3)}_{1},\alpha^{(1\times 2\times 3)}_{2},\cdots,\alpha^{(1\times 2\times
3)}_{n}]^{T}$. The previous studies established that the first-order linear
system is equivalent to the normal equation of the linear mixed-effects model
(e.g., [38, 37, 24, 7, 8]):

$\displaystyle\mathbf{y}=\mathbf{X}\beta+\mathbf{h}_{\mathbf{M}^{(1)}}+\mathbf{h}_{\mathbf{M}^{(2)}}+\mathbf{h}_{\mathbf{M}^{(3)}}+\mathbf{h}_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}}+\mathbf{h}_{\mathbf{M}^{(1)}\times\mathbf{M}^{(3)}}$
$\displaystyle+\mathbf{h}_{\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}}+\mathbf{h}_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}}+\epsilon,$
(8)

where $\beta$ is the coefficient vector of fixed effects,
$\mathbf{h}_{\mathbf{M}^{(1)}}$, $\mathbf{h}_{\mathbf{M}^{(2)}}$,
$\mathbf{h}_{\mathbf{M}^{(3)}}$,
$\mathbf{h}_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}}$,
$\mathbf{h}_{\mathbf{M}^{(1)}\times\mathbf{M}^{(3)}}$,
$\mathbf{h}_{\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}}$ and
$\mathbf{h}_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}}$
are independent random effects with distribution as
$\mathbf{h}_{\mathbf{M}^{(1)}}\sim
N(0,\tau^{(1)}\mathbf{K}^{(1)}),\tau^{(1)}=\frac{\sigma^{2}}{\lambda^{(1)}}$,
$\mathbf{h}_{\mathbf{M}^{(2)}}\sim
N(0,\tau^{(2)}\mathbf{K}^{(2)}),\tau^{(2)}=\frac{\sigma^{2}}{\lambda^{(2)}}$,
$\mathbf{h}_{\mathbf{M}^{(3)}}\sim
N(0,\tau^{(3)}\mathbf{K}^{(3)}),\tau^{(3)}=\frac{\sigma^{2}}{\lambda^{(3)}}$,
$\mathbf{h}_{M^{(1\times 2)}}\sim N(0,\tau^{(1\times 2)}\mathbf{K}^{(1\times
2)}),\tau^{(1\times 2)}=\frac{\sigma^{2}}{\lambda^{(1\times 2)}}$,
$\mathbf{h}_{M^{(1\times 3)}}\sim N(0,\tau^{(1\times 3)}\mathbf{K}^{(1\times
3)}),\tau^{(1\times 3)}=\frac{\sigma^{2}}{\lambda^{(1\times 3)}}$,
$\mathbf{h}_{M^{(2\times 3)}}\sim N(0,\tau^{(2\times 3)}\mathbf{K}^{(2\times
3)}),\tau^{(2\times 3)}=\frac{\sigma^{2}}{\lambda^{(2\times 3)}}$,
$\mathbf{h}_{M^{(1\times 2\times 3)}}\sim N(0,\tau^{(1\times 2\times
3)}\mathbf{K}^{(1\times 2\times 3)}),\tau^{(1\times 2\times
3)}=\frac{\sigma^{2}}{\lambda^{(1\times 2\times 3)}}$. $\epsilon$ is also an
independent random variable with the distribution $\epsilon\sim
N(0,\sigma^{2}\mathbf{I})$, where $\mathbf{I}$ is an identity matrix. The
effects obtained by minimizing the loss function is the same as the best
linear unbiased predictors (BLUPs) of the linear mixed-effects model in Eq.
(2.5).

By restricted maximum likelihood (ReML) approach we can estimate the variance
components [7, 8]. The solution gives the coefficients of the fixed effect,
$\beta$, and the random effect, $\alpha$.

###  2.6 Statistical testing

This section addresses the test statistic of the overall effect, marginal
effects, interaction effects, and composite effects.

####  2.6.1 Overall testing

According to our model, the testing overall effect

$H_{0}:h_{\mathbf{M}^{(1)}}=\cdots=h_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}}=\cdots=\cdots$
$=h_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}\times\mathbf{M}^{(4)}\times\mathbf{M}^{(5)}}=0$

is equivalent to test the variance components in Eq.(2.5),

$H_{0}:\tau^{(1)}=\tau^{(2)}=\tau^{(3)}=\tau^{(1\times 2)}=\tau^{(1\times
3)}=\tau^{2\times 3}=$ $\cdots=\tau^{(1\times 2\times 3\times 5\times 5)}=0.$

We know that kernel matrices are not block-diagonal. The parameter in variance
component analysis is placed on the edge of the parameter space when the null
hypothesis is true. Since the asymptotic distribution of a likelihood ratio
test (LRT) statistic in favor of the null hypothesis is neither a chi-square
distribution nor a mixture chi-square distribution, we can use a score test
statistic on the restricted likelihood [8, 38]. The score test statistic is
defined as

$\displaystyle
S(\sigma_{0}^{2})=\frac{1}{2\sigma^{2}_{0}}(\mathbf{y}-\mathbf{X}\hat{\beta})^{T}\mathbf{K}(\mathbf{y}-\mathbf{X}\hat{\beta}),$
(9)

where $\mathbf{K}=\mathbf{K}^{(1)}+\cdots+\mathbf{K}^{(1\times
2)}+\cdots+\mathbf{K}^{(1\times 2\times 3)}$ and $\hat{\beta}$ is the maximum
likelihood estimator (MLE) of the regression coefficients. Under the null
model, $\mathbf{y}=\mathbf{X}\beta+\epsilon_{0}$, $\sigma_{0}^{2}$ is the
variance of $\epsilon_{0}$ and the quadratic function $S(\sigma_{0}^{2})$ of
the variable $\mathbf{y}$ follows a weighted mixture of the chi-square
distribution. Using the Satterthwaite method, we are able to approximate the
distribution of $S(\sigma_{0}^{2})$ to a scaled chi-square distribution,
($S(\sigma_{0}^{2})\sim\gamma\chi^{2}_{\nu}$). For estimating the the scale
parameter $\gamma$ and the degrees of freedom $\nu$, we use the method of
moments on the mean ($\gamma\nu$) and variance ($2\gamma^{2}\nu$) of the test
statistic that obey:
$\hat{\gamma}=\frac{\rm{Var}[S(\sigma_{0}^{2})]}{2\rm{E}[S(\sigma_{0}^{2})]}$
and
$\hat{\nu}=\frac{2\rm{E}[S(\sigma_{0}^{2})^{2}]}{\rm{Var}[S(\sigma_{0}^{2})]}$.
After all, we use the scaled chi-square distribution
$\hat{\gamma}\chi^{2}_{\hat{\nu}}$ to compute the $p-$ value of the score
statistic $S(\hat{\sigma}_{0}^{2})$.

####  2.6.2 Testing composite effects

Unlike an interaction hypothesis testing, a composite hypothesis testing is
assumed that all lower order effects are statistically significant. Likewise
the over all testing, to test the 3rd order composite effect, testing the null
hypothesis
$H_{0}:h_{\mathbf{M}^{(1)}\times\mathbf{M}^{(2)}\times\mathbf{M}^{(3)}\times\mathbf{M}^{(4)}\times\mathbf{M}^{(5)}}(\cdot)=0$
is equivalent to testing the variance component: $H_{0}:\tau^{1\times 2\times
3\times 4\times 5}=0$. Let
$\Sigma=\sigma^{2}\mathbf{I}+\tau^{(1)}\mathbf{K}^{(1)}+\cdots+\tau^{1\times
2}\mathbf{K}^{(1\times 2)}+\cdots+\tau^{2\times 3\times 4\times
5}\mathbf{K}^{(2\times 3\times 4\times 5)}$, and all $\tau$, and $\sigma^{2}$
are model parameters in favour of the null model. The test statistic is
defined as:

$\displaystyle
S_{I}(\tilde{\theta})=\frac{1}{2\sigma^{2}_{0}}\mathbf{y}^{T}\mathbf{B}_{I}\mathbf{K}^{(1\times
2\times 3\times 4\times 5)}\mathbf{B}_{I}\mathbf{y},$ (10)

where
$\tilde{\theta}=(\sigma^{2},\tau^{(1)},\tau^{(2)},\tau^{(3)},\tau^{(1\times
2)},\tau^{(2\times 3)})$, and
$\mathbf{B}_{I}=\Sigma^{-1}-\Sigma^{-1}\mathbf{X}(\mathbf{X}^{T}\Sigma^{-1}\mathbf{X})^{-1}\mathbf{X}^{T}\Sigma^{-1}$
is the projection matrix under the null hypothesis.

We apply the Satterthwaite method to approximate the distribution of higher
order composite test statistic $S_{I}(\tilde{\theta})$ by a scaled chi-square
distribution with scaled $\gamma_{I}$ and degree of freedom $\nu_{I}$ i.e.,
$S(\tilde{\theta})\sim\gamma_{I}\chi^{2}_{\nu_{I}}$. Using MOM, we can compute
the scaled parameter and degree of freedom,
$\hat{\gamma}_{I}=\frac{\rm{Var}[S_{I}(\tilde{\theta})]}{2\rm{E}[S_{I}(\tilde{\theta})]}$
and
$\hat{\nu}_{I}=\frac{2\rm{E}[S_{I}(\tilde{\theta})]}{\rm{Var}[S_{I}(\tilde{\theta})]}$,
receptively. In practice, the unknown model parameters are estimated by their
respective ReML estimates favoring the null model. Finally, the $p-$ value of
an observed higher-order composite effect test score statistic
$S_{I}(\tilde{\theta})$ is computed by the scaled chi-square distribution
$\hat{\gamma}_{I}\chi^{2}_{\hat{\nu_{I}}}$.

##  3 Experiments

In this section, we analyze the performance of RobKAM on synthetic data and a
multi-omics dataset from osteoporosis studies. We compare and evaluate the
performance of RobKAM against the performance of standard baselines kernel
machine approach, principal component SKAT, and principal component regression
(pPCAR). For the genome data we use the IBS kernel [25, 7] and for all others
dataset we consider the Gaussian kernel (the median of the pairwise distance
as the bandwidth [25, 3]). For Fisher’s scoring algorithm (the ReML
algorithm), we follow the parameters setting as in [8] to optimize the
proposed and standard kernel machine approaches. We conduct a set of initial
points in (0, 1) and pick the point which maximized the ReML algorithm to
overcome the potential of a local minima.

###  3.1 Simulation studies

In simulation studies, we consider a similar setting as in previous work [7,
8] and simulate data in different values of three parameters
$(\alpha_{1},\alpha_{2},\alpha_{3})$ for evaluating the performance of the
test. For example, $\alpha_{1}=\alpha_{2}=\alpha_{3}=0$ means that all effects
have vanished. Then, we analyze the false positive rate for the score test in
favor of the overall impact. For the main effects and 2nd order interaction
effects (but no higher-order composite effect), we report evaluating the power
of the score test. We repeat $1000$ simulations for each parameters setting to
get consistent results.

To evaluate the power of the composite hypothesis test, we compute the
composite score test in different parameter settings. We state the power of
the higher-order composite score test of RobKM and state-of-the-arts methods
in Table 1 . In general, we have two observations. (i) the test’s false
positive rate for the higher-order composite effects score test is controlled
by fixing the minimal $p-$value threshold to $0.05$ as state-of-the-arts
methods ( the similar observation for the false positive rate) . (ii) by
considering the power analysis ($\alpha_{5}\geq 0$), we observed that the
proposed method performs better than other methods, and its power exceeds
$0.80$ (Table 1). However, we note that the state-of-the-art methods (pPCAR,
fPCAR and SKAT) can significantly overstate the false positive rates and lose
substantial statistical power.

We further visualize the receiver operating characteristic (ROC) for three
sample sizes, $n\in\\{100,500,1000\\}$. Figure 3 shows the visual of the the
ROC with related random and the same parameter values
($\alpha_{2}=\alpha_{3}$) but the linear parameter is fixed to $\alpha_{1}=1$.
We assign each number with a probability of $0.5$. Also, a random number is
uniformly distributed either in a range or at $0$. By taking a step size
$0.0001$, we plot the sensitivity against (1- specificity) for each $p$-value
in the range of $0-1$. In all scenarios, the power gain of RobRKM and the
standard kernel machine approach relative to the alternative ones is apparent.
This result support that RobRKM has a similar power as the standard kernel
Machine approach.

Table 1: The power of higher-order composite score test of the proposed
approach (RobKAM), and state-of-the-art methods, the stamdard kernel machine
(sKMA), using dimension reduction regression (pPCAR) and sequence kernel
association test (SKAT)

Parameters | Simulation. n=300
---|---
| RobKMR | State-of-the-art methods
|  | sKMA | pPCAR | SKAT
($\alpha_{1}$, $\alpha_{2}$, $\alpha_{3}$) |  |  |  |
(0.1, 0, 0) | $0.061$ | $0.057$ | $0.051$ | $0.044$
(0, 0, 0.1) | $0.852$ | $0.842$ | $0.056$ | $0.044$
(0, 0, 0.5) | $0.863$ | $0.846$ | $0.051$ | $0.042$
(0, 0, 1) | $0.812$ | $0.824$ | $0.052$ | $0.041$
(1, 1, 0.1) | $0.840$ | $0.851$ | $0.048$ | $0.044$
(1, 1, 0.5) | $0.836$ | $0.862$ | $0.051$ | $0.041$
(1, 1, 1) | $0.844$ | $0.859$ | $0.049$ | $0.042$

Figure 3: The receiver operating characteristics of the proposed approaches
(RobKMA), the standard kernel machine approach (stKMA), dimension reduction
regression (pPCAR, fPCAR), and sequence kernel association test (SKAT) with
four sample sizes, $n\in\\{100,300,500,1000\\}$ for all interrelated
parameters ($\alpha_{2}=\alpha_{3}$) values are random but the linear
parameter is fixed to $\alpha_{1}=1$.

###  3.2  Osteoporosis data analysis

We apply the proposed method to our generated multi-omics dataset from
Louisiana osteoporosis studies data (as stated in Section $3.1.2$) [58, 42].
The dataset is available at our lab, and some of the data has already been
deposited in dbGaP (phs001960.v1.p1). This dataset included genome (3,997,535
SNPs which annotated to 25, 442 genes), epigenome (46,690 CpG methylations
which annotated to 4,676 genes.), and transcriptome (22,682 genes expression
profiles), from 57 Caucasian females with high BMD and $51$ with low BMD.

We fix each feature (gene) of the genome, epigenome, and transcriptome data as
an individual testing unit. We reduced the dimensionality of all these
datasets to make doable. To that end, we apply different methods (the t-test,
canonical correlation analysis based gene shaving (CCAOut), kernel canonical
correlation analysis based gene shaving (KCCAOut), and Linear Models for
Microarray-based gene shaving (LIMMA)) methods for three datasets [44, 8].
Finally, for our experiment, we consider $100$, $53$, and $61$ genes for the
genome, epigenome, and transcriptome data, respectively. Hence, we have
$323300$ ($100\times 53\times 61)$ triplets to test the overall and higher-
order composite effect. The overall test of the RobKMR approach and stKMR
approach provides us with $102973$ and $75525$ significant triples ($p\leq
0.05$), respectively. In Figure 4, we exhibit the plot of $-\rm{log}_{10}(p)$
for the both RobKMR and stKMA approaches. The vertical solid, dotted, and
double dotted lines correspond to the p-values of $0.05$, $0.01$, and $0.001$,
respectively. Table 3 shows the number of significant genes selected at
different p-values by the stKMA and the RobKMR methods. This table clearly
shows that the proposed method (RobKMR) can identify a small set of genes at
all p-values.

Table 2 presents the ReML estimates of all parameters, $\sigma^{2}$,
$\tau^{(1)}$, $\tau^{(2)}$, $\tau^{(3)}$, $\tau^{(1\times 2)}$,
$\tau^{(1\times 3)}$, $\tau^{(2\times 3)}$, $\tau^{(1\times 2\times 3)}$ and
the $p$-values for both the proposed and SKAT methods for each of the $10$
triplets. By the proposed method, these $10$ triplets were identified to have
significant interactions at a level of $p\leq 0.0000624$. At this $p$-value,
we observe that the unique $8$ genes (DKK1, MTND5, WNT3, MPP7, ANAPC1, FUBP3,
YWHAE, LRP5) $5$ genes (SMTN, FASTKD2, COG3, DNMBP, NMBP ), and $4$ genes
(DRGX, CSMD3, SOX1, USP17L1), are selected from genome, epigenome and
transcriptome data, respectively. For the stKMA, at P-value $0.00001$ and
$0.000001$, the selected genes for three datasets are GPLD1, LINC00461, PANK1,
SOST, ATP2B2,, KNDC1, NUP214, CAPN3 PRPF38B and (LINC00461, SOST, KNDC1, CAPN3
PRPF38B, respectively. On the other hand, for the RobKMR, at same p-values,
the selected genes are DDK1, MTND5, COG3, FASTKD2, SMTN, CSMD3, DRGX and DDK1,
MTND5, FASTKD2, SMTN, CSMD3, DRGX, respectively. At these low p-values, we
observe the different set of genes for the stKMA and RobKMR methods. In
addition, a list of genes (genome, epigenome, and transcriptome data) for the
RobKMR is tabulated in Table 3 at a p-value of $0.001$. Figure 5 presents Venn
diagrams of the selected genes of each three datasets (genome, epigenome, and
transcriptome) using the proposed and other three methods. This figure also
shows the proposed method is able to identify a least set of genes than other
methods.

Figure 4:  Manhattan plots of negative log of the p-values $-log_{10}(p)$ in
the y-axis against triplets in the x-axis based on all significant overall
score tests of the standard (a) and the proposed robust (b) methods. Table 2:
The selected significant genes using the proposed method and state-of-the-art
methods. The $p-$value was set to be $0.0000624$ for top 10 triples (OV:
Overall effect and HOC: Higher-order composite effect)

|  |  | RobKMR | State-of-the-Art Methods
---|---|---|---|---
Genome | Epigenome | Transcriptome | $\sigma^{2}$ | $\tau^{(1)}$ | $\tau^{(2)}$ | $\tau^{(3)}$ | $\tau^{1\times 2}$ | $\tau^{1\times 3}$ | $\tau^{2\times 3}$ | $\tau^{1\times 2\times 3}$ | OV | HOC | StKMA | SKAT | pPCAR
DKK1 | SMTN | DRGX | $0.0179$ | $4e-04$ | $0.0019$ | $0.0073$ | $0.0044$ | $0.0001$ | $0.0383$ | $0.0000$ | $0.0911$ | $2.71e-10$ | $1.0000$ | $0.00012$ | $0.0393$
MTND5 | FASTKD2 | CSMD3 | $0.0343$ | $0e+00$ | $0.0001$ | $0.0001$ | $0.0001$ | $0.0001$ | $0.0000$ | $0.0100$ | $0.0433$ | $2.45e-07$ | $0.0004$ | $0.9853$ | $0.2354$
MTND5 | COG3 | CSMD3 | $0.0342$ | $5e-04$ | $0.0010$ | $0.0009$ | $0.0010$ | $0.0012$ | $0.0000$ | $0.0100$ | $0.0167$ | $8.67e-06$ | $1.0000$ | $0.0003$ | $0.1903$
WNT3 | DNMBP | SOX1 | $0.0343$ | $1e-03$ | $0.0010$ | $0.0010$ | $0.0010$ | $0.0010$ | $0.0000$ | $0.0000$ | $0.0000$ | $1.63e-05$ | $1.0000$ | $0.0124$ | $0.3547$
MPP7 | DNMBP | DRGX | $0.0340$ | $1e-03$ | $0.0010$ | $0.0010$ | $0.0005$ | $0.0010$ | $0.0000$ | $0.0000$ | $0.0000$ | $1.96e-05$ | $1.0000$ | $0.8562$ | $0.1529$
MPP7 | DNMBP | SOX1 | $0.0343$ | $1e-03$ | $0.0010$ | $0.0010$ | $0.0009$ | $0.0010$ | $0.0000$ | $0.0000$ | $0.0000$ | $2.10e-05$ | $1.0000$ | $0.6842$ | $0.0265$
ANAPC1 | DNMBP | USP17L1 | $0.0342$ | $1e-03$ | $0.0010$ | $0.0012$ | $0.0010$ | $0.0000$ | $0.0004$ | $0.0100$ | $0.0000$ | $5.31e-05$ | $1.0000$ | $1.0000$ | $0.2564$
FUBP3 | NMBP | SOX1 | $0.0343$ | $1e-03$ | $0.0010$ | $0.0010$ | $0.0010$ | $0.0010$ | $0.0000$ | $0.0000$ | $0.0001$ | $5.97e-05$ | $1.0000$ | $0.0125$ | $0.3602$
YWHAE | DNMBP | SOX1 | $0.0343$ | $1e-03$ | $0.0010$ | $0.0010$ | $0.0010$ | $0.0010$ | $0.0000$ | $0.0000$ | $0.0001$ | $6.24e-06$ | $0.0235$ | $0.0285$ | $0.4369$
LRP5 | DNMBP | DRGX | $0.0341$ | $8e-04$ | $0.0010$ | $0.0011$ | $0.0000$ | $0.0010$ | $0.0000$ | $0.0032$ | $0.0001$ | $6.54e-05$ | $0.0091$ | $0.0002$ | $0.3602$

Table 3: The number of significant genes selected at differect p-values by the
stKMA and RobKMR methods

| KMRHCE | RKMRHCE
---|---|---
P-values | Genome | Epigenome | Transcriptome | Genome | Epigenome | Transcriptome
$0.05$ | $98$ | $45$ | $36$ | $90$ | $40$ | $44$
$0.01$ | $90$ | $30$ | $21$ | $58$ | $23$ | $26$
$0.001$ | $45$ | $20$ | $13$ | $25$ | $7$ | $5$
$0.0001$ | $18$ | $6$ | $4$ | $11$ | $4$ | $4$
$0.00001$ | $4$ | $3$ | $2$ | $2$ | $3$ | $2$
$0.000001$ | $2$ | $1$ | $2$ | $2$ | $2$ | $2$

Figure 5:  Venn diagrams of the selected genes of each three datasets: (a)
genome, (b) epigenome, and (c) transcriptome using the proposed and other
three methods. Table 4: The selected significant genes by the RobKMA method at
p-value $=0.001$. The highlighted genes identified by RobKMA only.

Method |  |  |  | Genome |  |  |  |  |
---|---|---|---|---|---|---|---|---|---
| ANAPC1 | CLCN7 | CPN1 | CSF1 | DKK1 | FUBP3 | GPATCH1 | GPLD1 | INSIG2
Genome | JAG1 | LIN7C | LRP5 | MBL2 | MPP7 | MTND5P19 | NIPAL1 | NTAN1 | RGCC
| RUNX2 | SIDT1 | SLC25A13 | SOST | SOX6 | WNT3 | YWHAE |  |
|  |  |  |  |  |  |  |  |
Epigenome | COG3 | DNER | DNMBP | FASTKD2 | PLD5 | RBM38 | SMTN |  |
Transcriptome | CSMD3 | DRGX | SOX1 | UCKL1 | USP17L1 |  |  |  |

We also construct functional protein association networks for available
enrichment analysis using STRING 222https://string-db.org/. Figure 6 shows the
gene-gene networks based on the protein interactions among the selected genes
of all three datasets. ($37=25+07+05$, at p-values $\leq$ 0.001) by the
RobKMR. The network analysis demonstrates that the number of nodes, edges,
expected edges, average node degree, clustering coefficient, protein-protein
interaction enrichment p-values are $36$, $62$, $32$, $3.44$, $0.518$ and
$0.0000012$, respectively. In addition, this network has $2$ GO-terms, $525$
publications, $5$ pathways, and $6$ diseases significantly enriched. In this
figure, the color saturation of the edges represents the confidence score of a
functional association. Therefore, this network analysis confirms that the
selected genes Have significant interactions than expected. It also indicates
that the most selected genes may function collaboratively.

We access the biomedical and genomic information using the Human gene database
(GeneCards) 333 https://www.genecards.org/ and the GeneMANIA 444
https://genemania.org/  to confirm the biological roles of the selected genes
[22]. As we know, the critical goal of the GeneCards database is the
unequivocal identification of enhancer elements and uncovering their
connections to genes for understanding gene regulation and molecular pathways.
On the other hand, GeneMANIA finds other genes related to a gene or a set of
input genes using an extensive set of functional association data. This
organized data includes protein and genetic interactions, pathways, co-
expression, co-localization, and protein domain similarity. We used the
GeneCards database to provide insight into the gene regulatory elements
(promoters and enhancers) for $7$ selected genes (DKK1, MTND5, COG3, FASTKD2,
SMTN, CSMD3, and DRGX) at p-value $\leq 0.00001$ by the RobKMR. Table 5 shows
GeneHancer identifier, GeneHancer score, gene association score, total score,
major-related diseases, and PubMed database. This table shows that the
selected genes have a remarkable GeneHancer score, gene association score,
total score, and literature review in the past studies. According to the
disease annotation, the selected 7 genes are highly associated with complex
diseases, including the higher risk of developing osteoporosis.

As a comparison between the stKMA and RobKMR, stKMA extract $78(45+20+13)$
significant genes, while the RobKMR extract $37(25+7+5)$ significant genes at
a p-values $\leq 0.00001$. The RobKMR extracts $20$ unique genes out of 37 (i.
e., 17 genes are common for both methods). We intend to show the performance
that the RobKMR could find undiscovered genes and could exclude not
significant genes (in the meaning of robustness). To that end, we also conduct
the network analysis of genes related to an input gene or a set of input genes
using an extensive set of functional association data with the GeneMANIA.
Table 6 presents the network analysis of each gene along with the number of
edges, average node degree, average local clustering coefficient, expected
number of edges, and interaction enrichment p-values. The table report that
the RobKMR method has discovered different significant genes. To show the
performance of these genes, we also consider $7$ genes for an extensive set of
functional association data with the GeneMANIA. Figure 7 shows the genes
related networks to the input genes COG3, SMTN, DRGX, MTMD5, CSMD3, DKK1 and
FSTKD2 (at p-value $\leq 0.00001$). This figure suspects that the selected
genes have strong physical interactions, co-expression, predicted, co-
localization, genetic interactions, pathway, and shared protein domains
networks. Thus, the proposed robust methods (RobKMR) can find undiscovered
genes in addition to significant gene triplets.

Table 5: The GeneHancer identifier, GeneHancer score, gene association score,
total score, major-related diseases, and PubMed database (PMID) of 6 selected
genes

Gene | GeneHancer | GeneHancer | Association | Total | Transcription Factor | Related
---|---|---|---|---|---|---
ID | ID | Score | Score | Score | Binding Sites (TFs) | Disease
DKK1 | GH10J052312 | $2.0$ | $5256.8$ | $520.05$ | $140$ | bone formation and bone disease, Osteoporosis, cancer (Leukemia) and Alzheimer disease
MTND5 | GH21J043724 | $2.4$ | $9.7$ | $23.76$ | $340$ | brain and muscles
COG3 | GH20J056524 | $0.4$ | $500.7$ | $221.23$ |  $(2)$: $25367360$; $25130324$ | Bone marrow; thyroid
FASTKD2 | GGH02J206764 | $1.7$ | $257.7$ | $447.37$ | $179$ | Mitochondrial disease
SMTN | GH22J031077 | $2$ | $271.3$ | $549.79$ | $184$ | actin binding and structural constituent of muscle, Glomuvenous Malformations, and Viral Gastritis
CSMD3 |  GH08J113434 | $1.7$ | $254.5$ | $424.63$ | $38$ | Familial Adult Myoclonic Epileps and Trichorhinophalangeal Syndrome
DRGX |  GH10J049396 | $1$ | $5258.5$ | $270.84$ | $6$ | Chromosome 16P12.2-P11.2 Deletion Syndrome, 7.1- To 8.7-Mb

Figure 6: The functional protein association networks of the selected genes of
genome, epigenome, and transcriptome datasets by the proposed method (at
p-values $leq0.001$). Table 6: Network analysis of selected $20$ unique genes
by the proposed method at p-value $\leq$ $0.001$

Gene | Number of | Average node | Average local | Expected number | Interaction enrichment
---|---|---|---|---|---
| edges | degree | clustering coefficient | of edges | p-value
|  |  | Genome |  |
CLCN7 | $29$ | $5.27$ | $0.84$ | $10$ | $1.06e-06$
DKK1 | $55$ | $10.00$ | $1$ | $11$ | $1.0e-16$
INSIG2 | $30$ | $5.45$ | $0.775$ | $11$ | $7.37e-07$
JAG1 | $40$ | $7.27$ | $0.882$ | $11$ | $3.42e-11$
LIN7C | $32$ | $5.82$ | $0.874$ | $10$ | $4.25e-08$
MBL2 | $55$ | $10.00$ | $1$ | $10$ | $1.0e-16$
NTAN1 | $24$ | $4.36$ | $0.88$ | $10$ | $0.000179$
RGCC | $40$ | $7.27$ | $0.856$ | $25$ | $0.003$
SIDT1 | $19$ | $3.45$ | $0.891$ | $10$ | $0.0102$
SLC25A13 | $54$ | $9.82$ | $0.982$ | $11$ | $1.0e-16$
WNT3 | $45$ | $8.18$ | $8.18$ | $10$ | $2.66e-15$
|  |  | Epigenome |  |
COG3 | $53$ | $9.64$ | $0.968$ | $11$ | $1.0e-16$
DNER | $55$ | $10.00$ | $1$ | $11$ | $1.0e-16$
DNMBP | $25$ | $4.55$ | $0.838$ | $11$ | $0.000262$
RBM38 | $22$ | $4.00$ | $0.77$ | $13$ | $0.0205$
SMTN | $29$ | $5.27$ | $0.913$ | $11$ | $5.21e-06$
|  |  | Transcriptome |  |
DRGX | $16$ | $2.91$ | $0.921$ | $10$ | $0.541$
SOX1 | $49$ | $8.91$ | $0.901$ | $12$ | $7.77e-16$
UCKL1 | $40$ | $7.27$ | $0.872$ | $11$ | $5.05e-12$

Figure 7: The network of predicted related genes of each unique selected genes
by the proposed method at p-value $\leq$ 0.00001.

We infer a causal relationship to explore further whether the selected genes
will have unique associations with the BMD. Hill climbing, a greedy search
engine, obtain optimal solutions for convex problems instead of local optima
for a causal relationship. This approach selects the best successor node under
the evaluation function and commits the search [49]. We apply a Hill-Climbing
approach among genome and epigenome datasets (the select genes at p-values
$\leq 0.00001$). The causal relationship of BMD with (a) epigenome and (b)
genome and epigenome datasets are illustrated in Figure 8. We observed that
these two genes, FASTKD2 ( with one methylation profile) and COG3 (with three
methylation profiles), are directly related to BMD. Thus, This observation
concluded that the selected biomarkers may have significant impact on BMD but
are not general.

Figure 8: The causal relationship of BMD drives with the epigenome only (a)
and the epigenome and genome (b) datasets. Two genes of the epigenome, FASTKD2
(one methylation profile) and COG3 (three methylation profiles) and one gene
of genome (7 SNPs) data, are causal related to BMD.

###  3.3  Drug repositioning for osteoporosis

We conducted a literature review of OP disease for exploring candidate drugs
that have host transcriptome-guided 35 meta-drug agents. Thus we considered 35
drug agents to explore candidate drugs by molecular docking with our selected
$3$ genes (proposed receptors proteins). To offer in-silico validated
efficient candidate drugs for the treatment against OP, we employed a
molecular docking study of our proposed receptor proteins with the drug
agents.

We downloaded the 3D structure of the 20 proteins (MBL2, DNMBP, YWHAE, RBM38,
SMTN, CPN1, JAG1, LIN7C, MPP7, DKK1, CSF1, SLC25A13, ANAPC1, MTND5, NTAN1,
WNT3, SOST, INSIG2, RUNX2, CLCN7) from Protein Data Bank (PDB) with source
codes 1hup, 1ug1, 2br9, 2cqd, 2d87, 2nsm, 2vj2, 3lra, 3o46, 3s2k, 3uez, 4p5w,
4ui9, 5xtc, 6a0e, 6ahy, 6l6r, 6m49, 6vg8, 7bxu, respectively [10]. On the
other hand, the 3D structure of 17 proteins (DRGX, USP17L17, SOX1, LRP5, SOX6,
GPLD1, NIPAL1, PLD5, DNER, FUBP3, COG3, GPATCH1, RGCC, UCKL1, SIDT1, FASTKD2,
CSMD3) were downloaded from AlphaFold source using UniProt ID of A6NNA5,
D6RBQ6, O00570, O75197, P35712, P80108, Q6NVV3, Q8N7P1, Q8NFT8, Q96I24,
Q96JB2, Q9BRR8, Q9H4X1, Q9NWZ5, Q9NXL6, Q9NYY8, Q7Z407 [33]. We downloaded the
3D structures of 35 drugs from the PubChem database [32]. Then molecular
docking was carried out between $37$ proteins and $35$ meta-drug agents to
calculate the binding affinity scores (kcal/mol) for each pair of proteins and
drugs. Then we organized the proteins in descending order of row sums of the
binding affinity matrix and drug agents according to the column sums of the
scoring matrix to select a few drug agents as the candidate drugs. Figure 9
presents the binding affinity matrix. Thus we set eight top-ranked drug agents
(Tacrolimus, Ibandronate, Alendronate, Bazedoxifene, Goserelin, Raloxifene,
Buserelin, Prednisolone) as candidate drugs with average binding affinity
scores $-7.5$ kcal/mol $\leq$ against the $37$ proteins.

The docked complexes of the top three virtual hits from AutoDock-Vina docking
are further considered for protein-ligand interaction profiling. As shown in
Figure 10(a), the SIDT1_Tacrolimus complex showed one hydrogen bond with
Tyr590 residues. Although the ligand formed significant (key) hydrophobic
interactions with Leu587, Tyr590, Phe717 residues, and Tyr724 deposition
showed additional electrostatic interactions with the drug. On the other hand,
MTND5_Ibandronate (Figure 10(b)) complex showed five hydrogen bonds with
Leu429, Thr432, Arg436, Asn505, Asn509 residues and the significant
hydrophobic interactions with Arg357, Leu429residues. In the case of the
DKK1_Alendronate complex, Alendronate formed six hydrogen bonds with Thr221,
His229, Arg236, Cys237, Tyr238, Cys239 residues (see Figure 10(b)). As shown
in Figure 10(b), the FASTKD2_Bazedoxifene complex showed five hydrogen bonds
with Leu61, Asn62, Glu568, His603, Asp605 residues, and the significant
hydrophobic interactions with Phe64, Leu637, Val639, Ala643, Phe654, Leu655,
Lys658 residues.

Figure 9: Molecular docking results computed with autodock vina. Image of
binding affinities based on the 35 ordered drugs in the X-axis and ordered
$37$ proteins (selected genes by the RobKMR at p-values $\leq 0.001$) in the
Y-axis. The strong and weak binding affinities between proteins and drugs are
denoted by red and green colors, respectively. Figure 10:  The 3D structure
(left) and 2D Schematic diagram (right) of hub protein of the top four
potential targets (DKK1, MTND5, FASTKD2, and SIDT1) and top four lead drugs.
We calculated drug interacting amino acid within 4 amstrong. Lead four drugs,
Tacrolimus, Ibandronate, Alendronate, and Bazedoxifene, are selected by
investigating the binding affinity score.

##  4 Concluding remarks

This paper developed a novel robust approach to identify inter-related risk
factors in multi-omics data for drug repurposing in OP. Benchmarking
experiments based on simulation and real datasets analysis have shown that our
proposed approach provides a competitive performance compared with the
existing ones to derive a statistic for testing the inter-related risk
factors. The RobKMR method’s power, biological validations, and drug
repurposing are further demonstrated by its application to synthesized and a
real multi-omics dataset in osteoporosis studies.

Simulation studies show that the test’s false positive rate and the power
analysis for higher-order composite effect are mitigated by fixing the nominal
p-value threshold along with other state-of- the-art-methods. The ROC curves
also present the power gain by the RobKMR and stKMA methods in all scenarios.

The real datasets analysis confirmed that the RobKMR could select the inter-
related risk factors of OP for drug repurposing. The network analysis shows
that the selected genes of each omics also have significantly more
interactions than expected. The gene’s function is collaborative and
biologically relevant to OP. We observe that the selected genes are directly
related to the BMD by a causal analysis. Additionally, the top four genes
DKK1, MTND5, FASTKD2 (at p-value $\leq$ 0.00001), and SIDT1 (at p-value $\leq$
0.001) identify four lead drugs: Tacrolimus, Ibandronate, Alendronate, and
Bazedoxifene from all $30$ experimented candidates for drug repurposing in OP.

While our proposed method can identify stable biomarkers (DKK1, MTND5,
FASTKD2, and SIDT1) for OP studies, we acknowledge that further investigation
of these four genes is essential to improve OP studies for its drug discovery.
Further, the proposed approach can be applied to any disease model where
multi-omics datasets are available.

### Acknowledgments

This work is benefited by the support of U19AG05537301 and R01AR069055.

## References

  * [1] Md Ashad Alam, V. Calhoun, and Y. P. Wang.  Influence function of multiple kernel canonical analysis to identify outliers in imaging genetics data.  Proceedings of 7th ACM Conference on Bioinformatics, Computational Biology, and Health Informatics (ACM BCB),Seattle, WA, USA, pages 210–2198, 2016.
  * [2] Md Ashad Alam, V. Calhoun, and Y. P. Wang.  Identifying outliers using multiple kernel canonical correlation analysis with application to imaging genetics.  Computational Statistics & Data Analysis, 125:70–85, 2018.
  * [3] Md Ashad Alam and K. Fukumizu.  Higher-order regularized kernel canonical correlation analysis.  International Journal of Pattern Recognition and Artificial Intelligence, 29(4):1551005(1–24), 2015.
  * [4] Md Ashad Alam, K. Fukumizu, and Y.-P. Wang.  Robust Kernel (Cross-) Covariance Operators in Reproducing Kernel Hilbert Space toward Kernel Methods.  ArXiv e-prints.
  * [5] Md Ashad Alam, Kenji Fukumizu, and Yu-Ping Wang.  Influence function and robust variant of kernel canonical correlation analysis.  Neurocomputing, 304:12–29, 2018.
  * [6] Md Ashad Alam, Osamu Komori, Hong-Wen Deng, Vince D Calhoun, and Yu-Ping Wang.  Robust kernel canonical correlation analysis to detect gene-gene co-associations: A case study in genetics.  Journal of bioinformatics and computational biology, 17 (4):1950028, 2019.
  * [7] Md Ashad Alam, Hui-Yi Lin, Hong-Wen Dengc, V. Calhoun, and Y. P. Wang.  A kernel machine method for detecting higher order interactions in multimodal datasets: Application to schizophrenia.  Journal of Neuroscience Methods, 309:161–174, 2018.
  * [8] Md Ashad Alam, Chuan Qiu, Hui Shen, Yu-Ping Wang, and Hong-Wen Deng.  A generalized kernel machine approach to identify higher-order composite effects in multi-view datasets, with application to adolescent brain development and osteoporosis.  Journal of Biomedical Informatics, 120:103854, 2021.
  * [9] N. Aronszajn.  Theory of reproducing kernels.  Transactions of the American Mathematical Society, 68:337–404, 1950\.
  * [10] H M Berman, J Westbrook, Z. Feng, G. Gilliland, and etc.  The protein data bank.  Nucleic Acids Research, 28(1):235–42, 2000.
  * [11] K. M. Borgwardt and et al.  Protein function prediction via graph kernels,.  Bioinformatics, 21:i47–i56, 2005.
  * [12] G. Camps-Valls, J. L. Rojo-Alvarex, and M. Martinez-Romon.  Kernel Methods in Bioengineering, Signal and Image.  Idea Group publishing, London, 2007.
  * [13] Kai Cheng, Zhenzhou Lu, Yuhao Wei, Yan Shi, and Yicheng Zhou.  Mixed kernel function support vector regression for global sensitivity analysis.  Mechanical Systems and Signal Processing, 96:201–2014, 2017.
  * [14] A. Christmann and I. Steinwart.  On robustness properties of convex risk minimization methods for pattern recognition.  Journal of Machine Learning Research, 5:1007–1034, 2004.
  * [15] A. Christmann and I. Steinwart.  Consistency and robustness of kernel-based regression in convex risk minimization.  Bernoulli, 13(3):799–819, 2007.
  * [16] Graeme R. Clark and Emma L. Duncan.  The genetics of osteoporosis.  British Medical Bulletin, 113:73–81, 2015.
  * [17] Pietro Coretto, Angela Serra, and Roberto Tagliaferri.  Robust clustering of noisy high-dimensional gene expression data for patients subtyping.  Bioinformatics, 34 (23):4064–4072, 2018.
  * [18] M. Debruyne, M. Hubert, and J.V. Horebeek.  Model selection in kernel based regression using the influence function.  Journal of Machine Learning Research, 9:2377–2400, 2008.
  * [19] Lei Du, Kefei Liu, Lei Zhu, Xiaohui Yao, and etc.  Identifying progressive imaging genetic patterns via multi-task sparse canonical correlation analysis: a longitudinal study of the adni cohort.  Bioinformatics, 35(14):i474–i483, 2019.
  * [20] Aurelien Dugourd, Christoph Kuppe, Marco Sciacovelli, and etc.  Causal integration of multi-omics data with prior knowledge to generate mechanistic hypotheses.  Molecular Systems Biology, 17(1):1–17, 2021.
  * [21] Ruzong Fan, Paul S. Albert, and Enrique F. Schisterman.  A discussion of gene-gene and gene- environment interactions and longitudinal genetics analysis of complex traits.  Statistical Medicine, 31(22):2565–2572, 2012.
  * [22] Simon Fishilevich, Ron Nudel, and et al.  Genehancer: genome-wide integration of enhancers and target genes in genecards.  Database, 2017:1–17, 2017.
  * [23] Vittorio Fortino, Pia Kinaret, Nanna Fyhrquist, Harri Alenius, and Dario Greco.  A robust and accurate method for feature selection and prioritization from multi-class omics data.  PLoS One, 9(9):e107801, 2014.
  * [24] T. Ge, T. E. Nichols, D. Ghoshd, E. C. Morminoe, am M. R. Sabuncu J. W.Smoller, and the Alzheimer’s Disease Neuroimaging Initiative.  A kernel machine method for detecting effects of interaction between multidimensional variable sets: An imaging genetics application.  NeuroImage, 109:505–514, 2015.
  * [25] A. Gretton, K. Fukumizu, C. H. Teo, L. Song, B. Schölkopf, and A. Smola.  A kernel statistical test of independence.  In Advances in Neural Information Processing Systems, 20:585–592, 2008.
  * [26] F. R. Hampel, P. J. Rousseeuw E. M. Ronchetti, and W. A. Stahel.  Robust Statistics: The Approach Based on Influence Functions.  John Wiley & Sons, New York, 2011.
  * [27] F. R. Hampel, E. M. Ronchetti, and W. A. Stahel.  Robust Statistics.  John Wiley & Sons, New York, 1986.
  * [28] Y. Hasin, M. Seldin, and A. Lusis.  Multi-omics approaches to disease.  Genome Biol, 18 (83):1:15, 2017.
  * [29] P. J. Huber and E. M. Ronchetti.  Robust Statistics.  John Wiley & Sons, England, 2009.
  * [30] and Srikant Verma Indhupriya Subramanian, Shiva Kumar, Abhay Jere, and Krishanpal Anamika.  Multi-omics data integration, interpretation, and its application.  Bioinformatics and Biology Insights, 14:1–24, 2020.
  * [31] J. Kim and C. D. Scott.  Robust kernel density estimation.  Journal of Machine Learning Research, 13:2529–2565, 2012.
  * [32] So Yeon Kim, Hyun-Hwan Jeong, Jaesik Kim, Jeong-Hyeon Moon, and Kyung-Ah Sohn.  Robust pathway-based multi-omics data integration using directed random walks for survival prediction in multiple cancer studies.  Biology Direct, 14(8):1–13, 2019.
  * [33] Sunghwan Kim, Jie Chen, Tiejun Cheng, Asta Gindulyte, Jia He, and etc.  Pubchem 2019 update: improved access to chemical data.  Nucleic Acids Research, 8:47(D1):D1102–D1109, 2019.
  * [34] G. Kimeldorf and G. Wahhba.  Some results on tchebycheffian spline functions.  Journal of mathematical analysis and applications, 33(1):82–95, 1971.
  * [35] S. Y. Kung.  Kernel Methods and Machine Learning.  Cambridge University Press, New York, 2014.
  * [36] G. R. G. Lanckriet, T. De Bie, N. Cristianini, M. I. Jordan, and W. S. Noble.  A statistical framework for genomic data fusion.  Bioinformatics, 20:2626–2635, 2004.
  * [37] S. Li and Y. Cui.  Gene-centric gene-gene interaction: a model-based kernel machine method.  The Annals of Applied Statistics, 6(3):1134–1161, 2012.
  * [38] D. Liu, X. Lin, and D. Ghosh.  Semiparametric regression of multidimensional genetics pathway data: least squares kernel machines and linear mixed model,.  Biometrics, 630(4):1079–1088, 2007.
  * [39] A. C. A. Nascimento, R. B. C. Prudêncio, and Ivan G. Costa.  A multiple kernel learning algorithm for drug-target interaction prediction.  BMC Genomics, 17:46:1–16, 2016.
  * [40] Daniel E. Neafsey, Aimee R. Taylor, and Bronwyn L. MacInnis.  Advances and opportunities in malaria population genomics.  Nature Reviews Genetics, 22(8):502–517, 2021.
  * [41] Q. N. Peng, J. Zhao, and F. Xue.  A gene-based method for detecting gene-gene co-association in a case-control association study.  European Journal of Human Genetics, 18:582–587, 2010.
  * [42] Chuan Qiu, Fangtang Yu, Kuanjui Su, Lan Zhang, and et al.  Multi-omics data integration for identifying osteoporosis biomarkers and their biological interaction and causal mechanisms.  iScience, 23(2):100847, 2020.
  * [43] Anthony Rios, Ramakanth Kavuluru, and Zhiyong Lu.  Generalizing biomedical relation classification with neural adversarial domain adaptation.  Bioinformatics, 34(17):2973–2981, 2018.
  * [44] Manuela G M Rocha-Braz and Bruno Ferraz de Souza.  Genetics of osteoporosis: searching for candidate genes for bone fragility.  Archives of Endocrinology and Metabolism, 60(40):391–401, 2016\.
  * [45] B. D. Moor S. Yu, L-C. Tranchevent and Y. Moreau.  Kernel-based Data Fusion for Machine Learning.  Springer, Verlag Berlin Heidelberg, 2011.
  * [46] B. Schölkopf and A. J. Smola.  Learning with Kernels.  MIT Press, Cambridge MA, 2002.
  * [47] T L Stewart and S H Ralston.  Role of genetic factors in the pathogenesis of osteoporosis.  Journal of Endocrinology, 166:235–245, 2000.
  * [48] Lei Tang.  Multiomics sequencing goes spatial.  Nature Methods, 18 (31):31–35, 2021.
  * [49] S. Taskinen, C. Croux, A. Kankainen, E. Ollila, and H. Oja.  Canonical analysis based on scatter matrices.  Journal of Multivariate Analysis, 97(2):359–384, 2006.
  * [50] J. W. Tukey.  Exploratory Data Analysis.  Addison-Wesley, Reading, Massachusetts, 1977.
  * [51] Janani Venugopalan, Li Tong, Hamid Reza Hassanzadeh, and May D. Wang.  Multimodal deep learning models for early detection ofalzheimer’s disease stage.  Scientific Reports, 11(1):1–13, 2021.
  * [52] Xiaoqian Wang, Kamran Ghasedi Dizaji, and Heng Huang.  Conditional generative adversarial network for gene expression inference.  Bioinformatics, 34:i603–i611, 2018.
  * [53] Xun Wang, Lihong Diao, Dezhi Sun, Dan Wang, and etc.  Osteoporos atlas: a human osteoporosis-related gene database.  PeerJ, 7:e6778, 2019.
  * [54] Shaoping Wanga, Ang Li, Kuangyu Wena, and Ximing Wub.  Robust kernels for kernel density estimation.  Economics Letters, 191:109138, 2020.
  * [55] Kai Wu, Yasuyuki Taki, Kazunori Sato, Haochen Qi, Ryuta Kawashima, and Hiroshi Fukuda.  A longitudinal study of structural brain network changes with normal aging.  Forntiers in Human Neuroscience, 7(APR 2013):1–12, 2013.
  * [56] C. Xu, D. Tao, and C. Xu.  A survey of multi-view machine learning.  Neural Computation and Applications, 23(7-8):2031–2038, 2013.
  * [57] K. K. Yan, H. Zhao, and H. Pang.  A comparison of graph- and kernel-based – omics data integration algorithms for classifying complex traits.  BMC Bioinformatics, 18:539:1–13, 2017.
  * [58] Fangtang Yu, Chuan Qiu, Chao Xu, , and et al.  Mendelian randomization identifies cpg methylation sites with mediation effects for genetic influences on bmd in peripheral blood monocytes.  frontiers in Genetics, 11(60):1–14, 2020.
  * [59] Wen Zhong, Fredrik Edfors, Anders Gummesson, Göran Bergström, Linn Fagerberg, and Mathias Uhlen.  Next generation plasma proteome profiling to monitor health and disease.  Nature Communications, 12 (2493):4064–4072, 2021.