A framework for detecting noncoding rare-variant associations of large-scale whole-genome sequencing studies

Abstract

Large-scale whole-genome sequencing studies have enabled analysis of noncoding rare-variant (RV) associations with complex human diseases and traits. Variant-set analysis is a powerful approach to study RV association. However, existing methods have limited ability in analyzing the noncoding genome. We propose a computationally efficient and robust noncoding RV association detection framework, STAARpipeline, to automatically annotate a whole-genome sequencing study and perform flexible noncoding RV association analysis, including gene-centric analysis and fixed window-based and dynamic window-based non-gene-centric analysis by incorporating variant functional annotations. In gene-centric analysis, STAARpipeline uses STAAR to group noncoding variants based on functional categories of genes and incorporate multiple functional annotations. In non-gene-centric analysis, STAARpipeline uses SCANG-STAAR to incorporate dynamic window sizes and multiple functional annotations. We apply STAARpipeline to identify noncoding RV sets associated with four lipid traits in 21,015 discovery samples from the Trans-Omics for Precision Medicine (TOPMed) program and replicate several of them in an additional 9,123 TOPMed samples. We also analyze five non-lipid TOPMed traits.

Extended Data Fig. 1|. Rare variant (MAF < 0.01) distribution in the discovery phase using TOPMed cohorts (n=21,015).

Variant categories are defined by GENCODE VEP categories.

Extended Data Fig. 2|. Manhattan plots and Q-Q plots for unconditional gene-centric noncoding analysis and sliding window analysis of high-density lipoprotein cholesterol (HDL-C) in the discovery phase (n=21,015).

a, Manhattan plots for unconditional gene-centric noncoding analysis of protein-coding gene. The horizontal line indicates a genome-wide STAAR-O P value threshold of $3.57\times {10}^{-7}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/\left(20,000\times 7\right)=3.57\times {10}^{-7}\right)$. Different symbols represent the STAAR-O P value of the protein-coding gene using different functional categories (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS). Promoter_CAGE and promoter_DHS are the promoters with overlap of Cap Analysis of Gene Expression (CAGE) sites and DNase hypersensitivity (DHS) sites for a given gene, respectively. Enhancer_CAGE and enhancer_DHS are the enhancers in GeneHancer predicted regions with the overlap of CAGE sites and DHS sites for a given gene, respectively. b, Quantile-quantile plots for unconditional gene-centric noncoding analysis of protein-coding gene. Different symbols represent the STAAR-O P-value of the gene using different functional categories (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS). c, Manhattan plots for unconditional gene-centric noncoding analysis of ncRNA gene. The horizontal line indicates a genome-wide STAAR-O P value threshold of $2.50\times {10}^{-6}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/20,000=2.50\times {10}^{-6}\right)$. d, Quantile-quantile plots for unconditional gene-centric noncoding analysis of ncRNA gene. e, Manhattan plot for 2-kb sliding windows. The horizontal line indicates a genome-wide P value threshold of $1.88\times {10}^{-8}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/\left(2.66\times {10}^6\right)=1.88\times {10}^{-8}\right)$. f, Quantile-quantile plot for 2-kb sliding windows. In panels, a, c and e, the chromosome number are indicated by the colors of dots. In all panels, STAAR-O is a two-sided test.

Extended Data Fig. 3|. Manhattan plots and Q-Q plots for unconditional gene-centric noncoding analysis and sliding window analysis of low-density lipoprotein cholesterol (LDL-C) in the discovery phase (n=21,015).

a, Manhattan plots for unconditional gene-centric noncoding analysis of protein-coding gene. The horizontal line indicates a genome-wide STAAR-O P-value threshold of $3.57\times {10}^{-7}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/\left(20,000\times 7\right)=3.57\times {10}^{-7}\right)$. Different symbols represent the STAAR-O P-value of the protein-coding gene using different functional categories (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS). Promoter_CAGE and promoter_DHS are the promoters with overlap of Cap Analysis of Gene Expression (CAGE) sites and DNase hypersensitivity (DHS) sites for a given gene, respectively. Enhancer_CAGE and enhancer_DHS are the enhancers in GeneHancer predicted regions with the overlap of CAGE sites and DHS sites for a given gene, respectively. b, Quantile-quantile plots for unconditional gene-centric noncoding analysis of protein-coding gene. Different symbols represent the STAAR-O P-value of the gene using different functional categories (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS). c, Manhattan plots for unconditional gene-centric noncoding analysis of ncRNA gene. The horizontal line indicates a genome-wide STAAR-O P-value threshold of $2.50\times {10}^{-6}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/20,000=2.50\times {10}^{-6}\right)$. d, Quantile-quantile plots for unconditional gene-centric noncoding analysis of ncRNA gene. e, Manhattan plot for 2-kb sliding windows. The horizontal line indicates a genome-wide P-value threshold of $1.88\times {10}^{-8}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/\left(2.66\times {10}^6\right)=1.88\times {10}^{-8}\right)$. f, Quantile-quantile plot for 2-kb sliding windows. In panels, a, c and e, the chromosome number are indicated by the colors of dots. In all panels, STAAR-O is a two-sided test.

Extended Data Fig. 4|. Manhattan plots and Q-Q plots for unconditional gene-centric noncoding analysis and sliding window analysis of triglycerides (TG) in the discovery phase (n=21,015).

a, Manhattan plots for unconditional gene-centric noncoding analysis of protein-coding gene. The horizontal line indicates a genome-wide STAAR-O P-value threshold of $3.57\times {10}^{-7}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/\left(20,000\times 7\right)=3.57\times {10}^{-7}\right)$. Different symbols represent the STAAR-O P-value of the protein-coding gene using different functional categories (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS). Promoter_CAGE and promoter_DHS are the promoters with overlap of Cap Analysis of Gene Expression (CAGE) sites and DNase hypersensitivity (DHS) sites for a given gene, respectively. Enhancer_CAGE and enhancer_DHS are the enhancers in GeneHancer predicted regions with the overlap of CAGE sites and DHS sites for a given gene, respectively. b, Quantile-quantile plots for unconditional gene-centric noncoding analysis of protein-coding gene. Different symbols represent the STAAR-O P-value of the gene using different functional categories (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS). c, Manhattan plots for unconditional gene-centric noncoding analysis of ncRNA gene. The horizontal line indicates a genome-wide STAAR-O P-value threshold of $2.50\times {10}^{-6}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/20,000=2.50\times {10}^{-6}\right)$. d, Quantile-quantile plots for unconditional gene-centric noncoding analysis of ncRNA gene. e, Manhattan plot for 2-kb sliding windows. The horizontal line indicates a genome-wide P-value threshold of $1.88\times {10}^{-8}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/\left(2.66\times {10}^6\right)=1.88\times {10}^{-8}\right)$. f, Quantile-quantile plot for 2-kb sliding windows. In panels, a, c and e, the chromosome number are indicated by the colors of dots. In all panels, STAAR-O is a two-sided test.

Extended Data Fig. 5|. Manhattan plots and Q-Q plots for unconditional gene-centric noncoding analysis and sliding window analysis of total cholesterol (TC) in the discovery phase (n=21,015).

a, Manhattan plots for unconditional gene-centric noncoding analysis of protein-coding gene. The horizontal line indicates a genome-wide STAAR-O P-value threshold of $3.57\times {10}^{-7}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/\left(20,000\times 7\right)=3.57\times {10}^{-7}\right)$. Different symbols represent the STAAR-O P-value of the protein-coding gene using different functional categories (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS). Promoter_CAGE and promoter_DHS are the promoters with overlap of Cap Analysis of Gene Expression (CAGE) sites and DNase hypersensitivity (DHS) sites for a given gene, respectively. Enhancer_CAGE and enhancer_DHS are the enhancers in GeneHancer predicted regions with the overlap of CAGE sites and DHS sites for a given gene, respectively. b, Quantile-quantile plots for unconditional gene-centric noncoding analysis of protein-coding gene. Different symbols represent the STAAR-O P-value of the gene using different functional categories (upstream, downstream, UTR, promoter_CAGE, promoter_DHS, enhancer_CAGE, enhancer_DHS). c, Manhattan plots for unconditional gene-centric noncoding analysis of ncRNA gene. The horizontal line indicates a genome-wide STAAR-O P-value threshold of $2.50\times {10}^{-6}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/20,000=2.50\times {10}^{-6}\right)$. d, Quantile-quantile plots for unconditional gene-centric noncoding analysis of ncRNA gene. e, Manhattan plot for 2-kb sliding windows. The horizontal line indicates a genome-wide P-value threshold of $1.88\times {10}^{-8}$. The significant threshold is defined by multiple comparisons using the Bonferroni correction $\left(0.05/\left(2.66\times {10}^6\right)=1.88\times {10}^{-8}\right)$. f, Quantile-quantile plot for 2-kb sliding windows. In panels, a, c and e, the chromosome number are indicated by the colors of dots. In all panels, STAAR-O is a two-sided test.

Fig. 1 |. Workflow of STAARpipeline.

(a) Prepare the input data of STAARpipeline, including genotypes, phenotypes and covariates. (b) Annotate all variants in the genome using FAVORannotator through FAVOR database and calculate the (sparse) genetic relatedness matrix. (c) Define analysis units in the noncoding genome: eight functional categories of regulatory regions, sliding windows and dynamic windows using SCANG. (d) Obtain genome-wide significant associations and perform analytical follow-up via conditional analysis.