Last updated: 2016-04-13

Code version: 6738dfffa82af813517c5f8e732635ed66d8de3d

Test Model Performance with Real Data

The purpose of this file is to test the performance of the Gibbs sampler with real data. Presumably gene epxression data for site \(g\) is \(G_g \approx \theta_{cis} X_g + \epsilon_g\), meaning the mean standardized expression will approximately be \(G_g - \bar{G}_g \approx \theta_{cis} \left(X_g - \bar{X}_g\right) + \epsilon_g\). If \(\theta_{cis}\) is small then I should be able to simulate my own effects without much influence from the genotype. I will simulate data using one of the raw expression files Brian sent by: 1.) Center the true gene epxression matrix \(G \in \mathbb{R}^{n_{\text{ind}} \times n_{\text{genes}}}\) 2.) Chose an anchor gene \(g\) and chose the \(t\) most correlated genes using \(\Sigma_g = \frac{1}{n_{\text{ind}}} G^T G\) as a metric. 3.) Simulate genotypes \(X \in \mathbb{R}^{n_{\text{ind}} \times 1}\) for the anchor gene \(g\). 4.) Fix the cis-eQTL effect \(\theta_g > 0\) and partition the \(t\) nodes into D, I, U nodes. 5.) Simulate \(Y \in \mathbb{R}^{(t+1) \times 1}\) using \(\Sigma_g\) as the correlation matrix.

For these simulations, I will shrink the indirect effects by a constant \(\lambda = 1\).

Upload Source Functions

library('gtools')
source("../R/directionality_CGM.R")
source("../R/SimulateNetworkExpression.R")
source("../R/Model2.R")

Load Data

path.data <- "../data/raw_expression/skinsunexposedlowerleg_hcp_corrected_factors_genes_certain_biotypes_TPM.txt"
G.raw <- data.frame(read.table(path.data, sep="\t", dec=".", header=T, check.names=T))
n.ind <- ncol(G.raw) - 1    #Number of individuals
p <- nrow(G.raw)        #Number of genes
G.centered <- as.matrix(G.raw[,2:(n.ind+1)]) - cbind( apply(as.matrix(G.raw[,2:(n.ind+1)]), 1, mean) ) %*% rbind( rep(1, n.ind) )

Create Networks

Create networks of up to \(t\) highly correlated genes with at most \(m\) neighbors

m <- 4    #Minimum number of neighbors
t <- 9    #Maximum number of neighbors
cor.thresh <- 0.8   #Correlation threshold to be considered a neighbor

Network.genes <- list()     #The first index is the anchor gene
Network.Exp <- list()       #Raw gene expression of genes in the network
count <- 1
for (g in 1:p) {
  G.g <- G.centered[g,]
  corr.g <- G.g %*% t(G.centered) / n.ind
  corr.g <- corr.g / corr.g[g]
  ind.g <- which(abs(corr.g) >= cor.thresh)
  if (length(ind.g) >= m+1) {
    ind.g.sorted <- order(-corr.g)[1:min(t+1, length(ind.g))]
    Network.genes[[count]] <- ind.g.sorted
    Network.Exp[[count]] <- G.centered[ind.g.sorted,]   #an n.nei + 1 x n.ind matrix
    count <- count + 1
  }
}
n.networks <- count - 1

Simulate and Analyze Data

Define parameters for simulation

theta <- 0.4    #effect size = theta * standard deviation
maf <- 0.1      #minor allele frequency
D <- c(1)    #Direct effect indices
I <- c(2:m)    #Indirect effect indices
lambda <- 1        #Shrinkage parameter for indirect effects

n.D <- length(D)
n.I <- length(I)

###Some Parameters for the Gibbs Sampler###
n.iter <- 2000
n.burn <- 1000
theta.gibbs <- c(0.5,0.5,0.5)  #If Dirichlet = T, this is the prior on alpha. Otherwise, this is proportional to the probability of unaffected, indirectly affected, directly affected by the cis-eQTL
Dirichlet <- T

sigma.a <- c(0.2, 0.4, 0.6)         #Prior sd's on effect; used in BF calculation
weights.sigma <- c(1,1,1)                   #Relative weights of each element of sigma.a

Simulate and analyze data

SimResults.cgm <- list()     #posterior probability unaffected from CGM's method (condition on Sigma)
SimResults.stephens <- list()      #posterior probability unaffected from Matthew's method (Wishart prior on Sigma)
for (ind.net in 1:n.networks) {
  net.genes <- Network.genes[[ind.net]]
  raw.expr <- Network.Exp[[ind.net]]
  
  n.nei <- nrow(raw.expr) - 1
  m.wishart <- n.nei + 1    ##DOF used in Wishart prior in Bayes Factor analysis
  U <- ((m+1):(n.nei+1))
  n.U <- length(U)
  
  X.genotype <- rbinom(n.ind, 1, maf) + rbinom(n.ind, 1, maf)   #Genotype vector
  
  Sigma.net <- raw.expr %*% t(raw.expr) / n.ind   #Covariance matrix for the simulated network (needed for indirect effects)
  Sigma.12 <- Sigma.net[I, c(D, U)]
  Sigma.22 <- Sigma.net[c(D, U), c(D, U)]
  tmp.vec <- lambda * Sigma.12 %*% cbind( solve(Sigma.22, c(rep(1, n.D), rep(0, n.U))) )
  
  effect.size <- theta * Sigma.net[1,1]    #Effect of genotype at source gene
  
  mu.D <- effect.size * matrix( rep(X.genotype, n.D), nrow=n.D, ncol=n.ind, byrow=T )   #Mean for direct effect
  mu.I <- effect.size * lambda * cbind(tmp.vec) %*% rbind(X.genotype)    #Mean for indirect effects
  
  Y.sim <- raw.expr
  Y.sim[D,] <- Y.sim[D,] + mu.D
  Y.sim[I,] <-  Y.sim[I,] + mu.I   #A (n.nei + 1) x n.ind matrix
  
  suff.stat <- Suff.stat(t(Y.sim), X.genotype)
  
  gibbs.cgm <- Gibbs.dir.cgm(n.iter, n.burn, suff.stat, D[1], n.ind, sigma.a, weights.sigma, theta=theta.gibbs, dirichlet=Dirichlet, update.lambda=F)
  gibbs.stephens <- Gibbs.dir.1(n.iter, n.burn, suff.stat, D[1], n.ind, sigma.a, weights.sigma, m.wishart, theta.gibbs, Dirichlet)
  
  SimResults.cgm[[ind.net]] <- gibbs.cgm$post.mean.U
  SimResults.stephens[[ind.net]] <- gibbs.stephens$post.mean.U
}

Visualize Simulation Results

Look at: 1.) estimated conditional FDR vs. true FDR for both CGM’s and Matthew’s models 2.) True positive rate as a function of true FDR for both CGM’s and Matthew’s models

post.U.D.cgm <- c()
post.U.I.cgm <- c()
post.U.U.cgm <- c()
post.U.D.stephens <- c()
post.U.I.stephens <- c()
post.U.U.stephens <- c()

for (ind.net in 1:n.networks) {
  post.U.cgm <- SimResults.cgm[[ind.net]]
  post.U.stephens <- SimResults.stephens[[ind.net]]
  
  n.nei <- length(post.U.cgm) - 1
  U <- ((m+1):(n.nei+1))
  n.U <- length(U)
  
  if (n.D > 1) {
    post.U.D.cgm <- c(post.U.D.cgm, post.U.cgm[2:n.D])
  }
  post.U.I.cgm <- c(post.U.I.cgm, post.U.cgm[I])
  post.U.U.cgm <- c(post.U.U.cgm, post.U.cgm[U])
  
  if (n.D > 1) {
    post.U.D.stephens <- c(post.U.D.stephens, post.U.stephens[2:n.D])
  }
  post.U.I.stephens <- c(post.U.I.stephens, post.U.stephens[I])
  post.U.U.stephens <- c(post.U.U.stephens, post.U.stephens[U])
}

all.data.cgm <- c(post.U.D.cgm, post.U.I.cgm, post.U.U.cgm); order.data.cgm <- order(all.data.cgm)   #Posterior probabilities genes are UNAFFECTED by SNP
gene.label.cgm <- c(rep(2, (n.D-1)*n.networks), rep(1, n.I*n.networks), rep(0, length(all.data.cgm) - (n.D-1)*n.networks - n.I*n.networks))[order.data.cgm]
all.data.cgm <- sort(all.data.cgm)

all.data.stephens <- c(post.U.D.stephens, post.U.I.stephens, post.U.U.stephens); order.data.stephens <- order(all.data.stephens)   #Posterior probabilities genes are UNAFFECTED by SNP
gene.label.stephens <- c(rep(2, (n.D-1)*n.networks), rep(1, n.I*n.networks), rep(0, length(all.data.stephens) - (n.D-1)*n.networks - n.I*n.networks))[order.data.stephens]
all.data.stephens <- sort(all.data.stephens)

fdr.all.cgm <- Cond.FDR(all.data.cgm)
fdr.labels.all.cgm <- FDR.labels(gene.label.cgm)
sens.all.cgm <- Sens.labels(gene.label.cgm)
plot(fdr.all.cgm, fdr.labels.all.cgm, xlab="Conditional False Discovery Rate From Gibbs Sampler", ylab="True False Discovery Rate", main="False Discovery Rate in CGM's Method", type="l")
abline(a=0,b=1, col="red")

plot(fdr.labels.all.cgm, sens.all.cgm, xlab="True False Discovery Rate", ylab="Sensitivity = Fraction of Indirect and Directly Affected Genes Captured", main="ROC Plot for CGM's Method", type="l")

hist(1 - post.U.I.cgm, xlab="Posterior Probability Gene is Affected by SNP", main="Indirect Effect Sensitivity in CGM's Method")

fdr.all.stephens <- Cond.FDR(all.data.stephens)
fdr.labels.all.stephens <- FDR.labels(gene.label.stephens)
sens.all.stephens <- Sens.labels(gene.label.stephens)
plot(fdr.all.stephens, fdr.labels.all.cgm, xlab="Conditional False Discovery Rate From Gibbs Sampler", ylab="True False Discovery Rate", main="False Discovery Rate in Matthew's Method", type="l")
abline(a=0,b=1, col="red")

plot(fdr.labels.all.stephens, sens.all.stephens, xlab="True False Discovery Rate", ylab="Sensitivity = Fraction of Indirect and Directly Affected Genes Captured", main="ROC Plot for Matthews's Method", type="l")

hist(1 - post.U.I.stephens, xlab="Posterior Probability Gene is Affected by SNP", main="Indirect Effect Sensitivity in Matthews's Method")

Conclusion

We see that if we use real data as residuals and no indirect effect shrinkage (i.e. \(\lambda = 1\)), we get the right answer using both models assuming moderate effect sizes (standardize effect size \(=\) 0.4 \(\times\) sd of expression for anchor gene).

Session information

sessionInfo()
R version 3.1.3 (2015-03-09)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: OS X 10.10.5 (Yosemite)

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] gtools_3.4.2 knitr_1.12.3

loaded via a namespace (and not attached):
 [1] digest_0.6.8  evaluate_0.8  formatR_1.2.1 htmltools_0.3 magrittr_1.5 
 [6] rmarkdown_0.7 stringi_1.0-1 stringr_1.0.0 tools_3.1.3   yaml_2.1.13