BiocStyle 2.8.0

- 1 Introduction
- 2 Reading data into R
- 3 Subset-quantile within array normalization (SWAN)
- 4 Filter out poor quality probes
- 5 Extracting Beta and M-values
- 6 Testing for differential methylation using
- 7 Removing unwanted variation when testing for differential methylation
- 8 Testing for differential variability (DiffVar)
- 9 Gene ontology analysis
- 10 Session information
- References

The *missMethyl* package contains functions to analyse methylation
data from Illumina’s HumanMethylation450 and MethylationEPIC beadchip.
These arrays are a cost-effective alternative to whole genome bisulphite
sequencing, and as such are widely used to profile DNA methylation.
Specifically, *missMethyl* contains functions to perform SWAN
normalisation (Maksimovic, Gordon, and Oshlack 2012),
perform differential methylation analysis using **RUVm** (Maksimovic et al. 2015),
differential variability analysis (Phipson and Oshlack 2014) and gene set analysis
(Phipson, Maksimovic, and Oshlack 2016). As our lab’s research into specialised analyses of
these arrays continues we anticipate that the package will be
continuously updated with new functions.

Raw data files are in IDAT format, which can be read into R using the
*minfi* package (Aryee et al. 2014). Statistical analyses are usually
performed on M-values, and \(\beta\) values are used for visualisation,
both of which can be extracted from objects, which is a class of object created
by *minfi*. For detecting differentially variable CpGs we recommend that the
analysis is performed on M-values. All analyses described here are
performed at the CpG site level.

We will use the data in the *minfiData* package to demonstrate the
functions in *missMethyl*.
The example dataset has 6 samples across two slides. The sample
information is in the targets file. An essential column in the targets
file is the `Basename`

column which tells where the idat files to be
read in are located. The R commands to read in the data are taken from
the *minfi* User’s Guide. For additional details on how to read
the IDAT files into R, as well as information regarding quality control please
refer to the *minfi* User’s Guide.

```
library(missMethyl)
library(limma)
library(minfi)
```

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

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

`targets[,1:9]`

```
## Sample_Name Sample_Well Sample_Plate Sample_Group Pool_ID person age sex
## 1 GroupA_3 H5 <NA> GroupA <NA> id3 83 M
## 2 GroupA_2 D5 <NA> GroupA <NA> id2 58 F
## 3 GroupB_3 C6 <NA> GroupB <NA> id3 83 M
## 4 GroupB_1 F7 <NA> GroupB <NA> id1 75 F
## 5 GroupA_1 G7 <NA> GroupA <NA> id1 75 F
## 6 GroupB_2 H7 <NA> GroupB <NA> id2 58 F
## status
## 1 normal
## 2 normal
## 3 cancer
## 4 cancer
## 5 normal
## 6 cancer
```

`targets[,10:12]`

```
## Array Slide
## 1 R02C02 5723646052
## 2 R04C01 5723646052
## 3 R05C02 5723646052
## 4 R04C02 5723646053
## 5 R05C02 5723646053
## 6 R06C02 5723646053
## Basename
## 1 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646052/5723646052_R02C02
## 2 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646052/5723646052_R04C01
## 3 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646052/5723646052_R05C02
## 4 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646053/5723646053_R04C02
## 5 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646053/5723646053_R05C02
## 6 /home/biocbuild/bbs-3.7-bioc/R/library/minfiData/extdata/5723646053/5723646053_R06C02
```

`rgSet <- read.metharray.exp(targets = targets)`

The data is now an `RGChannelSet`

object and needs to be normalised and
converted to a `MethylSet`

object.

SWAN (subset-quantile within array normalization) is a within-array normalization method for Illumina 450k & EPIC BeadChips. Technical differencs have been demonstrated to exist between the Infinium I and Infinium II assays on a single Illumina HumanMethylation array (Bibikova et al. 2011, Dedeurwaerder, Defrance, and Calonne (2011)). Using the SWAN method substantially reduces the technical variability between the assay designs whilst maintaining important biological differences. The SWAN method makes the assumption that the number of CpGs within the 50bp probe sequence reflects the underlying biology of the region being interrogated. Hence, the overall distribution of intensities of probes with the same number of CpGs in the probe body should be the same regardless of assay type. The method then uses a subset quantile normalization approach to adjust the intensities of each array (Maksimovic, Gordon, and Oshlack 2012).

`SWAN`

can take a `MethylSet`

, `RGChannelSet`

or `MethyLumiSet`

as input. It
should be noted that, in order to create the normalization subset, `SWAN`

randomly selects Infinium I and II probes that have one, two and three
underlying CpGs; as such, we recommend using `set.seed`

before to ensure that
the normalized intensities will be
identical, if the normalization is repeated.

The technical differences between Infinium I and II assay designs can result in aberrant beta value distributions (Figure 1, panel “Raw”). Using SWAN corrects for the technical differences between the Infinium I and II assay designs and produces a smoother overall \(\beta\) value distribution (Figure 1, panel “SWAN”).

`mSet <- preprocessRaw(rgSet)`

`mSetSw <- SWAN(mSet,verbose=TRUE)`

```
## [SWAN] Preparing normalization subset
## 450k
## [SWAN] Normalizing methylated channel
## [SWAN] Normalizing array 1 of 6
## [SWAN] Normalizing array 2 of 6
## [SWAN] Normalizing array 3 of 6
## [SWAN] Normalizing array 4 of 6
## [SWAN] Normalizing array 5 of 6
## [SWAN] Normalizing array 6 of 6
## [SWAN] Normalizing unmethylated channel
## [SWAN] Normalizing array 1 of 6
## [SWAN] Normalizing array 2 of 6
## [SWAN] Normalizing array 3 of 6
## [SWAN] Normalizing array 4 of 6
## [SWAN] Normalizing array 5 of 6
## [SWAN] Normalizing array 6 of 6
```

```
par(mfrow=c(1,2), cex=1.25)
densityByProbeType(mSet[,1], main = "Raw")
densityByProbeType(mSetSw[,1], main = "SWAN")
```

Poor quality probes can be filtered out based on the detection p-value. For this example, to retain a CpG for further analysis, we require that the detection p-value is less than 0.01 in all samples.

```
detP <- detectionP(rgSet)
keep <- rowSums(detP < 0.01) == ncol(rgSet)
mSetSw <- mSetSw[keep,]
```

Now that the data has been `SWAN`

normalised we can extract \(\beta\) and
M-values from the object. We prefer to add an offset to the methylated
and unmethylated intensities when calculating M-values, hence we extract
the methylated and unmethylated channels separately and perform our own
calculation. For all subsequent analysis we use a random selection of
20000 CpGs to reduce computation time.

```
mset_reduced <- mSetSw[sample(1:nrow(mSetSw), 20000),]
meth <- getMeth(mset_reduced)
unmeth <- getUnmeth(mset_reduced)
Mval <- log2((meth + 100)/(unmeth + 100))
beta <- getBeta(mset_reduced)
dim(Mval)
```

`## [1] 20000 6`

```
par(mfrow=c(1,1))
plotMDS(Mval, labels=targets$Sample_Name, col=as.integer(factor(targets$status)))
legend("topleft",legend=c("Cancer","Normal"),pch=16,cex=1.2,col=1:2)
```

An MDS plot (Figure 2) is a good sanity check to make sure samples cluster together according to the main factor of interest, in this case, cancer and normal.

To test for differential methylation we use the *limma*
package (Smyth 2005), which employs an empirical Bayes framework based on
Guassian model theory. First we need to set up the design matrix.
There are a number of
ways to do this, the most straightforward is directly from the targets
file. There are a number of variables, with the `status`

column indicating
**cancer/normal** samples. From the `person`

column of the targets file, we
see that the **cancer/normal** samples are matched, with 3 individuals each
contributing both a **cancer** and **normal** sample. Since the
*limma* model framework can handle any experimental design which
can be summarised by
a design matrix, we can take into account the paired nature of the data
in the analysis. For more complicated experimental designs, please refer
to the *limma* User’s Guide.

```
group <- factor(targets$status,levels=c("normal","cancer"))
id <- factor(targets$person)
design <- model.matrix(~id + group)
design
```

```
## (Intercept) idid2 idid3 groupcancer
## 1 1 0 1 0
## 2 1 1 0 0
## 3 1 0 1 1
## 4 1 0 0 1
## 5 1 0 0 0
## 6 1 1 0 1
## attr(,"assign")
## [1] 0 1 1 2
## attr(,"contrasts")
## attr(,"contrasts")$id
## [1] "contr.treatment"
##
## attr(,"contrasts")$group
## [1] "contr.treatment"
```

Now we can test for differential methylation using the `lmFit`

and `eBayes`

functions from *limma*. As input data we use the matrix of
M-values.

```
fit.reduced <- lmFit(Mval,design)
fit.reduced <- eBayes(fit.reduced)
```

The numbers of hyper-methylated (1) and hypo-methylated (-1) can be
displayed using the `decideTests`

function in *limma* and the top
10 differentially methylated CpGs for *cancer* versus *normal* extracted using
`topTable`

.

`summary(decideTests(fit.reduced))`

```
## (Intercept) idid2 idid3 groupcancer
## Down 7056 0 121 568
## NotSig 3422 20000 19870 18996
## Up 9522 0 9 436
```

```
top<-topTable(fit.reduced,coef=4)
top
```

```
## logFC AveExpr t P.Value adj.P.Val B
## cg14034197 4.575822 -0.3196148 18.25222 7.004140e-06 0.03505176 4.171036
## cg13181745 4.519766 -0.9046829 17.41519 8.891317e-06 0.03505176 4.019872
## cg07952047 4.561251 -0.4189628 16.76257 1.079351e-05 0.03505176 3.892259
## cg15044384 4.251783 -0.3890717 16.12421 1.314170e-05 0.03505176 3.758349
## cg11946336 4.026042 -1.3425857 15.52598 1.591437e-05 0.03505176 3.623979
## cg25537217 4.800006 -0.2738630 15.09534 1.834732e-05 0.03505176 3.521512
## cg04398581 3.755296 0.1218274 15.08863 1.838857e-05 0.03505176 3.519877
## cg21159370 4.070400 0.4561686 14.96422 1.917434e-05 0.03505176 3.489309
## cg01952234 3.760289 -1.0274718 14.67533 2.115849e-05 0.03505176 3.416629
## cg09011883 4.097902 -1.3080656 14.59797 2.173067e-05 0.03505176 3.396755
```

Note that since we performed our analysis on M-values, the `logFC`

and
`AveExpr`

columns are computed on the M-value scale. For interpretability
and visualisation we can look at the \(\beta\) values. The beta values for
the top 4 differentially methylated CpGs shown in Figure 3.

```
cpgs <- rownames(top)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(beta[rownames(beta)==cpgs[i],]~design[,4],method="jitter",
group.names=c("Normal","Cancer"),pch=16,cex=1.5,col=c(4,2),ylab="Beta values",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(cpgs[i],cex.main=1.5)
}
```

Like other platforms, 450k array studies are subject to unwanted technical variation such as batch effects and other, often unknown, sources of variation. The adverse effects of unwanted variation have been extensively documented in gene expression array studies and have been shown to be able to both reduce power to detect true differences and to increase the number of false discoveries. As such, when it is apparent that data is significantly affected by unwanted variation, it is advisable to perform an adjustment to mitigate its effects.

*missMethyl* provides a *limma*-like interface to
functions from the CRAN package *ruv* that
enables the removal of unwanted variation when performing a differential
analysis (Maksimovic et al. 2015). All of the methods rely on negative control features
to accurately estimate the components of unwanted variation. Negative
control features are probes/genes/etc. that are known *a priori* to not
truly be associated with the biological factor of interest, but are
affected by unwanted variation. For example, in a microarray gene
expression study, these could be house-keeping genes or a set of
spike-in controls. Negative control features are extensively discussed
in Gagnon-Bartsch and Speed (2012) and Gagnon-Bartsch et al.
(2013). Once the unwanted factors are accurately estimated from
the data, they are adjusted for in the linear model that describes the
differential analysis.

If the negative control features are not known *a priori*, they can be
identified empirically. This can be achieved using a 2-stage approach,
**RUVm**, based on `RUV-inverse`

. Stage 1 involves performing a
differential methylation analysis using `RUV-inverse`

and the 613
Illumina negative controls (INCs) as negative control features. This
will produce a list of CpGs ranked by p-value according to their level
of association with the factor of interest. This list can then be used
to identify a set of empirical control probes (ECPs), which will capture
more of the unwanted variation than using the INCs alone. ECPs are
selected by designating a proportion of the CpGs least associated with
the factor of interest as negative control features; this can be done
based on either an FDR cut-off or by taking a fixed percentage of probes
from the bottom of the ranked list. Stage 2 involves performing a second
differential methylation analysis on the original data using
`RUV-inverse`

and the ECPs. For simplicity, we are ignoring the paired
nature of the **cancer** and **normal** samples in this example.

```
# get M-values for ALL probes
meth <- getMeth(mSet)
unmeth <- getUnmeth(mSet)
M <- log2((meth + 100)/(unmeth + 100))
grp <- factor(targets$status,levels=c("normal","cancer"))
des <- model.matrix(~grp)
des
```

```
## (Intercept) grpcancer
## 1 1 0
## 2 1 0
## 3 1 1
## 4 1 1
## 5 1 0
## 6 1 1
## attr(,"assign")
## [1] 0 1
## attr(,"contrasts")
## attr(,"contrasts")$grp
## [1] "contr.treatment"
```

```
INCs <- getINCs(rgSet)
head(INCs)
```

```
## 5723646052_R02C02 5723646052_R04C01 5723646052_R05C02
## 13792480 -0.3299654 -1.0955482 -0.5266103
## 69649505 -1.0354488 -1.4943396 -1.0067050
## 34772371 -1.1286422 -0.2995603 -0.8192636
## 28715352 -0.5553373 -0.7599489 -0.7186973
## 74737439 -1.1169178 -0.8656399 -0.6429681
## 33730459 -0.7714684 -0.5622424 -0.7724825
## 5723646053_R04C02 5723646053_R05C02 5723646053_R06C02
## 13792480 -0.6374299 -1.116598 -0.4332793
## 69649505 -0.8854881 -1.586679 -0.9217329
## 34772371 -0.6895514 -1.161155 -0.6186795
## 28715352 -1.7903619 -1.348105 -1.0067259
## 74737439 -0.8872082 -1.064986 -0.9841833
## 33730459 -1.5623138 -2.079184 -1.0445246
```

```
Mc <- rbind(M,INCs)
ctl <- rownames(Mc) %in% rownames(INCs)
table(ctl)
```

```
## ctl
## FALSE TRUE
## 485512 613
```

```
rfit1 <- RUVfit(data=Mc, design=des, coef=2, ctl=ctl) # Stage 1 analysis
rfit2 <- RUVadj(rfit1)
```

Now that we have performed an initial differential methylation analysis to rank the CpGs with respect to their association with the factor of interest, we can designate the CpGs that are least associated with the factor of interest based on FDR-adjusted p-value as ECPs.

```
top1 <- topRUV(rfit2, num=Inf)
head(top1)
```

```
## X1 X1 X1 p.BH p.ebayes p.ebayes.BH
## cg04743961 4.838190 26.74467 3.812882e-05 0.1401969 3.516091e-07 0.01017357
## cg07155336 5.887409 17.62103 1.608653e-04 0.1401969 3.583107e-07 0.01017357
## cg20925841 4.790211 26.69524 3.837354e-05 0.1401969 3.730375e-07 0.01017357
## cg03607359 4.394397 34.74068 1.542013e-05 0.1401969 4.721205e-07 0.01017357
## cg10566121 4.787914 21.80693 7.717708e-05 0.1401969 5.238865e-07 0.01017357
## cg07655636 4.571758 22.99708 6.424216e-05 0.1401969 6.080091e-07 0.01017357
```

```
ctl <- rownames(M) %in% rownames(top1[top1$p.ebayes.BH > 0.5,])
table(ctl)
```

```
## ctl
## FALSE TRUE
## 172540 312972
```

We can then use the ECPs to perform a second differential methylation
with `RUV-inverse`

, which is adjusted for the unwanted variation
estimated from the data.

```
# Perform RUV adjustment and fit
rfit1 <- RUVfit(data=M, design=des, coef=2, ctl=ctl) # Stage 2 analysis
rfit2 <- RUVadj(rfit1)
# Look at table of top results
topRUV(rfit2)
```

```
## X1 X1 X1 p.BH p.ebayes
## cg07155336 5.769286 15.345069 0.002005546 0.3431163 1.434834e-55
## cg06463958 5.733093 15.434797 0.001978272 0.3431163 6.749298e-55
## cg00024472 5.662959 15.946200 0.001832444 0.3431163 1.319390e-53
## cg02040433 5.651399 10.054445 0.005389436 0.3431163 2.146210e-53
## cg13355248 5.595396 9.963702 0.005504213 0.3431163 2.234891e-52
## cg02467990 5.592707 6.859614 0.013008521 0.3431163 2.499534e-52
## cg00817367 5.527501 13.070583 0.002921656 0.3431163 3.710480e-51
## cg11396157 5.487992 10.931263 0.004436178 0.3431163 1.873636e-50
## cg16306898 5.466780 5.573935 0.020790127 0.3431163 4.448085e-50
## cg03735888 5.396242 15.482605 0.001963955 0.3431163 7.700032e-49
## p.ebayes.BH
## cg07155336 6.966293e-50
## cg06463958 1.638433e-49
## cg00024472 2.135266e-48
## cg02040433 2.605027e-48
## cg13355248 2.022589e-47
## cg02467990 2.022589e-47
## cg00817367 2.573547e-46
## cg11396157 1.137091e-45
## cg16306898 2.399554e-45
## cg03735888 3.738458e-44
```

Note, at present does not support contrasts, so only one factor of interest can be interrogated at a time using a design matrix with an intercept term.

Rather than testing for differences in mean methylation, we may be interested in testing for differences between group variances. For example, it has been hypothesised that highly variable CpGs in cancer are important for tumour progression (Hansen et al. 2011). Hence we may be interested in CpG sites that are consistently methylated in the normal samples, but variably methylated in the cancer samples.

In general we recommend at least 10 samples in each group for accurate
variance estimation, however for the purpose of this vignette we perform
the analysis on 3 vs 3. In this example, we are interested in testing
for differential variability in the cancer versus normal group. Note
that when we specify the `coef`

parameter, which corresponds to the
columns of the design matrix to be used for testing differential
variability, we need to specify both the intercept and the fourth
column. The ID variable is a nuisance parameter and not used when
obtaining the absolute deviations, however it can be included in the
linear modelling step. For methylation data, the function will take
either a matrix of M-values, \(\beta\) values or a object as input. If
\(\beta\) values are supplied, a logit transformation is performed. Note
that as a default, `varFit`

uses the robust setting in the *limma*
framework, which requires the use of the *statmod* package.

`fitvar <- varFit(Mval, design = design, coef = c(1,4))`

The numbers of hyper-variable (1) and hypo-variable (-1) genes in **cancer**
vs **normal** can be obtained using `decideTests`

.

`summary(decideTests(fitvar))`

```
## (Intercept) idid2 idid3 groupcancer
## Down 0 2 4 1
## NotSig 19615 19996 19955 19996
## Up 385 2 41 3
```

```
topDV <- topVar(fitvar, coef=4)
topDV
```

```
## SampleVar LogVarRatio DiffLevene t P.Value
## cg05593887 13.742723 -0.4103645 -0.8161431 -5.255478 1.483965e-07
## cg23071808 5.551705 3.3795056 2.6509302 5.199805 2.004777e-07
## cg17942639 5.095744 4.0407816 2.6686110 4.745393 2.088204e-06
## cg17484671 4.590943 6.0095322 2.9026688 4.730653 2.245616e-06
## cg09770278 4.817211 2.2606740 2.0574707 4.261769 2.032784e-05
## cg21211367 4.499561 5.5939653 2.7003369 4.259559 2.052983e-05
## cg14479889 7.494275 3.2845218 2.3881498 4.242851 2.211936e-05
## cg18222083 5.066301 3.7112902 2.4550876 4.206153 2.603191e-05
## cg04186360 5.383690 2.2348168 2.0435622 4.200815 2.665310e-05
## cg23543123 4.617776 6.1347668 2.6889606 4.182078 2.894711e-05
## Adj.P.Value
## cg05593887 0.002004777
## cg23071808 0.002004777
## cg17942639 0.011228080
## cg17484671 0.011228080
## cg09770278 0.057894223
## cg21211367 0.057894223
## cg14479889 0.057894223
## cg18222083 0.057894223
## cg04186360 0.057894223
## cg23543123 0.057894223
```

An alternate parameterisation of the design matrix that does not include
an intercept term can also be used, and specific contrasts tested with
`contrasts.varFit`

.
Here we specify the design matrix such that the first two columns
correspond to the **normal** and **cancer** groups, respectively.

```
design2 <- model.matrix(~0+group+id)
fitvar.contr <- varFit(Mval, design=design2, coef=c(1,2))
contr <- makeContrasts(groupcancer-groupnormal,levels=colnames(design2))
fitvar.contr <- contrasts.varFit(fitvar.contr,contrasts=contr)
```

The results are identical to before.

`summary(decideTests(fitvar.contr))`

```
## groupcancer - groupnormal
## Down 1
## NotSig 19996
## Up 3
```

`topVar(fitvar.contr,coef=1)`

```
## SampleVar LogVarRatio DiffLevene t P.Value
## cg05593887 13.742723 -0.4103645 -0.8161431 -5.255478 1.483965e-07
## cg23071808 5.551705 3.3795056 2.6509302 5.199805 2.004777e-07
## cg17942639 5.095744 4.0407816 2.6686110 4.745393 2.088204e-06
## cg17484671 4.590943 6.0095322 2.9026688 4.730653 2.245616e-06
## cg09770278 4.817211 2.2606740 2.0574707 4.261769 2.032784e-05
## cg21211367 4.499561 5.5939653 2.7003369 4.259559 2.052983e-05
## cg14479889 7.494275 3.2845218 2.3881498 4.242851 2.211936e-05
## cg18222083 5.066301 3.7112902 2.4550876 4.206153 2.603191e-05
## cg04186360 5.383690 2.2348168 2.0435622 4.200815 2.665310e-05
## cg23543123 4.617776 6.1347668 2.6889606 4.182078 2.894711e-05
## Adj.P.Value
## cg05593887 0.002004777
## cg23071808 0.002004777
## cg17942639 0.011228080
## cg17484671 0.011228080
## cg09770278 0.057894223
## cg21211367 0.057894223
## cg14479889 0.057894223
## cg18222083 0.057894223
## cg04186360 0.057894223
## cg23543123 0.057894223
```

The \(\beta\) values for the top 4 differentially variable CpGs can be seen in Figure 4.

```
cpgsDV <- rownames(topDV)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(beta[rownames(beta)==cpgsDV[i],]~design[,4],method="jitter",
group.names=c("Normal","Cancer"),pch=16,cex=1.5,col=c(4,2),ylab="Beta values",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(cpgsDV[i],cex.main=1.5)
}
```

Testing for differential variability in expression data is
straightforward if the technology is gene expression microarrays. The
matrix of expression values can be supplied directly to the `varFit`

function.
For RNA-Seq data, the mean-variance relationship that occurs in count
data needs to be taken into account. In order to deal with this issue,
we apply a `voom`

transformation (Law et al. 2014) to obtain observation weights, which
are then used in the linear modelling step. For RNA-Seq data, the `varFit`

function will take a `DGElist`

object as input.

To demonstrate this, we use data from the *tweeDEseqCountData*
package. This data is part of the International HapMap project, consisting of
RNA-Seq profiles from 69 unrelated Nigerian individuals (Pickrell et al. 2010). The only
covariate is gender, so we can look at differentially variable expression between
males and females. We follow the code from the *limma* vignette to
read in and process the data before testing for differential variability.

First we load up the data and extract the relevant information.

```
library(tweeDEseqCountData)
data(pickrell1)
counts<-exprs(pickrell1.eset)
dim(counts)
```

`## [1] 38415 69`

```
gender <- pickrell1.eset$gender
table(gender)
```

```
## gender
## female male
## 40 29
```

```
rm(pickrell1.eset)
data(genderGenes)
data(annotEnsembl63)
annot <- annotEnsembl63[,c("Symbol","Chr")]
rm(annotEnsembl63)
```

We now have the counts, gender of each sample and annotation (gene
symbol and chromosome) for each Ensemble gene. We can form a `DGElist`

object
using the *edgeR* package.

```
library(edgeR)
y <- DGEList(counts=counts, genes=annot[rownames(counts),])
```

We filter out lowly expressed genes by keeping genes with at least 1 count per million reads in at least 20 samples, as well as genes that have defined annotation. Finally we perform scaling normalisation.

```
isexpr <- rowSums(cpm(y)>1) >= 20
hasannot <- rowSums(is.na(y$genes))==0
y <- y[isexpr & hasannot,,keep.lib.sizes=FALSE]
dim(y)
```

`## [1] 17310 69`

`y <- calcNormFactors(y)`

We set up the design matrix and test for differential variability. In
this case there are no nuisance parameters, so `coef`

does not need to
be explicitly specified.

```
design.hapmap <- model.matrix(~gender)
fitvar.hapmap <- varFit(y, design = design.hapmap)
```

`## Converting counts to log counts-per-million using voom.`

`fitvar.hapmap$genes <- y$genes`

We can display the results of the test:

`summary(decideTests(fitvar.hapmap))`

```
## (Intercept) gendermale
## Down 0 2
## NotSig 0 17308
## Up 17310 0
```

```
topDV.hapmap <- topVar(fitvar.hapmap,coef=ncol(design.hapmap))
topDV.hapmap
```

```
## Symbol Chr SampleVar LogVarRatio DiffLevene t
## ENSG00000213318 RP11-331F4.1 16 5.69839463 -2.562939 -0.9859943 -8.031243
## ENSG00000129824 RPS4Y1 Y 2.32497726 -2.087025 -0.4585620 -4.957005
## ENSG00000233864 TTTY15 Y 6.79004140 -2.245369 -0.6085233 -4.612934
## ENSG00000176171 BNIP3 10 0.41317384 1.199292 0.3632133 4.219404
## ENSG00000197358 BNIP3P1 14 0.39969125 1.149754 0.3353288 4.058147
## ENSG00000025039 RRAGD 6 0.91837213 1.091229 0.4926839 3.977022
## ENSG00000103671 TRIP4 15 0.07456448 -1.457139 -0.1520583 -3.911300
## ENSG00000171100 MTM1 X 0.44049558 -1.133295 -0.3334619 -3.896490
## ENSG00000149476 DAK 11 0.13289523 -1.470460 -0.1919880 -3.785893
## ENSG00000064886 CHI3L2 1 2.70234584 1.468059 0.8449434 3.782010
## P.Value Adj.P.Value
## ENSG00000213318 7.238039e-12 1.252905e-07
## ENSG00000129824 3.960855e-06 3.428120e-02
## ENSG00000233864 1.496237e-05 8.633290e-02
## ENSG00000176171 6.441668e-05 2.787632e-01
## ENSG00000197358 1.147886e-04 3.973982e-01
## ENSG00000025039 1.527695e-04 4.375736e-01
## ENSG00000103671 1.921104e-04 4.375736e-01
## ENSG00000171100 2.022293e-04 4.375736e-01
## ENSG00000149476 2.956364e-04 4.425050e-01
## ENSG00000064886 2.995692e-04 4.425050e-01
```

The log counts per million for the top 4 differentially variable genes can be seen in Figure 5.

```
genesDV <- rownames(topDV.hapmap)
par(mfrow=c(2,2))
for(i in 1:4){
stripchart(cpm(y,log=TRUE)[rownames(y)==genesDV[i],]~design.hapmap[,ncol(design.hapmap)],method="jitter",
group.names=c("Female","Male"),pch=16,cex=1.5,col=c(4,2),ylab="Log counts per million",
vertical=TRUE,cex.axis=1.5,cex.lab=1.5)
title(genesDV[i],cex.main=1.5)
}
```