####  Enrichment analysis course, SIB, June 26th 2020


# Enrichment analysis, SIB course, June 26th 2020

# load the packages needed for the R exercise
library(clusterProfiler)
library(pathview)

# 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


#### ---- 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.rds
# - Look for adatpive immune response gene set in results table
# - Count the number of negative NES values and the number of positive NES values


#### ---- 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.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.rds
# count the number of adj. p-values below 0.05


#### ---- Exercise 4

# Create figures for the GSEA results:
# - barplot of –log10(adj. p-value) of top 10 GO p-values
# - 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 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


#### ---- 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.rds"
# - count the number of significant adjusted p-values
# - use barplot() and gseaplot() for the visualization of the results