This is a demonstration of the effect of different normalisation methods available in the minfi package for the analysis of Illumina 450k data. In this case large global methylation changes are seen between samples, so can between-array normalisation approaches be applied without biasing the data? Further information and explanation can be see at this Biostars discussion.
Data is loaded into an RGset object.
RGset <- read.450k.exp( targets=sampleInfo.df )
RGset
## RGChannelSet (storageMode: lockedEnvironment)
## assayData: 622399 features, 24 samples
## element names: Green, Red
## An object of class 'AnnotatedDataFrame'
## sampleNames: 10006823075_R01C01 10006823075_R02C01 ...
## 10006823118_R06C02 (24 total)
## varLabels: Sample_Well Slide ... filenames (18 total)
## varMetadata: labelDescription
## Annotation
## array: IlluminaHumanMethylation450k
## annotation: ilmn12.hg19
Apply the normalisation methods
MSet.funnorm <- preprocessFunnorm(RGset) #using Functional Normalistion
MSet.norm <- preprocessIllumina(RGset, bg.correct = TRUE, normalize = "controls", reference = 1) #using Illumina's method
MSet.noob <- preprocessNoob(RGset) #using Noob
Compare the output by doing some mds plots
library(RColorBrewer)
library(gridExtra)
palette(brewer.pal(9, 'Set1'))
par(mfrow=c(2,2))
try(mdsPlot(getBeta(RGset), sampGroups = pData(RGset)$GROUP_ANON, legendPos = 'none', main='MDS plot using RGset', pch=16 ),silent=TRUE)
try(mdsPlot(getBeta(MSet.noob), sampGroups = pData(MSet.noob)$GROUP_ANON, legendPos = 'none', main='MDS plot using Noob', pch=16),silent=TRUE)
try(mdsPlot(getBeta(MSet.funnorm), sampGroups = pData(MSet.funnorm)$GROUP_ANON, legendPos = 'none', main='MDS plot using FunNorm', pch=16),silent=TRUE)
try(mdsPlot(getBeta(MSet.norm), sampGroups = pData(MSet.norm)$GROUP_ANON, legendPos = 'none', main='MDS plot using Illumina', pch=16), silent=TRUE)
And some density plots
par(mfrow=c(2,2))
densityPlot(getBeta(RGset), sampGroups = pData(RGset)$GROUP, main = "Density plot using RGset", xlab = "Beta", legend = FALSE)
densityPlot(getBeta(MSet.noob), sampGroups = pData(MSet.noob)$GROUP, main = "Density plot using Noob", xlab = "Beta", legend = FALSE)
densityPlot(getBeta(MSet.funnorm), sampGroups = pData(MSet.funnorm)$GROUP, main = "Density plot using FunNorm", xlab = "Beta", legend = FALSE)
densityPlot(getBeta(MSet.norm), sampGroups = pData(MSet.funnorm)$GROUP, main = "Density plot using Illumina", xlab = "Beta", legend = FALSE)
Would like to be able to distinguish between the density plots from different groups more easily, so define a custom density function:
customDensityMinfi
## function ( MSet.plot) {
##
## require(data.table)
## pd <- pData(MSet.plot) %>% as.data.frame() %>%
## mutate(SampleID=paste(Slide, Array, sep='_')) %>%
## dplyr::select(SampleID, Day, TREAT1_ANON, TREAT2_ANON)
## pd <- as.data.table(pd)
##
## ad <- getBeta(MSet.plot) %>% melt()
## colnames(ad) <- c('ProbeID', 'SampleID', 'beta')
## ad$SampleID <- as.character(ad$SampleID)
## ad <- as.data.table(ad)
##
## plot_data <- pd %>% inner_join(ad, by='SampleID')
##
## ggplot (plot_data , aes(x=beta, colour=as.factor(Day), linetype=as.factor(TREAT2_ANON))) +
## geom_density() +
## facet_wrap ( ~ TREAT1_ANON ) +
## theme_bw() + theme(legend.position='none') +
## xlab(label = 'beta value')
##
## }
Apply this to our four objects:
customDensityMinfi (RGset) + ggtitle('Density plot using RGset')
## Loading required package: data.table
## data.table 1.9.4 For help type: ?data.table
## *** NB: by=.EACHI is now explicit. See README to restore previous behaviour.
##
## Attaching package: 'data.table'
##
## The following objects are masked from 'package:dplyr':
##
## between, last
customDensityMinfi (MSet.noob) + ggtitle('Density plot using Noob')
customDensityMinfi (MSet.funnorm) + ggtitle('Density plot using Funnnorm')
customDensityMinfi (MSet.norm) + ggtitle('Density plot using Illumina')
