acyclic: Acyclic graph example

In qtlnet: Causal Inference of QTL Networks

Description

We generate synthetic data (sample size 300) according to a DAG composed by 100 nodes and 107 edges (exactly as in Figure 1). Each phenotype node is affected by three QTLs, and we allow only additive genetic effects. The QTLs for each phenotype are randomly selected among 200 markers, with 10 markers unevenly distributed on each of 20 autosomes. We allowed different phenotypes to potentially share common QTLs. For each phenotype, the regression coefficients with other phenotypes are chosen uniformly between 0.5 and 1; QTL effects are chosen between 0.2 to 0.6; and residual standard deviations are chosen from 0.1 to 0.5. For each realization we apply the QDG algorithm to infer causal directions for the edges of the skeleton obtained by the PC-skeleton algorithm.

Usage

data(acyclic)

Details

For cyclic graphs, the output of the qdg function computes the log-likelihood up to the normalization constant (un-normalized log-likelihood). We can use the un-normalized log-likelihood to compare cyclic graphs with reversed directions (since they have the same normalization constant). However we cannot compare cyclic and acyclic graphs.

References

Chaibub Neto et al. (2008) Inferring causal phenotype networks from segregating populations. Genetics 179: 1089-1100.

See Also

sim.cross, sim.geno, sim.map, skeleton, qdg, graph.qdg, generate.qtl.pheno

Examples

## Not run: 
## This reproduces Figure 1 exactly.
set.seed(3456789)

tmp <- options(warn=-1)
acyclic.DG <- randomDAG(n = 100, prob = 2 / 99)

options(tmp)

## Simulate cross object using R/qtl routines.
n.ind <- 300
mymap <- sim.map(len=rep(100,20), n.mar=10, eq.spacing=FALSE, include.x=FALSE)
mycross <- sim.cross(map=mymap, n.ind=n.ind, type="f2")
summary(mycross)
mycross <- sim.geno(mycross,n.draws=1)


## Produce 100 QTL at three markers apiece.
acyclic.qtl <- generate.qtl.markers(cross=mycross,n.phe=100)

## Generate data from directed graph.
bp <- runif(100,0.5,1)
stdev <- runif(100,0.1,0.5)
bq <- matrix(0,100,3)
bq[,1] <- runif(100,0.2,0.4)
bq[,2] <- bq[,1]+0.1
bq[,3] <- bq[,2]+0.1
## Generate phenotypes.
acyclic.data <- generate.qtl.pheno("acyclic", cross = mycross,
  bp = bp, bq = bq, stdev = stdev, allqtl = acyclic.qtl$allqtl)

acyclic.qdg <- qdg(cross=acyclic.data, 
		phenotype.names=paste("y",1:100,sep=""),
		marker.names=acyclic.qtl$markers, 
		QTL=acyclic.qtl$allqtl, 
		alpha=0.005, 
		n.qdg.random.starts=1,
		skel.method="pcskel")
save(acyclic.DG, acyclic.qtl, acyclic.data, acyclic.qdg,
  file = "acyclic.RData", compress = TRUE)

data(acyclic)

dims <- dim(acyclic.data$pheno)
SuffStat <- list(C = cor(acyclic.data$pheno), n = dims[1])
pc <- skeleton(SuffStat, gaussCItest, p = dims[2], alpha = 0.005)
summary(pc)

summary(graph.qdg(acyclic.qdg))
gr <- graph.qdg(acyclic.qdg, include.qtl = FALSE)
plot(gr)

## End(Not run)

qtlnet documentation built on April 14, 2020, 6:24 p.m.