# prep data
geneExpr = exprs(gset)
info = phenoData(gset)@data
colnames(info) = gsub(' ', '_', gsub(":ch1$", '', colnames(info)))

keep = grep("ESC", info$cell_line_of_origin, invert=TRUE)
geneExpr = geneExpr[,keep]
info = info[keep,]

keep = grep("fibro", info$differentiation_state, invert=TRUE)
geneExpr = geneExpr[,keep]
info = info[keep,]

info$Donor = gsub(".* (\\S+)$", "\\1", info$cell_line_of_origin)
info$cell_type[is.na(info$cell_type)] = 'neurospheres'

info$state = NA
info$state[info$differentiation_state == 'neurons at rest (day 45 of differentiation)'] = "rest"
info$state[info$differentiation_state == 'neurons kept in 67mM KCl for 9h'] = "KCl"

variancePartition

form <- ~ (1|differentiation_state) + (1|Donor) + (1|genotype)
vp = fitExtractVarPartModel( geneExpr, form, info)

plotVarPart(sortCols(vp))

dupCorList = foreach(ds = unique(info$differentiation_state) ) %do% {
  cat("\r", ds, '        ')  
  idx = info$differentiation_state == ds

  design = model.matrix( ~ genotype, info[idx,])
  dupcor <- duplicateCorrelation(geneExpr[,idx], design, block=info$Donor[idx])

  # dupcor = list(consensus=0.047)
  fitDupCor <- lmFit(geneExpr[,idx], design, block=info$Donor[idx], correlation=dupcor$consensus)
  fitDupCor = eBayes( fitDupCor )
  list(fitDupCor = fitDupCor, consensus=dupcor$consensus)
} 

neurons at rest (day 45 of differentiation) iPSC neurospheres (at day 7 in suspension) neurons kept in 67mM KCl for 9h

names(dupCorList) = unique(info$differentiation_state) 

dreamList = foreach(ds = unique( info$differentiation_state) ) %do% {
    
  cat("\r", ds, '        ')  
  idx = info$differentiation_state == ds

  # new syntax
  form <- ~ genotype + (1|Donor)     
  fitmm = dream( geneExpr[,idx], form, info[idx,])
  fitmm
}  

neurons at rest (day 45 of differentiation) Dividing work into 100 chunks... iPSC Dividing work into 100 chunks... neurospheres (at day 7 in suspension) Dividing work into 100 chunks... neurons kept in 67mM KCl for 9h Dividing work into 100 chunks...

names(dreamList) = unique(info$differentiation_state) 

figList = foreach( ds = names(dupCorList) ) %do% {

  # Compare p-values and make plot
  p1 = topTable(dupCorList[[ds]]$fitDupCor, coef="genotypeTS", number=Inf, sort.by="none")$P.Value
  p2 = topTable(dreamList[[ds]], coef="genotypeTS", number=Inf, sort.by="none")$P.Value
       
  plotCompareP( p1, p2, with(vpList[[ds]], Donor), dupCorList[[ds]]$consensus) + ggtitle( ds )
}
do.call("grid.arrange", c(figList, ncol=2))

Examine top difference between methods

library(lmerTest)
ds = names(dupCorList)[4]
         
# Find gene with greatest difference in -log10(p)
p1 = topTable(dupCorList[[ds]]$fitDupCor, coef="genotypeTS", number=Inf, sort.by="none")$P.Value
p2 = topTable(dreamList[[ds]], coef="genotypeTS", number=Inf, sort.by="none")$P.Value
j = which.max(abs(-log10(p1) + log10(p2)))

pvals = list()

# dream p-values
pvals[['dream']] = topTable(dreamList[[ds]], number=Inf, sort.by="none")[j,]$P.Value

# duplicateCorrelation results
pvals[['duplicateCorrelation']] = topTable(dupCorList[[ds]]$fitDupCor, number=Inf, sort.by="none")[j,]$P.Value

# p-value from lmer
idx = info$differentiation_state == ds
y = t(geneExpr[j,idx,drop=FALSE])
fitm = lmer( y ~ genotype + (1|Donor), info[idx,], REML=TRUE  )
# coef(summary(fitm))
pvals[['lmer()']] = coef(summary(fitm))[2,5]


# p-value from linear model
fit = lm( y ~ genotype, info[idx,]  )
# coef(summary(fit))
pvals[['lm()']] = coef(summary(fit))[2,4]

# linear model on individual averages
GE = data.table( Expression = as.numeric(y), Indiv =  info[idx,]$Donor, Genotype =  info[idx,]$genotype)
data = GE[,data.frame(Expression=mean(Expression)), by="Indiv"]
data = merge(data, unique(GE[,2:3]), by="Indiv", all.x=TRUE, all=FALSE)

pvals[['lm() on Individual averages (n=4)']] = coef(summary(lm( Expression ~ Genotype, data)))[2,4]

library(knitr)
kable(data.frame( p.value = sapply(pvals, function(x) format(x, digits=3, scientific=3))))

	p.value
dream	0.0684
duplicateCorrelation	7.96e-09
lmer()	0.0684
lm()	1.81e-08
lm() on Individual averages (n=4)	0.0044

We see that duplicateCorrelation and a simple linear model give very significant p-values for this probe (ILMN_2209027). However, dream and lmer() give this probe much larger p-values. Indeed, summarizing the data at the individual level and applying a linear model on the mean expression per individual for 4 individuals, gives a p-value more similar to dream/lmer. Thus, the dream/lmer method is approviately conservative and duplicateCorrelation/lm gives a very liberal p-value in this case.

Enrichments

Compare enrichment

Publication plot