Interpretable network-guided epistasis detection

  1. Diane Duroux
  2. Héctor Climente-González
  3. Chloé-Agathe Azencott
  4. Kristel Van Steen

Reviewed by

Read the full article See related articles

Listed in

Log in to save this article

Abstract

Background

Detecting epistatic interactions at the gene level is essential to understanding the biological mechanisms of complex diseases. Unfortunately, genome-wide interaction association studies involve many statistical challenges that make such detection hard. We propose a multi-step protocol for epistasis detection along the edges of a gene-gene co-function network. Such an approach reduces the number of tests performed and provides interpretable interactions while keeping type I error controlled. Yet, mapping gene interactions into testable single-nucleotide polymorphism (SNP)-interaction hypotheses, as well as computing gene pair association scores from SNP pair ones, is not trivial.

Results

Here we compare 3 SNP-gene mappings (positional overlap, expression quantitative trait loci, and proximity in 3D structure) and use the adaptive truncated product method to compute gene pair scores. This method is non-parametric, does not require a known null distribution, and is fast to compute. We apply multiple variants of this protocol to a genome-wide association study dataset on inflammatory bowel disease. Different configurations produced different results, highlighting that various mechanisms are implicated in inflammatory bowel disease, while at the same time, results overlapped with known disease characteristics. Importantly, the proposed pipeline also differs from a conventional approach where no network is used, showing the potential for additional discoveries when prior biological knowledge is incorporated into epistasis detection.

Article activity feed

  1. GigaScience

    This work has been peer reviewed in GigaScience (see paper https://doi.org/10.1093/gigascience/giab093), which carries out open, named peer-review.

    These reviews are published under a CC-BY 4.0 license and were as follows:

    Reviewer 2: Elisabetta Manduchi

    This manuscript presents a workflow for gene-gene epistasis detection which leverages functional annotation resources such as Biofilter (to reduce the search space) and FUMA (to map SNPs to genes) and investigates the results obtained via different SNP-gene mapping criteria (positional, eQTL, Chromatin contacts, and some combinations of these). Moreover, these results are compared with those obtained via a 'standard' analysis where no filtering is applied to pair selection and positional SNPgene mapping is used. Due to the challenges presented by GWAIS, leveraging functional genomics to focus the search is a valid strategy which has been employed in other recent works in the field. This is a nice work and the paper is generally well written, with sufficiently detailed methodological information. Below are some comments and questions.

    1. As indicated in recent GWAS works aimed at 'solving GWAS loci' (i.e. determining the genes affected by significant SNPs), it is not always the case that the gene affected by a SNP is that positionally closest to the SNP. Indeed, a SNP may not only affect a gene when it resides in its coding or promoter regions, but it may also affect a gene when it resides in a far away enhancer. This is why epigenetic information such as chromatin loops (referred by the authors as 'Chromatin') can be useful for SNP-gene mapping. In the presence of chromatin contacts or eQTL information, typically one would use the derived mapping to augment the positional mapping, which is always available. That is, if one had chromatin contacts data, they would use positional + Chromatin to map SNPs to genes. If one had eQTL data, they would use positional + eQTL to map to genes. If one had both, they would use positional + eQTL + Chromatin. From a biological interpretability perspective, there is no reason to exclude the positional information. For example, a SNP in the promoter of gene could interact with a SNP in a distal enhancer of another gene, affecting a specific trait. In view of this, the statement (lines 326-328) "Since the main objective of this protocol is to increase the biological interpretability of epistasis findings, we have excluded other combinations that mix functional and non-functional information (Positional + eQTL and Positional + Chromatin)" is not quite valid, as positional information is also functional. On the other hand, using eQTL only, Chromatin only, or eQTL + Chromatin, albeit interesting in terms of looking at how this type of reduction in the search space affects results, do not quite reflect a biologically guided approach.

    2. I wonder on whether the authors have considered filtering also by markers of relevant chromatin states. Information about open chromatin and other epigenetic marks could help further filtering SNPs, both in enhancers and promoters. This would be particularly useful for SNPs mapped via chromatin contacts, which are likely to contain many irrelevant signals.

    3. The eQTL and chromatin contact data used in this work were from all available tissues. Typically, GWAS related functional filtering is done using data from tissues relevant to the trait under investigation, when available. For IBD, it may help to restrict to intestinal tissues, immune cells (like T- cells, macrophages, dendritic cells), and possibly also nervous system cells (which, at least according to some, could also be among the potential 'culprit' IBD tissues).

    4. To adjust for population structure the authors regressed out the first 7 PCs from the phenotype. Given that the PCs are confounders, it would be good to discuss the impact of doing this as opposed to also regress the confounders out of the SNPs, i.e. testing the response residuals vs the SNP residuals. In the same spirit, it would be good to discuss the impact of the PC-SNP association on the p-value and type-I error results obtained by permuting the response residuals.

    5. Section 2.1 is somewhat too concise and may result unclear to the reader. Later in the Discussion (lines 229-239) the authors explain how their procedure corrects for multiple testing at the SNP model without additional corrections for multiple testing at the gene model (this is also implicitly described in Fig 7CD), but yet their procedure keeps type I error under control. However, it may be beneficial, for ease of reading, to expand section 2.1 (via text and/or figure) so to clarify better, at the onset, where and when multiple testing corrections are applied.

    6. In the absence of a replication data set, the authors assess the robustness of the gene pair results via 10 repetitions of the workflow using 80% of the discovery data set. It would be useful to include some discussion of how their results could be further assessed in other GWAS data sets (e.g. from UK biobank, etc.), in view of the fact that it is typically hard to reproduce epistasis findings, at least at the SNP level. Certainly one could first check whether the discovered SNP-SNP interactions are reproduced and limiting the analyses to those pairs would require a less severe multiple testing correction. But another approach may be to start with the discovered gene pairs and then analyze all pairs of SNPs mapping to these genes (not necessarily those discovered in this study), etc. Do the authors plan future follow-up studies on this?

    7. In section 2.7 the results of pathway analyses on 3 (eQTL, Positional, and Standard) of the 5 networks presented in Figure 3 are provided. What about the other 2?

    8. For these two points I defer to the editor:

    (i) The format of the manuscript is close to but does not exactly match the specifications at https://academic.oup.com/gigascience//pages/research. I do not know how strict these specifications are and I have no objections to the current format.

    (ii) Data availability is not discussed (as per Data and materials in https://academic.oup.com/gigascience/pages/instructions_to_authors). I imagine that the IIBDGC only makes publicly available the summary statistics. This is, however, common in the GWAS field.

    1. Some minor notes follow:

    (i) In the Author Summary the 'ATPM' acronym is used for the first time without explanation.

    (ii) In section 4.2 it would be helpful to re-iterate that the SNP-gene mapping for the Standard analysis was genomic proximity (this is only mentioned briefly at line 206).

    (iii) Typo at line 168 "the same than" should be "the same as".

    (iv) It should be specified which of the MigSigDB collections was used. Later in this section gene sets are referred to as 'pathways' but there is more than one pathway collection in MigSigDB.

    (v) In the formula at line 397 doesn't "tested gene sets" refer to "tested gene neighborhoods"? If so, it would be better to use the latter for clarity.

    (vi) There appear to be some typos in the caption for Supplementary Figure 1: "we computed three linear models using the different residuals as response variable and SNP interactions as dependent variables". I guess should be "SNP interactions as independent variables". Also, weren't the two individual SNPs also included as independent variables in these models?

  2. GigaScience

    This work has been peer reviewed in GigaScience (see paper https://doi.org/10.1093/gigascience/giab093), which carries out open, named peer-review.

    These reviews are published under a CC-BY 4.0 license and were as follows:

    Reviewer 1: Shing Wan Choi

    Here, the authors presented a pipeline for the analysis of epistasis effects in GWAS data (GWAIS). This is traditionally a difficult problem due to the large search space involved, which makes the analysis computationally intensive and suffers from multiple testing, a new method that restrict search space can definitely help making GWAIS more feasible and reproducible. I have some questions after reading the current paper, please excuse me if the information is already presented within the paper:

    1. I am not sure what the Standard model comprised of. According to the methodology section, the Standard analyzed all SNP pairs without prior filtering, does that mean all 14,501,130,150 SNP pairs (C(170301, 2)) were tested? Or was it not all SNP pairs were considered?

    2. When presenting the number of SNPs linked to each gene based on different criteria (e.g. position, eQTL or Chromatin contact), wouldn't the gene size be a major predictor of the number of SNPs link? It would most likely be the case for positional mapping, right?

    3. I am curious to see if restricting the eQTL and Chromatin information to disease specific tissue will help improving the performance of the current model.

    4. Very little information was provided for the PRS analysis. What genome wide association summary statistics were used? Did the authors perform high-resolution scoring? With shrinkage / thresholding done in normal PRS analysis, some SNPs' effect might be excluded or "shrinked" away from the PRS model, would that affect the interpretation of the PRS covariate analysis? E.g. maybe SNPs not included in the PRS model were those unaffected? (With PRSice, can use --print-snp to obtain list of SNPs that are included in the model)

    5. It seems like Biofilter provide a SNP-SNP interaction prediction model, how does that compare to what was presented here?

    6. Figure 3, the results from eQTL + Chromatin and Position + eQTL + Chromatin is identical. Together, it seems like the positional mapping does not contribute to the result at all, which is a bit surprising. Are there any explanation of that? Would it be due to mapping of the Immunochip array, or a characteristic of IBD?

    7. Given the sample size of the current data, a HWE threshold of 0.001 seems rather stringent. Will the result improve if a less stringent threshold is used (e.g. 1e-6?)

  3. Version published to 10.1093/gigascience/giab093
  4. Version published to 10.1101/2020.09.24.310136v2 on bioRxiv
  5. Version published to 10.1101/2020.09.24.310136v1 on bioRxiv