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).
##
## 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......
We default to plotting the first two principal components, but you can also choose others.
Differential expression analysis
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.
##
## 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.
## Intercept Autoimmune.retinopathy Retinitis.pigmentosa
## HLA-DQB1 7.422597 -0.9145722 -1.563955039
## KIT 3.799585 0.1385032 -0.009840124
## LAG3 6.266360 -0.3694319 -0.064424116
## SOCS3 7.227935 0.1626320 0.089525085
## TCF7 10.494227 -0.1130549 0.052214067
## IKBKB 7.820259 0.4208591 0.111493933
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")
Results of Limma analyses can be visualized using simple volcano plots.
##
## 'plotContrast' not provided, setting it to Autoimmune.retinopathy
You can add additional ggplot layers as well.
They can also be converted to a simple table or exported to a text file.
|
|
Log2FC (Autoimmune.retinopathy)
|
t (Autoimmune.retinopathy)
|
p-val (Autoimmune.retinopathy)
|
q-val (Autoimmune.retinopathy)
|
Log2FC (Retinitis.pigmentosa)
|
t (Retinitis.pigmentosa)
|
p-val (Retinitis.pigmentosa)
|
q-val (Retinitis.pigmentosa)
|
|
IKBKB
|
0.42
|
4.9
|
1.4e-05
|
0.0065
|
0.11
|
1.3
|
0.2000
|
0.73
|
|
MIF
|
-0.54
|
-4.3
|
8.5e-05
|
0.0180
|
-0.24
|
-1.9
|
0.0620
|
0.66
|
|
IL8
|
-1.20
|
-4.2
|
1.2e-04
|
0.0180
|
-0.59
|
-2.0
|
0.0480
|
0.66
|
|
CR1
|
0.70
|
4.1
|
1.8e-04
|
0.0210
|
0.43
|
2.5
|
0.0150
|
0.62
|
|
PLAU
|
0.92
|
4.0
|
2.6e-04
|
0.0240
|
0.63
|
2.7
|
0.0093
|
0.62
|
|
ITGAM
|
0.65
|
3.9
|
3.0e-04
|
0.0240
|
0.37
|
2.2
|
0.0310
|
0.62
|
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.
##
## Reading in .RCC files......
## Checking codeset contents......
## Checked gene name consistency in .RCC files.
## 15 genes labeled as 'housekeeping'.
## Warning in read_sampleData(dat, file.name = sampleTab, idCol = idCol, groupCol = groupCol, :
## idCol not provided. Assuming the first column of 'GSE117751_sample_data.csv' contains sample ID's.
## Warning in read_sampleData(dat, file.name = sampleTab, idCol = idCol, groupCol = groupCol, :
## Sample names in the two files don't match. NanoTube is
## assuming that samples are in the same order. Please confirm
## with your data.
##
## Averaging technical replicates.....
## [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).
Gene set analysis
Gene Set Enrichment Analysis can be conducted on the Limma results using the fgsea package (Sergushichev 2016). Before running it, it is useful to consider whether Gene Set Enrichment Analysis is appropriate for your specific data set, based on the number and types of genes represented in your chip, and whether any of them were actually differentially expressed.
Gene set databases can be loaded as a list object (either directly or in an .rds file), or as a .gmt (MSigDB or similar) or .tab (CPDB) file. We have included a list of pathways from WikiPathways (Martens et al. 2021).
The min.set option is important, as pathways containing only a few genes present in your data set will probably not provide informative enrichment statistics. We will discard all gene sets where fewer than min.set genes from that set are present in the analysis. You also have the option to rank genes by ‘coefficients’ (frequently, the log2FC) or the ‘t’ statistics. skip.first will skip the first column of the limma design if TRUE (default). Generally, this column is the Intercept of the regression, which is not useful for gene set analysis.
The limmaToFGSEA function will then conduct preranked analysis using fgsea for each column in the Limma coefficients or t-statistic matrix (possibly skipping the first). The output will be a list object containing the results for each analysis.
## [1] "Autoimmune.retinopathy" "Retinitis.pigmentosa"
We can order the pathways by p-value and view the top results. As previously reported in (Lundy et al. 2018), we see that immune pathways are significantly altered in autoimmune retinopathy patients. We also identify EGFR Signaling as a potential pathway of interest.
|
|
pathway
|
pval
|
padj
|
log2err
|
ES
|
NES
|
size
|
leadingEdge
|
|
1
|
EGF/EGFR Signaling Pathway
|
0.0008897
|
0.0115659
|
0.4772708
|
0.6726876
|
2.045655
|
15
|
STAT5A, STAT3 , RAF1 , JAK1 , MAPK1 , MAPK14, STAT5B, SRC , ABL1 , PRKCD
|
|
2
|
IL-1 signaling pathway
|
0.0061725
|
0.0341494
|
0.4070179
|
0.5511928
|
1.770841
|
19
|
IKBKB , IL1B , MAPK1 , MAPK14 , IL1R1 , IKBKG , MAPKAPK2, TOLLIP , IRAK3 , IRAK2 , NFKBIA , NFKB1 , TRAF6 , MYD88
|
|
3
|
White fat cell differentiation
|
0.0078806
|
0.0341494
|
0.3807304
|
0.7747081
|
1.720120
|
6
|
STAT5A, STAT5B, CEBPB , IRF3
|
|
4
|
Angiopoietin Like Protein 8 Regulatory Pathway
|
0.0203272
|
0.0601384
|
0.3524879
|
0.6602543
|
1.683155
|
9
|
RAF1 , MAPK1 , MAPK14, MAP4K4, MAPK11, MAP4K2
|
|
5
|
Kit receptor signaling pathway
|
0.0231301
|
0.0601384
|
0.3524879
|
0.5260238
|
1.616896
|
16
|
STAT5A, STAT3 , RAF1 , MAPK1 , PTPN6 , MAPK14, STAT5B
|
Another similar function, nsdiffToFGSEA, is provided to conduct fgsea on NanoStringDiff results. This one conducts analysis on a single preranked list.
After analysis, the leading edge genes can be extracted for gene sets (with some cutoff for enrichment statistics) using fgseaToLEdge.
The nominal p-value or NES (normalized enrichment score) can also be used as a cutoff. If NES is used, you can either select all gene sets with abs(NES) > cutoff, if nes.abs.cutoff == TRUE. Otherwise, you can select gene sets with NES > cutoff (if cutoff > 0) or NES < cutoff (if cutoff < 0).
A basic leading edge heatmap can then be drawn using the pheatmap package.
You can further cluster pathways by their leading edge genes. This is particularly useful when you have lots (tens to hundreds) of significantly enriched pathways, as you can prioritize certain ones.
| EGF/EGFR Signaling Pathway |
0.0009 |
0.0116 |
0.4773 |
0.6727 |
2.0457 |
15 |
1 |
x |
1 |
1 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| Kit receptor signaling pathway |
0.0231 |
0.0601 |
0.3525 |
0.5260 |
1.6169 |
16 |
1 |
|
1 |
1 |
1 |
0 |
0 |
0 |
0 |
0 |
1 |
0 |
0 |
1 |
1 |
0 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| IL-1 signaling pathway |
0.0062 |
0.0341 |
0.4070 |
0.5512 |
1.7708 |
19 |
2 |
x |
0 |
1 |
1 |
0 |
0 |
0 |
1 |
0 |
0 |
0 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
1 |
0 |
0 |
| White fat cell differentiation |
0.0079 |
0.0341 |
0.3807 |
0.7747 |
1.7201 |
6 |
3 |
x |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
0 |
0 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
1 |
1 |
| Angiopoietin Like Protein 8 Regulatory Pathway |
0.0203 |
0.0601 |
0.3525 |
0.6603 |
1.6832 |
9 |
4 |
x |
1 |
1 |
1 |
1 |
1 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
| B Cell Receptor Signaling Pathway |
0.0556 |
0.1205 |
0.3218 |
0.3990 |
1.4361 |
31 |
5 |
x |
1 |
1 |
1 |
0 |
0 |
0 |
1 |
1 |
1 |
1 |
1 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
0 |
In this example, only “EGF/EGFR Signaling Pathway” and “Kit receptor signaling pathway” are sufficiently similar to be clustered together (“Cluster 1”), and the other four pathways are deemed unique. In the heatmap above, we see that the Kit receptor pathway contains only one leading edge gene unique from the EGF/EGFR pathway (PTPN6). Based on a lower p-value and higher NES, you would consider EGFR Signaling to be more important in this analysis, while the Kit receptor pathway is largely redundant. This is denoted by the “Cluster.Max” variable, which identifies maximum enrichment in each cluster with an “x”.
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.
You can also use groupedGSEAtoStackedReport to generate the gene-level report for one comparison.
## Symbol Cluster Log2FC (Autoimmune.retinopathy) t (Autoimmune.retinopathy)
## 1 ABL1 1 0.1530830 1.503434
## 11 JAK1 1 0.3410859 2.938236
## 14 MAPK1 1 0.2430387 2.570981
## 16 MAPK14 1 0.2725272 2.243526
## 21 PRKCD 1 0.2232819 1.333450
## 22 PTPN6 1 0.2835481 2.433887
## 23 RAF1 1 0.4468993 3.289308
## 24 SRC 1 0.2795798 1.540417
## 25 STAT3 1 0.5492925 3.432543
## 26 STAT5A 1 0.3752109 3.658707
## 27 STAT5B 1 0.2944106 2.220156
## p-val (Autoimmune.retinopathy) q-val (Autoimmune.retinopathy)
## 1 0.1400740841 0.36434571
## 11 0.0053001150 0.07645014
## 14 0.0137036754 0.11953694
## 16 0.0300945473 0.15767315
## 21 0.1894411684 0.42939998
## 22 0.0191886640 0.13445041
## 23 0.0020151593 0.05642446
## 24 0.1308235984 0.34984763
## 25 0.0013383902 0.04550527
## 26 0.0006906003 0.03652508
## 27 0.0317580056 0.15767315
