Example gene set analysis: The case of SARS-CoV2 infection
Mark Ziemann & Kaumadi Wijesooriya
2025-04-02
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
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 svgGet expression data
I’m using some RNA-seq data looking at the difference in transcriptome caused by SARS-CoV-2 infection.
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)| SRR_accession | QC_summary | SRX_accession | SRS_accession | SRP_accession | Experiment_title | GEO_series | trt | |
|---|---|---|---|---|---|---|---|---|
| 149273 | SRR11517674 | PASS | SRX8089264 | SRS6456133 | SRP253951 | GSM4462336: Series5_A549_Mock_1; Homo sapiens; RNA-Seq | GSE147507 | 0 |
| 149274 | SRR11517675 | PASS | SRX8089265 | SRS6456134 | SRP253951 | GSM4462337: Series5_A549_Mock_2; Homo sapiens; RNA-Seq | GSE147507 | 0 |
| 149275 | SRR11517676 | PASS | SRX8089266 | SRS6456135 | SRP253951 | GSM4462338: Series5_A549_Mock_3; Homo sapiens; RNA-Seq | GSE147507 | 0 |
| 149276 | SRR11517677 | PASS | SRX8089267 | SRS6456136 | SRP253951 | GSM4462339: Series5_A549_SARS-CoV-2_1; Homo sapiens; RNA-Seq | GSE147507 | 1 |
| 149277 | SRR11517678 | PASS | SRX8089268 | SRS6456137 | SRP253951 | GSM4462340: Series5_A549_SARS-CoV-2_2; Homo sapiens; RNA-Seq | GSE147507 | 1 |
| 149278 | SRR11517679 | PASS | SRX8089269 | SRS6456139 | SRP253951 | GSM4462341: Series5_A549_SARS-CoV-2_3; Homo sapiens; RNA-Seq | GSE147507 | 1 |
w<-getDEE2("hsapiens",samplesheet$SRR_accession,metadata=mdat,legacy = TRUE)## For more information about DEE2 QC metrics, visit
## https://github.com/markziemann/dee2/blob/master/qc/qc_metrics.mdx<-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.
# filter out lowly expressed genes
x1<-x1[which(rowSums(x1)/ncol(x1)>=(10)),]
nrow(x)## [1] 39297nrow(x1)## [1] 15182Now 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.
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.
y <- DESeqDataSetFromMatrix(countData = round(x1), colData = s1, design = ~ trt)## converting counts to integer modey <- DESeq(y)## estimating size factors## estimating dispersions## gene-wise dispersion estimates## mean-dispersion relationship## final dispersion estimates## fitting model and testingde <- 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)| baseMean | log2FoldChange | lfcSE | stat | pvalue | padj | |
|---|---|---|---|---|---|---|
| ENSG00000113739 | 3523.734 | 4.004593 | 0.1012568 | 39.54888 | 0 | 0 |
| ENSG00000176153 | 8851.309 | -2.785312 | 0.0820808 | -33.93378 | 0 | 0 |
| ENSG00000058085 | 5767.111 | 3.699041 | 0.1148334 | 32.21223 | 0 | 0 |
| ENSG00000169710 | 6917.854 | -2.259158 | 0.0703230 | -32.12547 | 0 | 0 |
| ENSG00000170421 | 41561.919 | -2.233482 | 0.0742019 | -30.10007 | 0 | 0 |
| ENSG00000100297 | 2476.082 | -2.364383 | 0.0802583 | -29.45965 | 0 | 0 |
Now let’s have a look at some of the charts showing differential expression. In particular, an MA plot and volcano plot.
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)