2.1.1 Running Harman

Okay, so let’s run Harman.

olf.harman <- harman(olf.data, expt=olf.info$Treatment, batch=olf.info$Batch, limit=0.95)

This creates a harmanresults S3 object which has a number of slots. These include the centre and rotation slots of a prcomp object which is returned after calling prcomp(t(x)), where x is the datamatrix supplied. These two slots are required for reconstructing the data later. The other slots are the factors supplied for the analysis (specified as expt and batch), the runtime parameters and some summary stats for Harman. Finally, there are the original and corrected PC scores.

str(olf.harman)

## List of 7
##  $ factors   :'data.frame':  28 obs. of  4 variables:
##   ..$ expt         : Factor w/ 7 levels "1","2","3","4",..: 1 2 3 4 5 6 7 1 2 3 ...
##   ..$ batch        : Factor w/ 4 levels "1","2","3","4": 1 1 1 1 1 1 1 2 2 2 ...
##   ..$ expt.numeric : int [1:28] 1 2 3 4 5 6 7 1 2 3 ...
##   ..$ batch.numeric: int [1:28] 1 1 1 1 1 1 1 2 2 2 ...
##  $ parameters:List of 4
##   ..$ limit     : num 0.95
##   ..$ numrepeats: int 100000
##   ..$ randseed  : num 0
##   ..$ forceRand : logi FALSE
##  $ stats     :'data.frame':  28 obs. of  3 variables:
##   ..$ dimension : Factor w/ 28 levels "PC1","PC10","PC11",..: 1 12 22 23 24 25 26 27 28 2 ...
##   ..$ confidence: num [1:28] 0.95 0.949 0.95 0.95 0.951 ...
##   ..$ correction: num [1:28] 0.25 0.33 0.5 0.9 0.44 0.85 0.74 1 1 1 ...
##  $ center    : Named num [1:33297] 4.13 3.95 3.19 8.92 3.19 ...
##   ..- attr(*, "names")= chr [1:33297] "p1" "p2" "p3" "p4" ...
##  $ rotation  : num [1:33297, 1:28] -0.020637 -0.011028 -0.006488 -0.000785 -0.00551 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:33297] "p1" "p2" "p3" "p4" ...
##   .. ..$ : chr [1:28] "PC1" "PC2" "PC3" "PC4" ...
##  $ original  : num [1:28, 1:28] -26.72 3.79 -19.76 25.66 -3.14 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:28] "c1" "c2" "c3" "c4" ...
##   .. ..$ : chr [1:28] "PC1" "PC2" "PC3" "PC4" ...
##  $ corrected : num [1:28, 1:28] -27.17 3.34 -20.21 25.21 -3.59 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : chr [1:28] "c1" "c2" "c3" "c4" ...
##   .. ..$ : chr [1:28] "PC1" "PC2" "PC3" "PC4" ...
##  - attr(*, "class")= chr "harmanresults"

2.1.2 Inspecting results

Harman objects can be inspected with methods such as pcaPlot and arrowPlot, as well as the generic functions plot and summary.

plot(olf.harman)

arrowPlot(olf.harman, main='Arrowplot of correction')

Using summary or inspecting the stats slot we can see Harman corrected the first 7 PCs and left the rest uncorrected.

summary(olf.harman)

## Total factor levels:
## 
##  expt batch 
##     7     4 
## 
## Experiment x Batch Design:
## 
##     batch
## expt 1 2 3 4
##    1 1 1 1 1
##    2 1 1 1 1
##    3 1 1 1 1
##    4 1 1 1 1
##    5 1 1 1 1
##    6 1 1 1 1
##    7 1 1 1 1
## 
## Simulation parameters:
## 100000 simulations (with seed of 0). ForceRand is FALSE.
## 
## Harman results with confidence limit of 0.95:
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 
## 0.25 0.33 0.50 0.90 0.44 0.85 0.74 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 
## PC16 PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 
## 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 
## 
## Top batch-effected PCs:
##  PC1  PC2  PC5  PC3  PC7 
## 0.25 0.33 0.44 0.50 0.74

As the confidence (limit) was 0.95, Harman will look for batch effects until this limit is met. It does this by reducing the batch means in an iterative way using a binary search algorithm until the chance that biological variance is being removed (the factor given in expt) is too high. Consider the values specified in the confidence column as the likelihood that separation of samples in this PC is due to batch alone and not due to the experimental variable. If this confidence value is less than limit, Harman will not make another iterative correction. For example, on PC8 the confidence is about 0.835 which is below the limit of 0.95, so Harman did not alter data on this PC. Let’s check this on a plot.

arrowPlot(olf.harman, 1, 8)

The horizontal arrows show us that, on PC1 the scores were shifted, but on PC8 they were not.

If both dimensions were not shifted, the arrowPlot will default to printing x instead of arrows (can’t have 0 length arrows!)

arrowPlot(olf.harman, 8, 9)

We can also easily colour the PCA plot by the experimental variable to compare:

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')

2.1.3 Reconstruct the corrected data

So, if corrections have been made and we’re happy with the value of limit, then we reconstruct corrected data from the Harman adjusted PC scores. To do this, a harmanresults object is handed to the function reconstructData(). This function produces a matrix of the same dimensions as the input matrix filled with corrected data.

olf.data.corrected <- reconstructData(olf.harman)

This new data can be over-written into something like an ExpressionSet object with a command like exprs(eset) <- data.corrected. An example of this is below.

To convince you that Harman is working as it should in reconstructing the data back from the PCA domain, let’s test this principle on the original data.

olf.data.remade <- reconstructData(olf.harman, this='original')
all.equal(as.matrix(olf.data), olf.data.remade)

## [1] TRUE

So, the data is indeed the same (to within the machine epsilon limit for floating point error)!

2.1.4 How does the new data differ from the old data?

Let’s try graphically plotting the differences using image(). First, the rows (assays) are ordered by the variation in the difference between the original and corrected data.

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'))

We can see that most Harman adjusted probes were from batch 3 (samples 15 to 21), while batch 1 is relatively unchanged. In batches 2 and 3, often the same probes were adjusted to a larger extent in batch 3. This suggests some probes on this array design are prone to a batch effect.

2.2.1 Evidence of batch effects

While there isn’t glaring batch effects in PC1 and PC2, they are more apparent when plotting PC1 and PC3.

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')

2.2.2 Running Harman

imr90.hm <- harman(as.matrix(imr90.data), expt=imr90.info$Treatment,
                   batch=imr90.info$Batch)

We can see the most batch correction was actually on the 3rd principal component and the data was also corrected on the 1st and 4th component.

Plotting PC1 v. PC2 and PC1 v. PC3…

plot(imr90.hm, pc_x=1, pc_y=2)

plot(imr90.hm, pc_x=1, pc_y=3)

While no batch effect was found on PC2, there are batch effects on PC1, 3 and 4. An arrowPlot() shows the extent of the correction:

arrowPlot(imr90.hm, pc_x=1, pc_y=3)

On PC1 the data has been shifted only a little, while on PC3 a large batch-effect signature seems to be present. The corrected data on the PC1 v. PC3 plot is still a little separated by batch. If we wish to aggressively stomp on the batch effect, (with increased risk of removing some experimental variance), we may specify something like limit=0.9.

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)

This time the samples across batches are less separated on the plot and there is also a batch effect found on PC2.

2.2.3 Reconstruct the corrected data

Again, as before we reconstuct the data. Let’s also use this data to do a comparison: a PCA plot of the original data and shiny new reconstructed data.

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')

2.3.1 Running Harman

Let’s fire up harman and try to remove this batch effect.

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)

If we inspect the stats on the Harman runs,

Loess-normalised data

Quantile-normalised data

The correction statistic is 1.0 for all dimensions, so Harman didn’t find any batch effect!

2.3.2 More aggressive settings

While our intuition tells us there is a batch effect, with the default settings (limit=0.95), harman fails to find one. This is due to the fact there are only two replicates in each batch! Let’s step up the limit now, to a confidence interval of 0.75. So, we are 75% sure only technical (batch) variation and not biological variation is being removed.

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)

This time a batch effect is found on PC1.

Loess-normalised data

Quantile-normalised data

2.3.3 Harman plots

So what corrections did Harman make with limit=0.75? Let’s take a look…

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)

2.3.4 Reconstruct the corrected data

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')

3.1.1 Batch structure

The experimental design is very poor in controlling for a batch effect. Ideally, the five factors should be distributed across the five run dates (batches). Instead, four out of the five experimental factors are distributed across two batches. The sTCC+CIS samples are all within batch 5. In this instance, batch is completely confounded with the experimental variable, so any variation in the data will be a sum of the biological variance and technical variance. Any inference about the difference of sTCC+CIS to other groups needs to be made very guardedly. It is also worth noting the spread of the mTCC samples across batches is suboptimal. While there are 11 samples in batch 1, there is only a single sample in batch 2. Finally, batches 3 and 4 have only one sample type.

Let’s explore the multivariate grouping of the data by generating some PCA plots, first with colouring by batch and then by experimental variable.

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')

The batch structures are highly unbalanced and so a batch effect is not immediately obvious to casual observation. However, a closer examination of samples that are represented well across two batches (Biopsy and Normal), shows a batch effect is certainly present. Certainly, with an older microarray technology like the Affymetrix HG-U133A, batch effects will be present.

3.1.2 Running Harman

Even though the batch structures are highly unbalanced, this dataset is quite large - so, we expect Harman to be well powered to find batch effects. Let’s take a look.

expt <- pData(bladderEset)$outcome
batch <- pData(bladderEset)$batch
bladder_harman <- harman(exprs(bladderEset), expt=expt, batch=batch)
sum(bladder_harman$stats$correction < 1)

## [1] 46

Harman found 46 of the 57 PCs had a batch effect.

summary(bladder_harman)

## Total factor levels:
## 
##  expt batch 
##     5     5 
## 
## Experiment x Batch Design:
## 
##           batch
## expt        1  2  3  4  5
##   Biopsy    0  0  0  5  4
##   mTCC     11  1  0  0  0
##   Normal    0  4  4  0  0
##   sTCC-CIS  0 13  0  0  3
##   sTCC+CIS  0  0  0  0 12
## 
## Simulation parameters:
## 100000 simulations (with seed of 0). ForceRand is FALSE.
## 
## Harman results with confidence limit of 0.95:
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 PC15 
## 0.00 0.56 0.07 0.11 0.27 0.56 0.38 0.10 0.22 0.25 0.25 0.42 0.61 0.28 0.36 
## PC16 PC17 PC18 PC19 PC20 PC21 PC22 PC23 PC24 PC25 PC26 PC27 PC28 PC29 PC30 
## 0.46 0.24 1.00 0.63 0.15 0.24 0.51 0.70 0.41 0.35 0.60 0.44 0.37 0.81 0.38 
## PC31 PC32 PC33 PC34 PC35 PC36 PC37 PC38 PC39 PC40 PC41 PC42 PC43 PC44 PC45 
## 0.33 0.32 0.47 0.53 0.25 0.46 0.60 1.00 0.34 0.88 1.00 1.00 1.00 0.81 1.00 
## PC46 PC47 PC48 PC49 PC50 PC51 PC52 PC53 PC54 PC55 PC56 PC57 
## 0.48 1.00 0.75 1.00 0.72 0.89 1.00 0.97 0.94 1.00 1.00 0.03 
## 
## Top batch-effected PCs:
##  PC1 PC57  PC3  PC8  PC4 
## 0.00 0.03 0.07 0.10 0.11

The summary() shows PC3 to be highly affected by batch. Now plotting the original and corrected data.

par(mfrow=c(1,1))
plot(bladder_harman)

This time for PCs 1 and 3.

plot(bladder_harman, 1, 3)

Now the PC scores changes displayed on an arrowPlot.

arrowPlot(bladder_harman, 1, 3)

3.1.3 Reconstruct the corrected data

First create a new object and then fill the exprs slot of the ExpressionSet object with Harman corrected data. Alternatively, the existing object can have the exprs slot overwritten.

CorrectedBladderEset <- bladderEset
exprs(CorrectedBladderEset) <- reconstructData(bladder_harman)
comment(bladderEset) <- 'Original'
comment(CorrectedBladderEset) <- 'Corrected'

4.1.1 Loading 450K data

So, let’s get underway and analyse the toy dataset supplied with minfi.

library(minfiData, quietly = TRUE)
baseDir <- system.file("extdata", package="minfiData")
targets <- read.metharray.sheet(baseDir)

## [1] "/home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/SampleSheet.csv"

Read in the files, normalise using preprocessSWAN() and preprocessIllumina() and filter out poorly called CpG sites.

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,]

This dataset has quite a few phenotype variables of interest. The samples are paired cancer-normal samples from three people, 1 male and 2 females.

kable(pData(mSetSw)[,c('Sample_Group', 'person', 'sex', 'status', 'Array',
                       'Slide')])

Doing some exploratory data analysis using MDS (considering it is methylation data), we can see that status is separated on PC1 and sex on PC2. This makes biological sense as we know cancer and gender both have large effects on DNA methylation. In cancer samples, there is typically global hypomethylation and focal hypermethylation at some CpG islands. While in the case of gender, female and male samples have very different X chromosome methylation due to imprinting, as well as some autosomal CpG site differences.

If the plot is recoloured by “Slide”“, there is some suggestion a batch effect is intermingled with these two experimental variables. The tricky question is how to use Harman to remove the influence of Slide on the data, without also removing variance due to the experimental variables of interest?

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)

The experiment is very small, with only 6 (paired) samples. The MDS plots suggests there are two main phenotypes of influence, status and sex. The batch variable in this case is “Slide”. The sex phenotype variable, which is highly influencing the data, is spread unevenly across the two slides (batches).

table(gender=pData(mSetSw)$sex, slide=pData(mSetSw)$Slide)

##       slide
## gender 5723646052 5723646053
##      F          1          3
##      M          2          0

From the table, we can see that only slide 5723646052 has samples of male origin. If we specify a simple cancer-based model with expt=pData(mSetSw)$status and batch=pData(mSetSw)$Slide, Harman will attribute the male-specific methylation signature of the two paired male samples on slide 5723646052 as batch effect and will try and eliminate it.

There are two strategies here: 1. Form a compound factor by joining the status and sex factors together. 2. Generate two Harman corrections, one for status and one for sex.

As this experiment is so small, the first strategy won’t be very effective. Only the factor level “F normal” is shared by both batches and for only one replicate. It’s very hard for Harman to get an idea of the batch distributions.

cancer_by_gender_factor <- as.factor(paste(pData(mSetSw)$sex,
                                           pData(mSetSw)$status)
                                     )
table(expt=cancer_by_gender_factor, batch=pData(mSetSw)$Slide)

##           batch
## expt       5723646052 5723646053
##   F cancer          0          2
##   F normal          1          1
##   M cancer          1          0
##   M normal          1          0

In the second strategy, both levels of expt are shared across both levels of batch. This is a far more ideal way to find batch effects. Of course though, there is confounding between sex and Slide! So, in these cases, do not consider a batch confounded factor in downstream differential analysis.

table(expt=pData(mSetSw)$status, batch=pData(mSetSw)$Slide)

##         batch
## expt     5723646052 5723646053
##   cancer          1          2
##   normal          2          1

4.1.2 Appropriate normalisation

If the Beta-values are shifted away from 0 or 1 by a very small amount, (say 1e-7), they will generate extreme M-values. Instead a moderate correction (such as 1e-4) seems the preferred option.

par(mfrow=c(2,2))
library(lumi, quietly = TRUE)

## 
## Attaching package: 'lumi'

## The following objects are masked from 'package:affy':
## 
##     MAplot, plotDensity

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)

It can be seen that M-values have far more central-tendency (more Gaussian-like) than Beta-values.

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)

A comparison of M-values produced by the GenomeStudio-like preprocessIllumina() function and by preprocessSWAN().

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)

From here on, we will work with SWAN normalised M-values.

4.1.3 Harman correction of M-values

For this example, status is declared as the experimental variable. When comparing conditions with very large DNA methylation differences, such as cancer and non-neoplastic samples, the batch effect will not be so easy to observe amongst all that biological variation. This scenario changes with a more subtle phenotype of interest. More arrays in each batch will allow Harman to more easily find batch effects. Considering this experiment is so small, an aggressive setting (limit=0.65) is needed to find a batch effect across a number of PCs (PCs 2, 5 and 1 in particular). A plot shows the correction made was relatively minor.

methHarman <- harman(getM(mSetSw), expt=pData(mSetSw)$status,
                     batch=pData(mSetSw)$Slide, limit=0.65)
summary(methHarman)

## Total factor levels:
## 
##  expt batch 
##     2     2 
## 
## Experiment x Batch Design:
## 
##         batch
## expt     5723646052 5723646053
##   cancer          1          2
##   normal          2          1
## 
## Simulation parameters:
## 100000 simulations (with seed of 0). ForceRand is FALSE.
## 
## Harman results with confidence limit of 0.65:
##  PC1  PC2  PC3  PC4  PC5  PC6 
## 0.96 0.61 1.00 1.00 0.92 0.81 
## 
## Top batch-effected PCs:
##  PC2  PC6  PC5  PC1 
## 0.61 0.81 0.92 0.96

plot(methHarman, 2, 5)

4.1.4 Reconstruct the corrected data

It is difficult to tease out the batch effect with such a small experimental group. Rather than further reduce the limit, let’s just convert the data back.

ms_hm <- reconstructData(methHarman)
betas_hm <- m2beta(ms_hm)

5.1.1 Loading the mass-spec data

# Call dependencies
library(msmsEDA)
library(Harman)

data(msms.dataset)
msms.dataset

## MSnSet (storageMode: lockedEnvironment)
## assayData: 697 features, 14 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: U2.2502.1 U2.2502.2 ... U6.0302.3 (14 total)
##   varLabels: treat batch
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
##   pubMedIds: http://www.ncbi.nlm.nih.gov/pubmed/22588121 
## Annotation:  
## - - - Processing information - - -
##  MSnbase version: 1.8.0

5.1.2 Preprocessing

The data matrix in an MSnSet object from package msmsEDA may have all zero rows (if it’s a subset of a larger object) and some samples may have NA values, which correspond to proteins not identified in that particular sample. Principle components analysis cannot be undertaken on matrices with such features. So first we can use the wrapper function pp.msms.data(), which removes all zero rows and replaces NA with 0.

# Preprocess to remove rows which are all 0 and replace NA values with 0.
msms_pp <- pp.msms.data(msms.dataset)

5.1.3 Batch structure

# Create experimental and batch variables
expt <- pData(msms_pp)$treat
batch <- pData(msms_pp)$batch
table(expt, batch)

##       batch
## expt   0302 2502
##   U200    3    4
##   U600    3    4

5.1.4 Running 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)

## Total factor levels:
## 
##  expt batch 
##     2     2 
## 
## Experiment x Batch Design:
## 
##       batch
## expt   0302 2502
##   U200    3    4
##   U600    3    4
## 
## Simulation parameters:
## 100000 simulations (with seed of 0). ForceRand is FALSE.
## 
## Harman results with confidence limit of 0.95:
##  PC1  PC2  PC3  PC4  PC5  PC6  PC7  PC8  PC9 PC10 PC11 PC12 PC13 PC14 
## 0.06 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 1.00 0.35 
## 
## Top batch-effected PCs:
##  PC1 PC14 
## 0.06 0.35

The Harman result is interesting. This MS data seems to have the batch effect contained within the 1st and 14th PCs only. A marked difference compared to transcriptome microarray data and the like.

A plot rather convincingly shows that Harman was able to remove the batch effect.

plot(hm)

5.1.5 Reconstruct the corrected data

As usual we now reconstruct corrected data, but we add an extra transformation step on the end. As the data was Log2 transformed we convert it back to the original format (and subtract 1 as this was added before during the transformation into Log2 space). The corrected and transformed back data is placed into a new ‘MSnSet’ instance.

# 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