In yunzhang813/yunr: Simulation study on deconvolution for network estimate

Introduction

## The package can be installed using the following code
# library(devtools)
# install_github("yunzhang813/???", build_vignettes=TRUE)

## Load the package
library(simDeNet)

Sample data

data("celltype")

One-step simulation

## get two cell types
mu.T <- expr[,ctab$Fastq_file_name[which(ctab$X3_letter_code=="ASM")]]
mu.N <- expr[,ctab$Fastq_file_name[which(ctab$X3_letter_code=="AEC")]]

## parameters
n.samp <- 5
rho <- c(0.9,0.8,0.7)
block.size <- c(5,10,15)
dd <- NULL
str.type <- "interchangeable"
multiplier <- 1
select.gene <- "random"

## mixing proportion of cell type T
prop.T <- seq(0, 1, length=n.samp)

## one-step simulation
set.seed(999)
out.oneStepSim <- oneStepSim(n.samp, mu.T, mu.N, Sigma.T=NULL, Sigma.N=NULL, prop.T=prop.T, 
                      # structure for Sigma.T
                      block.size=block.size, rho=rho, dd=dd, str.type=str.type, multiplier=1,
                      # selected genes to add structure
                      select.gene=select.gene)

Deconvolution

## deconvolution
tt <- system.time(out.deconv <- deconv(mixed=out.oneStepSim$expr.mixed, ref=out.oneStepSim$expr.pure.N))

Correlation-based network

One-step analysis for WGCNA

## simulated data
data.list <- list(out.oneStepSim$expr.pure.T, out.oneStepSim$expr.mixed, out.deconv$expr.deconv)
names(data.list) <- c("pure","mixed","deconvoluted")
## ture structure
true.str <- out.oneStepSim$true.str.T

## one step analysis
out.WGCNA <- oneStepAnalysis(data.list,true.str,method="WGCNA")

## plot ROC curves
plot(out.WGCNA$pure$perf, lwd=2)
plot(out.WGCNA$mixed$perf, add=TRUE, col=2, lwd=2)
plot(out.WGCNA$deconvoluted$perf, add=TRUE, col=3, lwd=2)
legend("bottomright", paste(c("pure","mixed","deconv"), ": AUC =", round(c(out.WGCNA$pure$AUC, out.WGCNA$mixed$AUC, out.WGCNA$deconvoluted$AUC),3)), col=1:3, lty=1, lwd=5)
title(paste("WGCNA"))

Step-by-step analysis for WGCNA

Network estimation

library(WGCNA)
## If softPower = 1, the adjacency matrix is simply the absolute correlation matrix
softPower <- 1
## pure samples
datExpr.pure <- as.data.frame(t(out.oneStepSim$expr.pure.T)) # NOTE: rows=samples, columns=genes
est.str.pure <- adjacency(datExpr.pure, power = softPower)
## mixed samples
datExpr.mixed <- as.data.frame(t(out.oneStepSim$expr.mixed)) # NOTE: rows=samples, columns=genes
est.str.mixed <- adjacency(datExpr.mixed, power = softPower)
## deconvoluted samples
datExpr.deconv <- as.data.frame(t(out.deconv$expr.deconv)) # NOTE: rows=samples, columns=genes
est.str.deconv <- adjacency(datExpr.deconv, power = softPower)

ROC analysis

## evaluate ROC and AUC
out.ROC.pure <- eval.ROC(est.str=est.str.pure, true.str=out.oneStepSim$true.str.T, 
                         plot.ROC=TRUE, main="Pure samples")
out.ROC.mixed <- eval.ROC(est.str=est.str.mixed, true.str=out.oneStepSim$true.str.T,
                          plot.ROC=TRUE, col=2, main="Mixed samples")
out.ROC.deconv <- eval.ROC(est.str=est.str.deconv, true.str=out.oneStepSim$true.str.T,
                           plot.ROC=TRUE, col=3, main="Deconvoluted samples")

## combined plot
library(ROCR)
plot(out.ROC.pure$perf, lwd=2)
plot(out.ROC.mixed$perf, add=TRUE, col=2, lwd=2)
plot(out.ROC.deconv$perf, add=TRUE, col=3, lwd=2)
legend("bottomright", paste(c("pure","mixed","deconv"), ": AUC =", round(c(out.ROC.pure$AUC, out.ROC.mixed$AUC, out.ROC.deconv$AUC),3)), col=1:3, lty=1, lwd=5)
title(paste("WGCNA", paste0("S",n.samp), str.type, select.gene, sep = ", "))
abline(0,1,col="grey")

Other networks

Mutual information network

Network estimation

library(minet)

## pure tumor samples
net.pure <- minet(datExpr.pure, method = "aracne")
## mixed samples
net.mixed <- minet(datExpr.mixed, method = "aracne")
## ISOpure tumor samples
net.deconv <- minet(datExpr.deconv, method = "aracne")

ROC analysis

## evaluate ROC and AUC
out.ROC.pure <- eval.ROC(est.str=net.pure, true.str=out.oneStepSim$true.str.T, 
                         plot.ROC=TRUE, main="Pure samples")
out.ROC.mixed <- eval.ROC(est.str=net.mixed, true.str=out.oneStepSim$true.str.T, 
                          plot.ROC=TRUE, main="Mixed samples", col=2)
out.ROC.deconv <- eval.ROC(est.str=net.deconv, true.str=out.oneStepSim$true.str.T, 
                           plot.ROC=TRUE, main="Deconvoluted samples", col=3)

## combined plot
plot(out.ROC.pure$perf,
     # xlim=c(0,max(out.ROC.pure$xval,out.ROC.mixed$xval,out.ROC.deconv$xval)),
     # ylim=c(0,max(out.ROC.pure$yval,out.ROC.mixed$yval,out.ROC.deconv$yval)),
     lwd=2)
plot(out.ROC.mixed$perf, lwd=2, add=TRUE, col=2)
plot(out.ROC.deconv$perf, lwd=2, add=TRUE, col=3)
legend("bottomright", paste(c("pure","mixed","deconv"), ": AUC =", round(c(out.ROC.pure$AUC, out.ROC.mixed$AUC, out.ROC.deconv$AUC),3)), col=1:3, lwd=5)
title(paste("ARACNE", paste0("S",n.samp), str.type, select.gene, sep = ", "))

yunzhang813/yunr documentation built on May 4, 2019, 7:45 p.m.