{'fraction_non_alphanumeric': 0.0399476708872588, 'fraction_numerical': 0.02200626080456011, 'mean_word_length': 4.718407694362811, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Enhancing power of rare variant association test by Zoom-Focus Algorithm (ZFA) to locate optimal testing region

Motivation: Exome or targeted sequencing data exerts analytical challenge to test single nucleotide polymorphisms (SNPs) with extremely small minor allele frequency (MAF). Various rare variant tests were proposed to increase power by aggregating SNPs within a fixed genomic region, such as a gene or pathway. However, a gene could contain from several to thousands of markers, and not all of them may be related to the phenotype. Combining functional and non-functional SNPs in arbitrary genomic region could impair the testing power. Results: We propose a Zoom-Focus algorithm (ZFA) to locate the optimal testing region within a given genomic region, as a wrapper function to be applied in conjunction with rare variant association tests. In the first Zooming step, a given genomic region is partitioned by order of two, and the best partition is located within all partition levels. In the next Focusing step, boundaries of the zoomed region are refined. Simulation studies showed that ZFA substantially enhanced the statistical power of rare variant tests by over 10 folds, including the WSS, SKAT and W-test. The algorithm is applied on real exome sequencing data of hypertensive disorder, and identified biologically relevant genetic markers to metabolic disorder that are undiscoverable by testing using full gene. The proposed algorithm is an efficient and powerful tool to increase the effectiveness of rare variant association tests for exome sequencing datasets of complex disorder.


INTRODUCTION
The increasing depth and coverage of DNA sequencing technology have made possible the accumulation of human whole genome sequencing data. National-level efforts are undergone to collect big cohort deep sequencing data towards better knowledge of Mendelian and complex disorders (Ashley, 2015;Auffray, et al., 2016;Cyranoski, 2016;Jameson and Longo, 2015). These efforts press the demand of powerful and efficient methods for medical and clinical inference from the data. The challenge in analyzing the exome sequencing data set, besides the multiple-testing issue arising from the high dimensionality, centers on the extreme low allele frequencies -Classical statistical tests lose power on single nucleotide polymorphisms (SNPs) with small variances. In whole exome sequencing data, over 99% of the SNPs have minor allele frequency (MAF) below 1% (Consortium, 2015). A number of rare variant association tests have been proposed, which either pull adjacent SNPs together to increase minor allele frequencies and upweight the minor allele (burden tests) (Li and Leal, 2008;Liu and Leal, 2010;Madsen and Browning, 2009), or apply a linear mixed model on a certain genomic region to improve power (variance component tests) (Lee, et al., 2012;Neale, et al., 2011;Wu, et al., 2011). For most rare variant methods, a fixed genomic region for testing is assumed, such as a gene or a fixed window. However, in real data application, a gene may contain from several to thousands of SNPs. Directly applying rare variant association tests based on a gene or a fixed window may introduce a large number of un-necessary noises that could impair testing power (Santorico and Hendricks, 2016). Furthermore, many of the exome sequencing data have unknown gene functions and are difficult to pool without prior information. It is desirable to optimize the collapsing region such that minimum noise is included and power of the test can be enhanced.
In this paper, we propose a Zoom-Focus algorithm (ZFA) to optimize testing region for rare variant association tests. The algorithm can be applied as a wrapper function to existing tests and can increase the test power by over 10-fold in simulations studies. The algorithm takes two main steps, Zooming and Focusing. In the Zooming step, a fix genomic region is partitioned by an order of two, and search is conducted through all partition levels to identify the region with maximum information, evaluated by the smallest association p-value of that partition. Based on the best partition identified by Zooming, a next Focusing step is used to optimize the region's boundary, by adding or reducing adjacent SNPs around the boundaries. Simulation studies for various genetic scenarios demonstrated that the Zoom-Focus algorithm could substantially enhance the statistical power of different rare variant methods, including the WSS, SKAT and W-test. The ZFA is applied on real exome sequencing data of hypertensive disorder and identified biologically relevant genetic markers that are undiscoverable by un-optimized testing regions.

The Zoom-Focus Algorithm (ZFA)
The Zoom-Focus Algorithm is best described through an example. Suppose a gene or fixed windows contains 64 SNPs, among which 8 SNPs are causal. For the simplest scenario, the 8 causal markers cluster together, as shown in the red region in Diagram 1. (Other causal marker distributions are provided in the simulation study.) The ZFA first performs an exhaustive search in all possible binary partitions of the initial region to locate the best partition (Zooming), and then adjust the boundary of the zoomed region by considering increment or decrement of the bounds (Focusing) (Diagram 1). Diagram 1. The Zoom-Focus Algorithm (ZFA) Legend: The ZFA first performs an exhaustive search in all possible binary partitions in a given genomic region (Step1: Zooming); then adjust the boundary of the zoomed region by considering increment or decrement of the bounds (Step2: Focusing).
We will first define some notations, and then introduce the algorithm. Let P be the total number of SNPs in a fixed window; R is the maximum order of binary partition for P, and R= arg max{r: 2 r-1 ≤ P, r = 1, 2, …}; r is the order of partition. At a certain r, the number of partitioned regions is nr = 2 r-1 , and the size of partition d is the number of SNPs in a partitioned region, d =P/nr =P/2 r-1 . A higher order gives a smaller partitioned region size. c is the index for the c th partition, c = 1,…, nr. For the causal region in Diagram 1: P = 64, the causal region is located at 4 th order of binary division and the 2 nd partitioned region (c=2), thus r = 4, and d =8. We denote ϕ(d, c) = {xcj, j=1,…,d } as the region that contains SNPs in a partition size d and index c, and the causal region can be denoted as ϕ (8,2). Let f(•) be the Bonferroni corrected p-value calculated by rare variant method F(•) measured on region ϕ(d, c): where nr is the number of partitions for Bonferroni correction of multiple tests. It ensures p-values at the different partition order can be compared. Therefore, we have: Step 1: Zooming. Search all binary partitions of a given initial region, for an optimal partition of size d and index c: ϕ (d, c), such that f(d, c) is minimized: (d̂, ĉ) = arg min{ f(d, c), r = 0,…,R; c = 1 to nr } Eq 2.
Step 2: Focusing. Refine the boundary of (d̂, ĉ) by extending both lower and upper bounds, outwardly by d/2, and inwardly by d/4. Let LB denotes the lower bound of (d̂, ĉ), and UB denotes the upper bound of (d̂, ĉ). The refined lower bound (LBf) and upper bound (UBf) are: In sum, the Zooming step locates the best partition from global search, and the Focusing step refines the boundary of the zoomed region from local linear optimization. The terms are motivated from optical zooming lens.

Computation complexity
In the Zooming step, the maximum number of tests T(P) in a gene with P number of SNPs is: In the Focusing step, the worst-case scenario contains P/2 calculations. Thus the overall computation complexity of the ZFA is O(P). This is much efficient than searching all possible collapsing window size, by which the computation complexity is O(P!).

An alternative fast-Zoom method
The Zooming algorithm searches all binary partitions of a given region exhaustively. The computing speed of different statistical tests varies: some inherit probability distributions such as the W-test (Wang, et al., 2016), while some incorporates permutations to obtain the p-values such as the SKAT (Wu, et al., 2011). To assist the computation extensive methods to perform Zoom-Focus, we propose a fast-Zoom method: Instead of searching all possible partitions, the fast-Zoom performs a binary search, such that at each partition order r, the region is divided into two parts, only the part with smaller p-value is continued for the next level search, as illustrated in Diagram 2. The computation complexity of fast ZFA reduces to O(log2(P)).

Diagram 2. Fast-Zoom-Focus Algorithm (fast-ZFA)
Legend: In fast-Zoom, a binary search replaces the exhaustive search -only green colored regions are evaluated. Fast-ZFA's computation complexity is O(log2(P)), compared to ZFA's O(P).

Simulation study design
Rare variant association test's power is largely influenced by the distribution of causal variants in the evaluation region (Sham and Purcell, 2014). To fully explore all methods' performance with ZFA, three scenarios are considered, with varying distributions and different effect directions (Diagram 3): Scenario I: 8 causal SNPs cluster together in same effect direction Scenario II: 8 causal SNPs cluster in two groups in same effect direction Scenario III: 8 causal SNPs cluster together in opposite effect directions (5 SNPs show risk effect and 3 SNPs show protective effect)

Diagram 3. Causal variants distribution scenarios in simulation study
To evaluate the effect of Focusing, distribution of causal variants are simulated to carry the setting of Scenario I but with uneven boundaries. Receiver Operating Curves (ROCs) are used to compare the Zooming and full ZFA performances.
Each simulated data contains 15,360 SNPs, 1,000 cases and 1,000 controls. The data is divided into 120 regions of P=128 SNPs. In each dataset, 10 regions are randomly selected to include causal variants according to one of the three scenarios, and the rest 110 regions are noise. The proportion of functional markers in a causal gene is 6.25% (8/128), and the proportion of causal genes in entire dataset is 8.33% (10/120). Phenotypes are generated using a logistic regression model of the causal SNPs (Wu, et al., 2011): where X1 is a standard normal covariate, X2 is a dichotomous covariate that takes the value 0 with probability 0.5 and the value 1 otherwise. Gi are causal rare variants and βi are the effect size, βi=|log10MAF|×0.3, such that rarer variants have larger effects. We controlled prevalence by β0 and set it to 10% unless otherwise stated.

Statistical tests considered
The Zoom-Focus algorithm is applicable to all region-based rare variant association tests. Three representative methods are selected to combine with ZFA in the simulation studies, including the SKAT of variance component test; the WSS and the W-test of the burden test category. The SKAT is a quadratic score test that is composed of a weighted prediction error from a linear mixed model; it follows a mixture of Chi-squared distributions, and the p-value is obtained through Davies approximation (Wu, et al., 2011). The SKAT is advantageous when a few large effect variants are located in a genomic region, and when the effect directions are not identical (Lee, et al., 2014). The WSS uses Wilcoxon's rank sum statistic to measure a weighted count of the excess of mutations in the cases from the control, and use permutations to calculate p-values (Madsen and Browning, 2009). The WSS method is suitable for the scenario when the majority of causal variants have the same effect direction (Sha, et al., 2012). The W-test is a fast and model-free genetic association test (Sun, et al., 2016;Wang, et al., 2016). It tests the distributional differences of SNPs in the case group from that in the control group, and follows a Chi-squared distribution with dataset adaptive degrees of freedom. For rare variant application, the W-test is calculated on the summed contingency table of SNPs from a given region (Sun, et al., 2016). The WSS test and SKAT are computationally extensive and are applied with the fast-ZFA. The W-test is applied with both the fast and full ZFA.

Power, type I error calculation and ROC
One hundred replicated datasets are generated for each genetic scenario. In each replicate, 120 fixed regions are evaluated by alternative methods. A region is significant if it locates an optimal region with p-value smaller than 4.2×10 -4 (0.05/120), which is the Bonferroni corrected p-value at a family-wise error rate (FWER) of 0.05. An outcome is a true positive if the final optimal region overlaps the causal region. Power is the averaged true positive proportion from 100 simulated datasets, and the type I error rate is the average false positive rate of the 100 simulations with permuted phenotypes. Receiver Operating Characteristic (ROC) curves are used to assess the improvement of the Focusing step after the Zooming step.

Real data application
ZFA is applied on a real hypertensive disorder sequence data of the Genetics Analysis Workshop 19 (GAW19). The data consisted of 398 hypertensive patients and 1,453 healthy controls, and exome sequence of chromosome 3 is used. Variants with missing value percentage over 5%, MAF>1%, and inconsistent genotyping format are excluded (Laurie, et al., 2010). After quality control, 41,788 rare variants are remained. Full ZFA is applied on Chromosome 3 using an initial window size P = 256. The chromosome is then divided into 163 non-overlapping regions, and the remaining 60 SNPs are grouped as the last window. To improve the computation efficiency, zoomed region with p-value smaller than 6.1×10 -4 (0.1/164) are passed to the Focusing step. A genetic region will be regarded as significant if its final p-value reach the Bonferroni corrected significance level of 3.0×10 -4 (0.05/164). Initial window sizes P=128 and P=512 are also applied, with significance threshold adjusted accordingly. Genes that include or overlap with the returned region are reported as susceptible.

Simulation study
Fast-Zoom and full path Zoom are evaluated under three causal marker distribution scenarios ( The fast-Zoom loses certain power compared to full path zooming, however, it is efficient and useful for tests that require permutation to generate p-values. The results showed that Zooming could generally locate the causal region and greatly improve rare variant testing power, either using the fast or full path. The type I errors are reasonably controlled (Table 1). For fast-Zoom, the SKAT's type I error rate was 0.1%, WSS's 1.3%, and the W-test's 4.8%. For full path Zoom, the type I error of W-test was 6.6%, reasonable in view of its power gain. Figure 1 compared the ROC of full-path Zoom to the algorithm with Focusing step using W-test, simulated from scenario I. The ROCs showed that by applying the Zooming step alone, there was already considerable power, and the Focusing step further boosted the overall performance. The binary partition of Zooming resulted in discontinuous optimal region, and the true and false positive rates are returned in intervals, thus the ROC showed some form of a step function.

Computation time
On a PC cluster with 8 compute nodes each with 16GB RAM (2GB per core) and 2.66GHz E5430 CPU, the unit time elapsed for fast-Zoom in a region of 128 SNPs is 97.63 seconds(s) by WSS, 295.57s by SKAT, and 0.18s by W-test (Table 2). ZFA with Wtest is several hundred times faster than with WSS and SKAT. Full path Zoom by the Wtest took 2.28s. In real data, full ZFA with W-test took approximately 15 minutes to complete the whole exome evaluation of chromosome 3.

Real data application
On whole exome sequencing of chromosome 3 of hypertensive disorder, ZFA with W-test identified two significant regions, one contains 29 SNPs and the other has 27 SNPs (Table 3). Different initial window sizes are applied and gave the same optimal region (Supplementary Materials S1). This showed that the method is quite robust. The two regions are parts of the genes ITPR1 and MSL2/ PCCB (http://www.ncbi.nlm.nih.gov). Both genes were previously reported to associate with cardiovascular and metabolic traits (Consortium, 2013;Dehghan, et al., 2009;Kathiresan, et al., 2007;Sabater-Lleal, et al., 2013). The gene ITPR1 spans 251 SNPs, and MSL2 covers 85 SNPs (Table 3). We also calculated the p-values of these two genes using gene-based WSS, SKAT and W-test without region optimization, and none of them are significant) (Figure 2). The results showed that ZFA could greatly enhance the effective selection of rare variants with disease association on exome sequencing data. The phenotypes are also permuted for 100 times to evaluate possible false positives, and the average false positive rate is 1.7%, which indicated that there is no inflation of spurious association.

DISCUSSION
The search for rare variant's contribution of complex disease was triggered by the hypothesis that a few SNPs may account for large effect sizes (Manolio, et al., 2009). However, increasing evidences have shown that rare or low frequency genetic variants have modest-to-weak effect sizes (Auer and Lettre, 2015). These findings suggest that for rare variant association study, it is especially important to consider the effect of multiple loci jointly. Furthermore, functional annotation of rare variants is still evolving; there are multiple platforms to determine the actual regions of genes, some with overlapping boundaries Kircher, et al., 2014;Maurano, et al., 2012). There is also no consensus of whether the promoters, un-translated regions or introns should be included in gene-base tests (Auer and Lettre, 2015). Even for the coding variants, it is likely that only a part of them is functional and others represent random genetic variation (Ionita-Laza, et al., 2014). The problem is more serious with the increasing sequencing deptha common gene easily spans two to three thousand SNPs. The power of direct application of aggregation tests would be affected from the large proportion of possible noises, thus it is crucial to determine the region of testing as a starting point. The motivation of ZFA is to locate an optimal region with maximum association information and to exclude noise, based on the information of the data. In fact, the Zooming step is performing a feature selection within a given genomic region, and the features are SNP-clusters of varying sizes. Based on the best partition, boundaries are further refined by adding or subtracting adjacent SNPs surrounding them. In this way, noise can be discarded and statistical power can be improved. Simulation studies demonstrated that the ZFA can boost power of various rare variant tests by over ten-fold.
In real data analysis of one chromosome by ZFA, only two regions are significant and they belong to the protein coding genes ITPR1 and MSL2/PCCB. Both genes were previously found to have association with cardiovascular and metabolic disorders in independent studies: the ITPR1 was reported in the Framingham Heart Study to be associated with cholesterol levels (Kathiresan, et al., 2007). The MSL2/PCCB was found to have a strong association with fibrinogen levels, coronary artery disease, blood pressure, and body mass index (Consortium, 2013;Dehghan, et al., 2009;Sabater-Lleal, et al., 2013). In this dataset, neither of the two genes was significant by gene-based tests, but they are identified through optimizing the testing region. Specifically, the ITPR1 encodes an intracellular receptor for inositol 1,4,5-trisphosphate. It plays a role in mediating calcium release from the endoplasmic reticulum. The MSL2 is a subunit of a protein complex that functions in chromatin modification, and the protein encoded by PCCB is a subunit of the propionyl-CoA carboxylase (PCC) enzyme. Application of ZFA with rare variant aggregation tests increased the chance of identifying disease susceptible loci. More experiments are needed to confirm the discovered regions' underlying biological functions related to disease.
To conclude, we propose a Zoom-Focus Algorithm (ZFA) as an efficient, robust and powerful approach to enhance statistical power of rare variant association tests by locating the optimal testing region. The algorithm is flexible to be applied with various rare variant methods and scalable for whole genome analysis. It is a practical and powerful approach in exome sequencing data analysis towards elucidating the role of rare variant in complex disorder.

S1 Regions detected with different initial window sizes in real data application
When the initial window size P = 128, chromosome 3 is divided into 326 non-overlapping region, and the remaining 60 SNPs are grouped as the last window. Zoomed region with p-value smaller than 3.1×10-4 (0.1/327) are passed to the Focusing step. A genetic region will be regarded as significant if its final p-value reach the Bonferroni corrected significance level of 1.5×10-4 (0.05/327). When the initial window size P = 512, chromosome 3 is divided into 81 non-overlapping region, and the remaining 316 SNPs are grouped as the last window. Zoomed region with p-value smaller than 1.2×10-3 (0.1/82) are passed to the Focusing step. A genetic region will be regarded as significant if its final p-value reach the Bonferroni corrected significance level of 6.1×10-4 (0.05/82).