#### Enrichment analysis course, SIB, December 4th 2020 # load the packages needed for the R exercise library(clusterProfiler) library(pathview) library(org.Hs.eg.db) # Some reminders about the usage of R: # to look for help on a function and which arguments to change use ? ?t.test # to search for a keyword in all or in a particular package, use ?? ??adjust ??stats::adjust # Some reminders about the usage of R: # to look for help on a function and which arguments to change use ? ?t.test # to search for a keyword in all or in a particular package, use ?? ??adjust ??stats::adjust # to access a row in a data frame: iris[3, ] # to search for a word: which(iris$Species=="setosa") # to count events: table(iris$Species) # or summary(iris$Species) #### ---- Exercise 1 # Differential gene expression between two different # immune cell types, natural killer (NK) cells and CD4 T helper cells (Th)). # The data was generated by RNA sequencing and analyzed using limma. # Import file NK_vs_Th_diff_gene_exp.csv that contains gene ensembl ID, # gene symbols, log2 fold change, t-statistic and raw p-value # 1. Is CPS1 significantly differentially expressed at alpha=0.05? # 2. How many genes are up-regulated and how many genes are down-regulated # after BH p-value adjustment? # 3. Is CPS1 still differentially expressed after p-value adjustment at alpha=0.05? # Steps: # - import csv file NK_vs_Th_diff_gene_exercise_1.csv # - look for CPS1 gene in table # - create new column with adjusted p-values # - again look for CPS1 in table # To set the working directory (i.e. where your files are located) # for your own computer, change this path: setwd("/export/scratch/twyss/SIB_training") NK_vs_Th<-read.csv("NK_vs_Th_diff_gene_exercise_1.csv", header = T) NK_vs_Th[which(NK_vs_Th$symbol=="CPS1"),] NK_vs_Th$p.adj <- p.adjust(NK_vs_Th$P.Value, method = "BH") NK_vs_Th[which(NK_vs_Th$symbol=="CPS1"),] # ensembl_gene_id symbol logFC t P.Value p.adj # 348 ENSG00000021826 CPS1 -1.637697 -2.079127 0.04496309 0.1565113 #### ---- Exercise 2 # - Is the adaptive immune response gene set significantly enriched in genes up-regulated in NK vs Th? # - How many GO gene sets are significant after GSEA (use minGSSize=30, pvalueCutoff=0.1) ? # - Is the adaptive immune response gene set significant? Up-reg. or down-reg.? # - Are the majority of gene sets rather up-regulated or down-regulated? # Steps: # - Import the list of genes belonging to the adaptive immune # response: adaptive_immune_response_geneset_term2gene.csv" # - Run Fisher test of up-regulated genes # - Create named vector of t-statistics and run gsea of built-in GO gene sets: # set argument minGSSize=30. # NOTE: if the gsea fails, use readRDS() to import the pre-computed results # contained in file gseGO_Nk_vs_Th_results_122020.rds # - Look for adaptive immune response gene set in results table # - Count the number of negative NES values and the number of positive NES values adimresp_term2name<-read.csv("adaptive_immune_response_geneset_term2gene.csv", header=T) NK_up<-subset(NK_vs_Th, NK_vs_Th$p.adj<0.05&NK_vs_Th$logFC>0) # 1957 NK_up<-NK_up$symbol summary(NK_up %in% adimresp_term2name$V2) # Mode FALSE TRUE # logical 1882 75 NK_not_DE<-subset(NK_vs_Th, NK_vs_Th$p.adj>0.05) NK_not_DE<- NK_not_DE$symbol # 16633 summary(NK_not_DE %in% adimresp_term2name$V2) # Mode FALSE TRUE # logical 16363 270 cont.table<-matrix(c(75, 1882, 270, 16363), ncol=2, byrow=F) colnames(cont.table)<-c("up", "not_up") rownames(cont.table)<-c("in_set", "not_in_set") fisher.test(cont.table) # gseGO: gl<-NK_vs_Th$t names(gl)<-make.names(NK_vs_Th$symbol, unique = T) gl<-sort(gl, decreasing = T) gene.idtype.bods$hsa # set the random seed generator: set.seed(1234) GO_NK_Th<-gseGO(gl, ont="BP", OrgDb = org.Hs.eg.db, keyType = "SYMBOL", minGSSize = 30, seed=T) View(GO_NK_Th@result) GO_NK_Th@result[GO_NK_Th@result$Description=="adaptive immune response", ] # "Down-regulated" genesets: # 56/101 summary(GO_NK_Th@result$NES<0) summary(GO_NK_Th@result$NES>0) # up 101/56 # To order the GO results table according to decreasing NES value: ordered_GO <- GO_NK_Th@result[order(GO_NK_Th@result$NES, decreasing = T), ] # Second GSEA with same data to test "stability" GO_NK_Th_2<-gseGO(gl, ont="BP", OrgDb = org.Hs.eg.db, keyType = "SYMBOL", minGSSize = 30, seed=T) View(GO_NK_Th_2@result) #### ---- Exercise 3 # - First convert the gene symbols to EntrezID to perform a GSEA of KEGG # pathways (with argument minGSSize=30). # - Are the majority of gene sets rather up-regulated or down-regulated? # - Is there a KEGG immune-related gene set coming up? Is there a KEGG Natural # killer gene set coming up? # - If you want to see which genes are included in one of the built-in KEGG pathways, # where could you find this information? # - Import the hallmark gene sets and run a GSEA. How many significant gene sets are there? # Steps: # - Use this function to see the types of ID that can be converted: keytypes(org.Hs.eg.db) # Use bitr to convert Ensembl IDs to ENTREZID and create a named vector of t-statistics # to provide as gene list to the gseKEGG function, use minGSSize = 30 # NOTE: if the gsea fails, use readRDS to import the pre-computed results # contained in the file gseKEGG_Nk_vs_Th_results_122020.rds # - count the number of positive NES values and of negative NES values # - To search for a partial word in a table, eg. "immune", use grep() # - use str() to view the structure and content of a gsearesult object # - use read.gmt to import hallmark genesets in file h.all.v7.1.symbols.gmt # use function GSEA() providing a named vector of t-statistics # NOTE: if the gsea fails, use readRDS to import the pre-computed results # contained in the file hGSEA_Nk_vs_Th_results_122020.rds # count the number of adj. p-values below 0.05 # convert symbols EntrezID (=ncbi-geneid) keytypes(org.Hs.eg.db) gene_convert<-bitr(as.character(NK_vs_Th$ensembl_gene_id), fromType = "ENSEMBL", toType = c("SYMBOL", "ENTREZID"), OrgDb = "org.Hs.eg.db") gl_kegg<-cbind(ENSEMBL=as.character(NK_vs_Th$ensembl_gene_id), t=NK_vs_Th$t) gl_kegg<-merge(gl_kegg, gene_convert) gl_kegg_list<-as.numeric(gl_kegg$t) names(gl_kegg_list)<-gl_kegg$ENTREZID gl_kegg_list<-sort(gl_kegg_list, decreasing = T) # run the GSEA of KEGG pathways: KEGG_NK_Th<-gseKEGG(gl_kegg_list, organism = "hsa", keyType = "ncbi-geneid", minGSSize = 30) View(KEGG_NK_Th@result) summary(KEGG_NK_Th@result$NES<0) KEGG_NK_Th@result[grep("immune",KEGG_NK_Th@result$Description), ] KEGG_NK_Th@result[grep("killer",KEGG_NK_Th@result$Description), ] KEGG_NK_Th@geneSets$hsa03010 # Import hallmark ?GSEA() term2gene_h<-read.gmt("h.all.v7.1.symbols.gmt") h_NK_vs_Th<-GSEA(gl, TERM2GENE = term2gene_h, seed=T) View(h_NK_vs_Th@result) str(KEGG_NK_Th) strsplit(h_NK_vs_Th@result$core_enrichment, "/")[[1]] #### ---- Exercise 4 # Create figures for the GSEA results: # (- barplot of –log10(adj. p-value) of top 10 GO p-values) # - dotplot GO gene sets # - barcode plot of one of the significant hallmark gene sets (choose one) # - pathview map for KEGG Natural Killer mediated cytotoxicity (optional: with none-significant genes in grey) # Steps: # - use function barplot() and provide -log10(adj. p-value) of GO GSEA results and name the bars # - use function dotplot() with GO gseGO() object output # - use function gseaplot() and provide ID of gene set # - create a named vector of fold change values (optional: setting the logFC of non-significant # genes to 0), search the KEGG ID of the pathway (a hsa0000 type ID) to provide to pathview() function pdf("barplot_GO_top10.pdf") par(mar=c(5, 16, 3,3)) # bottom, left, top, right barplot(-log10(GO_NK_Th@result$p.adjust)[1:10], horiz = T, names=GO_NK_Th@result$Description[1:10], las=2, xlab = "-log10(adj.P.value", col="blue") dev.off() # dotplot(KEGG_NK_Th) gseaplot(KEGG_NK_Th, "hsa03010" ) # pathview: NK_vs_Th$logFC_0<-ifelse(NK_vs_Th$p.adj>0.05, 0, NK_vs_Th$logFC) genePW<-NK_vs_Th$logFC_0 names(genePW)<-NK_vs_Th$symbol pathview(gene.data = genePW, pathway.id = "hsa04650", species="hsa", gene.idtype = "SYMBOL") #### ---- Extra exercise for credits # - Perform GSEA of the NK vs Th data using the Reactome gene sets downloaded on the # MSigDB website. # - How many gene sets are significantly enriched? Generate an ordered barplot of # the NES of all genesets, and generate a barcode plot for the gene set with the lowest NES # Steps # - import the reactome gene sets in the file c2.cp.reactome.v7.1.symbols.gmt using read.gmt() # Use GSEA providing a sorted named vector of t-statistics and the reactome gene sets, # use minGSSize=30 # NOTE: if the GSEA fails, use readRDS to import the pre-computed results contained in the # file "reactomeGSEA_NK_vs_Th_results_122020.rds" # - count the number of significant adjusted p-values # - use barplot() and gseaplot() for the visualization of the results