nSolver Normalization

To normalize using standard nSolver steps, set normalization to “nSolver”. The three normalization steps under this method are:

1. Positive control normalization: The geometric mean of positive control probes is used to calculate a scaling factor for each sample.
2. Background assessment: Negative control probes are used to assess whether endogenous genes are truly expressed above the level of backgrounds. Two options are provided for this
   • bgType = "threshold" - The background threshold will be defined independently for each sample, using threshold = bgThreshold * sd(NegativeControls) + mean(NegativeControls), where bgThreshold is set by the user. For each gene, the proportion of samples with expression above the threshold is calculated, and only genes with a proportion greater than bgProportion (set by user) will be retained for analysis. For example, if bgThreshold is 2 and bgProportion is 0.5 (defaults), a gene will be retained for analysis if its expression is 2 standard deviations above the mean of negative control probes in at least half of the samples.
   • bgType = "t.test" - Expression of each gene is compared with all of the negative control probes (across samples) using a one-sided, two-sample t test. Genes will be retained for analysis if the endogenous gene expression is greater, with p < bgPVal.
   • bgSubtract - By setting this to TRUE, the calculated background level (mean + bgThreshold * sd) will be subtracted from the expression values for each gene in each sample. After subtraction, genes with negative values will be set to zero.
3. Housekeeping normalization: The geometric mean of housekeeping gene expression is used to calculate a scaling factor for each sample. Housekeeping genes (symbols or accessions) can be input using the housekeeping option. If not provided, genes marked as “Housekeeping” in the RCC files will be used. Alternatively, this can be skipped using skip.housekeeping = TRUE.

# Set housekeeping genes manually (optional)
hk.genes <- c("TUBB", "TBP", "POLR2A", "GUSB", "SDHA")

dat <- processNanostringData(nsFiles = example_data,
                             sampleTab = sample_info, 
                             idCol = "RCC_Name",
                             groupCol = "Sample_Diagnosis",
                             normalization = "nSolver",
                             bgType = "t.test", bgPVal = 0.01,
                             housekeeping = hk.genes, skip.housekeeping = FALSE)

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Averaging technical replicates.....
## Calculating positive scale factors......
## Checking endogenous genes against background threshold......
## Conducting housekeeping normalization......

# Automatically detect housekeeping genes

dat <- processNanostringData(nsFiles = example_data,
                             sampleTab = sample_info, 
                             idCol = "RCC_Name",
                             groupCol = "Sample_Diagnosis",
                             normalization = "nSolver",
                             bgType = "t.test", bgPVal = 0.01)

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Averaging technical replicates.....
## Calculating positive scale factors......
## Checking endogenous genes against background threshold......
## Conducting housekeeping normalization......

RUV-III Normalization

The RUV-III method normalizes data based on technical replicates. Replicate IDs must be included in the sample information table (‘sampleTab’ in processNanostringData), and the column containing these IDs is specified using the ‘replicateCol’ option. Additionally, the number of dimensions of unwanted data can be specified using the ‘n.unwanted’ option. If NULL (the default), the maximum possible number of factors will be identified and used. More information can be found in (Gagnon-Bartsch 2019).

Note that RUVIII returns log-transformed counts, in contrast to nSolver and RUVg normalization. This difference is accounted for in subsequent NanoTube analysis steps, but it should be considered by the user when using the normalized data for other purposes.

# This sample data table contains fake replicate identifiers (for use in example).
sample_info_reps <- system.file("extdata", 
                           "GSE117751_sample_data_replicates.csv", 
                           package = "NanoTube")

datRUVIII <- processNanostringData(nsFiles = example_data,
                             sampleTab = sample_info_reps, 
                             idCol = "RCC_Name",
                             replicateCol = "Replicate_ID",
                             groupCol = "Sample_Diagnosis",
                             normalization = "RUVIII",
                             bgType = "t.test", bgPVal = 0.01,
                             n.unwanted = 1)

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Conducting RUVIII normalization......

RUV-III can also be used without technical replicates, by generating PRPS (pseudo-replicates of pseudo-samples). This is described in (Molania et al. 2022), and it was demonstrated to provide superior normalization vs. other RNA-Seq normalization methods for TCGA RNA-Seq data sets. Software and vignettes for accomplishing this normalization can be found at https://github.com/RMolania/TCGA_PanCancer_UnwantedVariation.

RUVg Normalization

The RUVg method normalizes data based on housekeeping genes. These can be manually specified using the ‘housekeeping’ option in processNanostringData, or they can be automatically interpreted based on the CodeClass column of the RCC files. By default, RUVg will remove one dimension of unwanted data, but this can be adjusted using the ‘n.unwanted’ option. Additionally, the ‘RUVg.drop’ option can be used to drop some number of singular values. This defaults to zero, but it can be adjusted if the first singular value(s) identify the effect of interest. See (Bhattacharya et al. 2020) for more details.

hk.genes <- c("TUBB", "TBP", "POLR2A", "GUSB", "SDHA")

datRUVg <- processNanostringData(nsFiles = example_data,
                             sampleTab = sample_info, 
                             idCol = "RCC_Name",
                             groupCol = "Sample_Diagnosis",
                             normalization = "RUVg",
                             bgType = "t.test", bgPVal = 0.01,
                             n.unwanted = 1, RUVg.drop = 0,
                             housekeeping = hk.genes)

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Averaging technical replicates.....
## Conducting RUVg normalization......

## Warning: replacing previous import 'utils::findMatches' by
## 'S4Vectors::findMatches' when loading 'AnnotationDbi'

# Housekeeping genes are automatically identified
datRUVg <- processNanostringData(nsFiles = example_data,
                             sampleTab = sample_info, 
                             idCol = "RCC_Name",
                             groupCol = "Sample_Diagnosis",
                             normalization = "RUVg",
                             bgType = "t.test", bgPVal = 0.01,
                             n.unwanted = 1, RUVg.drop = 0)

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Averaging technical replicates.....
## Conducting RUVg normalization......

RLE Plots

Relative Log Expression (RLE) plots have been demonstrated as effective visualizations to assess whether NanoString data normalization has been successful (Gandolfo and Speed 2018). This method centers the log-expression of each gene to the median log-expression of that gene. The RUV implementation ruv_rle can be used to generate this plot, regardless of what normalization method was used (Gagnon-Bartsch 2019).

Note that for this function, you need to log-transform and transpose the expression data matrix. You’ll probably want to set the y-axis limits as well, using ‘ylim’.

# Here's how to do it with the NanoString ExpressionSet. If your data set is in
# list form, you would use dat$exprs instead of exprs(dat).

ruv::ruv_rle(t(log2(exprs(dat))), ylim = c(-2, 2))

Principal components analysis

PCA can be conducted after processing and normalization. We provide a standard plot using ggplot2 or an interactive plot using plotly (use the interactive.plot option).

dat <- processNanostringData(example_data, 
                             sampleTab = sample_info, 
                             idCol = "RCC_Name",
                             groupCol = "Sample_Diagnosis",
                             normalization = "nSolver", 
                             bgType = "t.test", 
                             bgPVal = 0.01)

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Averaging technical replicates.....
## Calculating positive scale factors......
## Checking endogenous genes against background threshold......
## Conducting housekeeping normalization......

library(plotly) 
nanostringPCA(dat, interactive.plot = TRUE)$plt

We default to plotting the first two principal components, but you can also choose others.

nanostringPCA(dat, pc1=3, pc2=4, interactive.plot = FALSE)$plt

Using Limma

After normalization, the data are likely to resemble a normal distribution, particularly with reasonable filtering of genes with expression below background levels. It is a good idea to check this before using limma (Ritchie et al. 2015). See limma vignette for full details and additional analysis options.

Assuming this information was provided in the processNanostringData step, your ExpressionSet will already contain the sample groups.

dat <- processNanostringData(example_data, 
                             sampleTab = sample_info, 
                             idCol = "RCC_Name",
                             groupCol = "Sample_Diagnosis",
                             normalization = "nSolver", 
                             bgType = "t.test", 
                             bgPVal = 0.01)

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Averaging technical replicates.....
## Calculating positive scale factors......
## Checking endogenous genes against background threshold......
## Conducting housekeeping normalization......

table(dat$groups)

## 
## Autoimmune retinopathy                   None   Retinitis pigmentosa 
##                     14                     14                     14

runLimmaAnalysis allows you to conduct group-vs-group comparisons by defining the base.group. In this example, setting base.group to “None” will cause Limma to build a linear model fitting the expression data, where the intercept is the average log2 expression of “None”, and the two other coefficients will correspond to the log2(FC) of Autoimmune retinopathy vs. None and Retinitis pigmentosa vs. None. This function will return standard Limma analysis results.

limmaResults <- runLimmaAnalysis(dat, base.group = "None")

head(signif(limmaResults$coefficients, digits = 2))

##          Intercept Autoimmune.retinopathy Retinitis.pigmentosa
## HLA-DQB1       7.4                  -0.91              -1.6000
## KIT            3.8                   0.14              -0.0098
## LAG3           6.3                  -0.37              -0.0640
## SOCS3          7.2                   0.16               0.0900
## TCF7          10.0                  -0.11               0.0520
## IKBKB          7.8                   0.42               0.1100

You can also directly define a design matrix, instead. For example, if this data were collected in two batches, you could include a batch term in the analysis.

# Generate fake batch labels
batch <- rep(c(0, 1), times = ncol(dat) / 2)

# Reorder groups ("None" first)
group <- factor(dat$groups, 
                levels = c("None", "Autoimmune retinopathy", 
                           "Retinitis pigmentosa"))

# Design matrix including sample group and batch
design <- model.matrix(~group + batch)

# Analyze data
limmaResults2 <- runLimmaAnalysis(dat, design = design)

# We can see there is no batch effect in this data, due to the somewhat uniform
# distribution of p-values. This is good, since the batches are fake.
hist(limmaResults2$p.value[,"batch"], main = "p-values of Batch terms")

Alternatively, a design matrix can be provided as ‘sampleData’ in the processNanostringData step. This one contains columns for the Intercept, the two diagnoses, and patient Age.

sample_info_design <- system.file("extdata", 
                           "GSE117751_design_matrix.csv", 
                           package = "NanoTube")

datDes <- processNanostringData(example_data, 
                             sampleTab = sample_info_design, 
                             idCol = "RCC_Name",
                             normalization = "nSolver", 
                             bgType = "t.test", 
                             bgPVal = 0.01)

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Averaging technical replicates.....
## Calculating positive scale factors......
## Checking endogenous genes against background threshold......
## Conducting housekeeping normalization......

head(pData(datDes))

##                                                            RCC_Name Intercept
## GSM3308226_20160715_5609WV-C.RCC   GSM3308226_20160715_5609WV-C.RCC         1
## GSM3308227_20170504_6090AA-A2.RCC GSM3308227_20170504_6090AA-A2.RCC         1
## GSM3308228_20170504_6121PF-A2.RCC GSM3308228_20170504_6121PF-A2.RCC         1
## GSM3308229_20170504_5993KP-A2.RCC GSM3308229_20170504_5993KP-A2.RCC         1
## GSM3308230_20170504_6536VE-A2.RCC GSM3308230_20170504_6536VE-A2.RCC         1
## GSM3308231_20170504_6568NC-A2.RCC GSM3308231_20170504_6568NC-A2.RCC         1
##                                   Autoimmune_retinopathy Retinitis_Pigmentosa
## GSM3308226_20160715_5609WV-C.RCC                       1                    0
## GSM3308227_20170504_6090AA-A2.RCC                      1                    0
## GSM3308228_20170504_6121PF-A2.RCC                      1                    0
## GSM3308229_20170504_5993KP-A2.RCC                      1                    0
## GSM3308230_20170504_6536VE-A2.RCC                      1                    0
## GSM3308231_20170504_6568NC-A2.RCC                      1                    0
##                                   Age normalization
## GSM3308226_20160715_5609WV-C.RCC   61       nSolver
## GSM3308227_20170504_6090AA-A2.RCC  71       nSolver
## GSM3308228_20170504_6121PF-A2.RCC  71       nSolver
## GSM3308229_20170504_5993KP-A2.RCC  66       nSolver
## GSM3308230_20170504_6536VE-A2.RCC  71       nSolver
## GSM3308231_20170504_6568NC-A2.RCC  59       nSolver

We can then use this pData object as our design matrix in limma analysis. Note that processNanostringData adds an additional column at the beginning and the end, so these should be removed for limma analysis. Based on this analysis, many genes appear to be correlated with patient age.

# Analyze data
limmaResults3 <- runLimmaAnalysis(datDes, 
                                  design = pData(datDes)[,2:(ncol(pData(datDes))-1)])


hist(limmaResults3$p.value[,"Age"], main = "p-values for Age term")

Results of Limma analyses can be visualized using simple volcano plots.

deVolcano(limmaResults, y.var = "p.value")

## 
## 'plotContrast' not provided, setting it to Autoimmune.retinopathy

You can add additional ggplot layers as well.

deVolcano(limmaResults, plotContrast = "Autoimmune.retinopathy", 
          y.var = "p.value") +
  geom_hline(yintercept = 2, linetype = "dashed", colour = "darkred") +
  geom_vline(xintercept = 0.5, linetype = "dashed", colour = "darkred") +
  geom_vline(xintercept = -0.5, linetype = "dashed", colour = "darkred")

They can also be converted to a simple table or exported to a text file.

limmaStats <- makeDiffExprFile(limmaResults, filename = NULL, returns = "stats")
limmaStats <- as.data.frame(limmaStats)

# Order by lowest to highest p-value for 'Autoimmune Retinopathy' vs. 'None'
knitr::kable(head(limmaStats[order(limmaStats$`p-val (Autoimmune.retinopathy)`, 
                                   decreasing = FALSE),]), 
      row.names = TRUE, format = "html", align = "c")

Using NanoStringDiff

NanoStringDiff models the data using a negative binomial approximation, which has been shown to be generally more accurate for gene expression count data (Wang et al. 2016). We have provided a function to convert processed, unnormalized, data to a NanoStringSet for use with NanoStringDiff.

# Remember to set normalization = "none"
datNoNorm <- processNanostringData(nsFiles = example_data,
                             sampleTab = sample_info, 
                             idCol = "RCC_Name",
                             groupCol = "Sample_Diagnosis",
                             normalization = "none")

## 
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Averaging technical replicates.....

nsDiffSet <- makeNanoStringSetFromEset(datNoNorm)

colnames(pData(nsDiffSet))

## [1] "Autoimmune retinopathy" "None"                   "Retinitis pigmentosa"

Then you can run as described in the NanoStringDiff vignette (see vignette for full details). Notably, the glm.LRT step is very slow for sample sizes above 10 (and may still be slow with fewer samples). Interested users could also consider using DESeq2 if a negative binomial method is desired, although that method would not explicitly handle the various control genes (Love, Huber, and Anders 2014).

### This block is not run! ###

# This step is fast
nsDiffSet <- NanoStringDiff::estNormalizationFactors(nsDiffSet)

# This step likely to take multiple hours on standard desktop computers.
result <- NanoStringDiff::glm.LRT(nsDiffSet, 
                                  design.full = as.matrix(pData(nsDiffSet)),
                                  contrast = c(1, -1, 0)) 
                                  #Contrast: Autoimmune retinopathy vs. None

Exporting GSEA results

fgseaPostprocessingXLSX allows you to output the results of gene set analyses to an Excel spreadsheet (fgsePostprocessing is similar, and provides .txt files). A summary sheet shows the overall GSEA results, while an additional table for each separate analysis (A vs. Control, B vs. Control, etc.) shows the differential expression statistics and expression profiles for leading edge genes. This step requires input of the FGSEA results, the leading edge results, and the Limma results. It will cluster the pathways if specified, prior to generating results tables.

fgseaPostprocessingXLSX(genesetResults = fgseaResults,
                        leadingEdge = leading.edge,
                        limmaResults = limmaResults,
                        join.threshold = 0.5,
                        filename = "analysis.xlsx")

You can also use groupedGSEAtoStackedReport to generate the gene-level report for one comparison.

results.AR <- groupedGSEAtoStackedReport( 
    fgsea.grouped,
    leadingEdge = leading.edge$Autoimmune.retinopathy,
    de.fit = limmaResults)

# View Cluster 1 gene statistics
resultsTab <- results.AR[results.AR$Cluster == 1, 1:6]
resultsTab[,2:6] <- lapply(resultsTab[,2:6], format, digits = 1, nsmall = 1)

knitr::kable(resultsTab, 
      row.names = FALSE, format = "html", align = "c")