workshop.rmd

---
title: "ISCB 2022 Gene set enrichment workshop"
author: "January Weiner"
date: "`r Sys.Date()`"
output: 
  html_document:
    code_folding: show  
toc: no
bibliography: bibliography.bib
link-citations: true
---

```{r,echo=FALSE}
## Set default options for the knitr RMD processing
knitr::opts_chunk$set(warning=FALSE,message=FALSE,fig.width=5,fig.height=5,cache=TRUE,autodep=TRUE, results="hide")
```

```{r libraries,cache=FALSE}
library(tidyverse)
library(colorDF)
library(ggplot2)
library(ggbeeswarm)
library(tmod)
library(colorDF)
library(cowplot)
theme_set(theme_minimal())
```


# A simple data set

## Basic analysis

This is a data set that I used in my poster. It is available from GEO and
contains RNA-Seq counts for healthy individuals and COVID-19 patients, described by @mick2020upper.

The data set has been already processed and downloaded; however, below you
will find the code to replicate this procedure (if you want).




```{r pheno_data_download,eval=FALSE}
#| class.source: fold-hide
## this sometimes fails. The RDS of the object is provided
Sys.setenv(VROOM_CONNECTION_SIZE=8*131072)

## with getGEO, we can only get the phenoData
geo <- getGEO("GSE156063")[[1]]
saveRDS(geo, file="GSE156063.geo")
```

```{r pheno_data_cleanup,eval=FALSE}
#| class.source: fold-hide
geo <- readRDS("GSE156063.geo")
covar <- as(phenoData(geo), "data.frame")

## remove columns with only one value
boring <- map_lgl(covar, ~ length(unique(.x)) == 1)
covar <- covar[ , !boring ] 

## clean up
covar <- covar %>% 
  dplyr::rename(gender = "gender:ch1") %>%
  mutate(disease = gsub(" .*", "", .data[["disease state:ch1"]])) %>%
  mutate(label = description) %>%
  mutate(group = disease) %>%
  arrange(description) %>%
  dplyr::select(all_of(c("title", "label", "gender", "disease", "group")))
```

```{r featuredata_download,eval=FALSE}
#| class.source: fold-hide
## the counts must be downloaded from GEO separately.
if(!file.exists("GSE156063_swab_gene_counts.csv.gz")) {
  download.file("https://ftp.ncbi.nlm.nih.gov/geo/series/GSE156nnn/GSE156063/suppl/GSE156063_swab_gene_counts.csv.gz",
                "GSE156063_swab_gene_counts.csv.gz")
}

counts <- read_csv("GSE156063_swab_gene_counts.csv.gz")
.tmp <- counts[[1]]
counts <- as.matrix(counts[,-1])
rownames(counts) <- .tmp
#counts <- counts[ , covar$description ]
counts <- counts[ , covar$label ]
lcpm   <- edgeR::cpm(counts, log=TRUE)
#stopifnot(all(colnames(counts) == covar$description))
stopifnot(all(colnames(counts) == covar$label))

annot <- data.frame(ENSEMBL = rownames(counts))

.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("ENTREZID"), keytype="ENSEMBL")
annot$ENTREZID <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]

.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("SYMBOL"), keytype="ENSEMBL")
annot$SYMBOL <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]

.tmp <- mapIds(org.Hs.eg.db, annot$ENSEMBL, column=c("GENENAME"), keytype="ENSEMBL")
annot$GENENAME <- .tmp[ match(annot$ENSEMBL, names(.tmp)) ]
```


```{r Preparation_of_the_covariates,cache=FALSE,eval=FALSE}
#| class.source: fold-hide
sel <- covar$group %in% c("no", "SC2")
counts <- counts[ , sel ]
covar  <- covar[ sel, ]
covar$group <- as.character(covar$group)

covar$group.disease <- paste0(covar$group, '_', covar$disease)
rownames(covar) <- covar$label
saveRDS(covar, file="covar.rds")
saveRDS(counts, file="counts.rds")

```

```{r DESeq2,cache=FALSE,eval=FALSE}
#| class.source: fold-hide
## DESeq2 calculations
## manual cache, since the operation takes a long time
sel <- !is.na(covar$group)

ds2 <- DESeqDataSetFromMatrix(counts[,sel], colData=covar[sel, ], design=~ disease)

ds2_file <- "ds2_cached.rds"
if(!file.exists(ds2_file)) {
  message("Running DESeq2")
  ds2 <- DESeq(ds2)
  saveRDS(ds2, file=ds2_file)
} else {
  message("Reading ds2 from manual cache")
  ds2 <- readRDS(ds2_file)
}
```

First, take a look at the structure of the data:


```{r results="markdown"}
covar  <- readRDS("data/covar.rds") # covariates
summary_colorDF(covar)
```





The results table is stored under "data/disease_SC2_vs_no.tsv". In
`data/annotation.tsv`, there is the annotation of the genes in this data
set.

```{r}
dsres <- read_table("data/disease_SC2_vs_no.tsv")
annot <- read_table("data/annotation.tsv")
dsres <- merge(annot, dsres, by="ENSEMBL")
```

Basic CERNO analysis:


```{r basic_tmod,results="markdown"}
#| fig.width=6,
#| fig.height=12
gl <- dsres$SYMBOL[ order(dsres$pvalue) ]
tres <- tmodCERNOtest(gl)
print(tres)
```

Which we can quickly visualize:


```{r basic_panelplot}
#| fig.width=6,
#| fig.height=12
ggPanelplot(tres)
```

We can also add information about whether the genes go up or down. Note
that we need to put `tres` inside a list, so `ggPanelplot` can match tmod
results with the `sgenes` object.


```{r basic_panelplot2}
#| fig.width=6,
#| fig.height=12
sgenes <- tmodDecideTests(dsres$SYMBOL, lfc = dsres$log2FoldChange,
                          pval=dsres$pvalue)
names(sgenes) <- "SC2_vs_no"
tres <- list("SC2_vs_no"=tres)
ggPanelplot(tres, sgenes=sgenes)
```


## Using another gene sets

Another built-in gene set are the cell surface markers.


```{r cellsignatures}
#| fig.width=6,
#| fig.height=12
data(cell_signatures)
tres2 <- tmodCERNOtest(gl, mset=cell_signatures)

tres2 <- list("SC2_vs_no"=tres2)
sgenes2 <- tmodDecideTests(dsres$SYMBOL, lfc = dsres$log2FoldChange,
                          pval=dsres$pvalue, mset=cell_signatures)
names(sgenes2) <- "SC2_vs_no"
ggPanelplot(tres2, sgenes=sgenes2, mset=cell_signatures)
```

The mset parameter is used in most of the functions to use another gene
sets.





## Inspecting individual enrichments

Let us consider the gene set "LI.M68" – RIG-1 like receptor signalling.
First, show the genes in this gene set:


```{r results="markdown"}
data(tmod) # make the built-in gene sets visible
getModuleMembers("LI.M68", mset=tmod)
```

Evidence plot – a ROC curve:


```{r evidence plots}
#| fig.width=12,
#| fig.height=12
par(mfrow=c(2,2))
evidencePlot(gl, m = "DC.M1.2", gene.labels=TRUE)
evidencePlot(gl, m = "LI.M68", gene.labels=TRUE)
evidencePlot(gl, m = "LI.M165")
evidencePlot(gl, m = "LI.M4.0")
```

## Eigengenes



The basic idea is as follos: run PCA only on a selected group of genes, use
PC1 instead of average gene expression of that gene set (mind the signs!).
Below, `counts.rds` contains the raw RNA-seq counts, which I first convert
to log~2~ counts per million.

```{r eigengenes_manual}
data(tmod)
counts <- readRDS("data/counts.rds")

lcpm <- edgeR::cpm(counts, log=TRUE)
mm <- getModuleMembers("LI.M68", mset=tmod)[[1]]
sel <- annot$SYMBOL %in% mm
covar$eig <- prcomp(t(lcpm[ sel, ]), scale.=TRUE)$x[,1]
```

We can use eigengenes to nicely plot the results:


```{r eigengenes}
ggplot(covar, aes(x=disease, y=eig)) + geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(cex=3)
```

(However, mind the sign!)

In tmod, eigengenes can be computed directly. Tmod makes sure that the
eigengene correlates positively with the majority of the genes in the gene
set and reverses the sign when necessary:


```{r eigengenes4}
#| fig.width=10,
#| fig.height=10
selm <- tres[[1]]$ID[1:4]
eig <- eigengene(lcpm, annot$SYMBOL)
eig <- as.data.frame(t(eig)) %>%
  select(all_of(selm)) %>%
  rownames_to_column("label") %>% 
  inner_join(covar, by="label")
plots <- map(selm, ~ ggplot(eig, aes_string(x="disease", y=.x)) +
             geom_boxplot(outlier.shape=NA) + 
             geom_beeswarm(cex=3) +
             ggtitle(tres[[1]]$Title[ match(.x, tres[[1]]$ID) ]))
plot_grid(plotlist=plots)
```



# Combining GSE with other tools

## Combining gene set enrichments with PCA

Using a method like tmodCERNOtest, tmodUtest or fsgsea makes it possible to
apply gene set enrichment to any single ranked list of genes.

(I have selected the components based on correlation with `covar$disease`)

```{r pca1}
#| fig.width=6,
#| fig.height=6
pca <- prcomp(t(lcpm), scale.=TRUE)
pca_df <- pca$x %>% as.data.frame %>% 
  rownames_to_column("label") %>% merge(covar)
ggplot(pca_df, aes(x=PC5, y=PC6, color=disease)) + geom_point()
```

We can sort all the genes by their decreasing absolute loadings in the pca object.



```{r}
#| fig.width=12,
#| fig.height=8
pcares <- map(1:10, ~ {
                gl <- annot$SYMBOL[ order(-abs(pca$rotation[,.x])) ]
                tmodCERNOtest(gl)
             })
names(pcares) <- colnames(pca$rotation)[1:10]
sgenes <- tmodDecideTests(annot$SYMBOL, 
                          lfc=pca$rotation[, 1:10],
                          pval=NULL,
                          lfc.thr = 0)
ggPanelplot(pcares, sgenes=sgenes, filter_row_auc = .85)
```




```{r}
ggplot(pca_df, aes(x=PC2, y=PC3, color=disease)) + geom_point()
```

Check the function `tmod_pca` for more details.


## Machine learning


```{r}
library(randomForest)
vars <- apply(lcpm, 1, var)
sel <- vars > quantile(vars, .75)
rf <- randomForest(t(lcpm[sel, ]), factor(covar$disease), importance=TRUE)
print(rf)
```


```{r}
if(require(myfuncs)) {
  rocplot(rf2pr(rf), confmat = F)
}

```

We can run gene set enrichment analysis on the genes ordered by variable
importance.



```{r}
rf$importance %>% as.data.frame() %>%
  rownames_to_column("ENSEMBL") %>%
  merge(annot) %>% arrange(-MeanDecreaseAccuracy) %>%
  pull(SYMBOL) %>% tmodCERNOtest()
```

## Correlation analysis

Let us consider the gene ANKRD22. Based on coexpression, what function is
it related to? We will first remove the disease effect so as we are
clustering only by variations in expression within the groups, and not by
changes between the groups.


```{r results="markdown"}
lcpm_c <- lm(t(lcpm) ~ disease, data=covar)$residuals
lcpm_c <- t(lcpm_c)

ankrd22_i <- which(annot$SYMBOL == "ANKRD22")
cor_ank <- cor(lcpm_c[ankrd22_i, ], t(lcpm_c)) %>% t() %>%
  as.data.frame() %>% rownames_to_column("ENSEMBL") %>%
  merge(annot) 
tcor <- cor_ank %>% arrange(-abs(V1)) %>% pull(SYMBOL) %>%
  tmodCERNOtest() %>% filter(N1 < 30 & adj.P.Val < 1e-3) %>%
  arrange(-AUC)
print(tcor)
```






# Using gene sets from msigdbr

Converting the data frame from the msigdbr package for use with tmod is straightforward:


```{r}
library(msigdbr)
msig_df <- msigdbr("Homo sapiens")

goset <- msig_df %>%
  filter(gs_subcat == "GO:BP") %>%
  makeTmodFromDataFrame(module_col="gs_name",
                        title_col="gs_name",
                        feature_col="gene_symbol")
```


```{r results="markdown"}
library(DT)
tmodCERNOtest(gl, mset=goset, qval = 1e-3) %>%
  filter(N1 > 5) %>%
  datatable() %>%
  formatSignif(c("adj.P.Val", "P.Value")) %>%
  formatRound(c("cES", "AUC"))
```

# Disentangling results

## Upset plots



```{r}
#| fig.width=14,
#| fig.height=8
tres[[1]] %>% filter(N1 < 20) %>% pull(ID) %>% upset()
```


## Overlaps between gene sets


```{r results="markdown"}
#| fig.width=12,
#| fig.height=8
foo <- tres[[1]] %>% 
  slice(1:20) %>%
  pull(ID) %>% 
  modCorPlot(stat = "overlap")
foo
```







## "Leading edge" analysis

Which genes are driving the enrichment?


```{r}
evidencePlot(gl, m = "LI.M165")
```

```{r}
evidencePlot(gl, m = "LI.M165", style="gsea")
```


Get the leading edge:


```{r results="markdown"}
lea <- tmodLEA(gl, "LI.M165")
lea
```


```{r}
clr <- ifelse(gl %in% lea[[1]], "red", "grey")
names(clr) <- gl
evidencePlot(gl, m = "LI.M165", 
               gene.labels=TRUE,
               gene.colors = clr)
```

And GSEA style:


```{r}
clr <- ifelse(gl %in% lea[[1]], "red", "grey")
names(clr) <- gl
evidencePlot(gl, m = "LI.M165", 
               gene.labels=TRUE,
               gene.colors = clr,
               style="gsea")
abline(v=attr(lea[[1]], "LEA"), col="red", lwd=2)
```





# The vaccination example


This example data set has been published by @weiner2019characterization.
RNA has been collected daily from individuals vaccinated with one of seven
vaccines over the course of several days. The data set has been collected
using microarrays. The pre-calculated results include two influenza
vaccines, one containing an adjuvant and the other not. For details how to
get this table, see the online tmod manual
[here](https://january3.github.io/tmod/articles/user_manual/tmod_user_manual.html),
section "Transcriptional responses to vaccination".

```{r vaccine_read_data}
tt <- readRDS("data/vaccination_toptable.rds")
```

The subset contains results for 5 time points for two vaccines (10
contrasts in total). The two vaccines are denoted with "A" (no adjuvant)
and "F" (adjuvant). We will run gene set enrichment on each of these time
points using `tmodCERNOtest`.



```{r vaccine_tmod}
contrasts_v <- paste0(rep(c("F", "A"), each=5), "_", 1:5)
sort_cols <- paste0("qval.", contrasts_v)
res <- map(sort_cols, ~ {
  tmodCERNOtest(tt$GENE_SYMBOL[ order(tt[[.x]]) ])
})
names(res) <- contrasts_v
```

The figure below shows the results.


```{r vaccine_panelplot0,fig.width=12,fig.height=8}
ggPanelplot(res, filter_row_q = 1e-3, filter_row_auc = .85, q_thr=0.01)
```

Do these gene sets contain genes that go up or down? For this, we need to
determine how many genes in each gene set are significantly going up or
down.



```{r vaccine_sgenes}
pval <- tt %>% select(starts_with("qval"))
lfc  <- tt %>% select(starts_with("logFC"))
sgenes <- tmodDecideTests(tt$GENE_SYMBOL, lfc=lfc, pval=pval)
names(sgenes) <- contrasts_v
```

```{r vaccine_panelplot,fig.width=12,fig.height=8}
ggPanelplot(res, sgenes=sgenes, filter_row_q = 1e-3, filter_row_auc = .85, q_thr=0.01)
```

## Concordance / discordance plots

Are responses on day 1 in Fluad similar to responses on day 1 in Agrippal?

 * test for interaction, or
 * use disco plots

$$d_s = \log_2 FC_A  \cdot \log_2 FC_B \cdot |\log_{10}p_A + \log_{10}p_B|$$


```{r}
disco <- function(lfc1, lfc2, pval1, pval2) {
	lfc1 * lfc2 * abs(log10(pval1) + log10(pval2))
}

tt <- tt %>% mutate(disco=disco(logFC.F_1, logFC.A_1, qval.F_1, qval.A_1)) %>%
  mutate(disco=ifelse(abs(disco) > .5, sign(disco), disco)) %>%
  arrange(abs(disco))



ggplot(tt, aes(x=logFC.F_1, y=logFC.A_1, color=disco)) +
  scale_color_gradient2(low="blue", mid="grey", high = "red") +
  geom_point(alpha=.1) +
  geom_abline(slope = 1, intercept = 0, color="grey")

```

How about F day 1 and F day 5?



```{r}
tt <- tt %>% mutate(disco=disco(logFC.F_1, logFC.F_5, qval.F_1, qval.F_5)) %>%
  mutate(disco=ifelse(abs(disco) > .5, sign(disco), disco)) %>%
  arrange(abs(disco))



ggplot(tt, aes(x=logFC.F_1, y=logFC.F_5, color=disco)) +
  scale_color_gradient2(low="blue", mid="grey", high = "red") +
  geom_point(alpha=.1) +
  geom_abline(slope = 1, intercept = 0, color="grey")
```

We can now make gene set enrichment of the *discordant* and *concordant*
gene sets.


```{r}
#| fig.width=6,
#| fig.height=12
tdisco <- list()
tdisco$discordant <- tt %>% arrange(disco) %>% pull(GENE_SYMBOL) %>%
  tmodCERNOtest() %>% filter(AUC > .55)
tdisco$concordant <- tt %>% arrange(-disco) %>% pull(GENE_SYMBOL) %>%
  tmodCERNOtest() %>% filter(AUC > .55)
sgenes <- tmodDecideTests(tt$GENE_SYMBOL, lfc=tt$disco, lfc.thr = 0)
sgenes <- list(discordant=sgenes[[1]], concordant=sgenes[[1]])
ggPanelplot(tdisco, sgenes=sgenes)
```

Note: there is more to the gene set enrichment analysis with disco. There
will be more on that in both the tmod package and the tmod manual.




# Bibliography