An efficient algorithm to perform multiple testing in epistasis screening

Abstract

Background: Research in epistasis or gene-gene interaction detection for human complex traits has grown over the last few years. It has been marked by promising methodological developments, improved translation efforts of statistical epistasis to biological epistasis and attempts to integrate different omics information sources into the epistasis screening to enhance power. The quest for gene-gene interactions poses severe multiple-testing problems. In this context, the maxT algorithm is one technique to control the false-positive rate. However, the memory needed by this algorithm rises linearly with the amount of hypothesis tests. Gene-gene interaction studies will require a memory proportional to the squared number of SNPs. A genome-wide epistasis search would therefore require terabytes of memory. Hence, cache problems are likely to occur, increasing the computation time. In this work we present a new version of maxT, requiring an amount of memory independent from the number of genetic effects to be investigated. This algorithm was implemented in C++ in our epistasis screening software MBMDR-3.0.3. We evaluate the new implementation in terms of memory efficiency and speed using simulated data. The software is illustrated on real-life data for Crohn's disease.

Results: In the case of a binary (affected/unaffected) trait, the parallel workflow of MBMDR-3.0.3 analyzes all gene-gene interactions with a dataset of 100,000 SNPs typed on 1000 individuals within 4 days and 9 hours, using 999 permutations of the trait to assess statistical significance, on a cluster composed of 10 blades, containing each four Quad-Core AMD Opteron(tm) Processor 2352 2.1 GHz. In the case of a continuous trait, a similar run takes 9 days. Our program found 14 SNP-SNP interactions with a multiple-testing corrected p-value of less than 0.05 on real-life Crohn's disease (CD) data.

Conclusions: Our software is the first implementation of the MB-MDR methodology able to solve large-scale SNP-SNP interactions problems within a few days, without using much memory, while adequately controlling the type I error rates. A new implementation to reach genome-wide epistasis screening is under construction. In the context of Crohn's disease, MBMDR-3.0.3 could identify epistasis involving regions that are well known in the field and could be explained from a biological point of view. This demonstrates the power of our software to find relevant phenotype-genotype higher-order associations.

Figure 1

Input/output formats of MBMDR-3.0.3. MBMDR-3.0.3 takes as argument a text file (possibly converted by our software from PLINK format) containing the trait and SNP values of the subjects under study and a set of command line parameters. If the a^th subject is a case (control), c_a=1(0) (a=1…s). SNP_b is a label referring to the b^th SNP (b=1,…M). The genotype of an individual a at locus b is denoted as g_ab (0 if homozygous for the first allele, 1 if heterozygous and 2 if homozygous for the second allele). The produced output is a text file containing the most significant SNP pairs in relation with the trait. (SNP_lj,SNP_rj) refers to the j^th best SNP pair, i.e. the pair with the j^th lowest p-value p_j. Our software has only one mandatory argument: the scale of the trait. Use either −−binary for a binary trait, or −−continuous for a continuous scale, or −−survival for a censored trait (in this case the trait column is replaced by two columns, one for the time variable and one for the censoring variable). We have developed an interactive help, accessible through −−help, describing all other options. For instance, -n sets the amount of p-values to compute (default: 1000), -p sets the amount of permutations to asses statistical significance (default: 999).

Figure 2

Classical versus Van Lishout’s implementation of maxT. In the classical maxT implementation, all T_i,j values are in memory. If only the x best p-values are envisaged then only the maximum M₁,…,M_B of the [T_1,n+1,…,T_1,m],…,[T_B,n+1,…,T_B,m] are needed, implying only temporary storage of the corresponding values.

Figure 3

MBMDR-3.0.3 parallel workflow. Step 1 of the maxT algorithm is first performed on the input file. This produce the file topfile.txt, containing the top pairs of SNPs and their corresponding test-statistics. Then, the computation of the permutations is split between the available machines. Finally, MBMDR-3.0.3 reads the produced permutation_x.txt files to create the final output file.

Figure 4

Decomposition of the different steps of the computation of T_i,j. c_a is 1 (0) if the a^th subject is a case (control) for the i^th permutation of the trait. g_alj and g_arj are 0, 1 or 2 depending on the genotype of the a^th subject for the j^th pair. A_mn and U_mn are respectively the number of affected/unaffected subjects, whose genotype g_kl= m and g_kr= n. R_mn is either “H” if the subjects whose genotype is m for SNP_lj and n for SNP_rj have a high statistical risk of disease, “L” if they have a low statistical risk and “O” if there is no statistical evidence.

Figure 5

SD plot. Synergy Disequilibrium (SD) plot of potential epistasis interactions between the loci indicated in Table 3