---
title: "Example gene set analysis: The case of SARS-CoV2 infection"
author: "Burnet Bioinformatics"
date: "`r Sys.Date()`"
output:
  html_document:
    toc: true
    toc_float: true
    code_folding: hide
    fig_width: 7
    fig_height: 5
theme: cosmo
---

Source: https://github.com/markziemann/SurveyEnrichmentMethods

## Intro

The dataset being used is SRP253951 and we are comparing A549 cells with and without infection with SARS-CoV-2.

Data are obtained from http://dee2.io/

Mitch devel was used (Apr 2, 2025). https://github.com/markziemann/mitch

```{r,begin}

suppressPackageStartupMessages({
library("mitch") #devel > install.packages("mitch",repos=NULL,type="source")
library("getDEE2") #BiocManager::install("getDEE2")
library("DESeq2") #BiocManager::install("DESeq2")
library("kableExtra")
})

knitr::opts_chunk$set(dev = 'svg') # set output device to svg

```

## Get expression data

I'm using some RNA-seq data looking at the difference in transcriptome caused by SARS-CoV-2 infection.

```{r,getdata,fig.width=7,fig.height=7}

name="SRP253951"
mdat<-getDEE2Metadata("hsapiens")
samplesheet <- mdat[grep("SRP253951",mdat$SRP_accession),]

samplesheet <- mdat[which(mdat$SRX_accession %in% c("SRX8089264","SRX8089265","SRX8089266","SRX8089267","SRX8089268","SRX8089269")),]

samplesheet$trt <- factor(c(0,0,0,1,1,1))
s1 <- samplesheet

s1 %>% kbl(caption = "sample sheet") %>% kable_paper("hover", full_width = F)

w<-getDEE2("hsapiens",samplesheet$SRR_accession,metadata=mdat,legacy = TRUE)
x<-Tx2Gene(w)
x <- x$Tx2Gene

# save the genetable for later
gt<-w$GeneInfo[,1,drop=FALSE]
gt$accession<-rownames(gt)

# counts 
x1<-x[,which(colnames(x) %in% samplesheet$SRR_accession)]

```

Here show the number of genes in the annotation set, and those detected above the detection threshold.

```{r,filter}

# filter out lowly expressed genes
x1<-x1[which(rowSums(x1)/ncol(x1)>=(10)),]
nrow(x)
nrow(x1)

```

Now multidimensional scaling (MDS) plot to show the correlation between the datasets.
If the control and case datasets are clustered separately, then it is likely that there will be many differentially expressed genes with FDR<0.05.

```{r,mds}

plot(cmdscale(dist(t(x1))), xlab="Coordinate 1", ylab="Coordinate 2", pch=19, col=s1$trt, main="MDS")

```

## Differential expression

Now run DESeq2 for control vs case.

```{r,deseq2}

y <- DESeqDataSetFromMatrix(countData = round(x1), colData = s1, design = ~ trt)
y <- DESeq(y)
de <- results(y)
de<-as.data.frame(de[order(de$pvalue),])
rownames(de)<-sapply(strsplit(rownames(de),"\\."),"[[",1)
head(de) %>% kbl() %>% kable_paper("hover", full_width = F)

```

Now let's have a look at some of the charts showing differential expression.
In particular, an MA plot and volcano plot.

```{r,deplots,fig.width=7,fig.height=7}

maplot <- function(de,contrast_name) {
  sig <-subset(de, padj < 0.05 )
  up <-rownames(subset(de, padj < 0.05 & log2FoldChange > 0))
  dn <-rownames(subset(de, padj < 0.05 & log2FoldChange < 0))
  GENESUP <- length(up)
  GENESDN <- length(dn)
  DET=nrow(de)
  SUBHEADER = paste(GENESUP, "up, ", GENESDN, "down", DET, "detected")
  ns <-subset(de, padj > 0.05 )
  plot(log2(de$baseMean),de$log2FoldChange, 
       xlab="log2 basemean", ylab="log2 foldchange",
       pch=19, cex=0.5, col="dark gray",
       main=contrast_name, cex.main=0.7)
  points(log2(sig$baseMean),sig$log2FoldChange,
         pch=19, cex=0.5, col="red")
  mtext(SUBHEADER,cex = 0.7)
}

make_volcano <- function(de,name) {
    sig <- subset(de,padj<0.05)
    N_SIG=nrow(sig)
    N_UP=nrow(subset(sig,log2FoldChange>0))
    N_DN=nrow(subset(sig,log2FoldChange<0))
    DET=nrow(de)
    HEADER=paste(N_SIG,"@5%FDR,", N_UP, "up", N_DN, "dn", DET, "detected")
    plot(de$log2FoldChange,-log10(de$padj),cex=0.5,pch=19,col="darkgray",
        main=name, xlab="log2 FC", ylab="-log10 pval", xlim=c(-6,6))
    mtext(HEADER)
    grid()
    points(sig$log2FoldChange,-log10(sig$padj),cex=0.5,pch=19,col="red")
}

maplot(de,name)
make_volcano(de,name)

```

## Gene sets from Reactome

In order to perform gene set analysis, we need some gene sets.

```{r,reactome}

if (! file.exists("ReactomePathways.gmt")) {
  download.file("https://reactome.org/download/current/ReactomePathways.gmt.zip", 
    destfile="ReactomePathways.gmt.zip")
  unzip("ReactomePathways.gmt.zip")
}
genesets<-gmt_import("ReactomePathways.gmt")

```

## FCS with Mitch

Mitch uses rank-ANOVA statistics for enrichment detection.

```{r,mitch1}

m <- mitch_import(de,DEtype = "DEseq2", geneTable = gt)
mres <- mitch_calc(m,genesets = genesets,priority="effect",cores=8)
m_up <- subset(mres$enrichment_result,p.adjustANOVA<0.05 & s.dist > 0)[,1]
m_dn <- subset(mres$enrichment_result,p.adjustANOVA<0.05 & s.dist < 0)[,1]
message(paste("Number of up-regulated pathways:",length(m_up) ))
message(paste("Number of down-regulated pathways:",length(m_dn) ))

up <- head(subset(mres$enrichment_result, s.dist > 0 & p.adjustANOVA < 0.05),20)
dn <- head(subset(mres$enrichment_result, s.dist < 0 & p.adjustANOVA < 0.05),20)

up %>% kbl(caption="Top up-regulated pathways") %>% kable_paper("hover", full_width = F)
dn %>% kbl(caption="Top down-regulated pathways") %>% kable_paper("hover", full_width = F)

```

A normal way to visualise the results is by barplot.

```{r,datavis,fig.height=7}

up <- head(subset(mres$enrichment_result, s.dist > 0 & p.adjustANOVA < 0.05),20)
dn <- head(subset(mres$enrichment_result, s.dist < 0 & p.adjustANOVA < 0.05),20)
top <- rbind(up,dn)
top <- top[order(top$s.dist),]

par(mar=c(5.1, 25.1, 4.1, 2.1))
gsn <- substr(top$set, start = 1, stop = 70 )
barplot(top$s.dist,horiz=TRUE,las=1,xlab="ES",names.arg=gsn,cex.names=0.6,main="Reactome")

```

There are some plots already provided by mitch.

```{r,mitch_plots}

mitch_plots(mres,"myplots.pdf")

```

Cool, but it doesn't show whether the gene sets have common driver genes.

```{r,network,fig.height=9,fig.width=14}

networkplot(mres,FDR=0.05,n_sets=20)

network_genes(mres,FDR=0.05,n_sets=20)

```

## Session information

```{r,session}

save.image("example6.Rdata")

sessionInfo()

```
