ScriptTiO2.Rmd

# Script and session info

Data pre-processing:
Purpose of script: pre-process gene-expression data that will be used in: https://github.com/laurent2207/TiO2-scripts.
Gene-expression data: GEO:GSE42069

Part 1:
Purpose of script: Takes one gene list for a process, finds relevant pathways and runs enrichment analysis (using transcriptomics data) to study how process is affected in dataset

Part 2:
Purpose of script: pre-process GO-term genelists that will be used in: https://github.com/laurent2207/TiO2-scripts.
GO-term genelists: Apoptopic process (GO:0006915), Inflammatory response (GO:0006954)
GO-term genelists: Cellular response to DNA damage stimulus (GO:0006974), response to oxidative stress (GO:0006979)

Part 3:
Purpose of script: pre-process GO-term genelists that will be used in: https://github.com/laurent2207/TiO2-scripts.

Author: Laurent Winckers, Martina Kutmon

Date Created: 2020-03-24

# Session info:
R version 3.6.3 (2020-03-24)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 17763)

Pacakges Data pre-processing: rstudioapi_0.11, biomaRt_2.42.0, dplyr_0.8.5, EnhancedVolcano_1.4.0 

Packages Part 1: rstudioapi_0.11, clusterProfiler_3.14.3, plyr_1.8.6, biomaRt_2.42.0, dplyr_0.8.5, data.table_1.12.8, pheatmap_1.0.12,
colorRamps_2.3, RColorBrewer_1.1-2, enrichplot_1.8.1, DOSE_3.14.0, org.Hs.eg.db_3.11.4, ggpubr_0.4., ggplot2_3.3.2

Packages Part 2: rstudioapi_0.11, biomaRt_2.42.0, clusterProfiler_3.14.3, tidyr_1.1.1, dplyr_0.8.5, ggplot2_3.3.2

Packages Part 3: rstudioapi_0.11, data.table_1.12.8, pheatmap_1.0.12, colorRamps_2.3, RColorBrewer_1.1-2, igraph_1.2.5, ggplot2_3.3.2


###############################
##### Data pre-processing #####
###############################

# Install required packages and set up BioMart
```{r}
library(rstudioapi)
library(biomaRt)
library(dplyr)
library(EnhancedVolcano)
```

# Set up environment
```{r}
#clear workspace and set string as factors to false
rm(list=ls())
options(stringsAsFactors = F)

#set working directroy
#setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) does not work in markdown as they are not linked to RStudio
```

# Load gene expression data (analysed with arrayanalysis.org)
```{r}
# Caco2
caco2_low <- read.table("./data/data_pre_processing/TiO2_24hrs_10_Caco2.txt", sep = "\t", header = T)
caco2_high <- read.table("./data/data_pre_processing/TiO2_24hrs_100_Caco2.txt", sep = "\t", header = T)
# SAE
SAE_low <- read.table("./data/data_pre_processing/TiO2_24hrs_10_SAE.txt", sep = "\t", header = T)
SAE_high <- read.table("./data/data_pre_processing/TiO2_24hrs_100_SAE.txt", sep = "\t", header = T)
# THP1
THP1_low <- read.table("./data/data_pre_processing/TiO2_24hrs_10_THP1.txt", sep = "\t", header = T)
THP1_high <- read.table("./data/data_pre_processing/TiO2_24hrs_100_THP1.txt", sep = "\t", header = T)

# Select necessary columns, remove columns that are not needed further down the process
# Caco2
caco2_low <- caco2_low[c(1,2,3,6)]
caco2_high <- caco2_high[c(1,2,3,6)]
# SAE
SAE_low <- SAE_low[c(1,2,3,6)]
SAE_high <- SAE_high[c(1,2,3,6)]
# THP1
THP1_low <- THP1_low[c(1,2,3,6)]
THP1_high <- THP1_high[c(1,2,3,6)]

# Change column names adressing respective cell line and either 10 ug/ml, low (L) or 100 ug/ml, high (H) concentration
# Caco2
colnames(caco2_low)[c(2,3,4)] <- c("caco2_L_logFC", "caco2_L_FC", "caco2_L_pval")
colnames(caco2_high)[c(2,3,4)] <- c("caco2_H_logFC", "caco2_H_FC", "caco2_H_pval")
# SAE
colnames(SAE_low)[c(2,3,4)] <- c("SAE_L_logFC", "SAE_L_FC", "SAE_L_pval")
colnames(SAE_high)[c(2,3,4)] <- c("SAE_H_logFC", "SAE_H_FC", "SAE_H_pval")
# THP1
colnames(THP1_low)[c(2,3,4)] <- c("THP1_L_logFC", "THP1_L_FC", "THP1_L_pval")
colnames(THP1_high)[c(2,3,4)] <- c("THP1_H_logFC", "THP1_H_FC", "THP1_H_pval")
```

# Identifier mapping
```{r}
### select ensembl IDs from one of the datasets
ids <- as.data.frame(caco2_low$ENSG_ID)
colnames(ids) <- "ENSG_ID"

### map Ensemble IDs to Entrez Gene and HGNC symbols
ensembl <- useEnsembl("ensembl", dataset = "hsapiens_gene_ensembl", mirror = "useast")
genes <- getBM(attributes = c('ensembl_gene_id', 'hgnc_symbol', 'entrezgene_id'), filters = 'ensembl_gene_id', values = ids$ENSG_ID, mart = ensembl)

### remove genes without Entrez Gene identifier from gene list 11791 -> 11324
genes <- genes[!(is.na(genes$entrezgene_id)),]
genes <- genes[!duplicated(genes$entrezgene_id),]
```

# Merge datasets with annotated gene identifiers
```{r}
colnames(genes)[1] <- "ENSG_ID"
data <- merge(genes, caco2_low, by = "ENSG_ID")
data <- merge(data, caco2_high, by = "ENSG_ID")
data <- merge(data, SAE_low, by = "ENSG_ID")
data <- merge(data, SAE_high, by = "ENSG_ID")
data <- merge(data, THP1_low, by = "ENSG_ID")
data <- merge(data, THP1_high, by = "ENSG_ID")
```

# Data visualization (volcano plots)
```{r}
cl <- EnhancedVolcano(data,title = "CACO2/HT29-MTX 10 μg/ml 24h", lab = data$hgnc_symbol, x = 'caco2_L_logFC',y = 'caco2_L_pval', xlim = c(-2, 2), FCcutoff = 0.26, pCutoff = 0.05, col = c("grey30", "orange", "royalblue", "darkorange4"), selectLab = '', ylim = c(0, 13.5))

ch <- EnhancedVolcano(data,title = "CACO2/HT29-MTX 100 μg/ml 24h", lab = data$hgnc_symbol, x = 'caco2_H_logFC',y = 'caco2_H_pval', xlim = c(-2, 2), FCcutoff = 0.26, pCutoff = 0.05, col = c("grey30", "orange", "royalblue", "darkorange4"), selectLab = '', ylim = c(0, 13.5))

sl <- EnhancedVolcano(data,title = "SAE 10 μg/ml 24h", lab = data$hgnc_symbol, x = 'SAE_L_logFC',y = 'SAE_L_pval', xlim = c(-2, 2), FCcutoff = 0.26, pCutoff = 0.05, col = c("grey30", "orange", "royalblue", "darkorange4"), selectLab = '', ylim = c(0, 13.5))

sh <- EnhancedVolcano(data,title = "SAE 100 μg/ml 24h", lab = data$hgnc_symbol, x = 'SAE_H_logFC',y = 'SAE_H_pval', xlim = c(-2, 2), FCcutoff = 0.26, pCutoff = 0.05, col = c("grey30", "orange", "royalblue", "darkorange4"), selectLab = '', ylim = c(0, 13.5))

tl <- EnhancedVolcano(data,title = "THP1 10 μg/ml 24h", lab = data$hgnc_symbol, x = 'THP1_L_logFC',y = 'THP1_L_pval', xlim = c(-2, 2), FCcutoff = 0.26, pCutoff = 0.05, col = c("grey30", "orange", "royalblue", "darkorange4"), selectLab = '', ylim = c(0, 13.5))

th <- EnhancedVolcano(data,title = "THP1 100 μg/ml 24h", lab = data$hgnc_symbol, x = 'THP1_H_logFC',y = 'THP1_H_pval', xlim = c(-2, 2), FCcutoff = 0.26, pCutoff = 0.05, col = c("grey30", "orange", "royalblue", "darkorange4"), selectLab = '', ylim = c(0, 13.5))

svg(file="./output/volcanoplots.svg", width = 17, height = 18)
par(mar=c(1, 1, 1, 1))
cowplot::plot_grid(cl, ch, sl, sh, tl, th, ncol=2, nrow=3)
dev.off()
```

# Save results and ranked value for GSEA 
```{r}
### save combined annotated gene-expression file
write.table(data, "./output/TiO2-dataset.txt", sep = "\t", quote = F, row.names = F)
```


##################
##### PART 1 #####
##################

# Install required packages and set up BioMart
```{r}
library(rstudioapi)
library(clusterProfiler)
library(plyr)
library(biomaRt)
library(dplyr)
library(data.table)
library(pheatmap)
library(colorRamps)
library(RColorBrewer)
library(enrichplot)
library(DOSE)
library(org.Hs.eg.db)
library(ggpubr)
library(ggplot2)
```

# Set up environment
```{r}
#clear workspace and set string as factors to false
rm(list=ls())
options(stringsAsFactors = F)

#set working directory
#setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) does not work in markdown as they are not linked to RStudio
```

# Enrichment analysis
```{r}
# load gene-expression file 
data <- read.table("output/TiO2-dataset.txt", header = T, sep ="\t")

### load geneset
wp2gene <- clusterProfiler::read.gmt("data/gmt_wp_Homo_sapiens.gmt")
wp2gene <- wp2gene %>% tidyr::separate(term, c("name","version","wpid","org"), "%")
wpid2gene <- wp2gene %>% dplyr::select(wpid,gene) #TERM2GENE
wpid2name <- wp2gene %>% dplyr::select(wpid,name) #TERM2NAME

comparisons <- c("caco2_L","caco2_H","SAE_L","SAE_H","THP1_L", "THP1_H")

source("functions/enrichment.R")

# set output destination
path = "output/ORApw_cb_"

enrichment(wpid2gene,wpid2name, data, comparisons, 0.26, path)
```


##################
##### PART 2 #####
##################

# Set up environment
```{r}
#clear workspace and set string as factors to false
rm(list=ls())
options(stringsAsFactors = F)

#set working directroy
#setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) does not work in markdown as they are not linked to RStudio
```

# Install required packages and set up BioMart
```{r}
library(rstudioapi)
library(biomaRt)
library(clusterProfiler)
library(tidyr)
library(dplyr)
library(ggplot2)

biomart <- biomaRt::useEnsembl("ensembl", dataset = "hsapiens_gene_ensembl", host="uswest.ensembl.org")
```

# Define GO terms of interest
```{r}
go.terms <- c("GO:0006915","GO:0006954","GO:0006974","GO:0006979")
```

# Load genesets
```{r}
### load genesets
wp2gene <- clusterProfiler::read.gmt("data/gmt_wp_Homo_sapiens.gmt")
wp2gene <- wp2gene %>% tidyr::separate(term, c("name","version","wpid","org"), "%")
wpid2gene <- wp2gene %>% dplyr::select(wpid,gene) #TERM2GENE
wpid2name <- wp2gene %>% dplyr::select(wpid,name) #TERM2NAME
```

# Create GO annotations for all GO terms of interest and run enrichment analysis
```{r}
source("functions/go_annotations.R")
source("functions/ora_enrichment.R")

# set output directory for ora_enrichtment function
path = "output/"

for(i in go.terms) {
  res <- go_annotations(go.term = i, biomart = biomart, output = paste("output/",gsub(":", "", i),".txt", sep=""))
  ora_enrichment(genes = res$entrezgene_id, wpid2gene = wpid2gene, wpid2name = wpid2name, path = path, prefix = paste0(gsub(":", "", i), "_ORA"))
  # update workflow data file
  file.copy(from = paste("output/", gsub(":", "", i),".txt",sep=""), to = paste("data/",gsub(":", "", i),".txt",sep=""), overwrite = TRUE)
}

```


##################
##### PART 3 #####
##################

# Set up environment
```{r}
#clear workspace and set string as factors to false
rm(list=ls())
options(stringsAsFactors = F)

#set working directroy
#setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) does not work in markdown as they are not linked to RStudio
```

# Install required packages
```{r}
library(rstudioapi)
library(data.table)
library(pheatmap)
library(colorRamps)
library(RColorBrewer)
library(igraph)
library(ggplot2)
```

# Load ORA results for GO-terms - only select rows that are significant (adjusted p-value < 0.05)
```{r}
path = paste0(getwd(), "/output")
ls_GO <- list.files(path = path, pattern = "_ORA.txt")
for (i in 1:length(ls_GO)){
  assign(paste0("GO", i), read.table(paste0(path,"/",ls_GO[i]), header = T, sep = "\t", quote = ""))
}

GO1 <- GO1[GO1$p.adjust < 0.05,]
GO2 <- GO2[GO2$p.adjust < 0.05,]
GO3 <- GO3[GO3$p.adjust < 0.05,]
GO4 <- GO4[GO4$p.adjust < 0.05,]

for (i in 1:nrow(GO1)){
GO1$Perc[i] <- (as.numeric(sub("\\/.*", "", GO1$GeneRatio[i])) / as.numeric(sub("\\/.*", "", GO1$BgRatio[i])))*100
}

for (i in 1:nrow(GO2)){
GO2$Perc[i] <- (as.numeric(sub("\\/.*", "", GO2$GeneRatio[i])) / as.numeric(sub("\\/.*", "", GO2$BgRatio[i])))*100
}

for (i in 1:nrow(GO3)){
GO3$Perc[i] <- (as.numeric(sub("\\/.*", "", GO3$GeneRatio[i])) / as.numeric(sub("\\/.*", "", GO3$BgRatio[i])))*100
}

for (i in 1:nrow(GO4)){
GO4$Perc[i] <- (as.numeric(sub("\\/.*", "", GO4$GeneRatio[i])) / as.numeric(sub("\\/.*", "", GO4$BgRatio[i])))*100
}

GO150 <- GO1[GO1$Perc > 50,]
GO160 <- GO1[GO1$Perc > 60,]
GO170 <- GO1[GO1$Perc > 70,]
GO180 <- GO1[GO1$Perc > 80,]

GO250 <- GO2[GO2$Perc > 50,]
GO260 <- GO2[GO2$Perc > 60,]
GO270 <- GO2[GO2$Perc > 70,]
GO280 <- GO2[GO2$Perc > 80,]

GO350 <- GO3[GO3$Perc > 50,]
GO360 <- GO3[GO3$Perc > 60,]
GO370 <- GO3[GO3$Perc > 70,]
GO380 <- GO3[GO3$Perc > 80,]

GO450 <- GO4[GO4$Perc > 50,]
GO460 <- GO4[GO4$Perc > 60,]
GO470 <- GO4[GO4$Perc > 70,]
GO480 <- GO4[GO4$Perc > 80,]


write.table(GO150, "output/GO1.txt", quote = F, sep = "\t", row.names = F)
write.table(GO250, "output/GO2.txt", quote = F, sep = "\t", row.names = F)
write.table(GO350, "output/GO3.txt", quote = F, sep = "\t", row.names = F)
write.table(GO450, "output/GO4.txt", quote = F, sep = "\t", row.names = F)
```

# Combine the significant results of all four GO-terms
```{r}
sigORA <- rbind(GO150, GO250, GO350, GO450)
sigORA <- as.data.frame(unique(sigORA[,c(1)]))
colnames(sigORA) <- c("ID")

for (i in 1:nrow(sigORA)) { 
  if (sigORA$ID[i] %in% GO1$ID) {
  sigORA$pid1[i] <- "1"
} else {
  sigORA$pid1[i] <- "0"
}
  if (sigORA$ID[i] %in% GO2$ID) {
  sigORA$pid2[i] <- "1"
} else {
  sigORA$pid2[i] <- "0"
} 
  if (sigORA$ID[i] %in% GO3$ID) {
  sigORA$pid3[i] <- "1"
} else {
  sigORA$pid3[i] <- "0"
}
  if (sigORA$ID[i] %in% GO4$ID) {
  sigORA$pid4[i] <- "1"
} else {
  sigORA$pid4[i] <- "0"
}}
```

# Read pathway ORA results and only select rows where pvalue < 0.05
```{r}
path = paste0(getwd(), "/output")
ls_pw <- list.files(path = path, pattern = "ORApw_")
for (i in 1:length(ls_pw)){
  assign(gsub(".txt", "", paste0(ls_pw[i])), read.table(paste0(path,"/",ls_pw[i]), header = T, sep = "\t", quote = ""))
}

ls_pw <- mget(ls(pattern = 'ORApw_'))
ls_pw <- lapply(ls_pw, function(x){x[x$pvalue<0.05,]}) 
#list2env(ls_pw, envir = .GlobalEnv)
```

# Filter ORA results of pathways - create plots
```{r}
source("functions/filteredPlot.R")

# set output directory
path = "output/"

dir.create(paste0(path, "nodesEdges"))

for (i in 1:length(ls_pw)){
filteredPlot(data = ls_pw[[i]], path, fileName = gsub("ORApw_", "", names(ls_pw[i])))
}
```