Biomarkers are often organized into networks, in which the strengths of network connections vary across subjects depending on subject-specific covariates. Variation of brain network connections, as subject-specific feature variables, has been found to predict disease clinical outcome in many studies. In this work, we develop a two-stage statistical method to estimate covariate-dependent brain networks to account for heterogeneity among network measures and evaluate their association with disease clinical manifestation. In the first stage, we propose a conditional Gaussian graphical model with mean and precision matrix depending on covariates to obtain subject-specific networks. In the second stage, we model the network connections estimated from the first step jointly with the biomarkers and covariates to identify important features of a clinical outcome via regularized regression. We assess the performance of our proposed method by extensive simulation studies and apply the method to a Huntington's disease (HD) study to investigate the effect of HD causal gene on the rate of change in motor symptom as mediated through brain subcortical and cortical gray matter atrophy connections.