Nothing
inst/doc/IntroductionToHarman.R
In Harman: The removal of batch effects from datasets using a PCA and constrained optimisation based technique
## ----style, eval=TRUE, echo = FALSE, results = 'asis'-------------------------
BiocStyle::markdown()
## ---- echo=FALSE, warning=FALSE-----------------------------------------------
library(knitr)
library(RColorBrewer)
## -----------------------------------------------------------------------------
library(Harman)
library(HarmanData)
data(OLF)
## -----------------------------------------------------------------------------
table(olf.info)
## -----------------------------------------------------------------------------
olf.info[1:7,]
## -----------------------------------------------------------------------------
head(olf.data)
## ----runharman----------------------------------------------------------------
olf.harman <- harman(olf.data, expt=olf.info$Treatment, batch=olf.info$Batch, limit=0.95)
## -----------------------------------------------------------------------------
str(olf.harman)
## ---- fig.height=4, fig.width=7-----------------------------------------------
plot(olf.harman)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(olf.harman, main='Arrowplot of correction')
## -----------------------------------------------------------------------------
summary(olf.harman)
## ---- echo=FALSE--------------------------------------------------------------
kable(olf.harman$stats)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(olf.harman, 1, 8)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(olf.harman, 8, 9)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2))
pcaPlot(olf.harman, colBy='batch', pchBy='expt', main='Col by Batch')
pcaPlot(olf.harman, colBy='expt', pchBy='batch', main='Col by Expt')
## -----------------------------------------------------------------------------
olf.data.corrected <- reconstructData(olf.harman)
## -----------------------------------------------------------------------------
olf.data.remade <- reconstructData(olf.harman, this='original')
all.equal(as.matrix(olf.data), olf.data.remade)
## ---- fig.height=9, fig.width=7-----------------------------------------------
olf.data.diff <- abs(as.matrix(olf.data) - olf.data.corrected)
diff_rowvar <- apply(olf.data.diff, 1, var)
by_rowvar <- order(diff_rowvar)
par(mfrow=c(1,1))
image(x=1:ncol(olf.data.diff),
y=1:nrow(olf.data.diff),
z=t(olf.data.diff[by_rowvar, ]),
xlab='Sample',
ylab='Probeset',
main='Harman probe adjustments (ordered by variance)',
cex.main=0.7,
col=brewer.pal(9, 'Reds'))
## ----cleanup1, echo=FALSE-----------------------------------------------------
rm(list=ls()[grep("olf", ls())])
## -----------------------------------------------------------------------------
require(Harman)
require(HarmanData)
data(IMR90)
## ---- echo=FALSE--------------------------------------------------------------
kable(imr90.info)
## -----------------------------------------------------------------------------
table(expt=imr90.info$Treatment, batch=imr90.info$Batch)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2), mar=c(4, 4, 4, 2) + 0.1)
imr90.pca <- prcomp(t(imr90.data), scale. = TRUE)
prcompPlot(imr90.pca, pc_x=1, pc_y=2, colFactor=imr90.info$Batch,
pchFactor=imr90.info$Treatment, main='IMR90, PC1 v PC2')
prcompPlot(imr90.pca, pc_x=1, pc_y=3, colFactor=imr90.info$Batch,
pchFactor=imr90.info$Treatment, main='IMR90, PC1 v PC3')
## -----------------------------------------------------------------------------
imr90.hm <- harman(as.matrix(imr90.data), expt=imr90.info$Treatment,
batch=imr90.info$Batch)
## ---- echo=FALSE--------------------------------------------------------------
kable(imr90.hm$stats)
## ---- fig.height=4, fig.width=7-----------------------------------------------
plot(imr90.hm, pc_x=1, pc_y=2)
plot(imr90.hm, pc_x=1, pc_y=3)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(imr90.hm, pc_x=1, pc_y=3)
## ---- fig.height=4, fig.width=7-----------------------------------------------
imr90.hm <- harman(as.matrix(imr90.data), expt=imr90.info$Treatment,
batch=imr90.info$Batch, limit=0.9)
plot(imr90.hm, pc_x=1, pc_y=3)
## ---- echo=FALSE--------------------------------------------------------------
kable(imr90.hm$stats)
## ---- fig.height=4, fig.width=7-----------------------------------------------
imr90.data.corrected <- reconstructData(imr90.hm)
par(mfrow=c(1,2))
prcompPlot(imr90.data, 1, 3, colFactor=imr90.hm$factors$batch, cex=1.5, main='PCA, Original')
prcompPlot(imr90.data.corrected, 1, 3, colFactor=imr90.hm$factors$batch, cex=1.5, main='PCA, Corrected')
## ----cleanup2, echo=FALSE-----------------------------------------------------
rm(list=ls()[grep("imr90", ls())])
## ----affydata, warning=FALSE, message=FALSE, comment=FALSE--------------------
library(affydata, quietly = TRUE, warn.conflicts = FALSE)
data(Dilution)
Dilution.log <- Dilution
intensity(Dilution.log) <- log(intensity(Dilution))
cols <- brewer.pal(nrow(pData(Dilution)), 'Paired')
## ---- echo=FALSE--------------------------------------------------------------
kable(pData(Dilution))
## ---- fig.height=4, fig.width=4-----------------------------------------------
boxplot(Dilution, col=cols)
## ----normalize, fig.height=7, fig.width=9-------------------------------------
Dilution.loess <- normalize(Dilution.log, method="loess")
Dilution.qnt <- normalize(Dilution.log, method="quantiles")
# define a quick PCA and plot function
pca <- function(exprs, pc_x=1, pc_y=2, cols, ...) {
pca <- prcomp(t(exprs), retx=TRUE, center=TRUE)
if(is.factor(cols)) {
factor_names <- levels(cols)
num_levels <- length(factor_names)
palette <- rainbow(num_levels)
mycols <- palette[cols]
} else {
mycols <- cols
}
plot(pca$x[, pc_x], pca$x[, pc_y],
xlab=paste('PC', pc_x, sep=''),
ylab=paste('PC', pc_y, sep=''),
col=mycols, ...)
if(is.factor(cols)) {
legend(x=min(pca$x[, pc_x]), y=max(pca$x[, pc_y]),
legend=factor_names, fill=palette, cex=0.5)
}
}
par(mfrow=c(1,2), mar=c(4, 4, 4, 2) + 0.1)
boxplot(Dilution.loess, col=cols, main='Loess')
boxplot(Dilution.qnt, col=cols, main='Quantile')
pca(exprs(Dilution.loess), cols=cols, cex=1.5, pch=16, main='PCA Loess')
pca(exprs(Dilution.qnt), cols=cols, cex=1.5, pch=16, main='PCA Quantile')
## ----affydata_harman1---------------------------------------------------------
library(Harman)
loess.hm <- harman(exprs(Dilution.loess),
expt=pData(Dilution.loess)$liver,
batch=pData(Dilution.loess)$scanner)
qnt.hm <- harman(exprs(Dilution.qnt),
expt=pData(Dilution.qnt)$liver,
batch=pData(Dilution.qnt)$scanner)
## ---- echo=FALSE, fig.align='left'--------------------------------------------
kable(loess.hm$stats)
## ---- echo=FALSE, fig.align='left'--------------------------------------------
kable(qnt.hm$stats)
## ----affydata_harman2---------------------------------------------------------
loess.hm <- harman(exprs(Dilution.loess),
expt=pData(Dilution.loess)$liver,
batch=pData(Dilution.loess)$scanner,
limit=0.75)
qnt.hm <- harman(exprs(Dilution.qnt),
expt=pData(Dilution.qnt)$liver,
batch=pData(Dilution.qnt)$scanner,
limit=0.75)
## ---- echo=FALSE--------------------------------------------------------------
kable(loess.hm$stats)
## ---- echo=FALSE--------------------------------------------------------------
kable(qnt.hm$stats)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(2,2), mar=c(4, 4, 4, 2) + 0.1)
plot(loess.hm, 1, 2, pch=17, col=cols)
plot(qnt.hm, 1, 2, pch=17, col=cols)
## ----affydata_reconstruct-----------------------------------------------------
Dilution.loess.hm <- Dilution.loess
Dilution.qnt.hm <- Dilution.qnt
exprs(Dilution.loess.hm) <- reconstructData(loess.hm)
exprs(Dilution.qnt.hm) <- reconstructData(qnt.hm)
## -----------------------------------------------------------------------------
detach('package:affydata')
## ---- message=FALSE, warning=FALSE--------------------------------------------
library(bladderbatch)
library(Harman)
# This loads an ExpressionSet object called bladderEset
data(bladderdata)
bladderEset
## ---- echo=FALSE--------------------------------------------------------------
kable(pData(bladderEset))
## -----------------------------------------------------------------------------
table(batch=pData(bladderEset)$batch, expt=pData(bladderEset)$outcome)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2))
prcompPlot(exprs(bladderEset), colFactor=pData(bladderEset)$batch,
pchFactor=pData(bladderEset)$outcome, main='Col by Batch')
prcompPlot(exprs(bladderEset), colFactor=pData(bladderEset)$outcome,
pchFactor=pData(bladderEset)$batch, main='Col by Expt')
## ---- cache=FALSE-------------------------------------------------------------
expt <- pData(bladderEset)$outcome
batch <- pData(bladderEset)$batch
bladder_harman <- harman(exprs(bladderEset), expt=expt, batch=batch)
sum(bladder_harman$stats$correction < 1)
## -----------------------------------------------------------------------------
summary(bladder_harman)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,1))
plot(bladder_harman)
## ---- fig.height=4, fig.width=7-----------------------------------------------
plot(bladder_harman, 1, 3)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(bladder_harman, 1, 3)
## -----------------------------------------------------------------------------
CorrectedBladderEset <- bladderEset
exprs(CorrectedBladderEset) <- reconstructData(bladder_harman)
comment(bladderEset) <- 'Original'
comment(CorrectedBladderEset) <- 'Corrected'
## ---- message=FALSE, warning=FALSE--------------------------------------------
library(limma, quietly=TRUE)
design <- model.matrix(~0 + expt)
colnames(design) <- c("Biopsy", "mTCC", "Normal", "sTCC", "sTCCwCIS")
contrast_matrix <- makeContrasts(sTCCwCIS - sTCC, sTCCwCIS - Normal,
Biopsy - Normal,
levels=design)
fit <- lmFit(exprs(bladderEset), design)
fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)
summary(decideTests(fit2))
## -----------------------------------------------------------------------------
fit_hm <- lmFit(exprs(CorrectedBladderEset), design)
fit2_hm <- contrasts.fit(fit_hm, contrast_matrix)
fit2_hm <- eBayes(fit2_hm)
summary(decideTests(fit2_hm))
## ----cleanup_bladderbatch, echo=FALSE-----------------------------------------
rm(list=ls()[!grepl("^pca$", ls())])
detach('package:bladderbatch')
## ---- message=FALSE, warning=FALSE--------------------------------------------
library(minfi, quietly = TRUE)
logit2
ilogit2
ilogit2(Inf)
ilogit2(-Inf)
## -----------------------------------------------------------------------------
betas <- seq(0, 1, by=0.05)
betas
newBetas <- shiftBetas(betas, shiftBy=1e-4)
newBetas
range(betas)
range(newBetas)
logit2(betas)
logit2(newBetas)
## ----loadminfidata, message=FALSE, warning=FALSE------------------------------
library(minfiData, quietly = TRUE)
baseDir <- system.file("extdata", package="minfiData")
targets <- read.metharray.sheet(baseDir)
## ----read450K-----------------------------------------------------------------
rgSet <- read.metharray.exp(targets=targets)
mSetSw <- preprocessSWAN(rgSet)
mSet <- preprocessIllumina(rgSet, bg.correct=TRUE, normalize="controls",
reference=2)
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mSetIl <- mSet[keep,]
mSetSw <- mSetSw[keep,]
## -----------------------------------------------------------------------------
kable(pData(mSetSw)[,c('Sample_Group', 'person', 'sex', 'status', 'Array',
'Slide')])
## ----mdsplots, fig.height=5, fig.width=7--------------------------------------
par(mfrow=c(1, 1))
cancer_by_gender_factor <- as.factor(paste(pData(mSetSw)$sex,
pData(mSetSw)$status)
)
mdsPlot(mSetSw, sampGroups=cancer_by_gender_factor, pch=16)
mdsPlot(mSetSw, sampGroups=as.factor(pData(mSet)$Slide), pch=16)
## -----------------------------------------------------------------------------
table(gender=pData(mSetSw)$sex, slide=pData(mSetSw)$Slide)
## -----------------------------------------------------------------------------
cancer_by_gender_factor <- as.factor(paste(pData(mSetSw)$sex,
pData(mSetSw)$status)
)
table(expt=cancer_by_gender_factor, batch=pData(mSetSw)$Slide)
## -----------------------------------------------------------------------------
table(expt=pData(mSetSw)$status, batch=pData(mSetSw)$Slide)
## ----shift, fig.height=7, fig.width=7-----------------------------------------
par(mfrow=c(2,2))
library(lumi, quietly = TRUE)
shifted_betas <- shiftBetas(betas=getBeta(mSetIl), shiftBy=1e-7)
shifted_ms <- beta2m(shifted_betas) # same as logit2() from package minfi
plot(density(shifted_ms, 0.05), main="Shifted M-values, shiftBy = 1e-7",
cex.main=0.7)
shifted_betas <- shiftBetas(betas=getBeta(mSetIl), shiftBy=1e-4)
shifted_ms <- beta2m(shifted_betas) # same as logit2() from package minfi
plot(density(shifted_ms, 0.05), main="Shifted M-values, shiftBy = 1e-4",
cex.main=0.7)
plot(density(beta2m(getBeta(mSetSw)), 0.05), main="SWAN normalised M-values",
cex.main=0.7)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2))
plot(density(shifted_betas, 0.1), main="Beta-values, shiftBy = 1e-4",
cex.main=0.7)
plot(density(shifted_ms, 0.1), main="M-values, shiftBy = 1e-4",
cex.main=0.7)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2))
plot(density(shifted_ms, 0.1),
main="GenomeStudio-like M-values, shiftBy = 1e-4", cex.main=0.7)
plot(density(getM(mSetSw), 0.1), main="SWAN M-values", cex.main=0.7)
## ---- fig.height=4, fig.width=7-----------------------------------------------
methHarman <- harman(getM(mSetSw), expt=pData(mSetSw)$status,
batch=pData(mSetSw)$Slide, limit=0.65)
summary(methHarman)
plot(methHarman, 2, 5)
## -----------------------------------------------------------------------------
ms_hm <- reconstructData(methHarman)
betas_hm <- m2beta(ms_hm)
## ----limma--------------------------------------------------------------------
library(limma)
group <- factor(pData(mSetSw)$status, levels=c("normal", "cancer"))
id <- factor(pData(mSetSw)$person)
design <- model.matrix(~id + group)
design
## ----fit----------------------------------------------------------------------
fit_hm <- lmFit(ms_hm, design)
fit_hm <- eBayes(fit_hm)
fit <- lmFit(getM(mSetSw), design)
fit <- eBayes(fit)
## ----limma_output-------------------------------------------------------------
summary_fit_hm <- summary(decideTests(fit_hm))
summary_fit <- summary(decideTests(fit))
summary_fit_hm
summary_fit
## ---- echo=FALSE, warning=FALSE, message=FALSE--------------------------------
rm(list=ls()[!grepl("^pca$", ls())])
#detach('package:minfiData', force = TRUE)
#detach('package:lumi', force = TRUE)
#detach('package:minfi', force = TRUE)
#detach('package:IlluminaHumanMethylation450kanno.ilmn12.hg19', force = TRUE)
## ----msms, message=FALSE, warning=FALSE---------------------------------------
# Call dependencies
library(msmsEDA)
library(Harman)
data(msms.dataset)
msms.dataset
## ----msms_pp------------------------------------------------------------------
# Preprocess to remove rows which are all 0 and replace NA values with 0.
msms_pp <- pp.msms.data(msms.dataset)
## ---- echo=FALSE--------------------------------------------------------------
kable(pData(msms_pp))
## -----------------------------------------------------------------------------
# Create experimental and batch variables
expt <- pData(msms_pp)$treat
batch <- pData(msms_pp)$batch
table(expt, batch)
## ----msms_harman--------------------------------------------------------------
# Log2 transform data, add 1 to avoid infinite values
log_ms_exprs <- log(exprs(msms_pp) + 1, 2)
# Correct data with Harman
hm <- harman(log_ms_exprs, expt=expt, batch=batch)
summary(hm)
## ---- fig.height=4, fig.width=7-----------------------------------------------
plot(hm)
## -----------------------------------------------------------------------------
# Reconstruct data and convert from Log2, removing 1 as we added it before.
log_ms_exprs_hm <- reconstructData(hm)
ms_exprs_hm <- 2^log_ms_exprs_hm - 1
# Place corrected data into a new 'MSnSet' instance
msms_pp_hm <- msms_pp
exprs(msms_pp_hm) <- ms_exprs_hm
## ---- echo=FALSE--------------------------------------------------------------
rm(list=ls()[!grepl("^pca$", ls())])
detach('package:msmsEDA')
## -----------------------------------------------------------------------------
table(imr90.info)
## ---- message=FALSE, warning=FALSE--------------------------------------------
library(HarmanData)
library(sva)
modcombat <- model.matrix(~1, data=imr90.info)
imr90.data.cb <- ComBat(dat=as.matrix(imr90.data), batch=imr90.info$Batch, mod=modcombat,
par.prior=TRUE, prior.plots=FALSE)
imr90.hr <- harman(imr90.data, expt = imr90.info$Treatment, batch = imr90.info$Batch)
imr90.data.hr <- reconstructData(imr90.hr)
prcompPlot(imr90.data, pc_x = 1, pc_y = 3,
colFactor = imr90.info$Batch,
pchFactor = imr90.info$Treatment,
main='PC1 versus PC3')
arrowPlot(imr90.hr, pc_x = 1, pc_y = 3, main='PC1 versus PC3')
## -----------------------------------------------------------------------------
library(RColorBrewer)
diffMap <- function(a, b, feature_order=1:nrow(a), batch, ...) {
image(t(abs(a[feature_order, ] - b[feature_order, ])),
col = brewer.pal(9, "Greys"),
xaxt='n',
yaxt='n',
xlab='Batch',
ylab='Ordered feature',
...)
axis(1, at=seq(0, 1, length.out=length(batch)), labels=batch)
}
## ---- fig.height=7, fig.width=7-----------------------------------------------
probe_diffs_hr <- order(rowSums(abs(as.matrix(imr90.data - imr90.data.hr))))
probe_diffs_cb <- order(rowSums(abs(as.matrix(imr90.data - imr90.data.cb))))
diffMap(a=imr90.data, b=imr90.data.hr, feature_order=probe_diffs_hr,
batch=imr90.info$Batch, main='Original - Harman')
## ---- fig.height=7, fig.width=7-----------------------------------------------
diffMap(a=imr90.data, b=imr90.data.cb, feature_order=probe_diffs_hr,
batch=imr90.info$Batch, main='Original - ComBat')
## ---- fig.height=7, fig.width=7-----------------------------------------------
diffMap(a=imr90.data.hr, b=imr90.data.cb, feature_order=probe_diffs_hr,
batch=imr90.info$Batch, main='Harman - ComBat')
Try the Harman package in your browser
Any scripts or data that you put into this service are public.
Harman documentation built on Nov. 8, 2020, 7:50 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----style, eval=TRUE, echo = FALSE, results = 'asis'-------------------------
BiocStyle::markdown()
## ---- echo=FALSE, warning=FALSE-----------------------------------------------
library(knitr)
library(RColorBrewer)
## -----------------------------------------------------------------------------
library(Harman)
library(HarmanData)
data(OLF)
## -----------------------------------------------------------------------------
table(olf.info)
## -----------------------------------------------------------------------------
olf.info[1:7,]
## -----------------------------------------------------------------------------
head(olf.data)
## ----runharman----------------------------------------------------------------
olf.harman <- harman(olf.data, expt=olf.info$Treatment, batch=olf.info$Batch, limit=0.95)
## -----------------------------------------------------------------------------
str(olf.harman)
## ---- fig.height=4, fig.width=7-----------------------------------------------
plot(olf.harman)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(olf.harman, main='Arrowplot of correction')
## -----------------------------------------------------------------------------
summary(olf.harman)
## ---- echo=FALSE--------------------------------------------------------------
kable(olf.harman$stats)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(olf.harman, 1, 8)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(olf.harman, 8, 9)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2))
pcaPlot(olf.harman, colBy='batch', pchBy='expt', main='Col by Batch')
pcaPlot(olf.harman, colBy='expt', pchBy='batch', main='Col by Expt')
## -----------------------------------------------------------------------------
olf.data.corrected <- reconstructData(olf.harman)
## -----------------------------------------------------------------------------
olf.data.remade <- reconstructData(olf.harman, this='original')
all.equal(as.matrix(olf.data), olf.data.remade)
## ---- fig.height=9, fig.width=7-----------------------------------------------
olf.data.diff <- abs(as.matrix(olf.data) - olf.data.corrected)
diff_rowvar <- apply(olf.data.diff, 1, var)
by_rowvar <- order(diff_rowvar)
par(mfrow=c(1,1))
image(x=1:ncol(olf.data.diff),
y=1:nrow(olf.data.diff),
z=t(olf.data.diff[by_rowvar, ]),
xlab='Sample',
ylab='Probeset',
main='Harman probe adjustments (ordered by variance)',
cex.main=0.7,
col=brewer.pal(9, 'Reds'))
## ----cleanup1, echo=FALSE-----------------------------------------------------
rm(list=ls()[grep("olf", ls())])
## -----------------------------------------------------------------------------
require(Harman)
require(HarmanData)
data(IMR90)
## ---- echo=FALSE--------------------------------------------------------------
kable(imr90.info)
## -----------------------------------------------------------------------------
table(expt=imr90.info$Treatment, batch=imr90.info$Batch)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2), mar=c(4, 4, 4, 2) + 0.1)
imr90.pca <- prcomp(t(imr90.data), scale. = TRUE)
prcompPlot(imr90.pca, pc_x=1, pc_y=2, colFactor=imr90.info$Batch,
pchFactor=imr90.info$Treatment, main='IMR90, PC1 v PC2')
prcompPlot(imr90.pca, pc_x=1, pc_y=3, colFactor=imr90.info$Batch,
pchFactor=imr90.info$Treatment, main='IMR90, PC1 v PC3')
## -----------------------------------------------------------------------------
imr90.hm <- harman(as.matrix(imr90.data), expt=imr90.info$Treatment,
batch=imr90.info$Batch)
## ---- echo=FALSE--------------------------------------------------------------
kable(imr90.hm$stats)
## ---- fig.height=4, fig.width=7-----------------------------------------------
plot(imr90.hm, pc_x=1, pc_y=2)
plot(imr90.hm, pc_x=1, pc_y=3)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(imr90.hm, pc_x=1, pc_y=3)
## ---- fig.height=4, fig.width=7-----------------------------------------------
imr90.hm <- harman(as.matrix(imr90.data), expt=imr90.info$Treatment,
batch=imr90.info$Batch, limit=0.9)
plot(imr90.hm, pc_x=1, pc_y=3)
## ---- echo=FALSE--------------------------------------------------------------
kable(imr90.hm$stats)
## ---- fig.height=4, fig.width=7-----------------------------------------------
imr90.data.corrected <- reconstructData(imr90.hm)
par(mfrow=c(1,2))
prcompPlot(imr90.data, 1, 3, colFactor=imr90.hm$factors$batch, cex=1.5, main='PCA, Original')
prcompPlot(imr90.data.corrected, 1, 3, colFactor=imr90.hm$factors$batch, cex=1.5, main='PCA, Corrected')
## ----cleanup2, echo=FALSE-----------------------------------------------------
rm(list=ls()[grep("imr90", ls())])
## ----affydata, warning=FALSE, message=FALSE, comment=FALSE--------------------
library(affydata, quietly = TRUE, warn.conflicts = FALSE)
data(Dilution)
Dilution.log <- Dilution
intensity(Dilution.log) <- log(intensity(Dilution))
cols <- brewer.pal(nrow(pData(Dilution)), 'Paired')
## ---- echo=FALSE--------------------------------------------------------------
kable(pData(Dilution))
## ---- fig.height=4, fig.width=4-----------------------------------------------
boxplot(Dilution, col=cols)
## ----normalize, fig.height=7, fig.width=9-------------------------------------
Dilution.loess <- normalize(Dilution.log, method="loess")
Dilution.qnt <- normalize(Dilution.log, method="quantiles")
# define a quick PCA and plot function
pca <- function(exprs, pc_x=1, pc_y=2, cols, ...) {
pca <- prcomp(t(exprs), retx=TRUE, center=TRUE)
if(is.factor(cols)) {
factor_names <- levels(cols)
num_levels <- length(factor_names)
palette <- rainbow(num_levels)
mycols <- palette[cols]
} else {
mycols <- cols
}
plot(pca$x[, pc_x], pca$x[, pc_y],
xlab=paste('PC', pc_x, sep=''),
ylab=paste('PC', pc_y, sep=''),
col=mycols, ...)
if(is.factor(cols)) {
legend(x=min(pca$x[, pc_x]), y=max(pca$x[, pc_y]),
legend=factor_names, fill=palette, cex=0.5)
}
}
par(mfrow=c(1,2), mar=c(4, 4, 4, 2) + 0.1)
boxplot(Dilution.loess, col=cols, main='Loess')
boxplot(Dilution.qnt, col=cols, main='Quantile')
pca(exprs(Dilution.loess), cols=cols, cex=1.5, pch=16, main='PCA Loess')
pca(exprs(Dilution.qnt), cols=cols, cex=1.5, pch=16, main='PCA Quantile')
## ----affydata_harman1---------------------------------------------------------
library(Harman)
loess.hm <- harman(exprs(Dilution.loess),
expt=pData(Dilution.loess)$liver,
batch=pData(Dilution.loess)$scanner)
qnt.hm <- harman(exprs(Dilution.qnt),
expt=pData(Dilution.qnt)$liver,
batch=pData(Dilution.qnt)$scanner)
## ---- echo=FALSE, fig.align='left'--------------------------------------------
kable(loess.hm$stats)
## ---- echo=FALSE, fig.align='left'--------------------------------------------
kable(qnt.hm$stats)
## ----affydata_harman2---------------------------------------------------------
loess.hm <- harman(exprs(Dilution.loess),
expt=pData(Dilution.loess)$liver,
batch=pData(Dilution.loess)$scanner,
limit=0.75)
qnt.hm <- harman(exprs(Dilution.qnt),
expt=pData(Dilution.qnt)$liver,
batch=pData(Dilution.qnt)$scanner,
limit=0.75)
## ---- echo=FALSE--------------------------------------------------------------
kable(loess.hm$stats)
## ---- echo=FALSE--------------------------------------------------------------
kable(qnt.hm$stats)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(2,2), mar=c(4, 4, 4, 2) + 0.1)
plot(loess.hm, 1, 2, pch=17, col=cols)
plot(qnt.hm, 1, 2, pch=17, col=cols)
## ----affydata_reconstruct-----------------------------------------------------
Dilution.loess.hm <- Dilution.loess
Dilution.qnt.hm <- Dilution.qnt
exprs(Dilution.loess.hm) <- reconstructData(loess.hm)
exprs(Dilution.qnt.hm) <- reconstructData(qnt.hm)
## -----------------------------------------------------------------------------
detach('package:affydata')
## ---- message=FALSE, warning=FALSE--------------------------------------------
library(bladderbatch)
library(Harman)
# This loads an ExpressionSet object called bladderEset
data(bladderdata)
bladderEset
## ---- echo=FALSE--------------------------------------------------------------
kable(pData(bladderEset))
## -----------------------------------------------------------------------------
table(batch=pData(bladderEset)$batch, expt=pData(bladderEset)$outcome)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2))
prcompPlot(exprs(bladderEset), colFactor=pData(bladderEset)$batch,
pchFactor=pData(bladderEset)$outcome, main='Col by Batch')
prcompPlot(exprs(bladderEset), colFactor=pData(bladderEset)$outcome,
pchFactor=pData(bladderEset)$batch, main='Col by Expt')
## ---- cache=FALSE-------------------------------------------------------------
expt <- pData(bladderEset)$outcome
batch <- pData(bladderEset)$batch
bladder_harman <- harman(exprs(bladderEset), expt=expt, batch=batch)
sum(bladder_harman$stats$correction < 1)
## -----------------------------------------------------------------------------
summary(bladder_harman)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,1))
plot(bladder_harman)
## ---- fig.height=4, fig.width=7-----------------------------------------------
plot(bladder_harman, 1, 3)
## ---- fig.height=5, fig.width=5-----------------------------------------------
arrowPlot(bladder_harman, 1, 3)
## -----------------------------------------------------------------------------
CorrectedBladderEset <- bladderEset
exprs(CorrectedBladderEset) <- reconstructData(bladder_harman)
comment(bladderEset) <- 'Original'
comment(CorrectedBladderEset) <- 'Corrected'
## ---- message=FALSE, warning=FALSE--------------------------------------------
library(limma, quietly=TRUE)
design <- model.matrix(~0 + expt)
colnames(design) <- c("Biopsy", "mTCC", "Normal", "sTCC", "sTCCwCIS")
contrast_matrix <- makeContrasts(sTCCwCIS - sTCC, sTCCwCIS - Normal,
Biopsy - Normal,
levels=design)
fit <- lmFit(exprs(bladderEset), design)
fit2 <- contrasts.fit(fit, contrast_matrix)
fit2 <- eBayes(fit2)
summary(decideTests(fit2))
## -----------------------------------------------------------------------------
fit_hm <- lmFit(exprs(CorrectedBladderEset), design)
fit2_hm <- contrasts.fit(fit_hm, contrast_matrix)
fit2_hm <- eBayes(fit2_hm)
summary(decideTests(fit2_hm))
## ----cleanup_bladderbatch, echo=FALSE-----------------------------------------
rm(list=ls()[!grepl("^pca$", ls())])
detach('package:bladderbatch')
## ---- message=FALSE, warning=FALSE--------------------------------------------
library(minfi, quietly = TRUE)
logit2
ilogit2
ilogit2(Inf)
ilogit2(-Inf)
## -----------------------------------------------------------------------------
betas <- seq(0, 1, by=0.05)
betas
newBetas <- shiftBetas(betas, shiftBy=1e-4)
newBetas
range(betas)
range(newBetas)
logit2(betas)
logit2(newBetas)
## ----loadminfidata, message=FALSE, warning=FALSE------------------------------
library(minfiData, quietly = TRUE)
baseDir <- system.file("extdata", package="minfiData")
targets <- read.metharray.sheet(baseDir)
## ----read450K-----------------------------------------------------------------
rgSet <- read.metharray.exp(targets=targets)
mSetSw <- preprocessSWAN(rgSet)
mSet <- preprocessIllumina(rgSet, bg.correct=TRUE, normalize="controls",
reference=2)
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mSetIl <- mSet[keep,]
mSetSw <- mSetSw[keep,]
## -----------------------------------------------------------------------------
kable(pData(mSetSw)[,c('Sample_Group', 'person', 'sex', 'status', 'Array',
'Slide')])
## ----mdsplots, fig.height=5, fig.width=7--------------------------------------
par(mfrow=c(1, 1))
cancer_by_gender_factor <- as.factor(paste(pData(mSetSw)$sex,
pData(mSetSw)$status)
)
mdsPlot(mSetSw, sampGroups=cancer_by_gender_factor, pch=16)
mdsPlot(mSetSw, sampGroups=as.factor(pData(mSet)$Slide), pch=16)
## -----------------------------------------------------------------------------
table(gender=pData(mSetSw)$sex, slide=pData(mSetSw)$Slide)
## -----------------------------------------------------------------------------
cancer_by_gender_factor <- as.factor(paste(pData(mSetSw)$sex,
pData(mSetSw)$status)
)
table(expt=cancer_by_gender_factor, batch=pData(mSetSw)$Slide)
## -----------------------------------------------------------------------------
table(expt=pData(mSetSw)$status, batch=pData(mSetSw)$Slide)
## ----shift, fig.height=7, fig.width=7-----------------------------------------
par(mfrow=c(2,2))
library(lumi, quietly = TRUE)
shifted_betas <- shiftBetas(betas=getBeta(mSetIl), shiftBy=1e-7)
shifted_ms <- beta2m(shifted_betas) # same as logit2() from package minfi
plot(density(shifted_ms, 0.05), main="Shifted M-values, shiftBy = 1e-7",
cex.main=0.7)
shifted_betas <- shiftBetas(betas=getBeta(mSetIl), shiftBy=1e-4)
shifted_ms <- beta2m(shifted_betas) # same as logit2() from package minfi
plot(density(shifted_ms, 0.05), main="Shifted M-values, shiftBy = 1e-4",
cex.main=0.7)
plot(density(beta2m(getBeta(mSetSw)), 0.05), main="SWAN normalised M-values",
cex.main=0.7)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2))
plot(density(shifted_betas, 0.1), main="Beta-values, shiftBy = 1e-4",
cex.main=0.7)
plot(density(shifted_ms, 0.1), main="M-values, shiftBy = 1e-4",
cex.main=0.7)
## ---- fig.height=4, fig.width=7-----------------------------------------------
par(mfrow=c(1,2))
plot(density(shifted_ms, 0.1),
main="GenomeStudio-like M-values, shiftBy = 1e-4", cex.main=0.7)
plot(density(getM(mSetSw), 0.1), main="SWAN M-values", cex.main=0.7)
## ---- fig.height=4, fig.width=7-----------------------------------------------
methHarman <- harman(getM(mSetSw), expt=pData(mSetSw)$status,
batch=pData(mSetSw)$Slide, limit=0.65)
summary(methHarman)
plot(methHarman, 2, 5)
## -----------------------------------------------------------------------------
ms_hm <- reconstructData(methHarman)
betas_hm <- m2beta(ms_hm)
## ----limma--------------------------------------------------------------------
library(limma)
group <- factor(pData(mSetSw)$status, levels=c("normal", "cancer"))
id <- factor(pData(mSetSw)$person)
design <- model.matrix(~id + group)
design
## ----fit----------------------------------------------------------------------
fit_hm <- lmFit(ms_hm, design)
fit_hm <- eBayes(fit_hm)
fit <- lmFit(getM(mSetSw), design)
fit <- eBayes(fit)
## ----limma_output-------------------------------------------------------------
summary_fit_hm <- summary(decideTests(fit_hm))
summary_fit <- summary(decideTests(fit))
summary_fit_hm
summary_fit
## ---- echo=FALSE, warning=FALSE, message=FALSE--------------------------------
rm(list=ls()[!grepl("^pca$", ls())])
#detach('package:minfiData', force = TRUE)
#detach('package:lumi', force = TRUE)
#detach('package:minfi', force = TRUE)
#detach('package:IlluminaHumanMethylation450kanno.ilmn12.hg19', force = TRUE)
## ----msms, message=FALSE, warning=FALSE---------------------------------------
# Call dependencies
library(msmsEDA)
library(Harman)
data(msms.dataset)
msms.dataset
## ----msms_pp------------------------------------------------------------------
# Preprocess to remove rows which are all 0 and replace NA values with 0.
msms_pp <- pp.msms.data(msms.dataset)
## ---- echo=FALSE--------------------------------------------------------------
kable(pData(msms_pp))
## -----------------------------------------------------------------------------
# Create experimental and batch variables
expt <- pData(msms_pp)$treat
batch <- pData(msms_pp)$batch
table(expt, batch)
## ----msms_harman--------------------------------------------------------------
# Log2 transform data, add 1 to avoid infinite values
log_ms_exprs <- log(exprs(msms_pp) + 1, 2)
# Correct data with Harman
hm <- harman(log_ms_exprs, expt=expt, batch=batch)
summary(hm)
## ---- fig.height=4, fig.width=7-----------------------------------------------
plot(hm)
## -----------------------------------------------------------------------------
# Reconstruct data and convert from Log2, removing 1 as we added it before.
log_ms_exprs_hm <- reconstructData(hm)
ms_exprs_hm <- 2^log_ms_exprs_hm - 1
# Place corrected data into a new 'MSnSet' instance
msms_pp_hm <- msms_pp
exprs(msms_pp_hm) <- ms_exprs_hm
## ---- echo=FALSE--------------------------------------------------------------
rm(list=ls()[!grepl("^pca$", ls())])
detach('package:msmsEDA')
## -----------------------------------------------------------------------------
table(imr90.info)
## ---- message=FALSE, warning=FALSE--------------------------------------------
library(HarmanData)
library(sva)
modcombat <- model.matrix(~1, data=imr90.info)
imr90.data.cb <- ComBat(dat=as.matrix(imr90.data), batch=imr90.info$Batch, mod=modcombat,
par.prior=TRUE, prior.plots=FALSE)
imr90.hr <- harman(imr90.data, expt = imr90.info$Treatment, batch = imr90.info$Batch)
imr90.data.hr <- reconstructData(imr90.hr)
prcompPlot(imr90.data, pc_x = 1, pc_y = 3,
colFactor = imr90.info$Batch,
pchFactor = imr90.info$Treatment,
main='PC1 versus PC3')
arrowPlot(imr90.hr, pc_x = 1, pc_y = 3, main='PC1 versus PC3')
## -----------------------------------------------------------------------------
library(RColorBrewer)
diffMap <- function(a, b, feature_order=1:nrow(a), batch, ...) {
image(t(abs(a[feature_order, ] - b[feature_order, ])),
col = brewer.pal(9, "Greys"),
xaxt='n',
yaxt='n',
xlab='Batch',
ylab='Ordered feature',
...)
axis(1, at=seq(0, 1, length.out=length(batch)), labels=batch)
}
## ---- fig.height=7, fig.width=7-----------------------------------------------
probe_diffs_hr <- order(rowSums(abs(as.matrix(imr90.data - imr90.data.hr))))
probe_diffs_cb <- order(rowSums(abs(as.matrix(imr90.data - imr90.data.cb))))
diffMap(a=imr90.data, b=imr90.data.hr, feature_order=probe_diffs_hr,
batch=imr90.info$Batch, main='Original - Harman')
## ---- fig.height=7, fig.width=7-----------------------------------------------
diffMap(a=imr90.data, b=imr90.data.cb, feature_order=probe_diffs_hr,
batch=imr90.info$Batch, main='Original - ComBat')
## ---- fig.height=7, fig.width=7-----------------------------------------------
diffMap(a=imr90.data.hr, b=imr90.data.cb, feature_order=probe_diffs_hr,
batch=imr90.info$Batch, main='Harman - ComBat')
Try the Harman package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.