dataesr/scientific-paragraphs-categorization huggingface dataset:
Indeed, due to the exponential complexity involved in the higher-order exhaustive search algorithms, they are not applicable to large datasets. To address the above issues, several parametric modelling approaches  9, 10  , machine learning algorithms  8, 11  , and combinatorial optimizations  12, 13  have been explored. But they are exclusively designed/used for detecting binary or higher-order interactions in case-control studies. There exist some approaches  [14] [15] [16] [17]  that can detect pairwise interactions in crosssectional studies with quantitative traits but they are simply not scalable to higher-order interactions. Moreover, current methods for detecting epistasis do not address the complexities and challenges involved in longitudinal datasets that allow studying the natural trajectory of traits and/or disease progression.
In this regard, we propose a novel approach of epistasis detection, called epiMEIF, using mixed-effect conditional inference forest (MEIF). The goal of our approach is to reveal higher-order epistatic interactions of genetic variants associated with the variation of complex quantitative phenotypes, primarily from cross-sectional and longitudinal GWAS datasets. The mixed effect conditional inference forest (MEIF) algorithm we propose is a combination of mixed effects model and conditional inference forest (see Figure  1  and Methods, for more details). Recursive partitioning approaches or tree-based algorithms like random forest have already proven to be effective for detecting the genetic loci and their interactions that impact the phenotypic outcome in case-control studies  7, 8, 18, 19. Nevertheless, MEIF is particularly beneficial over existing methods because it i) demonstrates how tree-based algorithms can be adapted for detecting higher-order interactions from cross-sectional and longitudinal datasets using conditional inference forest (cforest), ii) can handle 'large' datasets (as much as 10000 markers) and iii) simultaneously account for missing/censored genome data using cforest, and other confounding factors like population structure in the GWAS datasets using mixed effects model.
Subsequent to identifying interactions from MEIF, we propose to validate the epistatic interactions using two statistical testing approaches: Max T-test and anova test (see Methods, for more details). Overall, epiMEIF (epistasis detection from MEIF along with the additional testing strategies), provides a generalized way to obtain genetic variants and their higher order interactions from any GWAS data. To the best of our knowledge, alternative approach that can detect higher-order interactions in both cross-sectional and longitudinal datasets does not exist, therefore preventing us from evaluating our method via benchmarking. Hence, to evaluate our method, we applied it to an extensive simulation study and illustrated the power performance of the method under several practically relevant scenarios.
We then analysed its performances on a real dataset: natural variation of heart period in young flies (Drosophila) and heart period during aging in Drosophila. We evaluated the method's performance based on its ability to replicate in an independent cohort and to recover previously published genetic data associated with cardiac functions. We demonstrate that the obtained networks of statistical interactions provide insightful information with respect to cardiomyocytes structure and functions. These analyses illustrate the high performance of the method to identify higher-order interactions from both cross-sectional and longitudinal data. The proposed statistical methods are implemented in R and the source codes can be found in the Github repository (https://github. com/TAGC-NetworkBiology/epiMEIF).
Analysis on Simulated Data We evaluated the statistical power of the method with the help of synthetic datasets where the ground truth genetic architecture is known. We used the genotypes of the widely used Drosophila Melanogaster Genetic Reference Panel (DGRP  21  , see Methods) to design our simulation scenarios. The genotype data in our scenarios comprised 1,000 SNPs with a minor allele frequency (MAF) greater than 0. 1 subsampled from 2,456,752 genome-wide SNPs of the DGRP genome dataset. We considered three simulation scenarios where all are additive models containing four SNPs with marginal effects and 2 to 4 SNPs with interaction effects (see Figure  2 , to better understand the overall schematic of how the simulations are conducted). The scenarios differ only in the number of interaction pairs/sets and number of SNPs involved in these interactions.
Leveraging on this dataset, we investigated how Weighted epiMEIF would detect epistasis interactions among variants associated to natural variations of heart period (HP), which measures the duration of a complete cardiac contraction/relaxation cycle. To accommodate computational requirements of Weighted epiMEIF, interactions were tested on variants that show significant association with the quantitative trait from the single GWAS linear mixed model (LMM), with a nominal significance threshold of 10 -5 (3484 SNPs, see Methods). Finally, Weighted epiMEIF application on HP dataset led to a dense network of 35 interacting SNPs (Figure  3  and Figure  4a ). We have first illustrated the statistical interaction network obtained from Weighted epiMEIF (Figure  3 ) and then illustrated the biological significance of the corresponding network with the most connected nodes (Figure  4a ).
Validation of the results We then tested if the epiMEIF network for heart period was replicated in an independent cohort. Twenty DGRP lines -not included in the first dataset and therefore not used in building the statistical interaction network -were analysed for cardiac phenotypes at 1 week. We used Max-T and anova tests to test whether epiMEIF interactions were replicated in this independent population. Note that the interaction network comprised 45 binary interactions that are not part of any higher order interactions, 33 three-way interactions that are not part of any four-way interactions and 2 four-way interactions (see Supplementary Figure  2a  and Supplementary Table  2 ). For each interaction set of size n, anova/regression-based test detects if the n-snp interaction is statistically significant in the validation data (see Methods for more details).
Application to longitudinal data We have shown the efficacy of our approach with cross-sectional datasets. Noticeably, the epiMEIF method can be extended for detecting genetic interactions from longitudinal datasets (see Methods). We have rarely encountered genetic interactions being tested on longitudinal datasets in literature. Malzahn et al. ,  30  tested for gene-gene interaction using longitudinal nonparametric association test on Framingham Heart Study (FHS) cohorts, but they did not test for higher-order interactions. Our approach is particularly interesting because it also presents a way to explore higher-order interactions in longitudinal data. Similar to cross-section data above, we have tested the efficacy of our approach in longitudinal datasets via simulations and real data applications.
Analysis on Simulated Data We evaluated the statistical power of the method on longitudinal datasets using simulation scenarios similar to those used with cross-section data. Here as well, we have 3 simulation scenarios -AS1, AS2 and AS3. While AS1 considers the simplest scenario with 1 binary interaction, AS2 considers a comparatively complicated scenario with 2 binary SNP interactions. AS2 is designed as earlier; one set of SNPs involved in the interaction has a lower marginal effect and the other set of SNPs has a higher marginal effect. AS3 comprises 1 three-SNP interaction (see Methods for more details). We have conducted our simulations 100 times for each scenario, and we run each scenario with 30 or 50 additional variants, randomly chosen from the genome data (with 1000 SNPs), same as earlier.
Note that we have conducted the longitudinal simulations only with epiMEIF. Weighted epiMEIF cannot be easily adapted for longitudinal dataset as deciding a weight for the "Time" covariate (that captures the "longitudinal" component of the data) is challenging. Assigning a heavy weight to "Time" can lead to trees where "Time" is often selected at the root node, whereas assigning a low weight can end up with trees where "Time" is never selected for any node (see Methods for more details). Hence, to avoid complications we focus on epiMEIF application in longitudinal datasets. As in the previous simulations, we evaluated the performance of epiMEIF on longitudinal data based on "the overall power in capturing the true epistatic interactions" (see Table  2 ) and power to capture the true interactions as top-ranking (see Supplementary Figure  3 ).
Analysis on Real Data: epiMEIF network construction for natural variations of heart period aging in flies To further evaluate our method for identifying epistatic interactions from longitudinal data, we analyzed the aging of the heart function in the DGRP population. Several cardiac aging studies in Drosophila have revealed striking similarities with mammals, both in terms of heart physiology and transcriptional changes  22, 23, 31, 32. This highlights the conserved nature of cardiac aging across organisms. We analysed 165 of the DGRP fly lines included in the previous dataset for heart period at 4 weeks of age (12 individuals per line, see Methods). Overall, there is a marked increase of heart period in the DGRP from 1 week to 4 weeks, in agreement with previous observation (see Supplementary Figure  4 ).
Both 1 week and 4 weeks data on the 165 DGRP lines were used to identify epistatic interactions between variants associated to natural variations of heart period during aging. Similar to the previous 1-week study, epiMEIF interactions were identified on variants that show significant association with the aging of HP from the single GWAS longitudinal LMM, with a nominal significance threshold of 10 -5 (1682 ). Finally, epiMEIF application on HP aging dataset led to a dense network of 26 interacting SNPs; comprising 47 two-way that are not part of any 3-way interactions and 12 3-way interactions (Figure  4a , Supplementary Table  4 ). Strikingly, the epiMEIF network comprises tightly interconnected variants within genes involved in diverse biological processes, many of which have characterized cardiac function, either in flies or in mammals.
Validation of the results The twenty DGRP lines not included in the first dataset were also analyzed at 4 weeks of age, thus providing a validation set for cardiac aging which was used to test for replication of the interaction using Max-T/ Anova test in this independent cohort. 75% (35/47) two-way interactions and 58% (7/12) three-way interactions could be analysed in this validation cohort (see Supplementary Table  6 ). Of them, 43% two-way (15/35) and 43% three-way (3/7) were replicated (Figure  4b ). Several features of the network were replicated in the independent cohort, including the involvement of rdog, CG15611 and Astc-R2. This also confirmed the central positioning of DIP-kappa, whose implication in natural variations of cardiac aging therefore warrants further investigations.
Discussion Here, a new method is proposed for epistasis detection in large scale association studies with complex genetic traits using mixed effect conditional inference forest (epiMEIF). This method captures higher-order SNP interactions based on the tree structure in the cforest, and combined with mixed models, can handle a wide range of complex GWA studies. The effectiveness of epiMEIF is verified in extensive simulation scenarios reflecting a wide spectrum of complex models and with real datasets, illustrating its power for epistasis detection from both cross sectional and longitudinal data. The additional testing strategies applied a posteriori to the conditional inference forest in epiMEIF not only safeguards the method from detecting false positive interactions but also increases the reliability of the final selected interactions.
Moreover, by demonstrating the ability of the approach to replicate in an independent cohort we ensure that our approach does not suffer from winner's curse  37  which is quite common in GWAS. The ability of the approach to validate part of the higher-order interactions in an independent cohort also supports its' competency in detecting higher-order interactions. Additionally, for the cross-sectional datasets, we proposed an adaptation of epiMEIF, named Weighted epiMEIF that allows identifying genes associated with weak marginal effect variants. This is an important addition to the 'traditional' epistasis approaches that are biased towards variants with strong marginal effects  7, 8, 20. Similar to LMM, a major advantage of epiMEIF is their ability to effectively account for unwanted correlation between samples, thereby correcting for confounding factors such as population structure  38  or hidden covariates  39  and making it easily applicable for human GWAS studies as well.
However, unlike the standard LMM, epiMEIF can jointly model genetic effects of multiple loci or markers on the readout. This is particularly important because recent works have revealed that often, the single-locus association models are insufficient to explain the heritable component of complex traits  40, 41. The proposed approach has proven to overcome the caveats of existing GWAS approaches and though we have shown its applicability for identifying associations and epistasis detection, it can be also used for prediction  42, 43  and feature selection  20, 44. We however acknowledge that the method is not scalable to the entire genome dataset and do not provide the exhaustive list of interactions from the entire list of variants. Interestingly, the epiMEIF generate dense genetic interaction networks that are creating hubs around some focal genes.
Additionally, it has been often discovered that complex diseases are the results of interaction between a large number of units, and therefore, are more likely to be associated with genes that are well connected in the network  47. epiMEIF can comprehensively model the dynamics of longitudinal data and crosssectional data and offer to reveal the relationships between multiple variants, revealing the extensive networks of genetic interactions that are causing the complex diseases. The possibilities offered by epiMEIF will allow approaching questions that were previously difficult to address, as tools for the identification of epistatic interactions in complex datasets were lacking. In particular, we do not know when biomolecular interactions actually produce patterns of statistical epistasis, nor do we know how to biologically interpret statistical evidence of epistasis, since statistical interaction likely does not automatically entail interaction at a biological or mechanistic level.
METHODS We propose an approach for epistasis detection using mixed effect conditional inference forest (epiMEIF ) to identify the genetic variants and their epistatic interactions responsible for the complex quantitative traits. The novelty of this approach lies not only on the epistasis detection but also on the amalgamation of mixed effects modelling and conditional inference forest (cforest). We have divided the explanation of the epiMEIF method into four parts : 1) we explain the MEIF model, 2) we illustrate how the tree structure in the conditional inference forest (cforest) is utilized to detect high-order SNP interactions that impact the variation of the phenotype, 3) we explain how MEIF can be adapted to weighted MEIF, this is particularly useful for cross-sectional dataset, and 4) we show how the interactions from the MEIF can be validated using independent statistical tests thereby leading to a final Interaction Score Matrix ().
Mixed Effect Inference Forest (MEIF) Our method is primarily inspired by the Mixed effect random forest (MERF) proposed by Hajjem et al 50. It is an ensemble of mixed effects model and random forest where the application of random forest is extended to clustered data for a continuous outcome. The following model is proposed in MERF: 𝑦 𝑖 = 𝑍 𝑖 𝑏 𝑖 + 𝑓(𝑆) + 𝜖 𝑖 , 𝑤ℎ𝑒𝑟𝑒 𝑏 𝑖 ~𝑁𝑞 (0, 𝐷), 𝜖 𝑖 ~𝑁𝑛 𝑖 (0, 𝑅 𝑖 ), … … … … … … … …. 𝑒𝑞𝑢𝑎𝑡𝑖𝑜𝑛(1) In the equation above: •  is the cluster index. Assume that there are  clusters in the training data (e. g. , in the DGRP dataset, cluster denotes the DGRP strain). •   is the target variable. For cluster , there are   observations and   is the   × 1 vector of responses (e. g. , in the DGRP dataset,   is the phenotype of interest). •  is the set of covariates/markers and it is assumed to be  dimensional, e.
Note that we incorporate the cforest on the markers and time covariate only as we would like to capture the interactions of these variables. The covariate time ( ) can be omitted for cross-sectional data. Pooling data across all the clusters, we have the following MEIF model for cross-sectional dataset: 𝑌 = 𝑍𝑏 + 𝑎 ′ 𝑋 + 𝑓(𝑆 1 , 𝑆 2 , ⋯ , 𝑆 𝑁 ) + 𝜖 ⟹ ( 𝑦 1 : 𝑦 𝑖 : 𝑦 𝐾 )=( 𝑍 1 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑍 𝐾 ) ( 𝑏 1 : 𝑏 𝑖. 𝑏 𝐾 )+⟨𝑋 1 |. |𝑋 𝐿 ⟩ ( 𝑎 1 : 𝑎 𝑗. 𝑎 𝐿 ) + ( 1 ,  2 , ⋯ ,   ) +  For longitudinal dataset, the above models can be written as follows: ( 𝑌 1 : 𝑌 𝑡 ) = [ 𝑍 1 ⋯ 0 ⋮ ⋱ ⋮ 0 ⋯ 𝑍 𝑡 ] ( 𝑏 1 : 𝑏 𝑡 ) + ( 𝑎 1 ′ 𝑋 1 : 𝑎 𝑡 ′ 𝑋 𝑡 ) + 𝑓(𝑆 1 , 𝑆 2 , ⋯ , 𝑆 𝑁 , 𝑇) + 𝜖 where,  denote total number of time points in the longitudinal data and superscript  for each component denote the corresponding response and predictors at time point.
This can help in alleviating the computational burden of epiMEIF (illustrated in Figure  1 ). distribution of a group of individuals across the different genotype combinations arising from a two-SNP interaction in the Drosophila Cardiac Genetics dataset (see details on the dataset below) observed at two time points. In this article, 0 and 2 denote the major allele and minor allele, respectively. The highlighted combination 00 denote the major allele combinations of the binary interaction. We consider the population belonging to the major allele combinations (example: 00 in Supplementary Figure  7a ) as the reference population in our dataset. We are interested in identifying the SNP interactions where the longitudinal phenotypic trend of the alternative allele combinations (example: 02, 20, 22 in Supplementary Figure7a) is significantly different from the major allele combinations.
To address that, we remove the longitudinal trend of the reference population from the entire dataset and then explore the epistatic interactions using MEIF. This is illustrated in Supplementary Figure  7a. Supplementary Figure  7b  shows the step-by-step mechanism to compute the interactions sets and their strength of interaction from the forest in a longitudinal data. The different steps (5 steps) in Supplementary Figure  7b  can be summarized as follows: 1) An example of a tree of interest in MEIF is shown here (where Time is not in the root node). SNP A, SNP C and SNP D forms a tree with Time as the child node because the longitudinal phenotype trend of genotype combination 020 (SNP A=0, SNP C=2 and SNP D=0) is significantly different from 000 (SNP A=0, SNP C=0 and SNP D=0). 2) and 3) Whenever the SNPs like SNP A, SNP C and SNP D appears in combination with Time (multiple times), one can say that they are in epistasis, and they are attributing to the longitudinal variation of the phenotype.
Weighted Conditional Inference Forest A notable shortcoming of tree-based methods is that they are dependent on marginal effects. Random forests are essentially bagging approach which begins with the bootstrap sampling of individuals from the original dataset followed by feature selection for the construction of each node of the tree. To do so, instead of considering all variables from the initial GWAS dataset, a random subset of variables is picked without replacement. Hence at the tree learning step, the algorithm looks for a single SNP that well discriminates cases from controls or different groups of phenotypes. In practice, this is equivalent to looking for SNPs with high marginal effects. To overcome this, instead of randomly drawing the variables, we inject weights to each variable used in the random forest.
The weights assigned to each SNP is inversely proportional to the p-values and is plotted below. The rank of each SNP (based on the p-values) is shown in the x-axis. Note that in the cforest implementation we could control the weights of sampling only for the root node. For the subsequent child nodes each covariate has the same chance of getting sampled just like the traditional random forest. One hurdle for the longitudinal data is to decide on the weight for the "Time" covariate. Assigning a heavy weight can lead to trees where "Time" is often selected at the root node, whereas assigning a low weight can end up with trees where is "Time" is never selected for any node. So, we avoided the weighted cforest or Weighted epiMEIF application for longitudinal data and have applied it only for cross-sectional simulated and real datasets.
Max-T test: The proposed strategy for statistically validating the list of interactions using Max-T test is schematically shown in Supplementary Figure  9  and Supplementary Box 2. An example of a plot showing the impact of a 3 SNP epistatic interaction on the aging of the phenotype in the cardiac aging fly dataset is presented in Supplementary Figure  9  Note that the aging trend (Change in Phenotype from Age 1 to Age 4) for the combination of all 2's (minor allele combination) is significantly different from the aging trend for the combination of all zeros or major allele combination (highlighted by dotted blue curve). We are interested in highlighting such SNP combinations in this step. Supplementary Box 2 elaborates on the step-by-step process of validating an epistasis set of N SNPs (generated from MEIF).
Step 1: Randomly draw any N SNPs from the whole genome excluding the causal variants (n times, 1000 in our case). Step 2: For each N-way SNP set compute and evaluate the Max-T test using the following: a:  SNP interaction give rise to 2  genotype combinations denoted by the set  = {1, ⋯ , 2  }. b: Fit the regression: • For longitudinal dataset: (  () =    +  +  ), where  ∈ ,   () denote the phenotype of the population having genotype  at time point t and   capture the effect of the genotype  on the longitudinal variation of the phenotype. • For cross-sectional dataset: (  =   +  +  ), where  ∈ ,   denote the phenotype of the population having genotype  and   capture the effect of the genotype  on the variation of the phenotype. Note that we assume that the genotype combination 1 as the baseline here denoting the combination arising from the interactions of major alleles.
Anova test: To make the statistical validation test more robust we also perform the anova test apart from the Max-T test and claim that an interaction is validated if it is performing well in at least one of the two tests. We perform the anova test for the Age 1 and Aging dataset using the following: Age1 Dataset: For each N-way SNP set compute and evaluate the anova test using the following: • Model1: ~∑  + , where  ∈ {1,. , } ,   are the SNPs in the N-way SNP interaction set. • Model2: ~∑  + ∑  + , where  ∈ {1,. , 2  },   denote the different genotype combination arising from the interaction of N-SNPs. • Anova(Model1, Model2) tests if Model2 is significantly better than Model1. For Aging dataset, the test remains similar with the following modification: • Model1: ~∑:   + , where  ∈ {1,. , } ,   are the SNPs in the N-way SNP interaction set.