Briefings in Bioinformatics, 22(5), 2021, 1–11
https://doi.org/10.1093/bib/bbaa406
Problem Solving Protocol

via domain adaptation
Menglan Cai, Canh Hao Nguyen, Hiroshi Mamitsuka and Limin Li
Corresponding author: Limin Li, Xi’an Jiaotong University, 28 Xianning W Rd, Jiaoda Commerce Block, Beilin, Xi’an 710049, China.
Tel.: +86 029 82660967; E-mail: liminli@mail.xjtu.edu.cn

Abstract
Motivation: Gene set enrichment analysis (GSEA) has been widely used to identify gene sets with statistically significant
difference between cases and controls against a large gene set. GSEA needs both phenotype labels and expression of genes.
However, gene expression are assessed more often for model organisms than minor species. Also, importantly gene
expression are not measured well under specific conditions for human, due to high risk of direct experiments, such as
non-approved treatment or gene knockout, and then often substituted by mouse. Thus, predicting enrichment significance(on a phenotype) of a given gene set of a species (target, say human), by using gene expression measured under the same
phenotype of the other species (source, say mouse) is a vital and challenging problem, which we call CROSS-species gene set
enrichment problem (XGSEP). Results: For XGSEP, we propose the CROSS-species gene set enrichment analysis (XGSEA), withthree steps of: (1) running GSEA for a source species to obtain enrichment scores and p-values of source gene sets; (2)representing the relation between source and target gene sets by domain adaptation; and (3) using regression to predict
p-values of target gene sets, based on the representation in (2). We extensively validated the XGSEA by using five regression
and one classification measurements on four real data sets under various settings, proving that the XGSEA significantly
outperformed three baseline methods in most cases. A case study of identifying important human pathways for T -cell
dysfunction and reprogramming from mouse ATAC-Seq data further confirmed the reliability of the XGSEA. Availability:
Source code of the XGSEA is available through https://github.com/LiminLi-xjtu/XGSEA.
Key words: cross-species study; domain adaptation; gene set enrichment analysis

Introduction
Due to recent advancement of modern experimental technologies, currently we have a massive amount of basic biological
data. For example, next-generation sequencing technology has
made sequencing faster and lower-cost, generating an incredible
number of sequences. This situation makes bioinformatics tools
more promising in retrieving biological knowledge from the data.
For example, gene set enrichment analysis (GSEA) [1] has been
well used in biology and related areas, which can rank geneset(s) most relevant (precisely, statistically significant) to binarylabeled gene expression measurement. However, GSEA needs
gene expression data labeled binary, such as control and case,and is heavily affected by missing data.
Indeed gene expression are now measured by more speedy
and precise techniques like RNA-Seq than cDNA microarray,
while measuring gene expression is still costly both on money
and time. Existing expression data often have strong bias inmeasured organisms or species. Model organisms, such as

Menglan Cai is currently working toward the PhD degree in Xi’an Jiaotong University, China. Her research interest includes bioinformatics.
Canh Hao Nguyen has been working in machine learning and bioinformatics. His current interest includes machine learning for graph data withapplications in biological networks.
Hiroshi Mamitsuka is a professor of Kyoto University, Japan. His current research interests include mining from graphs and networks in biology andchemistry.
Limin Li is a professor of Xi’an Jiaotong University, China. Her research interests include machine learning and the applications in bioinformatics andbiostatistics.
Submitted: 18 September 2020; Received (in revised form): 12 December 2020
© The Author(s) 2021. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com

1

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

XGSEA: CROSS-species gene set enrichment analysis

2

Cai et al.

general domain adaption, a machine learning model, trained bya larger amount of labeled samples from a source domain, isapplied to a target domain with very few or no labeled samples[11]. This is exactly the same situation of XGSEP. A common way
of domain adaptation methods is to train a model so that themodel can reduce the probability gap between two domains. A
possible measure for the probability gap, i.e. the difference of
two data distributions, is maximum mean discrepancy (MMD)[12–15]. We will borrow the idea of domain adaptation and MMD
to solve XGSEP.
We propose a method, XGSEA, standing for Cross-species
Gene Set Enrichment Analysis (XGSEA). The XGSEA solves XGSEPby three steps: (1) we run GSEA over the source gene sets toobtain gene enrichment scores Es and gene enrichment significance vs . (2) We first define pairwise similarities among gene setsbased on M, and then propose a MMD-based domain adaptation
method to project Xs and Xt into a latent common space withaffine mappings Ps and Pt to obtain Zs and Zt , respectively, sothat (i) the probability gap between Zs and Zt in the latent spaceis minimized and (ii) Ps and Pt are smooth over the connection M
between source and target gene sets.
By solving this optimization problem, we can obtain the
optimal new representations Zs and Zt for source and target genesets, respectively. (3) We train a regression model by (Zs , Es ) andrun the trained model over Zt to predict enrichment scores Et
for target gene sets and finally p-values vt with the principle ofnull hypothesis. Schematically, we may be able to explain ourPs

Ptidea by using arrows: Xs → Zs and Xt → Zt , so that the adaptive
representations Zs and Zt for source and target gene sets shouldhave the smallest distribution divergence and preserve theirpairwise homology similarities.
The contribution of this work can be summarized into
three-fold: (1) we define a problem, XGSEP, which is helpful
for understanding a particular phenotype (label) of a species
with too limited data to run GSEA. (2) We propose a three-stepmethod called XGSEA for XGSEP through domain adaptationthat projects gene sets from two species into a common latent
space. This projection is formulated as a nonlinear optimization
problem, by which we can estimate the latent space and also
estimate the enrichment scores and p-values of target gene setsthrough the latent space. Furthermore, the computational complexity of the optimization problem is low enough so that thecomputation of the XGSEA becomes feasible over regular geneannotation matrices. (3) We empirically validated the XGSEA
by using four different real phenotypes with expression data.
The experimental results show that the XGSEA significantly
outperformed three baseline methods under various settings.
The advantage of the XGSEA was further confirmed by a casestudy of finding significant unknown human pathways for T-cell
dysfunction and reprogramming from a mouse ATAC-Seq data
set.

Method: CROSS-species Gene Set Enrichment
Analysis (XGSEA)
To the best of our knowledge, there is no existing work for
XGSEP. A similar problem setting might be cross-species gene set
analysis (XGSA) [16]. The goal of the XGSA is different with our
XGSEP. XGSA aims to compare a gene set from one species witha gene set from another species. That is, XGSA directly examines
if two gene sets (from two different species) are significantly
different or not, only through the homology between genes ingiven two gene sets.

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

Mus musculus, Caenorhabditis elegans, Arabidopsis thaliana, etc.,
are well measured, while data on minor species are relatively
insufficient. Additionally, human gene expression data areunable to be measured under some specific conditions, due tohigh risk of direct experiments on human, such as non-approved
treatment or gene knockout. On the other hand, mouse is usuallyused to study human disease [2, 3] because of lower cost, lower
risk and relatively strong homology relationship with human[4]. However, there exists essential differences between mouseand human [5–8]. Effective treatments developed by mouse
data often fail in human clinical trials [9, 10]. Thus it would
be strongly expected to develop a method to bridge the gap
between expression data of different species, such as humanand mouse.
We consider a problem of predicting enrichment significanceof given gene sets of one species (such as human) without
gene expression, by using sufficient gene expression data ofanother species (such as mouse). The assumption behind thisproblem is that both expression data are measured under thesame phenotype. We call this problem cross-species gene set
enrichment problem (XGSEP). Assume that we have enough data
behind XGSEP for human and mouse (more generally target andsource), except target expression data. A gene set, either frommouse or human, could be represented as a binary annotation
vector with dimension being the number of all genes in theexpression data, representing whether the corresponding geneis in the gene set. The enrichment significance (such as p-value)of a source gene set S with an annotation vector xs can becomputed by traditional GSEA. The goal of XGSEP is to predict
the enrichment significance for a target gene set T with an annotation vector xt , which might have a different dimension from xs
since the number the total genes for target (human) and source(mouse) are different. Note that the sequence homology between
target genes and source genes is assumed to be represented bybinary matrix M, which should be important information for theprediction.
A naive idea for XGSEP would be to first find a source geneset xs , most homologous to genes in a particular target geneset xt , by using M. Then, GSEA is run over source expression
data and xs . The resultant p-value for xs is considered as aprediction of the enrichment p-value for xt . The method is simple
and fast, but the homology relationship between source andtarget is often complex, and thus homologous source gene set
xs cannot be clearly defined. Also using M directly would be notrobust.
Our idea for XGSEP is, rather than focusing on only one geneset, to consider many gene sets at once and train a predictive
machine learning model by these gene sets. Suppose that wehave source gene sets S1 , . . . , Sm and target gene sets T1 , . . . , Tn ,1nwith annotation matrices Xs = [x1s , . . . , xms ] and Xt = [xt , . . . , xt ],respectively. Then the enrichment p-values for the source genesets can be computed beforehand (by traditional GSEA). The
goal of XGSEP is to predict enrichment p-values for target genesets x1t , . . . , xnt . Note that Xs (training data) and Xt (test data) are
different in size of rows (number of genes), and thus it is difficult
to compare the two matrices directly, meaning that a regularmachine learning model such as a classifier generated by Xs
cannot be run directly over test data Xt . Thus, a further idea isto transform both the target and source species into a common
space so that the target and source genes can be compared. However, this idea cannot be realized by regular machine learningmodels by the above problem of difference in size between training and test data. We solve this problem by domain adaptation,transfer learning between two domains: target and source. In

CROSS-species Gene Set Enrichment Analysis

Problem definition

Step 1: gene set enrichment analysis for source
Since gene expression Gs and phenotype ys are both available
for the source side, we can directly use regular GSEA to obtainp-values, vs,1 , · · · , vs,m for S1 , · · · , Sm , respectively. In fact, p-value
vs,i corresponds to null hypothesis H0s,i : gene set Si has no association with phenotype ys (against the entire set of genes) and canbe computed by the following procedure [1].
1a. Compute enrichment score E0s,i for gene set Si by using gene
expression Gs and phenotype ys .
1b. Permute the entries in ys and recompute the enrichmentscore for gene set Si . Repeat this step B times to generate
an empirical null distribution of the enrichment score: ENULL
with E1s,i , · · · , EBs,i .
1c. Compute empirical, nominal p-value vs,i for Si from null
distribution ENULL by using the positive (or negative) regionof the distribution corresponding to observed enrichmentscore E0s,i .
For source gene set Si , we can compute B+1 enrichmentscores E0s,i , · · · , EBs,i in 1a and 1b to compute p-value vs,i in 1c. Similarly for target gene set Tj , we can first predict B+1 enrichmentscores E0t,j , · · · , EBt,j for target gene set Tj and then p-value vt,j in 1c.

Step 2: domain adaptation for source and target genesets
We project the target and source genes into a common space, tomaximally use the information from the source gene side for thetarget gene sets. The objective function is formulated as follows.
We project Xs and Xt to a common subspace in Rd by using
affine mappings Ps ∈ Rp×d and Pt ∈ Rq×d , respectively, such thatTT
1n
Zs = [z1s , · · · , zms ] = Ps Xs and Zt = [zt , · · · , zt ] = Pt Xt .
In this process, we can set the following two reasonableobjectives:(1) Probability divergence between Zs and Zt should be small.

(2) Pairwise distances among the gene sets in Zs and Zt should
be preserved.
For the first objective, we use MMD [12, 14] to measure thedivergence. An empirical estimate of MMD can be defined asfollows:

D (Zs , Zt )

=



1 
1
φ(zis ) −φ(zit )2H ,
m
n
m

n

i=1

i=1

m


=

m,n
n
j
j
j

k(zis , zs )  k(zit , zt )k(zis , zt )
+
−2
2
2
m
n
mn
i,j=1
i,j=1
i,j=1

=

trace(KL),(1)

where φ(·) is a mapping to reproducible kernel Hilbert space H,k(·, ·) = (φ(·), φ(·)) is the kernel associated to this mapping, and

K=

Kss
Kts

Kst
Ktt


∈ R(m+n)×(m+n) ,(2)

where the (i, j)-element of Kab isj

Kab (i, j) = k(zia , zb ), a, b ∈ {s, t}, i = 1, · · · , m, j = 1, · · · , n,and the (i, j)-element of L is
⎧
2
⎪
⎨ 1/m
2
L(i, j) =1/n⎪⎩ −1/mni, j ∈ {1, · · · , m};i, j ∈ {m + 1, · · · , m + n}otherwise.

(3)

For the second objective, we can first define the pairwise
homologous similarity between source gene sets S1 , · · · , Sm andtarget gene sets T1 , · · · , Tn from given data directly as follows:

Wss (i, j) =ijmin{ |S,||S ∩S |

|Si ∩Sj |}|Sj |i, j ∈ {1, · · · , m};

Wtt (i, j) =|Ti ∩Tj | |Ti ∩Tj |min{ |T, |T | }i|ji, j ∈ {1, · · · , n};

Wst (i, j) =jimin{ |T,||T ∩S |i ∈ {1, · · · , m},ij|Si ∩Tj |}|Si |j ∈ {1, · · · , n},(4)where |A| is the number of genes in set A, Si = φM (Si ) ⊂ T is theset with the target genes homologous to the source genes in Si ,and Tj = φM (Ti ) ⊂ S is the set with the source genes homologous
to the target genes in Tj . The projections Ps and Pt should be
Wss Wstsmooth over homologous similarity matrix W =WTst Wttso that Zs and Zt could preserve the gene sets’ pairwise distances
captured by W.
Thus, the divergence D in Equation (1) should be minimized,
being regularized by the smoothness of the projections Ps and Pt
over similarity matrix W. Overall the objective function can be

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

We have two species, source and target. Let A = {a1 , · · · , ap } be asource (say mouse) gene set, and B = {b1 , · · · , bq } be a target (sayhuman) gene set. Let M ∈ Rp×q be a binary matrix of sequencehomology, where the (i, j)-element M(i, j) is one if source gene ai
is homologous to target gene bj ; otherwise zero. Suppose that
we have gene expression matrix Gs with phenotype vector ys forsource genes only, meaning that we can run GSEA over Gs andys to compute gene set enrichment significance for an arbitrary
source gene set.
Suppose further that we have multiple gene sets for both
source and target. Let S = {S1 , · · · , Sm } be m source gene sets andT = {T1 , · · · , Tn } be n target gene sets. We define an annotation
matrix for source gene sets S (for columns) by A (for rows) as
p×m
Xs = [x1s , · · · , xm
for source gene sets S1 , · · · , Sm , wheres ] ∈ {0, 1}jthe i-th element of xs is one if gene ai is in gene set Sj andzero otherwise. Similarly, let Xt = [x1t , · · · , xnt ] ∈ {0, 1}q×n be theannotation matrix for target gene sets T . Then the problem,XGSEP standing for CROSS-species Geneset Enrichment Problem,
is, given Gs , ys , Xs , Xt and M, to estimate the enrichment p-value ofeach gene set in T with respect to the same phenotype of ys . We
propose our method XGSEA, standing for CROSS-species GeneSet Enrichment Analysis, to solve XGSEP by using three steps.
Figure 1 shows a schematic picture of the three-step procedureof the XGSEA. Below, we will explain each of these three steps indetail.

3

4

Cai et al.

XGSEA-D: we train parameter α in logistic regression by source:

given as follows:

PTs Ps +PTt Pt =I

D (PTs Xs , PTt Xt ) + λ(

+

i,j=1

logit(vs,i ) := log(

l=1

i,j=1

m,n

j

Wst (i, j)zis − zt 22 +

1
2

n


j

Wtt (i, j)zit − zt 22 ),


vs,i
)=
α(l)zis (l) + i , fori = 1, · · · , m,
1 − vs,id1jWss (i, j)zis − zs 22
2
m

min

(5)

i,j=1

where the constraint PTs Ps + PTt Pt = I is used to avoid trivial
solutions.
The model (5) is rewritten in a simpler form and then can besolved by the optimization algorithm on Grassmann manifold.
The details of our optimization algorithm and our computational
complexity reduction strategy are given in Section 1.1 of the supplementary material. Algorithm 1 in the supplementary materialshows the pseudo code for the optimization algorithm.

Step 3: enrichment scores and p-values for target
In order to estimate p-values for target gene sets from adaptive
representations Zs and Zt (obtained in Step 2), we present thefollowing three methods:

and predict p-values for target using z1t , · · · , znt , directly.
XGSEA-E (enrichment scores and p-values): We already computeenrichment score Ebs,i of source gene set Si at the b-th permutation (b = 0 means no permutation) in Step 1. For any b ∈ {0, · · · , B},
we train a linear regression model for enrichment scores fromsource gene sets:

Ebs,i =dβ b (l)zis (l) + i , fori = 1, · · · , m,l=1j

and predict enrichment scores Ebt,j for target gene sets using zt .
We further estimate enrichment p-values for target gene sets byusing 1c of GSEA in Section ’Step 1: gene set enrichment analysis
for source’.
XGSEA-E± (two enrichment scores and p-values): we train two
linear regression models for positive and negative enrichment

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

Figure 1. Flow chart of the XGSEA: we (A) compute B+1 enrichment scores and p-values for each source gene set by GSEA, where B is the number of permutation, (B)obtain new representations for all source and target gene sets by domain adaptation and (C) predict enrichment p-values for target gene sets by a regression model
based on the new representations.

CROSS-species Gene Set Enrichment Analysis

scores, separately as follows:

Ebs,i =
Ebs,i =d+bil=1 γ b (l)zs (l) + i , ifEs,i ≥ 0
d
−
b
i
γ
(l)z
(l)
+

,
ifE
< 0,
i
s
l=1 b
s,i

Experiments
Comparison methods
We compared the XGSEA with three baseline methods, HM1 ,HMA and HMO , which all directly map each target gene to source
genes based on sequence homology, and estimate the enrichment p-value of target gene set T from enrichment p-valueof particular source gene set S. These three baseline methods
correspond to different strategies to generate S:
HM1 : S has a randomly chosen gene homologous to each genein T (i.e. |S| = |T|).
HMA : S has all genes homologous to each gene in T (i.e. |S| ≥|T|).
HMO : S has, out of gene sets predefined by biological pathways and GO terms, the set with genes most overlapped withthose in T.
Note that it is reasonable to compare our methods withthe three baseline methods due to two reasons. First, there areno other advanced methods proposed for XGSEP. Second, thethree baseline methods show fairly good performance in our
experiments later. Since we propose three methods, thus wecompared totally six methods : XGSEA-D, XGSEA-E, XGSEA-E±,
HM1 , HMA and HMO .

Data sets
To evaluate the performance of the XGSEA, we need target
expression data, so that we can compute ground truth enrichment p-values. We collected four gene expression data sets
as below, where each data set consists of human (target) andanother species (source: mouse or zebrafish) which share thesame phenotype. Table 1 shows the statistics of the four datasets.
1. Embryonic development (human and mouse): the two data
sets were collected from www.ncbi.nlm.nih.gov/geo withaccessing number GSE44183. Both gene expression data
sets were obtained from single cell RNA sequencing. In
the human data set, there are 29 samples with 14 766
genes and seven embryonic development stages, including
oocytes, pronucleus, zygote, 2-cell, 4-cell, 8-cell and morula.
For the mouse, there are 17 samples with gene expression levels of 13 879 genes at six embryonic development
stages, including oocytes, pronucleus,2-cell, 4-cell, 8-celland morula. These data sets were used in a cross-species
study [17] already, while this study is not on GSEA.

2. Brain cancer (human and mouse): the data sets of the two
species were downloaded from GEO with accession number
GSE45874 and GSE38591, respectively. Both data sets were
measured by microarray. The human data set has 44 030
genes with six disease and six control samples, while themouse data set has 9653 genes with four disease and four
control samples. These data sets were also used in anothercross-species study [18], while this study is also not on GSEA
at all.
3. Ovarian cancer (human and mouse): the two microarraygene expression data sets were downloaded from GEO withaccession number GSE6008 and GSE5987, respectively. The
human data set has 21 188 genes with 13 mucinous ovarian
tumors and four control samples, while the mouse data
set has 45 101 genes with seven disease and four control
samples. These data sets were also used in the cross-species
study [18].
4. Melanomas (human and zebrafish): the microarray gene
expression data sets of the two species were downloaded
from GEO with accession number GSE83343 and GSE83399,
respectively. The human data set has 42 346 genes witheight disease and four control samples, while the zebrafish
data set has 13 620 genes with five disease and three control
samples. These data sets were collected from two different
studies [19, 20].
We then accessed Ensembl BioMart through http://www.
ensembl.org/ [21] to retrieve homology relationships between
19 404 human and 19 614 mouse genes, and also 16 070 humanand 18 324 zebrafish genes. The homology data from Ensembl
is produced at the protein level rather than the DNA level bywhole-genome alignments of vertebrate species [22, 23].
We showed two homology matrices between human andmouse (left) and between human and zebrafish (right), respectively, in Figure S1 (see the supplementary material). We can seethat genes cannot be assigned in a simple manner of one-to-one
correspondence.
We collected 674 human gene sets (pathways) from Reactome
in Molecular Signatures Database (MSigDB), 2250 mouse genesets from http://baderlab.org/GeneSets and 1550 zebrafish genesets from http://bioinformatics.org/go2msig/. Note that for themelanomas data set, the genes in some gene sets do not havegene expression data, and thus GSEA could not be performed toobtain the ground truth p-values. After removing these gene sets,
only 664 human gene sets are left in the melanomas data set.

Experimental setting
In our experiments, we took human species as the target species,
and took mouse or zebrafish as the source species. We apply
our XGSEA approach to predict the enrichment p-values for the
674 human pathway gene sets T = {T1 , · · · , Tn }(n = 674), forembryonic development, brain, ovarian and melanomas, respectively. For the target gene sets T = {T1 , · · · , Tn }, we take thetraining source gene sets S = {S1 , · · · , Sn } in the XGSEA, where
Si corresponds to Ti , meaning that each gene in Si is homologous
to one or more genes in Ti .
To sufficiently evaluate our XGESA methods, we predictenrichment p-values for target gene sets with three experimental settings. Note that the homology between two genes can
be classified into four types: one-to-one, many-to-one, one-tomany and many-to-many, where one-to-one means only one
gene in one side is homologous to only one gene in the other(1)(1)
side. First level is for simple target gene sets T (1) = {T1 , · · · , Tn },Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022jfor i = 1, · · · , m. We then compute distances zt − z−t 2 andj+−zt − z+t 2 , where zt and zt are the centers for Zs with positive
and negative enrichment scores, respectively, and then assignjthe positive model to target gene set zt if the positive distanceis smaller; otherwise the negative model. We predict the enrichment scores of target gene sets for each b and estimate p-values
by using 1c of GSEA in Section ’Step 1: Gene set enrichment
analysis for source’.
More details of Step 3 can be found in Section 1.2 of the
supplementary material.

56

Cai et al.

Table 1. Statistics of expression data and human gene sets (T (3) , where the cutoff p-values were 0.01 and 0.05 for embryonic development andthe others, respectively)

Species

Gene expression
Number ofgenes

Number ofsamples

Test sets in human
Number of
Number of setslabelsEmbryonic

Human
Mouse
Human
Mouse
Human
Mouse
Human
Zebrafish

14 76613 87944 0309 653
21 18845 10142 346
13 620

29
17
12
8
17
11
12
8

7
6
2
2
2
2
2
2

Brain
Ovarian
Melanomas(1)

where each Ti⊂ Ti only includes the target genes in Ti with(1)label ’one-to-one’. For this case, each target gene g in set Ti
only has one homologous source gene, which does not have anyother homologous target gene except g. The second case is for(2)(2)
more complex target gene sets T (2) == {T1 , · · · , Tn } , where each(2)Ti ⊂ Ti only includes the target genes in Ti with label ’one-toone’ and ’one-to-many’. For this case, each target gene g in set(2)
Ti only has one homologous source gene, which may or may
not have other homologous target genes besides g. The thirdcase is the most complicated case with pathway target genesets T 3 = T = {T1 , · · · , Tn }, where the target genes may have anyof four labels.
In summary, we consider three levels for T , i.e. T (1) , T (2) and
(1)
(2)
(3)
T (3) , where Ti ⊂ Ti ⊂ Ti ( i = 1, . . . , n):(1)T(simple): each set in T (1) has one-to-one genes only. That
(1)
is, target gene g ∈ Ti has only one homologous source gene s,
which has no other homologous target genes except g.
T (2) (medium): each set in T (2) has one-to-one or many-to-one(2)
genes. That is, target gene g ∈ Ti has always only one homologous source gene s, which has one or more homologous target
genes including g.
T (3) (complex): each set in T (3) target gene g may have one ormore homologous source genes, and one of them s also may haveone or more homologous target genes, including g.
To evaluate the methods sufficiently, we conducted bootstrapping on 674 human gene sets in T (3) (or T (1) and T (2) )as follows. We first sampled 674 gene sets with replacementfrom T (3) for 20 times. For each sampling, the baseline
methods and our methods were then used to predict p-values
for the sampled gene sets that were used to calculate theevaluating measurements. Finally, we reported average resultsover the 20 trials. The parameters d and λ were chosenfrom {5, 10, 20, 30, 40, 50} and {0.01, 0.1, 1, 10, 100}, respectively,
to give the best performance under each experimental
setting.

Evaluating the XGSEA methods by regression andclassification measurements
The goal of the XGSEA methods is to predict p-values fortarget gene sets, and thus to further discover significant crossspecies gene sets. We first used five regression measurements
to evaluate the performance of our methods for p-value
prediction, and then used one classification measurement (AUC
area under ROC curve) to evaluate the performance of the
XGSEA methods for discovering significant cross-species genesets.

Number of
positive sets

674

24

674

24

674

13

664

15

Evaluating regression performance
In this section, we first predicted p-values for the target gene sets
in T (3) of the four data sets, respectively, by the XGSEA methods
and the baseline methods, and then compared the regression
performance by using five regression measurements, including
mean square errors (MSEs), mean absolute errors (MAEs), concordance index (CI), Pearson correlation and Cosine similarity. The
definition of the five regression measurements can be found in
Section 1.3 of the supplementary material. We first transformed
p-values to negative log p-values (− log p) and then computed themeasurements.
We reported the bootstrapped MSEs for four real data sets
under T (3) by the six methods in Table 2. We also reported theresults for the MAEs and the Cosine similarities in Table S1 andTable S2 in the supplementary material. Based on the resultsfor MSEs and MAEs, the three XGSEA methods performed better
than the three baseline methods for most cases, and the XGSEAD method tends to perform the best for all the four data sets.
As for the Cosine similarity, although the baseline methods areeven better than the XGSEA-D and XGSEA-E for most cases, the
XGSEA-E± performed the best among all the methods.
For the measurements of concordance index and Pearson correlation, our experiments showed that none of the six methods
could obtain reasonable values if we consider all the tested genesets. Thus, we took an alternative strategy to evaluate the performance by the measurements based on part of the gene sets.
Based on the ground truth p-values for the target gene sets, weselected the k most significant and the k most insignificant genesets. We then calculated the concordance indices and Pearson
correlations between the negative logarithm of the ground truth
p-values and the predicted p-values for the selected 2k gene sets.
We reported the bootstrapped CIs and Pearson correlations bychanging k from the set of {20 : 20 : 300} for all competing
methods on the four data sets with T (3) in Figure 2 and Figure
S2, respectively. The results show that the CIs or the Pearson
correlations basically decrease along with increasing k for mostof the methods and data sets, and the XGSEA methods could
obtain higher CIs and Pearson correlations than the baseline
methods for most cases.

Evaluating classification performance
The XGSEA methods predict p-values for cross-species gene sets.
One important application of these p-values is to discover thesignificant gene sets, by setting up a significant level. Thus, we
further evaluated the classification performance of the XGSEA
methods by comparing AUCs (area under ROC curves) with the

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

Data setCROSS-species Gene Set Enrichment Analysis

7

Table 2. Bootstrapped mean square errors of six competing methods on four data sets under gene set T (3)

HM1

HMA

HMO

XGSEA-DXGSEA-E

XGSEA-E±

Embryonic
Brain
Ovarian
Melanomas

0.33
0.18
0.20
0.23

0.32
0.18
0.19
0.22

0.41
0.16
0.16
0.14

0.30
0.15
0.12
0.13

0.34
0.18
0.13
0.13

0.32
0.15
0.12
0.13

Note: The best in each row are in bold.

Figure 2. Bootstrapped CIs of six methods on four data sets (T (3) )for the selected 2k human gene sets (target), including the k most significant and the k mostinsignificant gene sets based on the ground truth p-values of human gene sets.

baseline methods. We first labeled the target gene sets to positive and negative ones by setting a cutoff p-value (significancelevel) for the ground truth p-values so that a gene set is a positive
instance if the true p-value of this instance is lower than the
cutoff; otherwise a negative. This means that we can obtain
positive and negative gene sets by changing the cutoff p-value.
We then could compute the AUC values by comparing the predicted p-values and the ground truth labels for the target genesets. Note that the AUC values rely on the chosen cutoff p-values.
We first examined the performance of the competing methods by fixing the cutoff value for labeling. Table 3 shows bootstrapped AUCs under three different gene sets (T (1) , T (3) and T (3) )by all six methods, fixing the cutoff at 0.01 for embryonic development and 0.05 for the other data sets. This table shows that
the XGSEA significantly outperformed the baseline methods formost cases. For example, the XGSEA-E± achieved the best in 9

out of all 12 cases, followed by the XGSEA-E of three cases. Any
naive method could neither be the best nor the second best in allcases, the difference from the best being statistically significant
in t-test over 20 trials. Also the AUC of T (1) was not necessarily
higher than T (2) (also T (3) ), since each one-to-one homologous
gene pair between two species is not necessarily the same gene,which would be prediction-wise harder than the case that thetarget and source gene sets share the same gene.
We then checked how the AUCs change along with cutoff
p-values in the set {5e-1, 1e-1, 5e-2, 2.5e-2, 1e-2, 5e-3, 2.5e-3, 1e-3}.
Figure 3 shows the bootstrapped AUCs of all methods on all
four real data sets, under T (3) , by changing cutoff p-values. The
bootstrapped AUCs basically increase as the cutoff p-values
decrease, due to the decreasing number of the true positives. For
most cutoff p-values, the XGSEA methods could obtain higher
bootstrapped AUCs than the three baseline methods.

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

Data set8

Cai et al.

Table 3. Bootstrapped AUCs of six competing methods on four data sets and three target gene sets

Data set
Embryonic

Ovarian

Melanomas

HMA

HMO

XGSEA-DXGSEA-E

XGSEA-E±

0.81 (6.59e-06)
0.80 (2.12e-06)
0.79 (1.58e-09)
0.66 (3.14e-01)
0.59 (1.00e-04)
0.58 (1.75e-07)
0.45 (2.53e-12)
0.56 (6.72e-09)
0.57 (5.60e-12)0.72 (3.65e-12)
0.63 (6.14e-05)
0.44 (1.74e-16)

0.81 (6.59e-06)
0.80 (2.12e-06)
0.80 (4.43e-09)
0.66 (3.14e-01)
0.59 (1.00e-04)
0.60 (1.36e-05)
0.45 (2.53e-12)
0.56 (6.72e-09)
0.61 (1.50e-07)
0.72 (3.65e-12)
0.63 (6.14e-05)
0.44 (2.90e-15)

0.75 (1.48e-08)
0.74 (6.84e-09)
0.75 (3.73e-11)
0.58 (1.14e-05)
0.57 (5.59e-06)
0.55 (2.64e-07)
0.57 (1.65e-04)
0.50 (2.07e-08)
0.46 (6.60e-14)
0.47 (2.10e-16)
0.48 (8.01e-14)
0.59 (4.68e-06)

0.86
0.86
0.87
0.60
0.60
0.61
0.67
0.67
0.65
0.84
0.74
0.64

0.80
0.83
0.83
0.68
0.66
0.63
0.64
0.69
0.70
0.92
0.80
0.72

0.89
0.89
0.90
0.67
0.67
0.68
0.70
0.75
0.77
0.87
0.81
0.71

Notes: The best and second best in each row are in bold and underlined, respectively. The p-value by t-test between the best and each corresponding baseline method
is shown in brackets.

Figure 3. Bootstrapped AUCs on four data sets (T (3) ) with changed cutoff p-values (and the corresponding numbers of true positives in the brackets).

Discussion
Prediction results without bootstrapping
In the above section, we evaluated our XGSEA methods by bootstrapping gene sets in T (1) , T (2) or T (3) , and here we compared our
methods with the baseline methods on the original 674 gene sets
in T (3) without bootstrapping.

We changed cutoff p-values from the same set in Section’Evaluating classification performance’. Figure S3 (see thesupplementary material) shows the AUCs of all six methodson all four data sets under the original T (3) with varied
cutoff p-values. We can see that Figure S3 is similar withthe bootstrapped Figure 3, while the curves in Figure 3 seem
more stabilized and differentiable due to the bootstrapped

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

Brain

T (1)
T (2)
T (3)
T (1)
T (2)
T (3)
T (1)
T (2)
T (3)
T (1)
T (2)
T (3)

HM1

CROSS-species Gene Set Enrichment Analysis

Table 4. Eleven human pathways (with p-values) identified by the
XGSEA-E for T-cell dysfunction and reprogramming

p-value

Gene expression (Transcription)
A third proteolytic cleavage releases NICD
Signaling by NOTCH
Immune system
Signaling by NOTCH3Signaling by NOTCH4
NOTCH2 activation and transmission of signal to
the nucleus
Activated NOTCH1 transmits signal to the nucleus
Signaling by NOTCH2
Constitutive signaling by NOTCH1 HD+PEST
domain mutants
Signaling by NOTCH1

0.03
0.03
0.03
0.04
0.04
0.04
0.04
0.04
0.04
0.04
0.04

procedure. Both figures show the performance advantage ofthe XGSEA methods over the three baseline methods for mostcases.

Robustness against parameters
We first reported the most often best parameters in the bootstrapped results (Tables 2, 3, S1 and S2) in Table S3 (see the supplementary material), for the four data sets and four measurements including MSE, MAE, cosine similarity and AUC, respectively. We can see that different data sets tend to share similar
best parameters, while different measurements may prefer different parameters. Under the measurements MSE and AUC, four
data sets share the similar parameters d = 5 and λ = 0.01. MAE
prefers the same λ and a larger d = 50, while cosine similarity
prefers a larger λ = 100 and a larger d = 50.
We further examined the robustness of the XGSEA, regarding
the parameter λ variation. Figure S4 (see the supplementarymaterial) plotted the AUCs obtained by the XGSEA-E method
versus λ in the set {1e-4, 1e-3, 1e-2, 1e-1}, under three original
gene sets (T (1) , T (2) and T (3) ) of embryonic development andmelanomas. The AUCs for the best baseline methods were alsoplotted. This figure shows that AUC of the XGSEA-E was rather
stable within the given range, implying that the advantage overthe baseline methods will be kept constantly.

Effect of similarity and homology on predictive performance
We examined the contribution of three types of gene set similarity, i.e. Wss , Wst and Wtt , used in the XGSEA, by modifying theobjective function in the formulation of the XGSEA. The objectivefunction of the XGSEA is given by Equation (5), which has fourterms, where the first term is the divergence and the last threeterms are for Wss , Wst and Wtt . We then generated four different
variants of Equation (5), as follows:
MMD: only divergence, i.e. no terms on gene set similarity.
MMD+W: divergence and two terms on Wss and Wtt .
MMD+B: divergence and the term on Wst .
MMD+WB: original objective function, i.e. Equation (5).
We applied these four variants to embryonic development
data with target gene set T (3) . Table S4 (see the supplementarymaterial) shows AUCs obtained with the cutoff (for p-values) of
0.01. From Table S4 (see the supplementary material), MMD+WB(i.e. original Equation (5)) achieved the best result for the XGSEAE and the XGSEA-E±, and MMD was worst for them. This result
implies that all gene set similarity contribute to the performance
improvement.
We then evaluated the effect of sequence homology on predictive performance, by removing a certain amount of part insequence homology matrix M: being motivated by that lesshomology connectivity between two species would cause poorer
performance.
In more detail, we first randomly chose a certain numberof genes from the source and target gene sets, respectively,
and kept only the part corresponding to these genes in M.
Practically, we used 50 500 and 5000 for this number of
selected genes, resulting in three matrices: M50 , M500 and M5000 ,respectively. Using each of the four sequence homology matrices(including original M), we ran the XGSEA over embryonicdevelopment data under gene set T (3) to predict enrichment
p-values.
Table S4 (see the supplementary material) shows the performance results (AUC) of this experiment. The results showthat the AUC was reduced by decreasing the number of randomly selected genes, while if the selected number is 5000,
the performance was almost consistent with that of using theoriginal M, implying that interestingly 5000 genes might be goodenough.

Case study: identifying human pathways for T-cell
dysfunction and reprogramming from mouse ATAC-Seq
It is important for cancer immunotherapy to study the epigenetic regulation of T-cell dysfunction and therapeutic reprogrammability: a plastic dysfunctional state from which T-cells
can be rescued, and a fixed dysfunctional state in which cells
are resistant to reprogramming [24]. Identifying two (plastic or
fixed) dysfunctional chromatin states, through which T-cells intumours differentiate, would be very important to predict, forexample, if a patient will respond to a therapy. Using GSE89308
of GEO on ATAC-Seq data of mouse, with 22 samples and the twochromatin states [24], we ran the XGSEA-E (B = 100 000, λ=0.01
and d=5) to identify human pathways out of 1960 Reactome
pathways (downloaded from https://reactome.org/download-data).
Table 4 shows 11 human pathways identified by the XGSEA-E
at the cutoff of 0.05, where the top, ‘gene expression (transcription)’, and the fourth ‘immune system’ are large pathways with
1367 and 2296 genes, respectively. Obviously due to important
chromatin roles in transcription, ‘gene expression (transcription)’ is tightly related to the chromatin states. Also ‘immune
system’ definitely plays important roles in T-cell dysfunction
and reprogramming through a number of membrane proteins,
such as CD38, CD101, CD30L, CD5, TCF1, IRF4, BCL2, CD44, PD1,
LAG3 and CD62L [24].
The remaining nine pathways are all on Notch signaling
pathways, which affect T cells in various ways. Notch signaling
pathways play multiple essential roles in thymic T cell development and peripheral T cell differentiation [25]. For example, Delta-like ligand 4 (DLL4) interacts with Notch 1 to specify
thymic T cell commitment during lymphocyte development.
This Notch pathway regulates CD8+ T cells by directly upregulating mRNA expression of granzyme B and perforin to maintainmemory T cells [26].
Furthermore, the Notch pathway plays an important rolein antitumor immunity. CD8+ T cell-specific Notch2 deletion
impairs antitumor immunity, whereas the stimulation of the

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

Pathway910

Cai et al.

can be generalized or focused on another statistical problem.
Exploring more efficient, robust domain adaptation would beinteresting future work for Step 2. Reasonably in Step 3, we
can consider more sophisticated regression models. The most
key point of the XGSEA is Step 2, i.e. domain adaptation, which
would be useful for other problems between two species, suchas genome wide association studies between a well- and theother less-sequenced species. This direction of applying domain
adaptation to various problems would be also promising future
work. On the statistical side, we could also further consider the
problem of multiple testing and controlling the false discovery
rate or family wise error rate, which have been well studied inregular GSEA.

Data Availability
Conclusion
We have defined XGSEP for promoting GSEA on species withscarce expression data, and proposed the XGSEA with threesteps, which can be simply: (1) GSEA, (2) domain adaptation and(3) regression. Our empirical supervised validation, including
regression and classification, over four real data sets revealed
that the XGSEA outperformed three baseline approaches underthe measurements MSE, MAE, CI, Pearson correlation, cosine
similarity and AUC for most cases. Particularly, the advantage
was also proved statistically by bootstrapping and t-test. In
the case study, mouse ATAC-Seq expression data are used to
identify significant human pathways for T-cell dysfunction andreprogramming. The XGSEA found rather general two pathways
related with gene expression (transcription) and immune system, as well as nine Notch signal-related pathways, all being convincing, especially compared with pathways found by a baseline
approach. Our XGSEA methods also have potential applications
on minor species once the homology information between theminor species and one major species is available in the future.
There are two characteristics for the proposed XGSEA methods. On one hand, the XGSEA methods apply domain adaptation
to reduce the gap of samples from different species. It’s common that different species samples are drawn from distinguishdistributions. To shorten the distance of two sample distributions between different species for better downstream analysis,
the XGSEA methods take domain adaptation as an effective
strategy to search a latent subspace, in which two distributions
are as close as possible. When the subspace found and new
representations of samples computed for different species, traditional machine learning methods, such as linear regression,
can be used to train in one species and test in the other one.
Besides, compared with three naive methods, which project
linearly the target gene set to the corresponding source one
based on homology network directly, better new representationsof source and target gene sets are computed by the XGSEA methods with capturing the complex relationship among gene setsfrom different species. On the other hand, in the source species,
the XGSEA methods learn a regression model with enrichmentscores, instead of training a model with p-values directly. The
power of this strategy is shown by the performance of the XGSEA
methods on the results in Table 3, which further indicating thattraining on enrichment scores can capture more information
than training on p-values for cross-species gene sets enrichment
analysis.
Improvement of the XGSEA would be definitely interesting
future work. It would be worth working on exploring a better variation on each of the three steps of the XGSEA: Step 1

All data used in our experiments is available and the details ofdata source can be found in Data sets Section.

Supplementary Data
Supplementary data are available online at Briefings in Bioinformatics.

Funding
This work was supported by the National Natural Science Foundation of China (No.11631012 and No.11471256).

Conf lict of Interest
The authors declare no conflict interest.

References
1. Subramanian A, Tamayo P, Mootha VK, et al. Gene set
enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci
2005;102(43):15545–50.
2. Stuart JM. A gene-coexpression network for global discoveryof conserved genetic modules. Science 2003;302(5643):249–55.
3. Zheng-Bradley X, Rung J, Parkinson H, et al. Large scalecomparison of global gene expression patterns in humanand mouse. Genome Biol 2010;11(12):R124.
4. Debry RW, Seldin MF. Human/mouse homology relationships. Genomics 1996;33(3):0–351.
5. Liao BY, Zhang J. Null mutations in human and mouseorthologs frequently result in different phenotypes. Proc Natl
Acad Sci U S A 2008;105(19):6987–92.
6. Mestas J, Hughes CCW. Of mice and not men: differences between mouse and human immunology. J Immunol
2004;172(5):2731–8.
7. Geifman N, Rubin E. The mouse age phenome knowledgebase and disease-specific inter-species age mapping. PlosOne 2013;8.
8. Beura LK, Hamilton SE, Bi K, et al. Normalizing the environment recapitulates adult human immune traits in laboratory mice. Nature 2016;532:512–6.
9. Bugelski PJ, Martin PL. Concordance of preclinical and clinical pharmacology and toxicology of therapeutic monoclonal
antibodies and fusion proteins: cell surface targets. British
Journal of Pharmacol 2012;166.

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

Notch pathway can increase tumor suppression. Ezh2, a suppressor of the Notch pathway, regulates effector T-cell polyfunctionality and survival by targeting the Notch signaling pathway [27]. Down regulation of Ezh2 could elicit poor antitumor
immunity.
Besides, Delta-like 1-mediated Notch signaling enhances theconversion of human memory CD4 T cells into FOXP3-expressing
regulatory T cells [28]. These facts support the reliability of thepathways identified by the XGSEA.
On the other hand, we ran a naive approach, HMA , over thesame data, under the cutoff of 0.05, resulting in 20 pathways
showed in Table S5 (see the supplementary material). Although
the number of pathways is larger than Table 4, these 20 pathways
were diverse and less connected to the chromatin states, such as
only two being related to Notch signaling pathways.

CROSS-species Gene Set Enrichment Analysis

20. Venkatesan AM, Vyas R, Gramann AK, et al. Ligandactivated bmp signaling inhibits cell differentiation anddeath to promote melanoma. J Clin Invest 2018;128(1):294–308.
21. Durinck S, Spellman PT, Birney E, et al. Mapping identifiers
for the integration of genomic datasets with the r/bioconductor package biomart. Nat Protoc 2009;4(8):1184–91.
22. Clamp M, Andrews D, Barker D, et al. Ensembl 2002:
accommodating comparative genomics. Nucleic Acids Res
2002;31:38–54.
23. Flicek P. Ensembl 2014. Nucleic Acids Res 2014;42(Issue
D1):D749–55.
24. Philip M, Fairchild L, Sun L, et al. Chromatin states definetumour-specific t cell dysfunction and reprogramming.
Nature 2017;545(7655):452–6.
25. Freddy R, Robson MH, Fabienne TC. Regulation of innate
and adaptive immunity by notch. Nat Rev Immunol
2013;13(6):427–37.
26. Tsukumo S-i, Yasutomo K. Regulation of cd8+ t cells andantitumor immunity by notch signaling. Front Immunol
2018;9:101.
27. Ende Z, Tomasz M, Ilona K. Et. al. cancer mediates
effector t cell dysfunction by targeting micrornas andezh2 via glycolysis restriction. Nat Immunol 2016;17(1):95–103.
28. Mota C, Nunes-Silva V, Pires AR, et al. Delta-like 1-mediatednotch signaling enhances the in vitro conversion of humanmemory cd4 t cells into foxp3-expressing regulatory t cells.
J Immunol 2014;193(12):5854–62.

Downloaded from https://academic.oup.com/bib/article/22/5/bbaa406/6120324 by Library of Research Reactor Institute, Kyoto University user on 01 April 2022

10. Hünig and Thomas. The storm has cleared: lessons fromthe cd28 superagonist tgn1412 trial. Nat Rev Immunol
2012;12:317–8.
11. Pan SJ, Yang Q. A survey on transfer learning. IEEE Transaction
on Knowledge and Data Engineering 2010;22(10):1345–59.
12. Huang J, Gretton A, Borgwardt KM, et al. Correcting sample
selection bias by unlabeled data. In: Schölkopf B, Platt J,
Hoffman T (eds). Advances in Neural Information Processing
Systems 19. Cambridge, MA: MIT Press, 2006, 601–8.
13. Pan SJ, Kwok JT, Yang Q. Transfer learning via dimensionality
reduction. In: AAAI 2008, pages 677–682, 2008.
14. Arthur Gretton, Karsten M. Borgwardt, Malte J. Rasch, et al. A
kernel two-sample test. Journal of Machine Learning Research,
13:723–73, 2012.
15. Baktashmotlagh M, Harandi M, Salzmann M. Distributionmatching embedding for visual domain adaptation. Journalof Machine Learning Research 2016;17(108):1–30.
16. Djordjevic D, Kusumi K, Ho JWK. Xgsa: a statisticalmethod for cross-species gene set analysis. Bioinformatics
2016;32(17):i620–8.
17. Sun J, Jiang Z, Tian X, et al. A cross-species bi-clustering
approach to identifying conserved co-regulated genes. Bioinformatics 2016;32(12):i137–46.
18. Normand R, Du W. Found in translation: a machine learning model for mouse-to-human inference. Nat Methods
2018;15:1067–73.
19. FV Filipp CL, Boiko AD. Cd271 is a molecular switch withdivergent roles in melanoma and melanocyte development.
Sci Rep 2019;9(1):7696.

11