MicrobiomeGWAS: A Tool for Identifying Host Genetic Variants Associated with Microbiome Composition

Abstract

The microbiome is the collection of all microbial genes and can be investigated by sequencing highly variable regions of 16S ribosomal RNA (rRNA) genes. Evidence suggests that environmental factors and host genetics may interact to impact human microbiome composition. Identifying host genetic variants associated with human microbiome composition not only provides clues for characterizing microbiome variation but also helps to elucidate biological mechanisms of genetic associations, prioritize genetic variants, and improve genetic risk prediction. Since a microbiota functions as a community, it is best characterized by β diversity; that is, a pairwise distance matrix. We develop a statistical framework and a computationally efficient software package, microbiomeGWAS, for identifying host genetic variants associated with microbiome β diversity with or without interacting with an environmental factor. We show that the score statistics have positive skewness and kurtosis due to the dependent nature of the pairwise data, which makes p-value approximations based on asymptotic distributions unacceptably liberal. By correcting for skewness and kurtosis, we develop accurate p-value approximations, whose accuracy was verified by extensive simulations. We exemplify our methods by analyzing a set of 147 genotyped subjects with 16S rRNA microbiome profiles from non-malignant lung tissues. Correcting for skewness and kurtosis eliminated the dramatic deviation in the quantile-quantile plots. We provided preliminary evidence that six established lung cancer risk SNPs were collectively associated with microbiome composition for both unweighted (p = 0.0032) and weighted (p = 0.011) UniFrac distance matrices. In summary, our methods will facilitate analyzing large-scale genome-wide association studies of the human microbiome.

Keywords: gene–environment interaction; genome-wide association study; host genetics; microbiome; skewness and kurtosis; tail probabilities.

Figure 1

Microbiome distances are positively correlated with genetic distances at an associated SNP.

Figure 2

Define the joint test for testing $H_0: {\beta }_M={\beta }_I=0$ vs. ${\beta }_M>0 or {\beta }_I>0$. We assume that $Z_M~N\left(0,1\right)$, $Z_I~N\left(0,1\right)$ and $cor\left(Z_M,Z_I\right)=\rho$ under $H_0$. Details are in Appendix C.

Figure 3

Correcting tail probabilities for skewness and kurtosis. (A) The standard normal distribution $N\left(0,1\right)$ and an approximately normal distribution with positive skewness. The skewness has big impact when calculating the tail probability $P(Z>b)$ for a large value of $b$. (B) Numerical evaluation of tail probability approximation for $Z_M$. We used the unweighted UniFrac distance matrix of 500 samples from the American Gut Project (AGP). For each value of b (>0), we calculated p-values $P(Z_M>b)$ based on $N\left(0,1\right)$, skewness correction, both skewness and kurtosis correction, and 10⁸ simulations. (C) Skewness depends on minor allele frequency (MAF) of SNPs and the sample size of the study, calculated based on the weighted UniFrac distance matrix in AGP data. (D) Kurtosis depends on MAF of SNPs and the sample size, calculated based on the weighted UniFrac distance matrix in the AGP data.

Figure 4

Computation time for a microbiome GWAS with 500,000 SNPs. “Main”: computation time for testing main effect only. “All”: computation time for testing main effect, interaction and the joint null hypothesis $H_0: {\beta }_M=0,{\beta }_I=0$.

Figure 5

Results of analyzing the microbiome GWAS data of 147 adjacent normal lung tissues in the EAGLE study. (A) Skewness and kurtosis for the main effect test using the unweighted and the weighted UniFrac distance matrices. (B) Quantile–quantile (QQ) plot for association p-values using the unweighted UniFrac distance matrix. “Adjusted”: p-values were corrected for skewness and kurtosis. “Unadjusted”: p-values were approximated based on the asymptotic distribution $N\left(0,1\right)$. (C) Quantile–quantile (QQ) plot for association p-values using the weighted UniFrac distance matrix. (D) Manhattan plots based on the unweighted or the weighted UniFrac distance matrices. (E) Box plots for the top nine loci in microbiome GWAS analysis. Subject pairs are classified into three groups according to the genetic distance $\left|g_i-g_j\right|$ at the SNP. The y-coordinate is the microbiome distance.