An Efficient Sufficient Dimension Reduction Method for Identifying
      Genetic Variants of Clinical Significance
Authors: Momiao Xiong, Long Ma
ICLR 2013

Abstract: Fast and cheaper next generation sequencing technologies will generate unprecedentedly massive and highly-dimensional genomic and epigenomic variation data. In the near future, a routine part of medical record will include the sequenced genomes. A fundamental question is how to efficiently extract genomic and epigenomic variants of clinical utility which will provide information for optimal wellness and interference strategies. Traditional paradigm for identifying variants of clinical validity is to test association of the variants. However, significantly associated genetic variants may or may not be usefulness for diagnosis and prognosis of diseases. Alternative to association studies for finding genetic variants of predictive utility is to systematically search variants that contain sufficient information for phenotype prediction. To achieve this, we introduce concepts of sufficient dimension reduction and coordinate hypothesis which project the original high dimensional data to very low dimensional space while preserving all information on response phenotypes. We then formulate clinically significant genetic variant discovery problem into sparse SDR problem and develop algorithms that can select significant genetic variants from up to or even ten millions of predictors with the aid of dividing SDR for whole genome into a number of subSDR problems defined for genomic regions. The sparse SDR is in turn formulated as sparse optimal scoring problem, but with penalty which can remove row vectors from the basis matrix. To speed up computation, we develop the modified alternating direction method for multipliers to solve the sparse optimal scoring problem which can easily be implemented in parallel. To illustrate its application, the proposed method is applied to simulation data and the NHLBI's Exome Sequencing Project dataset

An Efficient Sufficient Dimension Reduction Method for Identifying 
Genetic Variants of Clinical Significance 
Momiao Xiong          Long Ma 
Division of Biostatistics                                                        Division of Biostatistics                        
The University of Texas School of Public Health                         The University of Texas School of Public Health 
Houston, TX 77030                                                                 Houston, TX 77030 
momiao.xiong@uth.tmc.edu                                                      marlone.zj@gmail.com 
Abstract 
Fast and cheaper next generation sequencing (NGS) technologies will generate unprecedentedly massive (thousands 
or even ten thousands of individuals) and highly-dimensional (up to hundreds of millions) genomic and epigenomic 
variation data. In the near future, a routine part of medical record will include the sequenced genomes.  A 
fundamental question is how to efficiently extract genomic and epigenomic variants of clinical utility which will 
provide information for optimal wellness and interference strategies. Traditional paradigm for identifying variants of 
clinical validity is to test association of the variants. However, significantly associated genetic variants may or may 
not be usefulness for diagnosis and prognosis of diseases. Alternative to association studies for finding genetic 
variants of predictive utility is to systematically search variants that contain sufficient information for phenotype 
prediction. To achieve this, we introduce concepts of sufficient dimension reduction (SDR) and coordinate 
hypothesis which project the original high dimensional data to very low dimensional space while preserving all 
information on response phenotypes. We then formulate clinically significant genetic variant discovery problem into 
sparse SDR problem and develop algorithms that can select significant genetic variants from up  to or even ten 
millions of predictors with the aid of dividing SDR for whole genome into a number of sub -SDR problems defined 
for genomic regions.  The sparse SDR is in turn formulated as sparse optimal scoring problem, but with penalty 
which can remove row vectors from the basis matrix. To speed up computation, we develop the modified alternating 
direction method for multipliers to solve the sparse optimal scoring problem which can easily be implemented in 
parallel. To illustrate its application, the proposed method is applied to simulation data and the NHLBI’s Exome 
Sequencing Project (ESP) dataset as well as the TCGA dataset. 
Introduction 
Purpose of this paper is to formulate clinically significant genetic variant discovery problem into sparse SDR 
problem and develop algorithms that can select significant genetic variants from up to millions of predictors. To 
achieve this, we first show that SDR for whole genome can be partitioned into a number of sub -SDR problems 
defined for divided genomic regions. Then, similar to Wang and Zhu’s approach, we formulate th e sparse SDR into 
sparse optimal scoring problem, but with penalty which can remove row vectors from the basis matrix. Since large -
scale discovery of genetic variants may involve millions of genetic variants, solving large sparse optimal scoring 
problem requires heavy computation. To speed up computation, we apply the alternating direction method for 
multipliers which can easily be implemented in parallel. To illustrate its application, the proposed method is applied 
to simulation data and the NHLBI’s Exome Sequencing Project (ESP) dataset and TCGA dataset. 
Methods 
Sufficient Dimension Reduction 
    Throughout the paper we consider continuous phenotype (response variable) and regression. In other words, we 
will focus on quantitative trait analysis. However, all discussed concepts can be extended to binary response variable 
and classification.  Let 
Y be a univariate response variable (phenotype) and 
X be a 
p dimensional vector of 
predictors (genotypes for genetic variants). Since dimension of genomic variation is extremely high , to reduce the 
impact of noise and irrelevant predictors, dimension reduction is a powerful tool for quantitative trait analysis and 
regression. Dimension reduction is to identify the best linear subspace that that best preserves information relevant 
to a regression (Nilsson et al. 2007). Dimension reduction consists of unsupervised dimension reduction and 
supervised dimension reduction. Principal component analysis (PCA) is a typical method for unsupervised 
dimension reduction which projects predictor data onto a linear space without response information.  Supervised 
dimension reduction is to discover the best subspace that maximally reduces the dimension of the input while 
preserving the information necessary to predict the response variable. The current popular supervised dimension 
reduction method is SDR which aims to find a linear subspace S such that the response Y is conditionally 
independent of the covariate vector X, given the projection of X on S: 
XPXY S|
,          (1) 
where 
 indicates independence and 
SP represents a projection on 
S . In other words, all information of 
X about 
Y
is contained in the space
S . The subspace S is referred to as a dimension reduction subspace. The subspace 
S
may not be unique. To uniquely describe dimension reduction subspace, we introduce central subspace (CS) that is 
defined the intersections of all reduction subspaces 
S satisfying conditional independence assertion. The CS is 
denoted by 
XYS | .  
Many methods have been developed for identifying CS. A popular sliced inverse regression for identifying the 
basis vector in the CS is to solve the following eigenequation:  
 xYXEXE ))|)((cov(
,        (2) 
where 
x are eigenvalues, and  
 is an eigenvector, respectively. Solutions to eigenequation (2) yields the basis 
matrices  
],...,[ 1 kB  for 
XYS | .  
Sparse SDR by Alternative Direction Method of Multipliers 
The eigenvalue problem can also been formulated as a constrained optimization problem (Chen and Li 1998; 
Wang and Zhu 2013): 
,,...,1,1,...,1,0,1
||||min 2
,
diijYYYY
XY
j
TT
ii
TT
i
iiRbR P
i
k
i



    t. s.
                 (3) 
where 
T
nyyY )](),...,([ 1  .  
To develop sparse SDR that can simultaneously reduce the dimension and the number of predictors, we first 
introduce a coordinated-independent penalty function. Let 
diB T
p
T
i ,...,1],,...,[],...,[ **
11  be a 
ip 
matrix which forms the basis matrix of the CS. We introduce the following penalty function to penalize the variable 
in all reduction directions toward zero (Chen et al. 2010): 



p
l
ll
p
l
l
T
llB
1
2
*
1
** |||)( . For simplicity of 
computation, we define 
r
ll
 2
* |||| , 
0r is a pre-specified parameter. We can chose 
5.0r as  Chen (2010) 
suggested.  
After introducing the penalty function, the sparse version of optimal scoring problem (3) for penalizing the 
variable can be defined as 
,,...,1,,0,1
||||||||min
1
1
2
*2
2,
diijYYYY
XY
j
TT
ii
TT
i
p
l
r
liiRR P
i
k
i

 



    t. s.
      
This is a bi-convex problem. It is convex in 
 for each 
 and convex in 
 for each 
 . It can be solved by a simple 
iterative algorithm.  The iterative process consists of two steps: (1) for fixed 
i we optimize with respect to 
i and 
for fixed 
i we optimize with respect to 
i . The algorithms are given bellow. 
Step 1: Initialization. Let 
nYYD T /  and  
TQ ]0,...,0,1[1  . We first initialize for 
dii ,...,1,)0(  : 
*
)0( )(~  DQQI T
iii
, 
,~~
~
)0()0(
)0(
)0(
i
T
i
i
i
D

],:[1 iii QQ    ,where 
*   is a random 
k -  vector. 
Step 2: Iterate between 
)()( ss   and  until convergence or until a specified maximum number of iterations 
(s=1,2,…)is reached: 
    Step A: For fixed 
,,...,1,)1( dis
i   we solve the following minimization problem: 





p
l
r
l
s
i
s
iR
X
Ps
i 1
1
2
*2
2
)()1( ||||||min
)(
Y||  
d
1i
,where 
)*()*(
1
)()(
1
)( ,...,[],...,[ s
p
ss
d
ssB  .  (4a)    
Step B: For fixed 
,,...,1,)( dis
i  we seek 
,,...,1,)( dis
i  which solve the following unconstrained optimization 
problem: 
  
.1,...,1,0,1
||||min
)()()()(
2
2
)()(
)(



ijYYYY
XY
s
j
TTs
i
s
i
TTs
i
s
i
s
iRks
i
    t. s.   
Solution to the above optimization leads to a nonlinear equation: 
.)(
)()(
)(1)()(
)(
s
i
TsT
i
s
i
TTs
i
s
is
i
XYn
XYDDQQI



     (4b) 
By Newton method, we obtain a solution
)(~ s
i  to equation (4b) . Set  
)()(
)(
)(
~~
~
s
i
sT
i
s
is
i
D
 .    
Results 
To illustrate  its application to selection of clinically useful genetic variants, the proposed method was  first  
applied to several simulation datasets.  Fourier series were used as the transformation function. The number of 
Fourier function in the simulations was 30.  We used the true positive rate (TPR), defined as the proportion of the 
correctly identified predictors, to measure how well the method selects the predictors. We considered two scenarios: 
50 SNPs and 100 SNPs  SNPs from chromosome 1 in the NHLBI’s Exome Sequencing Project (ESP) dataset with 
5,406 individuals and 1,779,016 SNPs.  The simulation models were given by 



p
l
ijiji bxy
1
1.0
, where 
i is distributed as a standard normal distribution 
).1,0(N  
The results were summarized in Table 1.  
Table 1. True positive rates for two simulated data. 
 Dataset Sample Size Number of SNPs Number of Causal SNPs TPR 
1 5406 50 10 100% 
2 5406 100 10 100% 
 
To further evaluate its performance, the proposed method was also applied to the real NHLBI’s Exome Sequencing 
Project (ESP) dataset with HDL phenotype. We discovered 863 SNPs that contribute the HDL variation. The top 10 
selected SNPs were listed in Table 2.  It is reported that the gene cholesteryl ester transfer protein, plasma (CETP) is 
associated with HDL (Braun  et al. 2012), CD36 is associated with atherosclerotic cardiovascular diseases (Yuasa-
Kawase et al. 2012),  
Table 2. Top 10 selected SNPs  
     CHR SNP Gene P-value CHR SNP Gene P-value 
chr16 rs34065661 CETP 9.624E-17 chr19 rs1052983 LILRA6 1.57313E-07 
chr7 rs3211938 CD36 4.629E-13 chr7 rs10085732 1.66681E-07 
chr6 rs17622 DDO 4.333E-08 chr19 rs1868953 3.45063E-07 
chr8 rs111855567 5.861E-08 chr19 rs117156027 3.45063E-07 
chr1 rs12088246 PTGFR 1.512E-07 chr1 rs79907831 SPOCD1 1.15888E-06 
 
References 
Nilsson J, Sha F, and Jordan MI. (2007). Regression on Manifolds Using Kernel Dimension Reduction. Proceedings 
of the 24th International Conference on Machine Learning, Corvallis, OR, 2007. 
Chen, C.H., Li, K.C., 1998. Can SIR be as popular as multiple linear regression? Statistica Sinica 8, 289–316. 
Wang T and Zhu L. (2013). Sparse sufficient dimension reduction using optimal scoring. Computational Statistics 
and Data Analysis. 57: 223-232. 


Peer review:
[Review]:
Summary of the paper:

This paper proposes a sparse extension of sufficient dimension reduction (the problem of findind a linear subspace so that the output and the input are conditionally independent given the projection of the input onto that subspace). 

The sparse extension is formulated through the eigenvalue formulation of sliced inverse regression.

The method is finally applied to identifying genetic variants of clinical significance.

Comments:

-Other sparse formulations of SIR have been proposed and the new method should be compared to it (see two below)

Lexin Li, Christopher J. Nachtsheim, Sparse Sliced Inverse Regression, Technometrics. Volume 48, Issue 4, 2006

Lexin Li. Sparse sufficient dimension reduction. Biometrika (2007) 94(3): 603-613 

-In experiments, it would have been nice to see another method run on these data.

-The paper appears out of scope for the conference