Add solution to exercise 3 for day 3

Day_3/Solutions/Exercise_3_Multicluster_R_code_2020.R

+# Let's learn how to perform comparative multi cluster analysis
+#
+# Matt Rogon
+# rogon@embl.de
+# 10.2020
+# ver. 0.5
+
+library(dplyr)
+library(STRINGdb)
+library(reshape2)
+library(ReactomePA)
+library(clusterProfiler)
+library(pathview)
+library(DOSE)
+library(org.Mm.eg.db)
+
+
+setwd("/Users/rogon/Work/01. Teaching/CBNA Courses/2020 CBNA-DeNBI Enrichment course - my materials/Part 2 - R-code session ClusterProfiler, ReactomePA, pathfindR/Exercise 4/")
+
+# load complete 5 cluster-dataset
+hit_list <- read.delim("550prot_clusters.txt", stringsAsFactors=FALSE)
+
+#--------------------------------------------------------------
+# let's annotate uniprot to entrez for ClusterProfiler
+#------------------------------------------------------------------------------------
+# Annotation with bitr in ClusterProfiler for enrichment check
+#------------------------------------------------------------------------------------
+keytypes(org.Mm.eg.db)
+entrez_map <- bitr(hit_list$UniProt, fromType="UNIPROT", toType=c("ENTREZID", "SYMBOL"), OrgDb="org.Mm.eg.db", drop= FALSE)
+write.table(entrez_map, file="output/entrez_map.txt", row.names = F, quote = FALSE, sep="\t")
+
+
+#--------------------------------------------------------------
+# Enrichment with ClusterProfiler
+#--------------------------------------------------------------
+# expression values:
+#data(geneList)
+#de <- names(geneList)[geneList > 1]
+
+# to do enrichment on a 'per cluster' basis I need the cluster column from the data frame: main_string_mapping_select
+input_for_enrichment <- merge(hit_list, entrez_map, by.x = "UniProt", by.y = "UNIPROT", all = TRUE)
+input_for_enrichment <- input_for_enrichment[complete.cases(input_for_enrichment),]
+
+
+
+# extract entrez id's from the annotation frame 'results' created earlier
+# for each cluster separately
+
+#c1 <- subset(input_for_enrichment, cluster==1)
+#c1_s <- dplyr::select(c1, ENTREZID)
+#c1_s <- unique(c1_s)
+#c1_list <- as.list(c1_s$ENTREZID)
+
+#same in 2 lines
+c1 <- dplyr::select(subset(input_for_enrichment, cluster ==1), ENTREZID)
+c1_list <- as.list(c1$ENTREZID)
+
+c2 <- dplyr::select(subset(input_for_enrichment, cluster ==2), ENTREZID)
+c2_list <- as.list(c2$ENTREZID)
+
+c3 <- dplyr::select(subset(input_for_enrichment, cluster ==3), ENTREZID)
+c3_list <- as.list(c3$ENTREZID)
+
+c4 <- dplyr::select(subset(input_for_enrichment, cluster ==4), ENTREZID)
+c4_list <- as.list(c4$ENTREZID)
+
+c5 <- dplyr::select(subset(input_for_enrichment, cluster ==5), ENTREZID)
+c5_list <- as.list(c5$ENTREZID)
+
+
+hit_list_Entrez_flat <- unlist(hit_list_Entrez)
+hit_list_Entrez2_flat <- unlist(hit_list_Entrez2)
+
+
+#---------------------------------------------------------------------------------------
+# Individual cluster enrichment
+#---------------------------------------------------------------------------------------
+
+go_1 <- enrichGO(c1_list, OrgDb = "org.Mm.eg.db", ont="BP", universe = input_for_enrichment$ENTREZID, pvalueCutoff=0.05, pAdjustMethod="BH", qvalueCutoff=0.1, minGSSize = 5)
+head(as.data.frame(go_1))
+
+
+#----
+# Multi-cluster analysis
+#----
+
+input_all = list()
+input_all$X1 <- c1_list
+input_all$X2 <- c2_list
+input_all$X3 <- c3_list
+input_all$X4 <- c4_list
+input_all$X5 <- c5_list
+summary(input_all)
+
+# Gene Ontology comparison
+res_go <- compareCluster(input_all, fun="enrichGO", OrgDb="org.Mm.eg.db", pvalueCutoff=0.1) 
+dotplot(res_go, title="GO Enrichment Comparison")
+
+res_reactome <- compareCluster(input_all, fun="enrichPathway", organism='mouse', pvalueCutoff=0.01, pAdjustMethod="BH", qvalueCutoff=0.1)
+dotplot(res_reactome, title="Reactome Enrichment Comparison")
+
+
+
+
+
+