forked from lmweber/PrinciplesSTA
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathdifferential-expression.qmd
66 lines (41 loc) · 2.15 KB
/
differential-expression.qmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
# Differential expression {#sec-differential-expression}
## Overview
In this chapter, we perform differential expression testing between clusters or spatial domains to identify representative marker genes for each cluster or spatial domain.
## Load data from previous steps
We start by loading the data object(s) saved after running the analysis steps from the previous chapters. Code to re-run the previous steps is shown in condensed form in @sec-analysis-steps.
```{r, message=FALSE, results='hide'}
library(SpatialExperiment)
library(here)
spe <- readRDS(here("outputs/spe_cluster.rds"))
```
## Differential expression testing
Identify representative marker genes for each cluster or spatial domain by testing for differential gene expression between clusters.
Here, we use the `findMarkers` implementation in `scran` [@Lun2016], using a binomial test, which tests for genes that differ in the proportion expressed vs. not expressed between clusters. This is a more stringent test than the default t-tests, and tends to select genes that are easier to interpret and validate experimentally.
```{r, message=FALSE}
library(scran)
library(scater)
library(pheatmap)
```
```{r}
# set gene names as row names for easier plotting
rownames(spe) <- rowData(spe)$gene_name
# test for marker genes
markers <- findMarkers(spe, test = "binom", direction = "up")
# returns a list with one DataFrame per cluster
markers
```
```{r, fig.width=5, fig.height=6}
# plot log-fold changes for one cluster over all other clusters
# selecting cluster 1
interesting <- markers[[1]]
best_set <- interesting[interesting$Top <= 5, ]
logFCs <- getMarkerEffects(best_set)
pheatmap(logFCs, breaks = seq(-5, 5, length.out = 101))
```
```{r, fig.width=7, fig.height=7}
# plot log-transformed normalized expression of top genes for one cluster
top_genes <- head(rownames(interesting))
plotExpression(spe, x = "label", features = top_genes)
```
## Pseudobulking
Alternatively, we can proceed by manually aggregating the counts per cluster or spatial domain, which is referred to as 'pseudobulking'. Then, we can perform differential expression testing between the pseudobulked clusters or spatial domains.