- 1 Version Information
- 2 Introduction
- 3 Package Content
- 4 Workflow steps
- 4.1 A. Reading Input Data
- 4.2 B. Creating an ExpressionSet Object
- 4.3 C. Principal Variance Component Analysis of the raw data
- 4.4 D. Surrogate Variable Analysis
- 4.5 E. Computing the correlation between the surrogate variables and the covariates
- 4.6 F. Principal Variance Component Analysis of the raw data with the surrogate variables included as covariates
- 4.7 G. Supervised normalization of Microarrays
- 4.8 H. Principal Variance Component Analysis on the normalized data

- 5 Discussion
- References

**R version**: R version 4.0.0 RC (2020-04-19 r78255)

**Bioconductor version**: 3.12

**Package**: 1.15.0

Every Gene Expression study has an underlying question which the experimenter tries to address.Before proceeding towards any downstream analysis, the researcher has to decide how to normalize the gene expression data to minimize the impact of technical and other confounding factors on estimation of the environmental or biological factors of interest. One approach is to simply remove all of the major principal components of variation or specified technical effects, but this can also remove the biological signal of interest.Supervised normalization strategies encourage the user to define the major parameters influencing the expression variation, and then systematically remove them while shielding the biological factors of interest.

Here we outline a dynamic workflow that recognizes that there is not any ‘one algorithm fits all data sets’ analysis approach that works for every experiment.The workflow employs three successive procedures available in Bioconductor: Principal Variance Component Analysis (PVCA) to explore how technical and biological factors correlate with the major components of variance in the data set; Surrogate Variable Analysis (SVA) to identify major unwanted sources of variation; and Supervised Normalization of Microarrays (SNM) to efficiently remove these sources of variation while retaining the biological factor(s) of interest.
The data is based on a study contrasting peripheral blood gene expression in acute myocardial infarction and coronary artery disease patients, described in Kim et al. (2014). Only a subset of the data is used, and the labels have been changed for privacy protection, but the full data set is accessible from the GEO data repository at
http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE49925.

The workflow package contains five functions:

1. expSetobj – function to create an Expression Set Object which encapsulates the expression values, covariate information, and the experimental metadata

2. pvcAnaly - Principal Variance Component Analysis functionality

3. surVarAnaly – function which identifies hidden or new surrogate variables that are potential covariates hidden in the data

4. conTocat – function to convert continuous variables to categorical variables (discretize)

5. snmAnaly – function which implements the Supervised Normalization of Microarrays (SNM) normalization technique

The workflow runs on these sample two data files:

1. CAD_Expression.csv - file containing gene expression values (8000 probes for 100 samples)

2. CAD_Exptdsgn.csv - file containing the phenotype information (100 samples with 10 covariates)

The covariates are Study (a large batch effect since samples were run 2 years apart), Rin (RNA integrity number, a measure of RNA quality), BMI, Age, Gender, CAD (disease status as Acute MI or Fine at the time of sampling of patients in a CAD cohort), Rin3, BMI3, Age3 (discrete categorizations of the three continuous measures, used in the PVCA), and Array (a sample number)
The concept of the analysis is to explore how adjusting for Study, Rin, and other covariates if desired, influences the interpretation of the impact of CAD on gene expression.That is, CAD is the biological variable of interest.

- File containing the gene expression values (log transformed)

Format - Features x Samples (where the features could be probes in the case of microarrays, or gene or exon names in the case of RNASeq data), preferably a comma separated file (.csv) - File containing the covariates

Covariates are the different phenotypes that describe sample attributes, and can be biological, technical, or environmental.In the workflow we will also identify hidden covariates

Format - Samples x Covariates, preferably a comma separated file (.csv)

You should start by examining the sample data to get an idea of how the files should be structured.

`## Bioconductor version 3.12 (BiocManager 1.30.10), R 4.0.0 RC (2020-04-19 r78255)`

```
## Old packages: 'GO.db', 'PFAM.db', 'TxDb.Dmelanogaster.UCSC.dm6.ensGene',
## 'TxDb.Ggallus.UCSC.galGal5.refGene', 'TxDb.Rnorvegicus.UCSC.rn5.refGene',
## 'human.db0', 'matrixTests', 'org.Ag.eg.db', 'org.At.tair.db', 'org.Bt.eg.db',
## 'org.Ce.eg.db', 'org.Cf.eg.db', 'org.Dm.eg.db', 'org.Dr.eg.db',
## 'org.EcK12.eg.db', 'org.EcSakai.eg.db', 'org.Gg.eg.db', 'org.Hs.eg.db',
## 'org.Mm.eg.db', 'org.Mmu.eg.db', 'org.Pf.plasmo.db', 'org.Pt.eg.db',
## 'org.Rn.eg.db', 'org.Sc.sgd.db', 'org.Ss.eg.db', 'org.Xl.eg.db', 'recipes',
## 'rematch2'
```

```
## Warning: replacing previous import 'Biobase::anyMissing' by
## 'matrixStats::anyMissing' when loading 'ExpressionNormalizationWorkflow'
```

```
## Warning: replacing previous import 'Biobase::rowMedians' by
## 'matrixStats::rowMedians' when loading 'ExpressionNormalizationWorkflow'
```

`## Loading required package: BiocGenerics`

`## Loading required package: parallel`

```
##
## Attaching package: 'BiocGenerics'
```

```
## The following objects are masked from 'package:parallel':
##
## clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
## clusterExport, clusterMap, parApply, parCapply, parLapply,
## parLapplyLB, parRapply, parSapply, parSapplyLB
```

```
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
```

```
## The following objects are masked from 'package:base':
##
## Filter, Find, Map, Position, Reduce, anyDuplicated, append,
## as.data.frame, basename, cbind, colnames, dirname, do.call,
## duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
## lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
## pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
## tapply, union, unique, unsplit, which, which.max, which.min
```

```
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
```

```
## Set your working directory (or use Control + Shift + H to choose it)
setwd("working directory")
## read in the file containing the gene expression values
expData_file_name <- system.file("folder containing the data", "CAD_Expression.csv", package="ExpressionNormalizationWorkflow")
exprs <- read.table(expData_file_name, header=TRUE, sep=",", row.names=1, as.is=TRUE)
## read in the file containing the covariates
expDesign_file_name <- system.file("folder containing the data", "CAD_ExptDsgn.csv", package="ExpressionNormalizationWorkflow")
covrts <- read.table(expDesign_file_name, header=TRUE, sep=",", row.names=1, as.is=TRUE)
```

```
## A001 A002 A003 A004 A005
## ILMN_2389211 14.091 13.892 13.747 14.056 14.243
## ILMN_1667796 13.872 13.903 14.273 14.276 14.125
## ILMN_1683271 13.419 13.665 13.180 13.487 14.252
## ILMN_2100437 14.237 13.638 14.160 14.004 13.945
## ILMN_2242491 13.692 13.830 13.667 13.756 14.015
```

```
## Study Rin3 Gender Age3 BMI3 CAD Array Rin Age BMI
## A001 A HIGH FEM ADULT OVER ACUTE 1 9.2 56 31.0
## A002 A HIGH FEM ADULT OVER ACUTE 2 9.1 54 26.9
## A003 A HIGH MAL ADULT NORM ACUTE 3 9.1 51 20.7
## A004 A HIGH MAL ADULT NORM ACUTE 4 9.3 52 24.4
## A005 A HIGH MAL MATURE OVER ACUTE 5 9.3 62 26.1
```

An ExpressionSet Class (Falcon, Morgan and Gentleman, 2006) is a Biobase data structure that is used to conveniently store experimental information and associated meta data, all in one place. This command creates an object of the the ExpressionSet Class, which stores the gene expression values and the covariate data.The function GEOquery (Davis, 2016) can be used to download ExpressionSet objects from directly from GEO (https://www.ncbi.nlm.nih.gov/gds)

`inpData <- expSetobj(exprs, covrts)`

PVCA (Bushel, 2013) estimates the proportion of the variance of the principal components of the gene expression data that can be attributed to each of the given covariates.The remaining fraction is assigned as “residual” variance.It efficiently combines principal component analysis (PCA) to reduce the feature space and variance component analysis (VCA) which fits a mixed linear model using factors of interest as random effects to estimate and partition the total variability. The variance explained is computed as a weighted average of the contributions of each factor to each PC, and you have the option of specifying how many principal components to include, by setting a threshold amount of the variance that needs to be explained by the identified PCs. Here PVCA is used to estimate the proportion of variance due to CAD as well as due to the covariates like Gender, BMI, Rin, and Study.Since BMI and Rin are continuous variables, we use a categorization of each into Obese, Overweight, Normal weight and High, Moderate, and Low quality RNA respectively, through the BMI3 and Rin3 columns.

```
## Set the covariates whose effect size on the data needs to be calculated
cvrts_eff_var <- c("CAD", "BMI3", "Rin3", "Gender", "Study")
## Set a PVCA Threshold Value between 0 & 1
## PVCA Threshold Value is the percentile value of the minimum amount of the variabilities that the selected principal components need to explain, here requiring 75% of the expression variance to be captured by the PCs
pct_thrsh <- 0.75
## Perform the PVCA analysis
pvcAnaly(inpData, pct_thrsh, cvrts_eff_var)
```

```
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
```

```
## $dat
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.02675754 0.03947988 0.2685198 0.02055874 0.03676186 0.6079222
##
## $label
## [1] "Rin3" "BMI3" "Study" "Gender" "CAD" "resid"
```

The analysis (also generated as a histogram) shows that the Study batch effect is a major technical artifact that explains ~27% of the expression variance.In comparison, CAD and BMI each only explain around 4% of the variance, RNA quality less than 3%, and Gender is a minor component.60% of the variance is residual.

Surrogate variables (Leek and Storey, 2007) are covariates constructed directly from high-dimensional data (gene expression or RNA-Seq data) that can be used in subsequent analyses to adjust for unknown/unmodeled covariates or latent sources of noise. The user provides one or more biological variable(s) that they are interested in estimating, and then estimates the surrogate variables (sv) given the presence of the biological signal. These are then appended as new covariates to the existing list of covariates, and we also add a step to convert them to categorical variables so as to estimate their contributions to the total variance.For a thorough introduction to SVA, refer to the paper by Leek and Storey (2007) .SVA is available from Bioconductor at https://bioconductor.org/packages/release/bioc/html/sva.html, and a version specifically for RNASeq data (svaseq) is also available from Github at https://github.com/jtleek/svaseq.

```
## Choose a biological variable that is to be used to calculate the surrogate variables
biol_var_sva <- "CAD"
## Perform SVA
sur_var_obj <- surVarAnaly(inpData, biol_var_sva)
```

```
## Number of significant surrogate variables is: 4
## Iteration (out of 5 ):1 2 3 4 5
```

```
## The newly generated surrogate variables sv1 through sv4 are appended to the ExpressionSet object
inpData_sv <- sur_var_obj$expSetobject
```

This step helps the researcher explore whether each newly identified surrogate variable is entirely independent of all the existing covariates or if they are correlated with either a technical factor or one of the biological covariates. If they are associated with a biological factor (especially the focus of the analysis, in this case CAD), then it would not be advisable to remove them during the normalization step, though they might be adjusted for.Otherwise, it will usually be advisable to remove them, or to remove the technical factor they capture.To see the correlations, we run a series of generalized linear models where the identified surrogate variables are modeled as a function of the existing covariates (Study, CAD, Rin, BMI; you could add Gender and Age if you wish – but they are not significant).

```
## Fit a generalized linear model for sv1
glm.sv1 <- glm(pData(inpData_sv)[,"sv1"]~pData(inpData_sv)[,"BMI"]+pData(inpData_sv)[,"Rin"]+pData(inpData_sv)[,"CAD"]
+pData(inpData_sv)[,"Study"])
summary(glm.sv1)
```

```
##
## Call:
## glm(formula = pData(inpData_sv)[, "sv1"] ~ pData(inpData_sv)[,
## "BMI"] + pData(inpData_sv)[, "Rin"] + pData(inpData_sv)[,
## "CAD"] + pData(inpData_sv)[, "Study"])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.127062 -0.040490 -0.009979 0.030200 0.139604
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.1832540 0.0711053 2.577 0.0115 *
## pData(inpData_sv)[, "BMI"] 0.0009422 0.0009208 1.023 0.3088
## pData(inpData_sv)[, "Rin"] -0.0159186 0.0074152 -2.147 0.0344 *
## pData(inpData_sv)[, "CAD"]FINE 0.0095449 0.0119928 0.796 0.4281
## pData(inpData_sv)[, "Study"]B -0.1699332 0.0129778 -13.094 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.00357263)
##
## Null deviance: 1.0000 on 99 degrees of freedom
## Residual deviance: 0.3394 on 95 degrees of freedom
## AIC: -272.79
##
## Number of Fisher Scoring iterations: 2
```

The analysis shows that the Study batch effect is very strongly associated with sv1.You can also see this by plotting the two measures from the pData(inpData_sv) file. Subsequent analyses of sv2 through sv4 with similar code shows that BMI is weakly associated with sv2, but strongly with sv3 (as is Study), while Rin is strongly associate with sv4.CAD is not associated with any of the surrogate variables, but Age contributes slightly to sv2 and sv3.

```
## Fit a generalized linear model for sv2
glm.sv2 <- glm(pData(inpData_sv)[,"sv2"]~pData(inpData_sv)[,"BMI"]+pData(inpData_sv)[,"Rin"]+pData(inpData_sv)[,"CAD"]
+pData(inpData_sv)[,"Study"])
summary(glm.sv2)
```

```
##
## Call:
## glm(formula = pData(inpData_sv)[, "sv2"] ~ pData(inpData_sv)[,
## "BMI"] + pData(inpData_sv)[, "Rin"] + pData(inpData_sv)[,
## "CAD"] + pData(inpData_sv)[, "Study"])
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -0.21940 -0.04787 -0.00900 0.04228 0.33138
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.261993 0.115073 -2.277 0.0250 *
## pData(inpData_sv)[, "BMI"] 0.003209 0.001490 2.154 0.0338 *
## pData(inpData_sv)[, "Rin"] 0.019394 0.012000 1.616 0.1094
## pData(inpData_sv)[, "CAD"]FINE 0.030581 0.019409 1.576 0.1184
## pData(inpData_sv)[, "Study"]B -0.016918 0.021003 -0.806 0.4225
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.009356872)
##
## Null deviance: 1.0000 on 99 degrees of freedom
## Residual deviance: 0.8889 on 95 degrees of freedom
## AIC: -176.51
##
## Number of Fisher Scoring iterations: 2
```

`## Output should be similar to that shown above for sv1`

The following PVCA step is performed to estimate the contributions of the newly identified surrogate variables to the overall expression variance, with or without the Study covariate in the model.

```
## First discretize the continuous surrogate variables
var_names <- c("sv1", "sv2", "sv3", "sv4")
pData(inpData_sv)<-conTocat(pData(inpData_sv), var_names)
## View them appended to the covariate matrix as additional covariate columns
#View(pData(inpData_sv))
## Include the surrogate variables as covariates in addition to BMI3, Rin3, CAD and Study (be sure to use categorical measures of BMI and Rin rather than the continuous measure)
cvrts_eff_var <- c("BMI3", "Rin3", "CAD", "Study", "sv1_cat", "sv2_cat", "sv3_cat", "sv4_cat")
## Again set the PVCA Threshold to explain 75% of the expression variation
pct_thrsh <- 0.75
## Perform PVCA
pvcAnaly(inpData_sv, pct_thrsh, cvrts_eff_var)
```

```
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
```

```
## $dat
## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 0.08128376 0.05778339 0.1442518 0.1313381 0.00861214 0.00638395 0.19577
## [,8] [,9]
## [1,] 0.01316618 0.3614106
##
## $label
## [1] "sv4_cat" "sv3_cat" "sv2_cat" "sv1_cat" "Rin3" "BMI3" "Study"
## [8] "CAD" "resid"
```

The analysis shows that each of the surrogate variables capture a large amount of the variance, while the Study batch effect remains important.The residual variance has reduced by almost half to 36%.For comparison, remove Study to see how much of it the surrogate variables capture

```
cvrts_eff_var <- c("BMI3", "Rin3", "CAD","sv1_cat", "sv2_cat", "sv3_cat", "sv4_cat")
## Again set the PVCA Threshold to explain 75% of the expression variation
pct_thrsh <- 0.75
## Perform PVCA
pvcAnaly(inpData_sv, pct_thrsh, cvrts_eff_var)
```

```
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
## boundary (singular) fit: see ?isSingular
```

```
## $dat
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 0.08063775 0.07554799 0.1537535 0.276012 0.007577982 0.007033493
## [,7] [,8]
## [1,] 0.01256891 0.3868684
##
## $label
## [1] "sv4_cat" "sv3_cat" "sv2_cat" "sv1_cat" "Rin3" "BMI3" "CAD"
## [8] "resid"
```

This implies that together the 4 SVs capture as much variance as Study alone.Both analyses have slightly reduced the CAD contribution relative to the analysis without surrogate variables.

SNM (Mecham, Nelson and Storey, 2010) is a study specific, customizable normalization approach that adjusts for specified biological and technical variables as far as possible without biasing the estimate of the biological source of interest. Much like SVA, the user chooses the biological variable(s), and then fits desired covariates, while also allowing for “intensity dependent” effects (here, just array/sample differences) to be adjusted for. The SNM software at http://www.bioconductor.org/packages/release/bioc/html/snm.html actually allows you to either remove the adjustment variables using the rm=TRUE option, or simply to fit them.You can also run both sequentially with different covariates (for example, removing a batch effect but adjusting for gender). SNM fits both categorical and continuous variables.In the workflow below we remove the Study or the surrogate variable effects, running 5 iterations.

```
## Choose the biological variableof interest
bv <- c("CAD")
## Choose your adjustment variable of interest, starting with just 'Study'
av <- c("Study")
## The intensity-dependent adjustment variables adjust for array effects
iv <- c("Array")
## Run SNM
sv_snmObj <- snmAnaly(exprs, pData(inpData_sv), bv, av, iv)
```

```
##
Iteration: 1
```

```
##
Iteration: 2
```