Since a total of n candidate genes were ranked from highest to lowest expression level, we denoted δ = ( δ 1 , δ 2 … δ n ) as the vector of differences, where each element δ i in δ represented the expression difference of a single gene g i between two groups. Accordingly, a corresponding gene vector g = ( g 1 , g 2 … g n ) was defined, where each element g i was the gene name using NCBI Gene Symbol. To calculate the ES in a KEGG pathway, it is necessary to determine whether the gene g i belongs to this pathway. We established the following definition: if g i belonged to the annotations of this pathway, r ( g i ) = 1 ; if g i did not belong to the annotations of this pathway, r ( g i ) = 0 , and a new vector r = ( r 1 , r 2 , … r n ) was generated, where each element r i consisted of 0 and 1. To facilitate the iteration in the computer program, the calculation of ES was rewritten into the following form according to the method described by Subramanian et al (Subramanian et al., 2005). e s 0 = 0 e s i = r i | δ i | p N R + ( r i − 1 ) 1 n − N H + e s i − 1 , i ≥ 1 where N H = ∑ r i , N R = ∑ r i | δ i | p , e s i represents the value of ES and p is a weighing factor with a range of [ 0 ,   1 ] . A value of p = 1 was set as a constant in the study. By the method of mathematical induction, it was easy to prove that e s n = 0 .

In particular, N R equal to 0 indicated that none of the genes in g belonged to the annotations in the KEGG pathway. To prevent a division by zero error (when N R = 0 ), we set any value of e s i equal to 0 under this condition. Thus, the maximum ES, e s m a x = s i g n ( e s k ) m a x ( | e s 1 | ,   | e s 2 | … | e s n | ) , was obtained, where | e s k | was the maximum element in the m a x ( ⋅ ) function, and s i g n ( ⋅ ) represented the signum function to indicate the positive and negative values of e s m a x .

To estimate the significance of the e s m a x generated from g in the KEGG pathway, a nominal p value was calculated using an empirical phenotype-based permutation method derived from Subramanian et al. (2005). Given the null hypothesis that the order of the elements in g was random, the way to reject the null hypothesis was to calculate the probability of e s m a x in its randomly generated distribution, which was derived from a set of randomly assigned g . When the order of elements in g was disrupted and the other parameters remained unchanged, the new corresponding value of e s m a x was related only to the random rearrangement of elements in g . After repeating the above random rearrangement process v times, a series of values of e s m a x ′ with a total number of v were obtained based on the null hypothesis, and the absolute value was taken: ( | e s m a x 1 ′ | , | e s m a x 2 ′ | , … | e s m a x v ′ | ) . The empirical cumulative distribution function (CDF) F ^ v   for e s m a x ′ (derived from the null hypothesis) can be described as follows: F v ^ ( x ) = 1 v ∑ i = 1 v I ( | e s m a x i ′ | ≤ x ) where I ( ⋅ ) is the indicator function of an event and v corresponds to the number of permutations in Subramanian et al (Subramanian et al., 2005). According to the empirical CDF, the p value for an e s m a x in a pathway was calculated as follows: p = 1 − F v ^ ( | e s m a x | )

In an enriched pathway, the leading-edge subset in g appears prior to the peak score for a positive e s m a x and appears subsequent to the peak score for a negative e s m a x . The leading-edge subset L can be described as follows: L = { { g j | r ( g j ) = 1 ,   j < i }   i f   e s m a x > 0 { g j | r ( g j ) = 1 ,   j ≥ i }   i f   e s m a x < 0   { ∅ }   i f   e s m a x = 0   where i represents the position number of   e s m a x .

According to prior knowledge in the Reactome database, there are three categories of gene effects on another: upregulation, downregulation, or no effect, where the symbols ‘ → ’ and ‘⇢’ represent upregulation and downregulation, respectively. For a vector containing n genes, its correlation matrix M was a square matrix of order n , whose elements were the interaction information among these genes. We defined its element a i j as follows: a i j = { 1   i f   g i → g j   and   i ≠ j − 1   i f   g i ⋯ > g j   and   i ≠ j 0   i f   g i   had   no   effect   on   g j   and   i ≠ j 1   i f   i = j where i and j represent the index numbers of g i and g j in the tuple, respectively.

For example, for a 4-dimensional vector of a gene set ( A ,  B ,  C ,  D ) with the following relationship: { ( A → B ) ,  ( A ⋯ > C ) ,  ( A → D ) ,  ( B ⋯ > C ) ,  ( C → D ) } (1) the correlation matrix was described as follows: M = [ 1 1 − 1 1 0 1 − 1 0 0 0 1 1 0 0 0 1 ]

Genes from the union of the leading-edge subsets were used to iteratively generate gene transactions. After obtaining the union of the leading-edge subsets with m genes, the vector l consisting of these genes and the correlation matrix M reflecting their interactions were created. Elements in l were considered basic gene interactions in the first iteration. For an iteration vector, a gene would be appended to a gene transaction if it had an interaction (upregulation or downregulation) with the gene at the end of the transaction, and all possible new interactions were used as the basic gene interactions for the next iteration. In the study, gene interactions with a maximum of 4 items were generated.

For example, let the 4-dimensional vector l = ( A ,  B ,  C ,  D ) be the union of the leading-edge subsets, and let the interactions be as described above in (1). For 4 iterations, the sequences of gene interactions I = ( I 1 , I 2 , I 3 , I 4 ) were generated as follows: I 1 = { ( A ) ( B ) ( C ) ( D ) }     I 2 = { ( A ,  B ) ( A ,  C ) ( A ,  D ) ( B ,  C ) ( C , D ) }    I 3 = { ( A ,  B ,  C ) ( A ,  C ,  D ) ( B ,  C ,  D ) }    I 4 = { ( A ,  B ,  C ,  D ) } (3)

To detect frequent gene sets from the gene transactions, the value s u p p o r t was used to evaluate the frequency of a gene set. For the vector a of a gene set, the s u p p o r t of a was defined as the total count of a in all transactions.

For example, the s u p p o r t of ( A ,  C ) in interactions I in (3) was 2 ( ( A ,  C ) in I 2 and ( A ,  C ,  D ) in I 3 ).

Given a threshold m i n for s u p p o r t (a threshold m i n = 3 was set in this study), gene interactions I = ( I 1 , I 2 … I t ) generated from t iterations and setting both the frequent gene set and the infrequent gene set, F and F ¯ , to ∅ initially, the algorithm in the study was described in two steps below:

Step 1

Start with i = 1 . Traverse all interactions in I . Append the interactions whose s u p p o r t ≥ m i n to F and the interactions whose s u p p o r t < m i n to F ¯ , respectively. Based on the Apriori principle, if a transaction is found to be infrequent, then all its super interactions are also infrequent. Remove the transaction in I i + 1 (if it exists) when a transaction in F ¯ belongs to the transaction in I i + 1 .

Step 2

Repeat Step 1 until no further successful extensions are found. F is the returned result.

For example, for the gene interactions I described above in (3) and a given threshold m i n = 3 , the process of the Apriori algorithm was as follows: i = 1 :   I 1 = { ( A ) ( B ) ( C ) ( D ) }     s u p p o r t = { ( A ) : 7 ( B ) : 6 ( C ) : 8 ( D ) : 6 }    F = { ( A ) ( B ) ( C ) ( D ) }    F ¯ = ∅

In i = 1 , no interactions in I 2 were removed. i = 2 :   I 2 = { ( A ,  B ) ( A ,  C ) ( A ,  D ) ( B ,  C ) ( C , D ) }     s u p p o r t = { ( A ,  B ) : 3 ( A ,  C ) : 2 ( A ,  D ) : 1 ( B ,  C ) : 2 ( C , D ) : 4 }     F = { { ( A ) ( B ) ( C ) ( D ) } { ( A , B ) ( C , D ) } }    F ¯ = { ( A , C ) ( A , D ) ( B , C ) }

In i = 2 , ( A ,  C ,  D ) and ( B ,  C ,  D ) were super interactions in F ¯ and were removed in I 3 . i = 3 :   I 3 = { ( A ,  B ,  C ) }     s u p p o r t = { ( A ,  B ,  C ) : 2 }     F = { { ( A ) ( B ) ( C ) ( D ) } { ( A , B ) ( C , D ) } }     F ¯ = { { ( A , C ) ( A , D ) ( B , C ) } { ( A , B , C ) } }

In i = 3 , ( A , B ,  C ,  D ) was a super transaction in F ¯ and was removed in I 4 . Iteration was terminated, and the frequent gene set F was obtained.

After the frequent gene set was obtained by the Apriori algorithm, a new statistic d was introduced to evaluate the real ‘weights’ of the frequent gene set consisting of expression levels obtained from RNA-Seq. Given an element vector g from the frequent gene set, d was defined as the multiplication of its expression vector δ on its correlation matrix M : d = δ M

For example, given an element vector g = ( A , B ) from the frequent gene set F , its expression vector δ = ( a , b ) and correlation matrix M = [ 1 1 0 1 ] [derived from the relationship (A → B)] were: d = [ a a + b ]

The modulus of d should be: ‖ d ‖ = a 2 + ( a + b ) 2

Since gene B was upregulated by gene A according to the prior knowledge, theoretically, if the expression value a of gene A was increased ( a > 0 ) , its expression value b should also be increased ( b > 0 ) . However, under the same condition where the expression value a was increased ( a > 0 ) , if the expression value b was decreased ( b < 0 ) , or if the increase in the expression value b was not as large as in the first case, ‖ d ‖ would be smaller in the second case than in the first case, indicating that the weights in the second case were smaller than those in the first case. Specifically, for a single gene  S , ‖ d ‖ is the absolute value of its expression level: ‖ d ‖ = ‖ ( s ) [ 1 ] ‖ = | s |

Any frequent gene set could be described by d no matter how complicated the interactions among genes, and its ‖ d ‖ represented the comprehensive weights related to their interactions and expression levels.

For instance, given the 4-dimensional vector of gene set l = ( A ,  B ,  C ,  D ) , its expression vector δ = ( a , b , c , d ) and its correlation matrix M given above, its comprehensive weights were described as follows: d = ( a , b , c , d ) [ 1 1 − 1 1 0 1 − 1 0 0 0 1 1 0 0 0 1 ] = [ a a + b − a − b + c a + c + d ] ‖ d ‖ = a 2 + ( a + b ) 2 + ( − a − b + c ) 2 + ( a + c + d ) 2

According to this method, the weights associated with any frequent gene set will be referred to as its value of ‖ d ‖ , and the frequent gene set with the highest weights could then be identified by ranking the value of ‖ d ‖ . Since the value of ‖ d ‖ will increase with the increase in the number of items, the value of ‖ d ‖ must be compared among frequent gene sets with the same number of genes.