BACKGROUND
the pairwise relationships between different clinical traits  and/or gene expression traits  have been successfully described with undirected gene co-expression networks  <cit> .
while gene expression traits  and clinical traits represent different quantities, both can be described in undirected trait networks.
by definition, these undirected networks cannot be used to describe causal relationships between the traits.
causal information can be encoded by directed networks where a → b if trait a causally influences trait b. we refer to the process of assigning a causal direction to at least some of the edges in a trait network as 'edge orienting'.
experimental edge orienting approaches include transgenic modifications, viral-mediated over-expression, and chemical perturbation of genes.
edge orienting methods can also be based on various approaches that involve multiple perturbations, such as genetic- and time series experiments  <cit> , or by integrating protein interaction and gene expression data  <cit> .
using genetic markers for orienting the edges of trait networks generated in genetic experiments provides significant statistical power and specificity for recovering directed edges  <cit> .
since randomization is the most convincing method for establishing causal relationships between two traits  <cit> , it is natural to make use of genetically randomized genotypes  to derive causality tests that are less susceptible to confounding by hidden variables  <cit> .
if a trait a is significantly associated with a genetic marker m, variation in m must be a cause of variation in a  since the randomization of marker alleles during meiosis precedes their effect on trait a. since the orientation of the edge between m and a is unambiguous, m is referred to as a causal anchor of a  <cit> .
we follow the convention of path analysis to represent a causal model by a directed graph.
for example, the directed graph m → a → b implies that the genetic marker m has a causal effect on trait a, which in turn has a causal effect on trait b. a causal graph encodes independencies between variables.
conditional independence can be determined by the graphical property of d-separation  <cit> .
if two traits a and b are d-separated in the graph by a set of variables s, then the two traits are independent given the variables in s. for example, m → a → b implies that m and b are independent after conditioning on a.
d-separation predicts the correlational consequences of conditioning in the causal graph  <cit> .
by testing the correlational predictions and assuming no false independencies , one can sometimes orient edges using observational data alone  <cit> .
RESULTS
correlation-based tests of causal models
for simplicity, we assume that the genetic markers are single nucleotide polymorphisms .
for a given sample , a bi-allelic snp can take on one of three possible genotypes.
by default, we assume an additive genetic effect and encode these genotypes as  <dig>   <dig>  or  <dig>  but alternative marker codings could also be considered.
to quantify the linear relationship between a snp marker m and a trait a, we use the correlation coefficient cor.
ordinal variables are routinely used in path analysis and structural equation modelling  <cit> .
to determine whether trait a mediates the effect of marker m on trait b  one can assess how conditioning on a affects the correlation between m and b. to quantify the linear relationship between m and b after conditioning on a, we use the partial correlation coefficient:
cor=cor−corcor2)2).
if the causal model m → a → b is correct, then the partial correlation coefficient cor is expected to be  <dig>
we use fisher's z transform to assess the statistical significance of a sample correlation coefficient r  <cit> :
zfisher= <dig> n−3 log,
where n denotes the sample size; zfisher asymptotically follows a normal distribution  with mean μ and variance  <dig>  under the null hypothesis of zero correlation, μ =  <dig> and zfisher follows a standard normal distribution.
for brevity, we denote the fisher transformations of the correlation coefficients cor and cor by z = zfisher) and z = zfisher), respectively.
if the causal graph m → a → b  is correct, z follows a standard normal distribution.
thus, if the p-value corresponding to z is high , the data fit the assumed causal graph.
using path analysis rules, the causal graph m → a → b  implies the following relationships between the correlation coefficients
cor = corcor
we refer to the marker m as a candidate common pleiotropic anchor  of a and b. if the expected values of cor and cor are non-zero and the causal model holds, eq.
implies that the genetic marker m will be significantly correlated with both a and b. thus, the marker m can be confirmed as a pleiotropic anchor of a and b by confirming the fit of the causal model m → a → b. we will now consider a situation where the correlation between a and b stems from a hidden confounder c, i.e.
m → a ← c → b. the graph implies that a and b are correlated due to the shared confounder c. the correlation cor is expected to be  <dig> since the arrows between m and b collide at a, i.e, m and b are d-separated without conditioning.
in this situation z = zfisher) follows a standard normal distribution.
if the p-value corresponding to z is high , the data fit a confounded model.
in contrast, the partial correlation cor is expected to be non-zero since conditioning on a 'activates' the causal flow through the collider node, i.e.
it induces conditional dependence  <cit> .
the opposite  causal graph m → a ← b also implies that the expected value of cor is zero since the causal paths collide at a. conditioning on a activates this collider node, and the partial correlation cor is expected to be non-zero.
similarly, one can show that the model a ← m → b implies that cor is expected to be zero.
under this causal model, z asymptotically follows a standard normal distribution.
in contrast, cor is expected to be non-zero.
these considerations illustrate that one can test the predicted correlational consequences of a causal model and thus evaluate its fit.
we will now consider the situation of multiple markers .
denote by
ma={ma,ma,...ma}
a set of candidate common pleiotropic anchors of a and b. analogous to eq.
, the causal model ma → a → b implies cor,b)=cor,a)cor.
the model implies that the partial correlations cor,b|a) are expected to be zero, i.e.
zfisher,b|a)) is predicted to follow a standard normal distributions.
frequently an additional set of markers mb is also available for trait b .
for example, when one marker is available for each trait, i.e.
ma→a→b←mb, the correlation cor) is expected to be  <dig> since the causal arrows 'collide' at b  <cit> .
geometrically speaking, the expected zero correlation between a and mb implies that the corresponding standardized vectors are orthogonal.
therefore, we refer to marker mb as an orthogonal causal anchor  with respect to the edge a → b. we will argue that the availability of orthogonal causal anchors significantly improves the recovery of the causal signal .
if the model ma→a→b←mb is correct, cor,b|a),cor) are expected to be zero and z,b|a) and z) asymptotically follow standard normal distributions.
local sem-based edge orienting scores
while sems can be used to study the fit of multi-trait causal models  <cit>  we only consider the local causal models depicted in figure  <dig> since the proposed neo method evaluates the orientation of each edge separately based on the best causal anchors available.
the fit of each single marker model in figure  <dig> can be tested using a chi-square test with  <dig> degree of freedom.
we refer to the resulting p-value as the model p-value.
in the methods section, we review and discuss the use of model p-values for quantifying the fit of a causal model.
the main point is that the higher the model p-value, the better the causal model fits the data.
to summarize the genetic evidence in favor of a given edge orientation a → b, we propose the use of edge orienting scores.
the higher the value of an edge orienting score for the orientation a → b, the stronger genetic evidence favors this causal model.
in the following, we propose local sem-based edge orienting  scores for orientation a → b. for a single genetic marker m and traits a and b, we consider the  <dig> different local causal models depicted in figure  <dig>  additional single marker models are possible.
however, under the constraint that the markers are causal anchors , then the five models pictured for nodes  in figure  <dig> exhaust all possible three node models that both  explain a – b and  associations and  can be tested.
the critical technical issue is having degrees of freedom  remaining after estimating the model parameters.
if the degrees of freedom are  <dig>  the model p-values cannot be calculated.
the  <dig> different local causal models depicted in figure  <dig> are used to compute the following model p-values: p, p, p, p, and p. while a detailed analysis should consider all model p-values, we find it useful to summarize the genetic evidence in favor of a given orientation a → b  using a single number: the local sem-based edge orienting next best  score.
the leo.nb score is defined by dividing the model p-value for a → b by the p-value of the best fitting alternative model, i.e.
the best of models 2– <dig> in figure  <dig>  the chi-square test p-value of the best fitting alternative model is the maximum p-value of the alternative causal models.
specifically, we define the single-marker leo.nb score as follows:
leo.nb.singlemarker=log⁡10max⁡,p,p,p)).
a positive leo.nb score indicates that the p-value in favor of model a → b is higher than that of any of the competing models in figure  <dig>  a negative leo.nb score indicates that the a → b model is inferior to at least one alternative model.
in our simulations, we use a threshold of  <dig> for leo.nb.singlemarker.
multi-marker leo.nb score
it is straightforward to generalize the single marker leo.nb score  to a set of genetic markers ma={ma,ma,...} .
we refer to the resulting edge orienting score as the leo.nb.cpa score since it is based on the set of candidate pleiotropic anchors ma :
leo.nb.cpa=log⁡10max⁡,p,p,p)).
note that the multi-marker models used in the definition of leo.nb.cpa correspond to the single marker models of figure  <dig> with m replaced by ma.
if an additional genetic marker set mb={mb,mb,...mb} associated with trait b is available , we propose to use another edge orienting score.
according to our consistency assumption, ma contains markers that are more strongly correlated with a than with b. similarly, mb holds markers more strongly correlated with b than with a. if the orientation a → b is correct, then each of the ma markers has a pleiotropic effect by impacting first a and subsequently b. furthermore, we refer to the markers in mb as candidate orthogonal causal anchors  since the model a → b implies that these markers impact b, but are independent of both a and ma.
we define the likelihood-based orthogonal causal anchor  score by assessing whether the model ma → a → b ← mb has a higher p-value than the alternative models depicted in figure  <dig>  specifically, we define
leo.nb.oca=log⁡10max⁡,p,p)).
note that model  <dig> in the denominator involves a hidden confounder c. the use of two independent genetic marker sets  alleviates the problem of model identifiability that may plague a cpa based edge orienting score.
model equivalence is also a key consideration in choosing which models to compare.
from the standpoint of model equivalence, we note that the multiple anchor models presented in figure  <dig> include a model with a hidden  variable connecting a and b, and that no such model is included in the single anchor model comparisons.
such a model was found to be indistinguishable from the models with a collider node, such as single anchor models  <dig> and model  <dig>  in the single marker case, both the collider node and the hidden variable models test for independence in the marginal relationship between the anchor and the more distal trait node.
future research may lead to an understanding of what type of data allow one to consider additional alternative models for the edge score computation.
it should be straightforward to adapt the proposed leo score to additional models as long as their model p-values can be calculated.
correlated markers, which are frequently encountered in practice such as in haplotype blocks, may compromise the performance of edge orienting scores.
leo scoring allows multiple parents of a node to be correctly accounted for within each model.
moreover, the parents  of a model are allowed to co-vary.
by contrast, the orthogonal causal anchor set is, by definition, penalized for any covariation with the pleiotropic anchors.
thresholds for the edge orienting scores
for the single marker score leo.nb.singlemarker, we use a threshold of  <dig>  which implies that the model p-value of the causal model is  <dig> =  <dig> fold higher than that of the next best model.
for the leo.nb.cpa and the leo.nb.oca, we use lower thresholds of  <dig>  and  <dig> , respectively.
using simulation studies presented in the additional file, we found that these thresholds lead to false positive rates that are often substantially below  <dig> .
similar to other statistical procedures, neo is susceptible to the pitfalls of multiple testing that may inflate the false positive rate.
permutation procedures and data dependent schemes  may inform the user on how to pick a threshold for a particular application.
further, we provide r software code for carrying out both single edge and multi-edge simulation studies.
simulation studies can be used to determine the power and false positive rates in different settings .
in practice, one often observes strong dependence relationships between genetic markers.
our simulations show that correlations between genetic markers can reduce the power of edge orienting scores.
further, we mention that the neo software implements an option for removing redundant markers that are highly correlated with each other.
the removal of redundant markers may alleviate the loss of power.
overview of network edge orienting with neo
we now provide a detailed step-by-step description of a typical neo analysis.
an overview is also provided in figure  <dig>
step 1: integrate traits  and snps
neo takes trait and genetic marker data as input.
traits can include microarray gene expression data, clinical phenotypes, or other quantitative variables.
each snp or trait is a node in the network, and the neo software evaluates and scores the edge between traits a and b if the absolute correlation |cor| lies above a user-specified threshold.
for each edge a – b, neo generates edge orienting scores for both possible orientations: a → b and b → a. if an erroneous edge exists between two traits, then it is meaningless to orient it.
the neo software can be used to orient any edge that the user chooses to consider.
to allow the user to judge whether the existence of an edge is supported by the data, the neo software outputs a wald test statistic of the path coefficient, the corresponding p-value, and the correlation between the two traits.
if the wald test p-value is insignificant, orienting the edge may be meaningless.
step 2: genetic marker selection and assignment to traits
edge orienting scores will only be generated for edges whose traits have been anchored to at least one genetic marker.
two basic approaches for anchoring traits to markers are implemented in the neo software: a manual selection by the user or an automatic selection by the software itself.
manual snp selection
neo provides great flexibility to the user on how to anchor traits to markers.
for example, the user can manually assign snps to the traits .
this flexibility entails that the user carefully studies what constitutes a significant relationship between traits and markers and between the traits in the data set.
the user may wish to anchor traits to snps that have been implicated by prior genetic analyses.
for example, results from previous quantitative trait locus studies may implicate genetic markers associated with a trait.
multiple comparison issues are just starting to be addressed in the sem literature  <cit> .
edge scores cannot be computed when an overly strict multiple testing control results in no causal anchors.
on the other hand, an overly lax multiple testing control may result in spurious causal anchors which may lead to erroneous edge scores.
we recommend that conservative measures of genome-wide qtl significance  <cit>  and false discovery rate be applied when selecting the initial causal anchor.
once a causal anchor has been established as obtaining genome-wide significance, neo can be used to evaluate the fit of different causal models.
automatic snp selection
neo can also be used to automatically relate  traits to snps.
the automated snp selection methods consider each trait a in isolation from the other traits when defining a preliminary genetic marker set .
toward this end, the user can choose 1) a greedy approach based on univariate linear regression models, 2) a forward-stepwise approach based on multivariate linear regression models, or 3) both.
the greedy snp selection approach defines m'a as the set of markers with the k highest absolute correlations with a. the greedy approach is equivalent to using univariate linear regression models to relate a to each marker separately and subsequently picking the k most significant markers.
for creating multivariate linear qtl models, neo also implements forward-stepwise marker selection.
the forward-stepwise marker selection method may avoid a pitfall that plagues the greedy snp selection: if several genetic markers are located very close to each other , the greedy snp selection may pick all of them before considering snps at other loci associated with the same trait.
for this reason, we recommend combining greedy and forward-stepwise snp selection methods.
once the preliminary sets of markers m'a and m'b are obtained, neo evaluates the consistency of each set.
we utilize a marker assignment consistency heuristic: a genetic marker can only serve as causal anchor for one trait.
to fulfill this heuristic, a snp is moved from m'a to m'b if its correlation with b is stronger than that with trait a. we denote the resulting consistent genetic marker sets by m"a and m"b. the resulting consistent genetic marker sets may be comprised of dozens of snps.
therefore, it can be useful to further filter the snps according to their joint predictive power for the trait.
toward this end, we use the akaike information criterion  in conjunction with multivariate regression models to select genetic markers from within the consistent genetic marker sets  <cit> .
specifically, to define the final genetic marker set ma for trait a, we use the aic criterion to find a parsimonious multivariate regression model of a using predictors from within m"a. the final sets of markers ma and mb are thus comprised of consistent genetic markers that according to the aic criterion best predict their respective traits; we use these final sets as causal anchors in computing the edge orienting scores.
the forward-stepwise approach based on multivariate linear regression models is akin to a legal courtroom where two cases are built, weighed, and judged.
broadly, the strongest genetic support  for the genetic influence on a and b are built independently, using aic-based halting criteria.
after consistency checks, these multivariate eqtl models are weighed by embedding them in causal models  and models are then compared using sem fitting indices.
when candidate cpa markers can be found for a and ocas for b, the neo method provides stringent consistency checks and balances against over-fitting.
we consider automated snp selection particularly useful when no prior evidence suggests causal anchors for the traits.
step 3: compute local edge orienting scores for aggregating the genetic evidence in favor of a causal orientation
both leo.nb.cpa and leo.nb.oca scores are computed for each edge orientation .
we recommend using the leo.nb.oca score  as the primary edge orienting score if markers affect both a and b. however, if the results of the leo.nb.cpa score strongly disagree with those of the leo.nb.oca score, the latter should not be trusted.
as described in the next step, all fitting indices should be considered before calling an edge causal.
step 4: for each edge, evaluate the fit of the underlying local sem models
edges with high edge orienting scores may not necessarily correspond to causal relationships.
although edge orienting scores flag interesting edges, they are no substitute for carefully evaluating the fit of the underlying local sems.
since a leo.nb score is defined as a ratio of two model p-values, it is advisable to check whether both p-values are small, as this would indicate poor fit of either model.
if the model p-value of the confounded model a ← c → b is high, the correlation between a and b may be largely due to a hidden confounder c. neo  also report a wald test statistic for the path coefficient from a → b. if the wald test for an edge is significant, the data support its existence.
apart from the model p-value, many other sem model fitting indices have been defined by contrasting the observed covariance matrix sm × m with the fitted covariance matrix Σ as detailed in the methods section.
the neo software reports the standard sem fitting indices  <cit>  that are implemented in the r package sem  <cit>  including the root mean square error of approximation , comparative fit index , standardized root mean square residual , bic.
since a single fitting index reflects only a particular aspect of model fit, a favorable value of that index does not by itself demonstrate good model fit; it is important to assess the model fit based on multiple indices.
we follow the following standard guidelines for interpreting these indices  <cit> .
before calling an edge a → b causal, we recommend verifying that the corresponding causal model has a high model p-value , a low rmsea score , a low srmsr , a high cfi , and a significant wald test p-value .
step 5: robustness analysis with respect to snp selection parameters
since the edge orienting scores for an edge a – b critically depend on the input genetic marker sets ma and mb, we also recommend carrying out a robustness analysis with respect to different marker sets.
in particular, the automated snp selection results should be carefully evaluated with regard to the threshold parameters that were used to define the marker sets.
for example, when using a greedy snp selection strategy, it is advisable to study how the leo.nb score is affected by altering the number of most highly correlated snps.
for a given edge and a given edge orienting score , neo implements a robustness analysis with respect to automatic marker selection .
a robustness plot shows how the leo.nb.oca score  depends on sets of automatically selected snp markers .
when using the default snp selection method , robustness step k corresponds to choosing the top k snps by greedy and forward selection for each trait.
since the greedy and forward snp selection may select the same snps, step k typically involves fewer than 2k snps per trait.
the edge orienting results should be relatively robust with respect to different choices of k.
step 6: repeat the analysis for the next a-b trait-trait edge and apply edge score thresholds to orient the network
neo orients each edge separately in an undirected input trait network.
the results are order-independent.
for each edge, neo repeats steps 1– <dig> until all edges have been assigned edge orienting scores.
once each edge has been scored, the user can generate a global, directed network by choosing an edge score  and a corresponding threshold .
neo output and r software
the primary output of neo is an excel spreadsheet which reports likelihood-based edge scores  and other edge scores that are described in the neo manual.
for each edge, the neo spreadsheet also contains hyperlinks that allows the user to access the log file for each edge.
the log file contains a host of information regarding computation of the edge orienting scores including sem model p-values, wald test statistics for each path coefficient, and the snp identifiers for the causal anchor sets ma and mb.
although the main output of neo are scores for every edge orientation, one can construct a global directed network by thresholding an edge orienting score.
neo uses the r software package sem  <cit>  to compute model p-values and other fitting indices.
the neo software is documented in a series of separate tutorials that illustrate real data applications and simulation studies.
these tutorials and the real data can be downloaded from our webpage.
applications
research goals that can be addressed with neo
neo can be used to address the following four research goals.
on the simplest level, neo can be used to assign edge orienting scores to a single edge using manually chosen genetic markers .
when dealing with a single edge and multiple genetic markers, the neo software can automatically select markers for edge orienting.
since the automatic marker selection entails certain parameter choices, we recommend carrying out a robustness analysis with respect to adding or removing genetic markers.
when dealing with a single trait a and manually selected genetic markers, the software can be used to screen for other traits that are causal or reactive to trait a. for example, in table  <dig> we screen for genes that are reactive to gene expression trait insig <dig>   when dealing with multiple edges, neo can be used to arrive at a global directed network.
this can be done by thresholding a chosen edge score.
if the resulting global network is acyclic  then d-separation  <cit>  and standard sem model fitting indices can be used to evaluate the fit of the global causal model to the data.
using the female mouse liver gene expression data, we report edge scores for the known causal relationships insig <dig> → dhcr <dig> and insig <dig> → fdft <dig> and the other edges depicted in figure  <dig>  the table represents a condensed summary of the neo software spreadsheet.
the high value of leo.nb.oca =  <dig>  suggests that this causal model is  <dig> ≈  <dig>  times more likely than the next best  local model.
similarly, leo.nb.oca =  <dig>  suggests that the causal model is  <dig> times more likely than the next best local model.
the fourth column reports the marginal pearson correlation coefficient, while the three path columns  give details for each individual edge in the sem models.
the last five columns summarize the fits of the two best fitting sem models shown in figures  <dig> and .
the model probability column  was computed using a central χ <dig> statistic with the  <dig> degrees of freedom.
the high, non-significant model p-values suggest good fit.
the root mean square error of approximation  is a standard sem fit evaluation index that, similar to the χ <dig> stastic, evaluates the overall fit of the sem model; a value smaller than  <dig>  is desirable.
here we used the female mouse liver bxh data to illustrate that neo can be used to identify genes that are reactive to a given trait .
the table reports the  <dig> genes with highest leo.nb.oca scores.
since  <dig> of the  <dig> genes are already known to be reactive to insig <dig>  these results represent a highly significant validation success of neo; using the  <dig> genes on the array and assuming that there are  <dig> known genes downstream of insig <dig> , the fisher exact p-value of validation success is p =  <dig>  × 10- <dig>  the neo analysis in female mice also implicates  <dig> novel genes.
pubmed searches on these genes did not turn up any information about a role of these genes in liver or sterol homeostasis.
the table also reports the analysis results using the male mice of the bxh cross.
the validation  column shows a star  if the initial finding in female liver was replicated in the independent test set of  <dig> f <dig> male mice; we defined validation success as leo.nb.oca score above  <dig>  using the default settings of automatic snp marker selection.
the male liver analysis confirms three of the nine novel genes suggested from the female analysis: tlcd <dig>  slc25a <dig>  and qdpr.
the fact that not all genes can be replicated in the male data may reflect known differences between female and male mouse liver tissue expression profiles  <cit> .
mouse data description
we illustrate our methods using data from a previously studied f <dig> mouse intercross   <cit>  involving two inbred mouse strains .
the strain c57bl/6j is susceptible to a variety of atherosclerosis, diabetes, obesity, and heart disease related traits to which c3h/hej is resistant.
the f <dig> offspring mice are expected to show a significant spectrum of atherosclerosis and metabolic syndrome responses to a high-fat diet.
the mice were genotyped at  <dig> genetic markers  across the mouse genome.
a variety of physiological traits were measured, including mouse body weight, fat mass, insulin, glucose, free fatty-acid levels in the blood, and cholesterol fractions .
here we focus on gene expression data in mouse liver tissue.
since significant differences in the gene expression profiles between male and female mice have been observed  <cit> , we analyzed each gender separately.
application i: studying the causal relationships between insig <dig>  fdft <dig> and dhcr7
here we use both manually and automatically selected snp markers to compute edge orienting scores to the known causal relationships between genes in the cholesterol biosynthesis pathways.
the gene expression levels of insig <dig> serve as a sensitive proxy for the activation level of the srebp transcription factors  <cit> , allowing us to study the known biology of those genes in the cholesterol biosynthesis pathway.
we used the mouse liver gene expression data of the bxh mouse cross to determine whether two known causal edge orientations  <cit> insig <dig> → fdft <dig> and insig <dig> → dhcr <dig> result in high leo.nb scores.
for the female mice of the bxh cross, qtl analysis of insig <dig> expression implicated two candidate pleiotropic anchors  on chromosomes  <dig> and  <dig> .
together these  <dig> snps explained  <dig>  percent  of the variation of insig <dig>  as candidate orthogonal causal anchor of fdft <dig>  we selected a highly significant snp on chromosome  <dig> as can be see from the single marker lod score curve in figure  <dig>  similarly, we found a candidate orthogonal causal anchor for dhcr <dig> chromosome  <dig>  figures  <dig> show the causal models used to compute the model p-value in the numerators of the leo.nb.oca score and the leo.nb.oca score, respectively.
in table  <dig>  we provide more details on the edge scores of the causal models in figure  <dig>  we find that leo.nb.oca =  <dig> , which lies above the recommended threshold of  <dig> .
further, we find that the wald test of the path coefficient is highly significant .
the model p-value of the causal model is p =  <dig>  and the rmsea is ≤  <dig> .
these results suggest that there is indeed a causal relationship insig <dig> → fdft <dig>  for the edge orientation insig <dig> → dhcr <dig>  leo.nb.oca =  <dig>  and the wald test is highly significant at z =  <dig> , and the rmsea is  <dig> .
these results confirm the known causal relationship: insig <dig> → dhcr <dig>
for the female bxh mice, we also used automatic snp selection to compute leo.nb scores.
for edge orientations insig <dig> → dhcr <dig> and insig <dig> → fdft <dig>  the results of a robustness analysis are presented in figures  <dig> and  <dig>  respectively.
the robustness analysis suggests that both edges are causal in female mice since the leo.nb.cpa scores remain above the recommended threshold of  <dig> .
however, the robustness analysis of leo.nb.oca for edge insig <dig> → fdft <dig>  shows that for a particular set of automatically selected markers, the score dips below the recommended threshold of  <dig>  for this score.
since automatic snp selection is particularly vulnerable to false positive causal anchors, it is advisable to replicate the neo analysis in an independent data set.
for example, we also used automatic snp selection to compute edge orienting scores in male mice of the bxh cross.
although causal relationships may differ between male and female mice, replication in male mice certainly provides evidence that the reported causal relationships are true.
figures  <dig> and  <dig> show the results of a robustness analysis for leo.nb.oca and leo.nb.oca, respectively.
overall, we find that automatic marker selection with the leo.nb.oca and leo.nb.cpa scores provide evidence of the reported causal relationships in both male and female mouse liver data.
application ii: screening for genes that are reactive to insig1
in this application, we illustrate that for a single trait  and manually selected genetic markers neo can be used to screen for other traits that are reactive to the trait in question.
we again used the above-mentioned genetic markers on chromosomes  <dig> and  <dig> as causal anchors for insig <dig>  for each gene expression trait b, we computed a leo.nb score for the edge insig <dig> → b. table  <dig> reports details for  <dig> highest ranking genes.
prior literature  <cit>  suggests that  <dig> out of the  <dig> genes are reactive to insig <dig> and are part of the well-studied sterol homeostasis pathway.
since so many known sterol regulated positive controls are recovered simultaneously, these findings are highly significant.
using the  <dig> array genes , and assuming that there are  <dig> known genes downstream of insig <dig> , we compute the fisher exact test for the set of  <dig> predicted versus  <dig> known downstream genes giving a p-value of  <dig>  × 10- <dig> for the predicted novel gene set.
moreover, our analysis also implicates nine novel genes as being affected by the same pathway in female liver.
a pubmed literature search on these genes did not suggest known relationships to liver or sterol impacted gene expression.
neo gene screening requires careful validation.
for example, we report the results of a neo analysis in male mouse liver data in table  <dig>  of the nine novel genes suggested from the female analysis, the male liver analysis confirms three of these: tlcd <dig>  slc25a <dig>  and qdpr.
the disparity between male and female mice may reflect the tissue-specific expression and regulation of sexually dimorphic genes  <cit> .
relationship to prior work
the above application describes the use of neo for finding reactive genes to insig <dig>  here we contrast this neo based gene screening method to a related gene screening method  <cit>  that utilizes a hybrid between the single anchor and the candidate pleiotropic anchor approach  and ).
similar to our computation of the leo.nb.cpa score, the authors use a forward-backward stepwise regression procedure to build the initial genetic model for the downstream trait.
for each locus retained in the genetic model for a given trait, their lcms  test evaluates genes for causality by comparing three of the five single anchor models  shown in our figure  <dig> for smallest aic.
taking just the genes for which the causal model – model  <dig> of figure  <dig> for a single gene a and trait b – fits best for at least two common pleiotropic markers, the candidate causal gene list is generated by ranking genes according to 'the amount of genetic variance of the trait that was causally explained by variation in their transcript abundance,' which amounts to comparing the p-values of the cpa models for all final candidate genes.
in contrast, neo makes use of orthogonal anchors and fits multiple orthogonal anchors simultaneously.
our simulations suggest power advantages of the resulting leo.nb.oca score .
application iii: fsp <dig> is upstream of a biologically interesting gene co-expression module in female bxh mice
here we illustrate how neo can be used to assign edge orienting scores to a single edge using manually and automatically chosen genetic markers.
specifically, we computed edge orienting scores for the edge fsp <dig> → meblue where fsp <dig>  corresponds to a pro-apoptotic gene that is related to metabolic syndrome: fsp27-null mice have been found to be resistant to obesity and diabetes; fsp <dig> expression is halved in obese humans after weight loss; and fsp <dig> regulates lipolysis in white human adipocytes  <cit> .
the other quantitative trait, meblue, represents the activation status of an entire pathway.
more specifically, meblue is the module eigengene  of the biologically important 'blue' gene co-expression module described in  <cit> .
this gene co-expression module was comprised of highly correlated genes and meblue is a summary gene expression trait that best represents the expressions of the blue module genes.
to study whether fsp <dig> causally affects meblue, we conducted both manual and automatic snp selection approaches.
we used a previously identified snp marker on chromosome  <dig>  that affected the expression of the blue module genes  <cit>  and of several physiologic traits as the manually chosen input snp.
this genetic marker was previously referred to as a module quantitative trait locus  since it was found to affect the gene expression profiles of most blue module genes.
using this snp, we found highly causal leo.nb scores between fsp <dig> and meblue.
the leo.nb.oca scores passed the threshold of  <dig>  and the leo.nb.cpa scores passed the threshold of  <dig> .
we also used the automatic snp selection strategies to assess the causal relationship between fsp <dig> and meblue.
the results of a robustness analysis can be found in figure  <dig>  we find that the causal relationship fsp <dig> → meblue is highly robust with respect to different automatic marker selection methods.
simulation studies
multi-edge simulation model that involves one hundred snps
neo analysis can orient the edges of a multi-trait network by automatically selecting markers for each trait separately.
the analyses proceed in a stepwise fashion: edges are oriented one at the time.
for each edge, neo computes edge orienting scores .
by thresholding these edge orienting scores, one can arrive at a globally oriented trait network.
the details of the simulation model and relevant r code is presented in an r software tutorial on our webpage.
briefly, we simulated a causal network between five gene expressions  and a trait .
each of the  <dig> traits was simulated to be under the causal influence of  <dig> snps.
we added  <dig> noise snps so that the data contained  <dig> snps.
we simulated the following causal relationships between the traits:
e <dig> → e <dig>
e <dig> → e <dig>
e <dig> ← hiddenconfounder → e <dig>
e <dig> → trait
trait → e <dig>
note that the correlation between traits e <dig> and e <dig> was entirely due a hidden confounder.
the heatmap plot in figure  <dig> depicts the true causal model.
note that a red square in the i-th row and j-th column indicates that trait i causally affects trait j. the rows and columns of the heatmap are ordered according to a hierarchical clustering tree, which was constructed using average linkage hierarchical clustering with the dissimilarity diss =  <dig> - |cor|.
figure  <dig> shows the corresponding heatmap of the observed network that was reconstructed using the leo.nb.oca score.
figure  <dig> shows an alternative output graph of neo.
blue edges indicate significant correlations  and a leo.nb.oca score is added to each edges whose leo.nb.oca score passes a user-supplied threshold.
we find that all true causal edges are correctly retrieved at the recommended leo.nb.oca threshold of  <dig> .
figure  <dig> shows the results of a robustness analysis for the leo.nb.oca and leo.nb.cpa scores for the edge orientation e <dig> → trait.
the leo.nb.oca and the leo.nb.cpa scores exceed their respective threshold of  <dig>  and  <dig>  for all steps of the robustness analysis, i.e., they retrieve the orientation correctly.
alternative simulation models can be explored using our online tutorial.
single edge simulation model parameterized with the heritability
in additional file  <dig>  we describe several simulation studies that use a single edge simulation model.
briefly, we simulated two traits a and b that are anchored to genetic marker sets ma and mb, respectively.
the correlation cor results from both a causal influence of b on a and from a hidden confounder c. this single edge model is used i) to study the choice of thresholds for the leo.nb scores, ii) to compare the leo.nb.cpa with leo.nb.oca scores, and iii) to evaluate automatic snp selection methods.
the results of these simulations are described in additional file  <dig> and in our online r software tutorials.
discussion
we propose methods for using multiple genetic markers to recover causal trait-trait relationships in systems genetic studies.
neo will be particularly useful for the analysis of experiments in which common genetic variations are leveraged to explore complex genetic traits.
we propose several edge orienting scores that measure the genetic evidence in favor of a given edge orientation a → b. while several methods exist for constructing undirected gene co-expression networks based on thousands of genes, we have evaluated the neo method for inferring directed networks involving relatively few genes .
future research could explore the use of the method for inferring directed networks involving thousands of genes.
our simulation studies show that orthogonal causal anchors lead to powerful edge scores that may outperform scores based only on candidate pleiotropic anchors .
to afford flexibility to the user, the neo software provides several options for anchoring the traits to genetic markers , computing local edge scores , and diagnosing poor model fit .
neo provides multiple options for automatically anchoring a trait to genetic markers: greedy, forward, and combined  snp marker selection.
while our simulation studies suggest that these three snp marker selection methods have similar performance, we find that the combined snp marker selection performs best when signal snps are in high linkage disequilibrium with noise snps .
neo's local, stepwise approach for orienting edges of a trait network allows one to orient networks involving hundreds or even thousands of traits.
since the calculation of edge orienting scores is based on local causal models, neo is relatively robust with regard to mistaken orientation of some edges in the global network.
although neo performs well in simulation studies and the reported real data applications, we note that it has several limitations.
the first limitation is that it requires the availability of genetic markers that are significantly associated with at least one trait per edge.
spurious associations between the markers and traits will result in meaningless edge orienting scores.
although the multi-marker score  is quite robust to noise snps in our simulations, false-positive input snps will result in unreliable edge scores.
the automatic snp selection is particularly vulnerable to false positives and its results should be carefully validated using biological experiments or causality analysis of independent data.
the second limitation is that the resulting global trait network may contain loops, i.e.
it may be cyclic.
in contrast, a directed acyclic graph  has no cycles.
dags appear in models where it does not make sense for a trait to have a path to itself.
while local dags are used for orienting individual edges, the reconstructed global trait network may no longer be acyclic.
acyclicity is theoretically desirable since it allows one to test causal predictions using pearl's formalism of d-separation  <cit> .
the constraint of acyclic graphs in many network learning algorithms is often more a mathematical convenience than reflective of biology; cycles may reflect feedback loops for maintaining homeostasis.
when more and more edges are oriented, as in the ic/ic*  <cit>  and pc/pc*  <cit>  algorithms, an error in one part of the network can propagate and cause erroneous orientations in unrelated portions of the network.
most often these errors arise due to confusion between confounded and truly causal flows.
to avoid being misled, neo deliberately discards the evidence from correlated trait neighbors in the undirected graph during leo scoring.
by computing local edge orienting scores without regard to a global acyclicity constraint, the analysis is relatively robust to mis-oriented neighboring edges.
neo uses causal anchors for each edge separately and thus allows the genetic data to speak for themselves.
the proposed leo.nb scores are local in that they orient one edge at a time without regard to the orientations of the other edges.
the reconstruction of the global network should be taken with a grain of salt.
while we report one simulation model where the global network was reconstructed correctly, future research should carefully evaluate the performance of the neo approach for inferring global networks.
a potential use of neo is to use it for initializing an iterative edge orienting algorithms for large networks that maximizes a global sem fitting index.
the third limitation is that the sem-based edge orienting scores assume linear relationships between traits and snp markers.
this is mathematically convenient but non-linear effects are common and have been reported in the literature  <cit> .
the neo approach works in the domain of linear graphical models since it is based on correlations and sems.
akin to the use of pearson versus spearman correlation, the software also offers the option of modelling monotone quantitative relationships in neo by converting all data to ranks before further processing.
neo will not work for traits that satisfy non-monotone relationships.
a fourth limitation is that the influence of genetic markers may be indirect.
neo may miss some relationships.
while snp changes must be upstream  of gene expression and phenotype manifestations, this does not preclude some snps from modifying the action of other snps, and the effect of such modifiers may become apparent only in particular contexts.
causal inference and structural equation modeling assume that relevant traits and causal anchors have been included in the causal model.
under-specified causal models, i.e.
models that omit important variables, may mislead the user to detect spurious causal relationships.
neo leads to relatively simple causal networks that do not incorporate dynamic or hierarchical properties .
given all these potential limitations it is reassuring that neo performs well at retrieving known causal relationships in the reported real data applications.
since neo focuses on individual edges, we expect that neo will be particularly useful for identifying traits that are causal for  a given trait.
for example, we illustrate that neo can be used to identify gene expression traits that are reactive to insig <dig>  the sem r package can be used to evaluate the global fit of an acyclic multi-trait network.
the neo algorithm computes an edge score for each edge without regard to the information gained from neighboring edges.
neo aims to harness the power of the established upstream causal anchors  as fully as possible; thus, it is appropriate when genetic variations are a major source of the variation in the traits.
to the extent the environment  is also varied, the neo approach may be less effective.
it is plausible that additional assumptions may allow one to use unshielded colliders to improve the causal inference  <cit> .
this is a promising avenue of future research.
we focused on the use of snp markers which capture only a limited amount of the sequence information of each individual.
in the not too distant future, it will be economically feasible to obtain the sequence information of each study subject.
since sequence information is likely to enhance the causal anchor assignment, sequence data may greatly improve the power of the neo method.
apart from the common genetic variation that perturbs gene expression in mouse crosses, neo can also be applied to orient edges on the basis of causal anchors from population-based allelic association studies, cell hybrids, or transfected cells.
CONCLUSIONS
natural randomization of alleles that occurs during meiosis can be used to study the causal information flow through trait networks.
for example, we use mouse cross data to retrieve known causal relationships in the sterol biosynthesis pathway.
we find that the proposed edge scores  are quite robust with respect to adding extraneous noise snps.
combined with the use of orthogonal causal anchors, the proposed edge orienting scores can provide a strong basis for further experimental evaluation of the predicted causal relationships.

 the pairwise relationships between different clinical traits and/or gene expression traits have been successfully described with undirected gene co-expression networks. while gene expression traits and clinical traits represent different quantities, both can be described in undirected trait networks. by definition, these undirected networks cannot be used to describe causal relationships between the traits. causal information can be encoded by directed networks where a → b if trait a causally influences trait b. we refer to the process of assigning a causal direction to at least some of the edges in a trait network as'edge orienting'. experimental edge orienting approaches include transgenic modifications, viral-mediated over-expression, and chemical perturbation of genes. edge orienting methods can also be based on various approaches that involve multiple perturbations, such as genetic- and time series experiments, or by integrating protein interaction and gene expression data. using genetic markers for orienting the edges of trait networks generated in genetic experiments provides significant statistical power and specificity for recovering directed edges. since randomization is the most convincing method for establishing causal relationships between two traits, it is natural to make use of genetically randomized genotypes to derive causality tests that are less susceptible to confounding by hidden variables. if a trait a is significantly associated with a genetic marker m, variation in m must be a cause of variation in a since the randomization of marker alleles during meiosis precedes their effect on trait a. since the orientation of the edge between m and a is unambiguous, m is referred to as a causal anchor of a. we follow the convention of path analysis to represent a causal model by a directed graph. for example, the directed graph m → a → b implies that the genetic marker m has a causal effect on trait a, which in turn has a causal effect on trait b. a causal graph encodes independencies between variables. conditional independence can be determined by the graphical property of d-separation. if two traits a and b are d-separated in the graph by a set of variables s, then the two traits are independent given the variables in s. for example, m → a → b implies that m and b are independent after conditioning on a. d-separation predicts the correlational consequences of conditioning in the causal graph. by testing the correlational predictions and assuming no false independencies, one can sometimes orient edges using observational data alone.  correlation-based tests of causal models for simplicity, we assume that the genetic markers are single nucleotide polymorphisms. for a given sample, a bi-allelic snp can take on one of three possible genotypes. by default, we assume an additive genetic effect and encode these genotypes as or but alternative marker codings could also be considered. to quantify the linear relationship between a snp marker m and a trait a, we use the correlation coefficient cor. ordinal variables are routinely used in path analysis and structural equation modelling. to determine whether trait a mediates the effect of marker m on trait b one can assess how conditioning on a affects the correlation between m and b. to quantify the linear relationship between m and b after conditioning on a, we use the partial correlation coefficient: cor=cor−corcor2)2). if the causal model m → a → b is correct, then the partial correlation coefficient cor is expected to be we use fisher's z transform to assess the statistical significance of a sample correlation coefficient r: zfisher= n−3 log, where n denotes the sample size; zfisher asymptotically follows a normal distribution with mean μ and variance under the null hypothesis of zero correlation, μ = and zfisher follows a standard normal distribution. for brevity, we denote the fisher transformations of the correlation coefficients cor and cor by z = zfisher) and z = zfisher), respectively. if the causal graph m → a → b is correct, z follows a standard normal distribution. thus, if the p-value corresponding to z is high, the data fit the assumed causal graph. using path analysis rules, the causal graph m → a → b implies the following relationships between the correlation coefficients cor = corcor we refer to the marker m as a candidate common pleiotropic anchor of a and b. if the expected values of cor and cor are non-zero and the causal model holds, eq. implies that the genetic marker m will be significantly correlated with both a and b. thus, the marker m can be confirmed as a pleiotropic anchor of a and b by confirming the fit of the causal model m → a → b. we will now consider a situation where the correlation between a and b stems from a hidden confounder c, i.e. m → a ← c → b. the graph implies that a and b are correlated due to the shared confounder c. the correlation cor is expected to be since the arrows between m and b collide at a, i.e, m and b are d-separated without conditioning. in this situation z = zfisher) follows a standard normal distribution. if the p-value corresponding to z is high, the data fit a confounded model. in contrast, the partial correlation cor is expected to be non-zero since conditioning on a'activates' the causal flow through the collider node, i.e. it induces conditional dependence. the opposite causal graph m → a ← b also implies that the expected value of cor is zero since the causal paths collide at a. conditioning on a activates this collider node, and the partial correlation cor is expected to be non-zero. similarly, one can show that the model a ← m → b implies that cor is expected to be zero. under this causal model, z asymptotically follows a standard normal distribution. in contrast, cor is expected to be non-zero. these considerations illustrate that one can test the predicted correlational consequences of a causal model and thus evaluate its fit. we will now consider the situation of multiple markers. denote by ma={ma,ma,...ma} a set of candidate common pleiotropic anchors of a and b. analogous to eq., the causal model ma → a → b implies cor,b)=cor,a)cor. the model implies that the partial correlations cor,b|a) are expected to be zero, i.e. zfisher,b|a)) is predicted to follow a standard normal distributions. frequently an additional set of markers mb is also available for trait b. for example, when one marker is available for each trait, i.e. ma→a→b←mb, the correlation cor) is expected to be since the causal arrows'collide' at b. geometrically speaking, the expected zero correlation between a and mb implies that the corresponding standardized vectors are orthogonal. therefore, we refer to marker mb as an orthogonal causal anchor with respect to the edge a → b. we will argue that the availability of orthogonal causal anchors significantly improves the recovery of the causal signal. if the model ma→a→b←mb is correct, cor,b|a),cor) are expected to be zero and z,b|a) and z) asymptotically follow standard normal distributions. local sem-based edge orienting scores while sems can be used to study the fit of multi-trait causal models we only consider the local causal models depicted in figure since the proposed neo method evaluates the orientation of each edge separately based on the best causal anchors available. the fit of each single marker model in figure can be tested using a chi-square test with degree of freedom. we refer to the resulting p-value as the model p-value. in the methods section, we review and discuss the use of model p-values for quantifying the fit of a causal model. the main point is that the higher the model p-value, the better the causal model fits the data. to summarize the genetic evidence in favor of a given edge orientation a → b, we propose the use of edge orienting scores. the higher the value of an edge orienting score for the orientation a → b, the stronger genetic evidence favors this causal model. in the following, we propose local sem-based edge orienting scores for orientation a → b. for a single genetic marker m and traits a and b, we consider the different local causal models depicted in figure additional single marker models are possible. however, under the constraint that the markers are causal anchors, then the five models pictured for nodes in figure exhaust all possible three node models that both explain a – b and associations and can be tested. the critical technical issue is having degrees of freedom remaining after estimating the model parameters. if the degrees of freedom are the model p-values cannot be calculated. the different local causal models depicted in figure are used to compute the following model p-values: p, p, p, p, and p. while a detailed analysis should consider all model p-values, we find it useful to summarize the genetic evidence in favor of a given orientation a → b using a single number: the local sem-based edge orienting next best score. the leo.nb score is defined by dividing the model p-value for a → b by the p-value of the best fitting alternative model, i.e. the best of models 2– in figure the chi-square test p-value of the best fitting alternative model is the maximum p-value of the alternative causal models. specifically, we define the single-marker leo.nb score as follows: leo.nb.singlemarker=log⁡10max⁡,p,p,p)). a positive leo.nb score indicates that the p-value in favor of model a → b is higher than that of any of the competing models in figure a negative leo.nb score indicates that the a → b model is inferior to at least one alternative model. in our simulations, we use a threshold of for leo.nb.singlemarker. multi-marker leo.nb score it is straightforward to generalize the single marker leo.nb score to a set of genetic markers ma={ma,ma,...}. we refer to the resulting edge orienting score as the leo.nb.cpa score since it is based on the set of candidate pleiotropic anchors ma: leo.nb.cpa=log⁡10max⁡,p,p,p)). note that the multi-marker models used in the definition of leo.nb.cpa correspond to the single marker models of figure with m replaced by ma. if an additional genetic marker set mb={mb,mb,...mb} associated with trait b is available, we propose to use another edge orienting score. according to our consistency assumption, ma contains markers that are more strongly correlated with a than with b. similarly, mb holds markers more strongly correlated with b than with a. if the orientation a → b is correct, then each of the ma markers has a pleiotropic effect by impacting first a and subsequently b. furthermore, we refer to the markers in mb as candidate orthogonal causal anchors since the model a → b implies that these markers impact b, but are independent of both a and ma. we define the likelihood-based orthogonal causal anchor score by assessing whether the model ma → a → b ← mb has a higher p-value than the alternative models depicted in figure specifically, we define leo.nb.oca=log⁡10max⁡,p,p)). note that model in the denominator involves a hidden confounder c. the use of two independent genetic marker sets alleviates the problem of model identifiability that may plague a cpa based edge orienting score. model equivalence is also a key consideration in choosing which models to compare. from the standpoint of model equivalence, we note that the multiple anchor models presented in figure include a model with a hidden variable connecting a and b, and that no such model is included in the single anchor model comparisons. such a model was found to be indistinguishable from the models with a collider node, such as single anchor models and model in the single marker case, both the collider node and the hidden variable models test for independence in the marginal relationship between the anchor and the more distal trait node. future research may lead to an understanding of what type of data allow one to consider additional alternative models for the edge score computation. it should be straightforward to adapt the proposed leo score to additional models as long as their model p-values can be calculated. correlated markers, which are frequently encountered in practice such as in haplotype blocks, may compromise the performance of edge orienting scores. leo scoring allows multiple parents of a node to be correctly accounted for within each model. moreover, the parents of a model are allowed to co-vary. by contrast, the orthogonal causal anchor set is, by definition, penalized for any covariation with the pleiotropic anchors. thresholds for the edge orienting scores for the single marker score leo.nb.singlemarker, we use a threshold of which implies that the model p-value of the causal model is = fold higher than that of the next best model. for the leo.nb.cpa and the leo.nb.oca, we use lower thresholds of and, respectively. using simulation studies presented in the additional file, we found that these thresholds lead to false positive rates that are often substantially below. similar to other statistical procedures, neo is susceptible to the pitfalls of multiple testing that may inflate the false positive rate. permutation procedures and data dependent schemes may inform the user on how to pick a threshold for a particular application. further, we provide r software code for carrying out both single edge and multi-edge simulation studies. simulation studies can be used to determine the power and false positive rates in different settings. in practice, one often observes strong dependence relationships between genetic markers. our simulations show that correlations between genetic markers can reduce the power of edge orienting scores. further, we mention that the neo software implements an option for removing redundant markers that are highly correlated with each other. the removal of redundant markers may alleviate the loss of power. overview of network edge orienting with neo we now provide a detailed step-by-step description of a typical neo analysis. an overview is also provided in figure step 1: integrate traits and snps neo takes trait and genetic marker data as input. traits can include microarray gene expression data, clinical phenotypes, or other quantitative variables. each snp or trait is a node in the network, and the neo software evaluates and scores the edge between traits a and b if the absolute correlation |cor| lies above a user-specified threshold. for each edge a – b, neo generates edge orienting scores for both possible orientations: a → b and b → a. if an erroneous edge exists between two traits, then it is meaningless to orient it. the neo software can be used to orient any edge that the user chooses to consider. to allow the user to judge whether the existence of an edge is supported by the data, the neo software outputs a wald test statistic of the path coefficient, the corresponding p-value, and the correlation between the two traits. if the wald test p-value is insignificant, orienting the edge may be meaningless. step 2: genetic marker selection and assignment to traits edge orienting scores will only be generated for edges whose traits have been anchored to at least one genetic marker. two basic approaches for anchoring traits to markers are implemented in the neo software: a manual selection by the user or an automatic selection by the software itself. manual snp selection neo provides great flexibility to the user on how to anchor traits to markers. for example, the user can manually assign snps to the traits. this flexibility entails that the user carefully studies what constitutes a significant relationship between traits and markers and between the traits in the data set. the user may wish to anchor traits to snps that have been implicated by prior genetic analyses. for example, results from previous quantitative trait locus studies may implicate genetic markers associated with a trait. multiple comparison issues are just starting to be addressed in the sem literature. edge scores cannot be computed when an overly strict multiple testing control results in no causal anchors. on the other hand, an overly lax multiple testing control may result in spurious causal anchors which may lead to erroneous edge scores. we recommend that conservative measures of genome-wide qtl significance and false discovery rate be applied when selecting the initial causal anchor. once a causal anchor has been established as obtaining genome-wide significance, neo can be used to evaluate the fit of different causal models. automatic snp selection neo can also be used to automatically relate traits to snps. the automated snp selection methods consider each trait a in isolation from the other traits when defining a preliminary genetic marker set. toward this end, the user can choose 1) a greedy approach based on univariate linear regression models, 2) a forward-stepwise approach based on multivariate linear regression models, or 3) both. the greedy snp selection approach defines m'a as the set of markers with the k highest absolute correlations with a. the greedy approach is equivalent to using univariate linear regression models to relate a to each marker separately and subsequently picking the k most significant markers. for creating multivariate linear qtl models, neo also implements forward-stepwise marker selection. the forward-stepwise marker selection method may avoid a pitfall that plagues the greedy snp selection: if several genetic markers are located very close to each other, the greedy snp selection may pick all of them before considering snps at other loci associated with the same trait. for this reason, we recommend combining greedy and forward-stepwise snp selection methods. once the preliminary sets of markers m'a and m'b are obtained, neo evaluates the consistency of each set. we utilize a marker assignment consistency heuristic: a genetic marker can only serve as causal anchor for one trait. to fulfill this heuristic, a snp is moved from m'a to m'b if its correlation with b is stronger than that with trait a. we denote the resulting consistent genetic marker sets by m"a and m"b. the resulting consistent genetic marker sets may be comprised of dozens of snps. therefore, it can be useful to further filter the snps according to their joint predictive power for the trait. toward this end, we use the akaike information criterion in conjunction with multivariate regression models to select genetic markers from within the consistent genetic marker sets. specifically, to define the final genetic marker set ma for trait a, we use the aic criterion to find a parsimonious multivariate regression model of a using predictors from within m"a. the final sets of markers ma and mb are thus comprised of consistent genetic markers that according to the aic criterion best predict their respective traits; we use these final sets as causal anchors in computing the edge orienting scores. the forward-stepwise approach based on multivariate linear regression models is akin to a legal courtroom where two cases are built, weighed, and judged. broadly, the strongest genetic support for the genetic influence on a and b are built independently, using aic-based halting criteria. after consistency checks, these multivariate eqtl models are weighed by embedding them in causal models and models are then compared using sem fitting indices. when candidate cpa markers can be found for a and ocas for b, the neo method provides stringent consistency checks and balances against over-fitting. we consider automated snp selection particularly useful when no prior evidence suggests causal anchors for the traits. step 3: compute local edge orienting scores for aggregating the genetic evidence in favor of a causal orientation both leo.nb.cpa and leo.nb.oca scores are computed for each edge orientation. we recommend using the leo.nb.oca score as the primary edge orienting score if markers affect both a and b. however, if the results of the leo.nb.cpa score strongly disagree with those of the leo.nb.oca score, the latter should not be trusted. as described in the next step, all fitting indices should be considered before calling an edge causal. step 4: for each edge, evaluate the fit of the underlying local sem models edges with high edge orienting scores may not necessarily correspond to causal relationships. although edge orienting scores flag interesting edges, they are no substitute for carefully evaluating the fit of the underlying local sems. since a leo.nb score is defined as a ratio of two model p-values, it is advisable to check whether both p-values are small, as this would indicate poor fit of either model. if the model p-value of the confounded model a ← c → b is high, the correlation between a and b may be largely due to a hidden confounder c. neo also report a wald test statistic for the path coefficient from a → b. if the wald test for an edge is significant, the data support its existence. apart from the model p-value, many other sem model fitting indices have been defined by contrasting the observed covariance matrix sm × m with the fitted covariance matrix Σ as detailed in the methods section. the neo software reports the standard sem fitting indices that are implemented in the r package sem including the root mean square error of approximation, comparative fit index, standardized root mean square residual, bic. since a single fitting index reflects only a particular aspect of model fit, a favorable value of that index does not by itself demonstrate good model fit; it is important to assess the model fit based on multiple indices. we follow the following standard guidelines for interpreting these indices. before calling an edge a → b causal, we recommend verifying that the corresponding causal model has a high model p-value, a low rmsea score, a low srmsr, a high cfi, and a significant wald test p-value. step 5: robustness analysis with respect to snp selection parameters since the edge orienting scores for an edge a – b critically depend on the input genetic marker sets ma and mb, we also recommend carrying out a robustness analysis with respect to different marker sets. in particular, the automated snp selection results should be carefully evaluated with regard to the threshold parameters that were used to define the marker sets. for example, when using a greedy snp selection strategy, it is advisable to study how the leo.nb score is affected by altering the number of most highly correlated snps. for a given edge and a given edge orienting score, neo implements a robustness analysis with respect to automatic marker selection. a robustness plot shows how the leo.nb.oca score depends on sets of automatically selected snp markers. when using the default snp selection method, robustness step k corresponds to choosing the top k snps by greedy and forward selection for each trait. since the greedy and forward snp selection may select the same snps, step k typically involves fewer than 2k snps per trait. the edge orienting results should be relatively robust with respect to different choices of k. step 6: repeat the analysis for the next a-b trait-trait edge and apply edge score thresholds to orient the network neo orients each edge separately in an undirected input trait network. the results are order-independent. for each edge, neo repeats steps 1– until all edges have been assigned edge orienting scores. once each edge has been scored, the user can generate a global, directed network by choosing an edge score and a corresponding threshold. neo output and r software the primary output of neo is an excel spreadsheet which reports likelihood-based edge scores and other edge scores that are described in the neo manual. for each edge, the neo spreadsheet also contains hyperlinks that allows the user to access the log file for each edge. the log file contains a host of information regarding computation of the edge orienting scores including sem model p-values, wald test statistics for each path coefficient, and the snp identifiers for the causal anchor sets ma and mb. although the main output of neo are scores for every edge orientation, one can construct a global directed network by thresholding an edge orienting score. neo uses the r software package sem to compute model p-values and other fitting indices. the neo software is documented in a series of separate tutorials that illustrate real data applications and simulation studies. these tutorials and the real data can be downloaded from our webpage. applications research goals that can be addressed with neo neo can be used to address the following four research goals. on the simplest level, neo can be used to assign edge orienting scores to a single edge using manually chosen genetic markers. when dealing with a single edge and multiple genetic markers, the neo software can automatically select markers for edge orienting. since the automatic marker selection entails certain parameter choices, we recommend carrying out a robustness analysis with respect to adding or removing genetic markers. when dealing with a single trait a and manually selected genetic markers, the software can be used to screen for other traits that are causal or reactive to trait a. for example, in table we screen for genes that are reactive to gene expression trait insig when dealing with multiple edges, neo can be used to arrive at a global directed network. this can be done by thresholding a chosen edge score. if the resulting global network is acyclic then d-separation and standard sem model fitting indices can be used to evaluate the fit of the global causal model to the data. using the female mouse liver gene expression data, we report edge scores for the known causal relationships insig → dhcr and insig → fdft and the other edges depicted in figure the table represents a condensed summary of the neo software spreadsheet. the high value of leo.nb.oca = suggests that this causal model is ≈ times more likely than the next best local model. similarly, leo.nb.oca = suggests that the causal model is times more likely than the next best local model. the fourth column reports the marginal pearson correlation coefficient, while the three path columns give details for each individual edge in the sem models. the last five columns summarize the fits of the two best fitting sem models shown in figures and. the model probability column was computed using a central χ statistic with the degrees of freedom. the high, non-significant model p-values suggest good fit. the root mean square error of approximation is a standard sem fit evaluation index that, similar to the χ stastic, evaluates the overall fit of the sem model; a value smaller than is desirable. here we used the female mouse liver bxh data to illustrate that neo can be used to identify genes that are reactive to a given trait. the table reports the genes with highest leo.nb.oca scores. since of the genes are already known to be reactive to insig these results represent a highly significant validation success of neo; using the genes on the array and assuming that there are known genes downstream of insig, the fisher exact p-value of validation success is p = × 10- the neo analysis in female mice also implicates novel genes. pubmed searches on these genes did not turn up any information about a role of these genes in liver or sterol homeostasis. the table also reports the analysis results using the male mice of the bxh cross. the validation column shows a star if the initial finding in female liver was replicated in the independent test set of f male mice; we defined validation success as leo.nb.oca score above using the default settings of automatic snp marker selection. the male liver analysis confirms three of the nine novel genes suggested from the female analysis: tlcd slc25a and qdpr. the fact that not all genes can be replicated in the male data may reflect known differences between female and male mouse liver tissue expression profiles. mouse data description we illustrate our methods using data from a previously studied f mouse intercross involving two inbred mouse strains. the strain c57bl/6j is susceptible to a variety of atherosclerosis, diabetes, obesity, and heart disease related traits to which c3h/hej is resistant. the f offspring mice are expected to show a significant spectrum of atherosclerosis and metabolic syndrome responses to a high-fat diet. the mice were genotyped at genetic markers across the mouse genome. a variety of physiological traits were measured, including mouse body weight, fat mass, insulin, glucose, free fatty-acid levels in the blood, and cholesterol fractions. here we focus on gene expression data in mouse liver tissue. since significant differences in the gene expression profiles between male and female mice have been observed, we analyzed each gender separately. application i: studying the causal relationships between insig fdft and dhcr7 here we use both manually and automatically selected snp markers to compute edge orienting scores to the known causal relationships between genes in the cholesterol biosynthesis pathways. the gene expression levels of insig serve as a sensitive proxy for the activation level of the srebp transcription factors, allowing us to study the known biology of those genes in the cholesterol biosynthesis pathway. we used the mouse liver gene expression data of the bxh mouse cross to determine whether two known causal edge orientations insig → fdft and insig → dhcr result in high leo.nb scores. for the female mice of the bxh cross, qtl analysis of insig expression implicated two candidate pleiotropic anchors on chromosomes and. together these snps explained percent of the variation of insig as candidate orthogonal causal anchor of fdft we selected a highly significant snp on chromosome as can be see from the single marker lod score curve in figure similarly, we found a candidate orthogonal causal anchor for dhcr chromosome figures show the causal models used to compute the model p-value in the numerators of the leo.nb.oca score and the leo.nb.oca score, respectively. in table we provide more details on the edge scores of the causal models in figure we find that leo.nb.oca =, which lies above the recommended threshold of. further, we find that the wald test of the path coefficient is highly significant. the model p-value of the causal model is p = and the rmsea is ≤. these results suggest that there is indeed a causal relationship insig → fdft for the edge orientation insig → dhcr leo.nb.oca = and the wald test is highly significant at z =, and the rmsea is. these results confirm the known causal relationship: insig → dhcr for the female bxh mice, we also used automatic snp selection to compute leo.nb scores. for edge orientations insig → dhcr and insig → fdft the results of a robustness analysis are presented in figures and respectively. the robustness analysis suggests that both edges are causal in female mice since the leo.nb.cpa scores remain above the recommended threshold of. however, the robustness analysis of leo.nb.oca for edge insig → fdft shows that for a particular set of automatically selected markers, the score dips below the recommended threshold of for this score. since automatic snp selection is particularly vulnerable to false positive causal anchors, it is advisable to replicate the neo analysis in an independent data set. for example, we also used automatic snp selection to compute edge orienting scores in male mice of the bxh cross. although causal relationships may differ between male and female mice, replication in male mice certainly provides evidence that the reported causal relationships are true. figures and show the results of a robustness analysis for leo.nb.oca and leo.nb.oca, respectively. overall, we find that automatic marker selection with the leo.nb.oca and leo.nb.cpa scores provide evidence of the reported causal relationships in both male and female mouse liver data. application ii: screening for genes that are reactive to insig1 in this application, we illustrate that for a single trait and manually selected genetic markers neo can be used to screen for other traits that are reactive to the trait in question. we again used the above-mentioned genetic markers on chromosomes and as causal anchors for insig for each gene expression trait b, we computed a leo.nb score for the edge insig → b. table reports details for highest ranking genes. prior literature suggests that out of the genes are reactive to insig and are part of the well-studied sterol homeostasis pathway. since so many known sterol regulated positive controls are recovered simultaneously, these findings are highly significant. using the array genes, and assuming that there are known genes downstream of insig, we compute the fisher exact test for the set of predicted versus known downstream genes giving a p-value of × 10- for the predicted novel gene set. moreover, our analysis also implicates nine novel genes as being affected by the same pathway in female liver. a pubmed literature search on these genes did not suggest known relationships to liver or sterol impacted gene expression. neo gene screening requires careful validation. for example, we report the results of a neo analysis in male mouse liver data in table of the nine novel genes suggested from the female analysis, the male liver analysis confirms three of these: tlcd slc25a and qdpr. the disparity between male and female mice may reflect the tissue-specific expression and regulation of sexually dimorphic genes. relationship to prior work the above application describes the use of neo for finding reactive genes to insig here we contrast this neo based gene screening method to a related gene screening method that utilizes a hybrid between the single anchor and the candidate pleiotropic anchor approach and ). similar to our computation of the leo.nb.cpa score, the authors use a forward-backward stepwise regression procedure to build the initial genetic model for the downstream trait. for each locus retained in the genetic model for a given trait, their lcms test evaluates genes for causality by comparing three of the five single anchor models shown in our figure for smallest aic. taking just the genes for which the causal model – model of figure for a single gene a and trait b – fits best for at least two common pleiotropic markers, the candidate causal gene list is generated by ranking genes according to'the amount of genetic variance of the trait that was causally explained by variation in their transcript abundance,' which amounts to comparing the p-values of the cpa models for all final candidate genes. in contrast, neo makes use of orthogonal anchors and fits multiple orthogonal anchors simultaneously. our simulations suggest power advantages of the resulting leo.nb.oca score. application iii: fsp is upstream of a biologically interesting gene co-expression module in female bxh mice here we illustrate how neo can be used to assign edge orienting scores to a single edge using manually and automatically chosen genetic markers. specifically, we computed edge orienting scores for the edge fsp → meblue where fsp corresponds to a pro-apoptotic gene that is related to metabolic syndrome: fsp27-null mice have been found to be resistant to obesity and diabetes; fsp expression is halved in obese humans after weight loss; and fsp regulates lipolysis in white human adipocytes. the other quantitative trait, meblue, represents the activation status of an entire pathway. more specifically, meblue is the module eigengene of the biologically important'blue' gene co-expression module described in. this gene co-expression module was comprised of highly correlated genes and meblue is a summary gene expression trait that best represents the expressions of the blue module genes. to study whether fsp causally affects meblue, we conducted both manual and automatic snp selection approaches. we used a previously identified snp marker on chromosome that affected the expression of the blue module genes and of several physiologic traits as the manually chosen input snp. this genetic marker was previously referred to as a module quantitative trait locus since it was found to affect the gene expression profiles of most blue module genes. using this snp, we found highly causal leo.nb scores between fsp and meblue. the leo.nb.oca scores passed the threshold of and the leo.nb.cpa scores passed the threshold of. we also used the automatic snp selection strategies to assess the causal relationship between fsp and meblue. the results of a robustness analysis can be found in figure we find that the causal relationship fsp → meblue is highly robust with respect to different automatic marker selection methods. simulation studies multi-edge simulation model that involves one hundred snps neo analysis can orient the edges of a multi-trait network by automatically selecting markers for each trait separately. the analyses proceed in a stepwise fashion: edges are oriented one at the time. for each edge, neo computes edge orienting scores. by thresholding these edge orienting scores, one can arrive at a globally oriented trait network. the details of the simulation model and relevant r code is presented in an r software tutorial on our webpage. briefly, we simulated a causal network between five gene expressions and a trait. each of the traits was simulated to be under the causal influence of snps. we added noise snps so that the data contained snps. we simulated the following causal relationships between the traits: e → e e → e e ← hiddenconfounder → e e → trait trait → e note that the correlation between traits e and e was entirely due a hidden confounder. the heatmap plot in figure depicts the true causal model. note that a red square in the i-th row and j-th column indicates that trait i causally affects trait j. the rows and columns of the heatmap are ordered according to a hierarchical clustering tree, which was constructed using average linkage hierarchical clustering with the dissimilarity diss = - |cor|. figure shows the corresponding heatmap of the observed network that was reconstructed using the leo.nb.oca score. figure shows an alternative output graph of neo. blue edges indicate significant correlations and a leo.nb.oca score is added to each edges whose leo.nb.oca score passes a user-supplied threshold. we find that all true causal edges are correctly retrieved at the recommended leo.nb.oca threshold of. figure shows the results of a robustness analysis for the leo.nb.oca and leo.nb.cpa scores for the edge orientation e → trait. the leo.nb.oca and the leo.nb.cpa scores exceed their respective threshold of and for all steps of the robustness analysis, i.e., they retrieve the orientation correctly. alternative simulation models can be explored using our online tutorial. single edge simulation model parameterized with the heritability in additional file we describe several simulation studies that use a single edge simulation model. briefly, we simulated two traits a and b that are anchored to genetic marker sets ma and mb, respectively. the correlation cor results from both a causal influence of b on a and from a hidden confounder c. this single edge model is used i) to study the choice of thresholds for the leo.nb scores, ii) to compare the leo.nb.cpa with leo.nb.oca scores, and iii) to evaluate automatic snp selection methods. the results of these simulations are described in additional file and in our online r software tutorials. discussion we propose methods for using multiple genetic markers to recover causal trait-trait relationships in systems genetic studies. neo will be particularly useful for the analysis of experiments in which common genetic variations are leveraged to explore complex genetic traits. we propose several edge orienting scores that measure the genetic evidence in favor of a given edge orientation a → b. while several methods exist for constructing undirected gene co-expression networks based on thousands of genes, we have evaluated the neo method for inferring directed networks involving relatively few genes. future research could explore the use of the method for inferring directed networks involving thousands of genes. our simulation studies show that orthogonal causal anchors lead to powerful edge scores that may outperform scores based only on candidate pleiotropic anchors. to afford flexibility to the user, the neo software provides several options for anchoring the traits to genetic markers, computing local edge scores, and diagnosing poor model fit. neo provides multiple options for automatically anchoring a trait to genetic markers: greedy, forward, and combined snp marker selection. while our simulation studies suggest that these three snp marker selection methods have similar performance, we find that the combined snp marker selection performs best when signal snps are in high linkage disequilibrium with noise snps. neo's local, stepwise approach for orienting edges of a trait network allows one to orient networks involving hundreds or even thousands of traits. since the calculation of edge orienting scores is based on local causal models, neo is relatively robust with regard to mistaken orientation of some edges in the global network. although neo performs well in simulation studies and the reported real data applications, we note that it has several limitations. the first limitation is that it requires the availability of genetic markers that are significantly associated with at least one trait per edge. spurious associations between the markers and traits will result in meaningless edge orienting scores. although the multi-marker score is quite robust to noise snps in our simulations, false-positive input snps will result in unreliable edge scores. the automatic snp selection is particularly vulnerable to false positives and its results should be carefully validated using biological experiments or causality analysis of independent data. the second limitation is that the resulting global trait network may contain loops, i.e. it may be cyclic. in contrast, a directed acyclic graph has no cycles. dags appear in models where it does not make sense for a trait to have a path to itself. while local dags are used for orienting individual edges, the reconstructed global trait network may no longer be acyclic. acyclicity is theoretically desirable since it allows one to test causal predictions using pearl's formalism of d-separation. the constraint of acyclic graphs in many network learning algorithms is often more a mathematical convenience than reflective of biology; cycles may reflect feedback loops for maintaining homeostasis. when more and more edges are oriented, as in the ic/ic* and pc/pc* algorithms, an error in one part of the network can propagate and cause erroneous orientations in unrelated portions of the network. most often these errors arise due to confusion between confounded and truly causal flows. to avoid being misled, neo deliberately discards the evidence from correlated trait neighbors in the undirected graph during leo scoring. by computing local edge orienting scores without regard to a global acyclicity constraint, the analysis is relatively robust to mis-oriented neighboring edges. neo uses causal anchors for each edge separately and thus allows the genetic data to speak for themselves. the proposed leo.nb scores are local in that they orient one edge at a time without regard to the orientations of the other edges. the reconstruction of the global network should be taken with a grain of salt. while we report one simulation model where the global network was reconstructed correctly, future research should carefully evaluate the performance of the neo approach for inferring global networks. a potential use of neo is to use it for initializing an iterative edge orienting algorithms for large networks that maximizes a global sem fitting index. the third limitation is that the sem-based edge orienting scores assume linear relationships between traits and snp markers. this is mathematically convenient but non-linear effects are common and have been reported in the literature. the neo approach works in the domain of linear graphical models since it is based on correlations and sems. akin to the use of pearson versus spearman correlation, the software also offers the option of modelling monotone quantitative relationships in neo by converting all data to ranks before further processing. neo will not work for traits that satisfy non-monotone relationships. a fourth limitation is that the influence of genetic markers may be indirect. neo may miss some relationships. while snp changes must be upstream of gene expression and phenotype manifestations, this does not preclude some snps from modifying the action of other snps, and the effect of such modifiers may become apparent only in particular contexts. causal inference and structural equation modeling assume that relevant traits and causal anchors have been included in the causal model. under-specified causal models, i.e. models that omit important variables, may mislead the user to detect spurious causal relationships. neo leads to relatively simple causal networks that do not incorporate dynamic or hierarchical properties. given all these potential limitations it is reassuring that neo performs well at retrieving known causal relationships in the reported real data applications. since neo focuses on individual edges, we expect that neo will be particularly useful for identifying traits that are causal for a given trait. for example, we illustrate that neo can be used to identify gene expression traits that are reactive to insig the sem r package can be used to evaluate the global fit of an acyclic multi-trait network. the neo algorithm computes an edge score for each edge without regard to the information gained from neighboring edges. neo aims to harness the power of the established upstream causal anchors as fully as possible; thus, it is appropriate when genetic variations are a major source of the variation in the traits. to the extent the environment is also varied, the neo approach may be less effective. it is plausible that additional assumptions may allow one to use unshielded colliders to improve the causal inference. this is a promising avenue of future research. we focused on the use of snp markers which capture only a limited amount of the sequence information of each individual. in the not too distant future, it will be economically feasible to obtain the sequence information of each study subject. since sequence information is likely to enhance the causal anchor assignment, sequence data may greatly improve the power of the neo method. apart from the common genetic variation that perturbs gene expression in mouse crosses, neo can also be applied to orient edges on the basis of causal anchors from population-based allelic association studies, cell hybrids, or transfected cells.  natural randomization of alleles that occurs during meiosis can be used to study the causal information flow through trait networks. for example, we use mouse cross data to retrieve known causal relationships in the sterol biosynthesis pathway. we find that the proposed edge scores are quite robust with respect to adding extraneous noise snps. combined with the use of orthogonal causal anchors, the proposed edge orienting scores can provide a strong basis for further experimental evaluation of the predicted causal relationships. 
the neo software can be used to orient the edges of gene co-expression networks or quantitative trait networks if the edges can be anchored to genetic marker data. we develop and implement network edge orienting methods and software that address the challenges of inferring unconfounded and directed gene networks from microarray-derived gene expression data by integrating mrna levels with genetic marker data and structural equation model comparisons. we demonstrate the utility of neo by recovering known causal relationships in the sterol homeostasis pathway using liver gene expression data from an f mouse cross. several authors have argued that genetic markers can serve as causal anchors for orienting the edges of a trait network. the availability of hundreds of thousands of genetic markers poses new challenges: how to relate traits to multiple genetic markers, how to score the genetic evidence in favor of an edge orientation, and how to weigh the information from multiple markers. the neo software implements several manual and automatic methods for incorporating genetic information to anchor traits. the neo software allows the user to carry out a robustness analysis with regard to genetic marker selection. systems genetic studies have been used to identify genetic loci that affect transcript abundances and clinical traits such as body weight. to summarize the genetic evidence in favor of a given edge orientation, we propose local sem-based edge orienting scores that compare the fit of several competing causal graphs. the pairwise correlations between gene expression traits and/or clinical traits can be used to define undirected trait networks.
 systems genetic studies have been used to identify genetic loci that affect transcript abundances and clinical traits such as body weight. the pairwise correlations between gene expression traits and/or clinical traits can be used to define undirected trait networks. several authors have argued that genetic markers can serve as causal anchors for orienting the edges of a trait network. the availability of hundreds of thousands of genetic markers poses new challenges: how to relate traits to multiple genetic markers, how to score the genetic evidence in favor of an edge orientation, and how to weigh the information from multiple markers.  we develop and implement network edge orienting methods and software that address the challenges of inferring unconfounded and directed gene networks from microarray-derived gene expression data by integrating mrna levels with genetic marker data and structural equation model comparisons. the neo software implements several manual and automatic methods for incorporating genetic information to anchor traits. the networks are oriented by considering each edge separately, thus reducing error propagation. to summarize the genetic evidence in favor of a given edge orientation, we propose local sem-based edge orienting scores that compare the fit of several competing causal graphs. sem fitting indices allow the user to assess local and overall model fit. the neo software allows the user to carry out a robustness analysis with regard to genetic marker selection. we demonstrate the utility of neo by recovering known causal relationships in the sterol homeostasis pathway using liver gene expression data from an f mouse cross. further, we use neo to study the relationship between a disease gene and a biologically important gene co-expression module in liver tissue.  the neo software can be used to orient the edges of gene co-expression networks or quantitative trait networks if the edges can be anchored to genetic marker data. r software tutorials, data, and supplementary material can be downloaded from:. 