1.1 Clustering on Heatmap

Now that you understand the concepts of hierarchical clustering both at the sample and at the gene level, we can use a heatmap function to explore the visual consequences of clustering. Here, we can make use of the pheatmap() function, which by default will do the clustering of the rows and columns.

R

cl <- pheatmap( logcpm[top_var,] , scale = "row" , color = colorRampPalette(c("navy","white","firebrick"))(90), border_color = NA, cluster_cols = F,cutree_rows = 2)

gene_clusters <- cutree(cl$tree_row,k = 2)

2.1 enrichR

You can list all available databases by using the command listEnrichrDbs() function.

library(enrichR)

genes_cluster1 <- names(gene_clusters)[gene_clusters == 1]
genes_cluster2 <- names(gene_clusters)[gene_clusters == 2]

head(logcpm[genes_cluster1,])
head(logcpm[genes_cluster2,])

2.1.1 GO enrichment (with enrichR)

R

go_cluster2 <- enrichr(genes = genes_cluster2,databases = "GO_Biological_Process_2018")
go_cluster2 <- go_cluster2$GO_Biological_Process_2018
go_cluster2 <- go_cluster2[order(go_cluster2$P.value),]
go_cluster2[1:5,]

mypar(1,1,mar=c(2,20,2,2))
barplot( -log10(go_cluster2$P.value[15:1]), horiz = T, names.arg = go_cluster2$Term[15:1],las=1)

Exploratory Questions

• Which is the most enriched GO term for cluster1?
• How many genes from your data set were detected in this most GO term? What is the percentage of genes out of the total list that defines the GO term?
• Which genes from your data set belong to this most enriched GO term?
• Some genes make part of several GO terms, can you list some?

2.1.2 KEGG enrichment (with enrichR)

R

kegg_cluster2 <- enrichr(genes = genes_cluster2,databases = "KEGG_2019_Human")
kegg_cluster2 <- kegg_cluster2$KEGG_2019_Human
kegg_cluster2 <- kegg_cluster2[order(kegg_cluster2$P.value),]
kegg_cluster2[1:5,]

mypar(1,1,mar=c(2,20,2,2))
barplot( -log10(kegg_cluster2$P.value[15:1]), horiz = T, names.arg = kegg_cluster2$Term[15:1],las=1)

Exploratory Questions

• Which is the most enriched KEGG pathway for cluster1?
• How many genes from your data set were detected in this most enriched KEGG pathway ?
• What is the percentage of genes out of the total list that defines the KEGG pathway ?
• Which genes from your data set belong to this most enriched KEGG pathway ?
• Some genes make part of several KEGG pathways, can you list some?
• Taking the results from GO term enrichment and the KEGG enrichment, what is a most likely function happening in your data set?

---

2.2 GSEA

MSigDB is one of the largest curated databases of signatures and pathways created by the BROAD INSTITUTE. The Molecular Signatures Database (MSigDB) is accessible here.

For this exercise, you can download a couple of data sets:

As you could already notice, the differences in gene expression between days t0 to t24 are very clear, but not the ones from t2. Here we will illustrate a gene set enrichment comparing t2 to t0 in order to explore which pathways are UP or DOWN regulated at the first 2 hours of the experiment.

# Create a gene rank based on the gene expression
mean_t0 <- rowMeans( logcpm[,metadata$Time == "t0"] )
mean_t2 <- rowMeans( logcpm[,metadata$Time == "t2"] )
gene_rank <- mean_t0 - mean_t2

# Load hallmark pathways
hallmark_pathways <- gmtPathways("data/MSigDB_files/h.all.v6.2.symbols.gmt.txt")

Once our list of genes are sorted (from highest to lowest expression in t2), we can proceed with the enrichment itself. Here, we will be using the fgsea() function. This will result in a table containing information for several pathways. We can then sort and filter those pathways to visualize only the top ones. You can select/filter them by either p-value or normalized enrichemnet score (NES).

R

# Perform enrichemnt analysis
fgseaRes <- fgsea( pathways=hallmark_pathways, stats=gene_rank, minSize=15, maxSize=500, nperm=10000)

# Filter the results table to show only the top 10 UP and DOWN regulated processes (optional)
topPathwaysUp <- fgseaRes[ES > 0][head(order(pval), n=10), pathway]
topPathwaysDown <- fgseaRes[ES < 0][head(order(pval), n=10), pathway]
topPathways <- c(topPathwaysUp, rev(topPathwaysDown))

# Nice summary table (shown as a plot)
plotGseaTable(hallmark_pathways[topPathways], gene_rank, fgseaRes, gseaParam = 0.5)

Exploratory Questions

• Which is the most significantly up-regulated pathway at time t2?
• What do the length of the black bars mean?
• Is there any pathway that is significantly down-regulated at t2? Which columns tells you what is Up of downregulated?
• Which is the most enriched pathway at time t6?
• What happens if you do t2-t0 instead of t0-t2 ? What changes and why?

Checking for individual enrichment:

R

topPathways

# Nice plot to see enrichment of one specific pathway
plotEnrichment(hallmark_pathways[["HALLMARK_MYC_TARGETS_V2"]],  gene_rank)

Exploratory Questions

• What do the black lines in the X-axis mean?
• By looking at the green line distribution, and the comparison made (t0 - t2), at which time point is this pathways enriched.

Exploratory Questions (optional)

• What if instead of the hallmark gene set, you use the GO? Try that.
• How different are the results from EnrichR and GSEA using GO annotation? Remember to use the same comparison (i.e. t2 vs t0).

2.3 Online enrichment tools:

Below you can find a list with the most commonly used online enrichment tools. Each one differs a bit on what they can do. The packages are sorted in order, we would like you to work on the workshop.

In case you want to test the online tools above, you can use the code below to copy the gene vector into memory, and then paste it to the webtools. Should work for both Mac and Windows users.

clip <- pipe("pbcopy", "w")
write.table(genes_cluster1, file=clip, sep = '\t', row.names = FALSE)         
close(clip)

2.4 Protein-protein interactions

In the course we will not have much time to work on gene expression networks. However, there are some nice databases that allow you to easily visualize protein-protein interactions. Please paste the list of genes using the command above into those sites to have a visual insight on how the Up- or Down- regulated genes interact with each other.

End of document