--- title: "Example gene set analysis: The case of SARS-CoV2 infection" author: "Mark Ziemann & Kaumadi Wijesooriya" 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 Here we are performing an analysis of some gene expression data to demonstrate the difference between ORA and FCS methods and to highlight the differences caused by improper background gene set use. 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/ ```{r,begin} suppressPackageStartupMessages({ library("getDEE2") library("DESeq2") library("clusterProfiler") library("mitch") library("kableExtra") library("eulerr") }) ``` ## 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) 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) )) head(mres$enrichment_result,10) %>% kbl() %>% kable_paper("hover", full_width = F) m_up_nom <- subset(mres$enrichment_result,pANOVA<0.05 & s.dist > 0)[,1] m_dn_nom <- subset(mres$enrichment_result,pANOVA<0.05 & s.dist < 0)[,1] ``` ## ORA with clusterprofiler Clusterprofiler uses a hypergeometric test. Firstly I will conduct the analysis separately for up and down regulated genes and with the correct backgound (as intended by the developers). ```{r,cp1} genesets2 <- read.gmt("ReactomePathways.gmt") de_up <- rownames(subset(de,log2FoldChange>0,padj<0.05)) de_up <- unique(gt[which(rownames(gt) %in% de_up),1]) de_dn <- rownames(subset(de,log2FoldChange<0,padj<0.05)) de_dn <- unique(gt[which(rownames(gt) %in% de_dn),1]) de_bg <- rownames(de) de_bg <- unique(gt[which(rownames(gt) %in% de_bg),1]) c_up <- as.data.frame(enricher(gene = de_up, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr")) c_up <- rownames(subset(c_up, p.adjust < 0.05)) c_dn <- as.data.frame(enricher(gene = de_dn, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr")) c_dn <- rownames(subset(c_dn, p.adjust < 0.05)) c_up_nom <- as.data.frame(enricher(gene = de_up, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="none" )) c_up_nom <- rownames(subset(c_up_nom, pvalue < 0.05)) c_dn_nom <- as.data.frame(enricher(gene = de_dn, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="none")) c_dn_nom <- rownames(subset(c_dn_nom, pvalue < 0.05)) ``` Now performing ORA with clusterprofiler with whole genome background list ```{r,cp3} wg_bg <- w$GeneInfo$GeneSymbol f_up <- as.data.frame(enricher(gene = de_up, universe = wg_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr")) f_up <- rownames(subset(f_up, p.adjust < 0.05)) f_dn <- as.data.frame(enricher(gene = de_dn, universe = wg_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr")) f_dn <- rownames(subset(f_dn, p.adjust < 0.05)) f_up_nom <- as.data.frame(enricher(gene = de_up, universe = wg_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="none")) f_up_nom <- rownames(subset(f_up_nom, pvalue < 0.05)) f_dn_nom <- as.data.frame(enricher(gene = de_dn, universe = wg_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="none")) f_dn_nom <- rownames(subset(f_dn_nom, pvalue < 0.05)) f_de_nom <- union(f_up_nom,f_dn_nom) ``` ## Bar chart of significance Here the idea is to classify the pathways as 1 not significant, 2 nominally significant, or 3 FDR significant ```{r,stacked1} c_up_df <- as.data.frame(enricher(gene = de_up, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr" ,qvalueCutoff=1,pvalueCutoff=1)) n_c_up_ns <- nrow( subset(c_up_df,pvalue>0.05) ) n_c_up_nom <- nrow( subset(c_up_df,pvalue<0.05 ) ) n_c_up_fdr <- nrow( subset(c_up_df,p.adjust<0.05) ) c_dn_df <- as.data.frame(enricher(gene = de_dn, universe = de_bg, maxGSSize = 5000, TERM2GENE = genesets2, pAdjustMethod="fdr" ,qvalueCutoff=1,pvalueCutoff=1)) n_c_dn_ns <- nrow( subset(c_dn_df,pvalue>0.05) ) n_c_dn_nom <- nrow( subset(c_dn_df,pvalue<0.05 ) ) n_c_dn_fdr <- nrow( subset(c_dn_df,p.adjust<0.05) ) n_m_up <- length(m_up) n_m_dn <- length(m_dn) n_m_up_nom <- length(m_up_nom) n_m_dn_nom <- length(m_dn_nom) par(mar=c(5,10,5,2)) ngenes <- c("ORA up fdr"=n_c_up_fdr,"ORA up nom"=n_c_up_nom, "ORA dn fdr"=n_c_dn_fdr,"ORA dn nom"=n_c_dn_nom, "FCS up fdr"=n_m_up,"FCS up nom"=n_m_up_nom, "FCS dn fdr"=n_m_dn,"FCS dn nom"=n_m_dn_nom ) barplot(ngenes, horiz=TRUE,las=1) c_up <- subset(c_up_df,p.adjust<0.05)$ID c_dn <- subset(c_dn_df,p.adjust<0.05)$ID c_up_nom <- subset(c_up_df,pvalue<0.05)$ID c_dn_nom <- subset(c_dn_df,pvalue<0.05)$ID n_c_fdr <- length(union(c_dn,c_up)) n_c_nom <- length(union(c_dn_nom,c_up_nom)) n_m_fdr <- length(subset(mres$enrichment_result,p.adjustANOVA<0.05 )[,1]) n_m_nom <- length(subset(mres$enrichment_result,pANOVA<0.05 )[,1]) par(mar=c(5,5,5,2)) ngenes <- c("ORA FDR<0.05"=n_c_fdr,"ORA p<0.05"=n_c_nom,"FCS FDR<0.05"=n_m_fdr,"FCS p<0.05"=n_m_nom) barplot(ngenes,ylab="no. gene sets") text((0:3*1.2)+0.7,ngenes-50,labels=ngenes,cex=1.1) ``` ## Venn diagram comparison The Venn (or Euler to be more correct) diagram is useful to visualise the overlaps between sets. ```{r,venn1} par(cex.main=0.5) par(mar=c(2,2,2,2)) v1 <- list("FCS up"=m_up, "FCS dn"=m_dn, "ORA up"=c_up,"ORA dn"=c_dn) plot(euler(v1),quantities = TRUE, edges = "gray", main="FCS compared to ORA") v0 <- list("FDR up"=m_up, "FDR dn"=m_dn, "Nom up"=m_up_nom,"Nom dn"=m_dn_nom) plot(euler(v0),quantities = TRUE, edges = "gray", main="Effect of FDR correction on FCS results") v0 <- list("FDR up"=c_up, "FDR dn"=c_dn, "Nom up"=c_up_nom,"Nom dn"=c_dn_nom) plot(euler(v0),quantities = TRUE, edges = "gray", main="Effect of FDR correction on ORA results") ora_nom <- union(c_up_nom,c_dn_nom) ora_fdr <- union(c_up,c_dn) fcs_nom <- union(m_up_nom,m_dn_nom) fcs_fdr <- union(m_up,m_dn) v3 <- list("ORA nom"=ora_nom, "ORA FDR"=ora_fdr, "FCS nom"=fcs_nom,"FCS FDR"=fcs_fdr) plot(euler(v3),quantities = TRUE, edges = "gray", main="Effect of FDR correction") v2 <- list("ORA up"=c_up,"ORA dn"=c_dn, "ORA* up"=f_up,"ORA* dn"=f_dn ) plot(euler(v2),quantities = TRUE, edges = "gray", main="Effect of inappropriate background* (whole genome)") ``` ## Jaccard calculation ```{r,jacc1} # FCS vs ORA cm <- length(intersect(c(c_up,c_dn), c(m_up,m_dn))) / length(union(c(c_up,c_dn), c(m_up,m_dn))) #FCS fdr vs nom fcs_fdr_nom <- length(intersect(c(fcs_nom), c(fcs_fdr))) / length(union(c(fcs_nom), c(fcs_fdr))) #ORA fdr vs nom ora_fdr_nom <- length(intersect(c(ora_nom), c(ora_fdr))) / length(union(c(ora_nom), c(ora_fdr))) m_up <- gsub("^","up ",m_up) m_dn <- gsub("^","dn ",m_dn) m_de <- union(m_up,m_dn) c_up <- gsub("^","up ",c_up) c_dn <- gsub("^","dn ",c_dn) c_de <- union(c_up,c_dn) f_up <- gsub("^","up ",f_up) f_dn <- gsub("^","dn ",f_dn) f_de <- union(f_up,f_dn) f_up_nom <- gsub("^","up ",f_up_nom) f_dn_nom <- gsub("^","dn ",f_dn_nom) f_de_nom <- union(f_up,f_dn_nom) # ORA vs ORA* cf <- length(intersect(c_de, f_de )) / length(union(c_de, f_de)) # ORA vs ORA*nom cfn <- length(intersect(c_de, f_de_nom )) / length(union(c_de, f_de_nom)) dat <- c( "FCS vs ORA"=cm, "FCS: FDR vs nominal"=fcs_fdr_nom, "ORA: FDR vs nominal"= ora_fdr_nom, "ORA vs ORA*"=cf, "ORA vs ORA*nom"=cfn) dat saveRDS(dat,file = "ex6dat.rds") par(mar=c(5,10,3,1)) barplot(rev(dat),xlab="Jaccard index",horiz = TRUE, las =1, xlim=c(0,.8) , main=name) text( x=rev(dat)-0.05 , y= 1:length(rev(dat))*1.2-0.5, labels = signif(rev(dat),2)) ``` ## Session information ```{r,session} sessionInfo() ```