edgeRによる遺伝子発現解析・GSEAとclusterProfierによるパスウェイ解析

# DEG
library(edgeR)

# enrichment analysis
library(clusterProfiler)

library(VennDiagram)
# library(ggVennDiagram)

library(EnhancedVolcano)

# basic
library(ggrepel)
library(beeswarm)
library(org.Mm.eg.db)
library(doParallel)
library(foreach)
library(pheatmap)
library(rtracklayer)
library(rvest) # this is for extraction of html table; results of leading edge analysis.
library(xml2) # this is for extraction of html table; results of leading edge analysis.
library(openxlsx) # to read xlsx file.
library(ComplexHeatmap)
library(tidyverse)
library(data.table)
library(gridExtra)

# Read GTF.
gtf <- readGFF("~/reference_seq/gencode.vM36.primary_assembly.annotation.gtf")
gene_id_gene_name <- gtf %>% select(all_of(c("gene_id", "gene_name")))%>% distinct(gene_id, .keep_all = TRUE)

# Read original data.
count_original <- read_tsv("/mnt/storage/Data/output_featurecounts/count.txt", skip = 1)
# Remove several columns that will not be used for further analysis
count_original_2 <- count_original %>% select(-c(2:6))
# Rename the columns into simple one as possible as can be.
colnames(count_original_2) <- c("gene_id", "T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9", "T10", "T11", "T12", "T13", "T14", "T15", "T16", "T17", "T18", "T19", "T20", "B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8", "B9", "B10", "B11", "B12", "B13", "B14", "B15", "B16", "B17", "B18", "B19", "B20")

# QC to understand situation of the libraries
check_lib_size <- colSums(count_original_2[,2:ncol(count_original_2)]) %>% data.frame() %>% rename("library_size" = ".")
# Calculate total read count regarding to each gene over all sample.
gene_count <- rowSums(count_original_2[,2:ncol(count_original_2)])%>% data.frame() %>% rename("gene_count_over_sample" = ".") 
rownames(gene_count) <- count_original_2$gene_id

# Check which sample is largely difference from other one.
## Prepare dataset for hclust
count_original_2_transpose <- count_original_2 %>% t() %>% data.frame()
colnames(count_original_2_transpose) <- count_original_2_transpose[1,]
count_original_2_transpose <- count_original_2_transpose[-1,]

# hierarchical clustering
check_cluster <- hclust(dist(count_original_2_transpose), method = "ward.D2", members = NULL)
plot(check_cluster)

# PCA for count_original_2
## prepare dataset
count_original_2_transpose <- apply(count_original_2_transpose, 2, as.numeric) %>% data.frame(row.names = colnames(count_original_2[,2:41])) 

# Calculate PCA
pca_count_original_2 <- prcomp(as.matrix(count_original_2_transpose))

# Check contribution
plot(pca_count_original_2)

# Check each PC
plot(pca_count_original_2$x[,2], pca_count_original_2$x[,3], pch = 16, col = c(rep("red", 5), rep("blue", 5), rep("lightblue", 5), rep("orange", 5)))

# Add label for each dot
text(pca_count_original_2$x[,2], pca_count_original_2$x[,3], labels = rownames(count_original_2_transpose), pos = 4, cex = 0.7, col =c(rep("red", 5), rep("blue", 5), rep("lightblue", 5), rep("orange", 5)))

# Plot # ggplot2 version for count_original_2
pca_count_original_2_ggplot <- data.frame(
  "PC2" = pca_count_original_2$x[,2],
  "PC3" = pca_count_original_2$x[,3],
  "color" = c(rep("red", 5), rep("blue", 5), rep("lightblue", 5), rep("orange", 5)),
  "label" = rownames(count_original_2_transpose)
)

ggplot(pca_count_original_2_ggplot, aes(x = PC2, y = PC3, color = color, label = label)) +
  geom_point(size = 3) +
  geom_text_repel(size = 3) +
  scale_color_identity() +
  labs(x = "PC2", y = "PC3") +
  theme_classic()

# Remove count zero genes over all samples.
count_original_3 <- count_original_2 %>% column_to_rownames("gene_id")
count_original_4 <- count_original_3[rowSums(count_original_3) != 0,]


# Check each gene using histgram.
hist(gene_count$gene_count_over_sample, breaks = 1e9, xlim =c(0, 1e5), ylim = c(0, 300))
hist(gene_count$gene_count_over_sample, breaks = 1e9, xlim =c(0, 1e5), ylim = c(0, 50))


# PCA, for count_original_4

count_original_4_transpose <- count_original_4 %>% t() %>% data.frame()
pca_count_original_4 <- prcomp(as.matrix(count_original_4_transpose), scale = TRUE)
plot(pca_count_original_4)
plot(pca_count_original_4$x[,2], pca_count_original_4$x[,3], pch = 16, col = c(rep("red", 5), rep("blue", 5), rep("lightblue", 5), rep("orange", 5)))
text(pca_count_original_4$x[,2], pca_count_original_4$x[,3], labels = rownames(count_original_2_transpose), pos = 4, cex = 0.7, col =c(rep("red", 5), rep("blue", 5), rep("lightblue", 5), rep("orange", 5)))


# Plot # ggplot2 version for count_original_4 
pca_count_original_4_ggplot <- data.frame(
  "PC2" = pca_count_original_4$x[,2],
  "PC3" = pca_count_original_4$x[,3],
  "color" = c(rep("red", 5), rep("blue", 5), rep("lightblue", 5), rep("orange", 5)),
  "label" = rownames(count_original_4_transpose)
)

ggplot(pca_count_original_4_ggplot, aes(x = PC2, y = PC3, color = color, label = label)) +
  geom_point(size = 3) +
  geom_text_repel(size = 3) +
  scale_color_identity() +
  labs(x = "PC2", y = "PC3") +
  theme_classic()

group <- data.frame(
  "tumor" = rep(x = c("FourT1", "EO771", "FourT1", "EO771"), each = 10),
  "tissue" = rep(x = c("tumor", "blood"), each = 20),
  "drug" = rep(x = c("CTRL", "CPD", "CTRL", "CPD"), each = 5, times = 2)
  )
rownames(group) <- colnames(count_original_4)

group <- group %>% tidyr::unite(c("tumor", "tissue", "drug"), col = "comparison", sep = "_", remove = FALSE) %>% data.frame()

group$comparison <- factor(group$comparison,
                           levels = c("FourT1_tumor_CTRL", "FourT1_tumor_CPD", "EO771_tumor_CTRL", "EO771_tumor_CPD", "FourT1_blood_CTRL", "FourT1_blood_CPD", "EO771_blood_CTRL", "EO771_blood_CPD"),
                           labels = c("FourT1_tumor_CTRL", "FourT1_tumor_CPD", "EO771_tumor_CTRL", "EO771_tumor_CPD", "FourT1_blood_CTRL", "FourT1_blood_CPD", "EO771_blood_CTRL", "EO771_blood_CPD"),
       exclude = NA, ordered = is.ordered(group$comparison), nmax = NA)

# Prepare DGEList object
y_raw <- edgeR::DGEList(count_original_4, group = group$comparison)
# Normalization 
y_raw <- calcNormFactors(y_raw)
# Set comparison pairs 
design <- model.matrix(~0+group, data=y_raw$samples)
# Calculate dispersion 
y_raw <- edgeR::estimateDisp(y_raw)
# check dispersion on graph
plotBCV(y_raw)
# Fitting to NB for analysis of raw count data.
fit_raw <- glmQLFit(y_raw, design = design, robust=FALSE)

# define comparison pair again.
FT1_tumor <- c(-1, 1, 0,0,0,0,0,0)
EO771_tumor <- c(0,0,-1,1,0,0,0,0)
FT1_blood <- c(0,0, 0,0,-1, 1, 0,0)
EO771_blood <- c(0,0,0,0,0,0,-1,1)

# Calculate DEG
qlf_FT1_tumor <- glmQLFTest(fit_raw, contrast = FT1_tumor)
qlf_EO771_tumor <- glmQLFTest(fit_raw, contrast = EO771_tumor)
qlf_FT1_blood <- glmQLFTest(fit_raw, contrast = FT1_blood)
qlf_EO771_blood <- glmQLFTest(fit_raw, contrast = EO771_blood)

# DEG results with FDR correction

# 1
deg_FT1_tumor <- topTags(qlf_FT1_tumor,
        n = nrow(qlf_FT1_tumor$table),
        adjust.method = "BH") %>% data.frame() %>% rownames_to_column(var = "gene_id")

deg_FT1_tumor <- inner_join(deg_FT1_tumor, gene_id_gene_name, by = "gene_id")
deg_FT1_tumor %>% readr::write_tsv(file = "~/output/qlf_FT1_tumor.tsv",  col_names = TRUE)

# 2
deg_EO771_tumor <- topTags(qlf_EO771_tumor,
        n = nrow(qlf_EO771_tumor$table),
        adjust.method = "BH") %>% data.frame() %>% rownames_to_column(var = "gene_id")
deg_EO771_tumor <- inner_join(deg_EO771_tumor, gene_id_gene_name, by = "gene_id")
deg_EO771_tumor %>% readr::write_tsv(file = "~/output/qlf_EO771_tumor",  col_names = TRUE)

# 3
deg_FT1_blood <- topTags(qlf_FT1_blood,
        n = nrow(qlf_FT1_blood$table),
        adjust.method = "BH") %>% data.frame() %>% rownames_to_column(var = "gene_id")
deg_FT1_blood <- inner_join(deg_FT1_blood, gene_id_gene_name, by = "gene_id")
deg_FT1_blood %>% readr::write_tsv(file = "~/output/qlf_FT1_blood.tsv",  col_names = TRUE)

# 4
deg_EO771_blood <- topTags(qlf_EO771_blood,
        n = nrow(qlf_EO771_blood$table),
        adjust.method = "BH") %>% data.frame() %>% rownames_to_column(var = "gene_id")
deg_EO771_blood <- inner_join(deg_EO771_blood, gene_id_gene_name, by = "gene_id")
deg_EO771_blood %>% readr::write_tsv(file = "~/output/qlf_EO771_blood.tsv",  col_names = TRUE)

#5
deg <- list()
deg[[1]] <- deg_FT1_tumor
deg[[2]] <- deg_EO771_tumor
deg[[3]] <- deg_FT1_blood
deg[[4]] <- deg_EO771_blood

deg_FT1_tumor_0.05less <- deg_FT1_tumor %>% filter(FDR < 0.05)
deg_EO771_tumor_0.05less <- deg_EO771_tumor %>% filter(FDR < 0.05)
deg_FT1_blood_0.05less <- deg_FT1_blood %>% filter(FDR < 0.05)
deg_EO771_blood_0.05less <- deg_EO771_blood %>% filter(FDR < 0.05)

deg_0.05less <- list()
deg_0.05less[[1]] <- deg_FT1_tumor_0.05less
deg_0.05less[[2]] <- deg_EO771_tumor_0.05less
deg_0.05less[[3]] <- deg_FT1_blood_0.05less
deg_0.05less[[4]] <- deg_EO771_blood_0.05less

# Output normalized count
normalized_raw_counts_df <- cpm(y_raw, normalized.lib.sizes = TRUE,
log = FALSE) %>% data.frame() %>% rownames_to_column(var = "gene_id") %>% 
  write_tsv("~/output/normalized_raw_counts_df.tsv",
  col_names = TRUE)
  
# Add gene_name
normalized_counts_gene_name <- inner_join(normalized_raw_counts_df, gene_id_gene_name, by = "gene_id")

# Store the dataframe of each comparison group to list
normalized_counts <- list()
normalized_counts[[1]] <- normalized_counts_gene_name[,1:11]
normalized_counts[[2]] <- normalized_counts_gene_name[, c(1, 12:21)]
normalized_counts[[3]] <- normalized_counts_gene_name[, c(1, 22:31)]
normalized_counts[[4]] <- normalized_counts_gene_name[, c(1, 32:41)]

# calculate mean, FC, and label DEGs with red, blue or black.
for(n in 1:length(normalized_counts)){
normalized_counts[[n]] <- normalized_counts[[n]] %>% mutate(
    mean_DMSO = apply(normalized_counts[[n]][,2:6], MARGIN = 1, FUN = mean),
    mean_CPD = apply(normalized_counts[[n]][,7:11], MARGIN = 1, FUN = mean),
    FC = apply(normalized_counts[[n]][,7:11], MARGIN = 1, FUN = mean)/apply(normalized_counts[[n]][,2:6], MARGIN = 1, FUN = mean),
    logFC = log2(apply(normalized_counts[[n]][,7:11], MARGIN = 1, FUN = mean)/apply(normalized_counts[[n]][,2:6], MARGIN = 1, FUN = mean)))
#
normalized_counts[[n]] <- inner_join(normalized_counts[[n]], deg[[n]], by = c("gene_id"), suffix = c("_count", "_deg"))
#
normalized_counts[[n]] <- normalized_counts[[n]] %>% mutate(
  color = case_when(
    logFC_count > log2(1.5) & FDR < 0.01 ~ "red",
    logFC_count < log2(0.75) & FDR < 0.01 ~ "blue",
    TRUE ~ "black"))
}

#
normalized_counts_volcano <- list()
temp_descending <- list()
temp_descending_gene <- list()
temp_ascending <- list()
temp_ascending_gene <- list()
label_volcano <- list()

for (i in 1:4){
  temp_descending[[i]] <- normalized_counts[[i]][-1*log10(normalized_counts[[i]]$FDR) >= -1*log10(0.05) & normalized_counts[[i]]$logFC_deg >= log2(1.5),]
  temp_descending[[i]]$intensity <- -1*log10(temp_descending[[i]]$FDR)*temp_descending[[i]]$logFC_deg
  temp_descending[[i]] <- temp_descending[[i]] %>% arrange(desc(intensity))
  temp_descending_gene[[i]] <- temp_descending[[i]]$gene_name[1:20]
#
  temp_ascending[[i]] <- normalized_counts[[i]][-1*log10(normalized_counts[[i]]$FDR) >= -1*log10(0.05) & normalized_counts[[i]]$logFC_deg <= log2(0.75),]
  temp_ascending[[i]]$intensity <- log10(temp_ascending[[i]]$FDR)*temp_ascending[[i]]$logFC_deg
  temp_ascending[[i]] <- temp_ascending[[i]]  %>% arrange(desc(intensity))
  temp_ascending_gene[[i]] <- temp_ascending[[i]]$gene_name[1:20]
#
  normalized_counts_volcano[[i]] <- normalized_counts[[i]]
  normalized_counts_volcano[[i]]$label_volcano <- rep("", nrow(normalized_counts_volcano[[i]]))
  normalized_counts_volcano[[i]]$label_volcano_check <- rep(NA, nrow(normalized_counts_volcano[[i]]))
#
  label_volcano[[i]] <- c(temp_descending_gene[[i]], temp_ascending_gene[[i]])
  normalized_counts_volcano[[i]]$label_volcano_check <- normalized_counts_volcano[[i]]$gene_name %in% label_volcano[[i]]
#
for(n in 1:nrow(normalized_counts_volcano[[i]])) {
  if(normalized_counts_volcano[[i]]$label_volcano_check[n]){
    normalized_counts_volcano[[i]]$label_volcano[n] <- normalized_counts_volcano[[i]]$gene_name[n]}
  else{normalized_counts_volcano[[i]]$label_volcano[n] <- NA}
  }
}

# plot
par(mfrow = c(1, 4), mar = c(1,2,1,1.5))
for(i in 1:4){
plot(normalized_counts_volcano[[i]]$logFC_deg, -1*log10(normalized_counts_volcano[[i]]$FDR), col = normalized_counts_volcano[[i]]$color, pch = 16, cex = 0.8)
text(
  normalized_counts_volcano[[i]]$logFC_deg, -1*log10(normalized_counts_volcano[[i]]$FDR),
  labels = normalized_counts_volcano[[i]]$label_volcano,
  pos = 4,  
  cex = 0.7, col = "black")
}

volcano <- list()
for(i in 1:4){
volcano[[i]] <- EnhancedVolcano(toptable = normalized_counts[[i]],
                x = "logFC_deg",
                y = "FDR",
                lab = normalized_counts[[i]]$gene_name,
                pCutoff = 0.05,
                FCcutoff = 1.5,
                title = NULL,
                subtitle = NULL,
                caption = NULL,
                xlab = NULL,
                ylab = NULL,
                axisLabSize = 11,
                labSize = 3,
                pointSize = 1) + theme(legend.position = "none")}
#
grid.arrange(grobs = volcano, ncol = 4)

par(mfrow = c(1, 4), mar = c(1,2,1,1.5))
plotSmear(qlf_FT1_tumor, de.tags = deg_FT1_tumor_0.05less$gene_id)
plotSmear(qlf_EO771_tumor, de.tags = deg_EO771_tumor_0.05less$gene_id)
plotSmear(qlf_FT1_blood, de.tags = deg_FT1_blood_0.05less$gene_id)
plotSmear(qlf_EO771_blood, de.tags = deg_EO771_blood_0.05less$gene_id)

par(mfrow = c(1, 4), mar = c(1,2,1,1.5))
for(i in 1:4){
plot(normalized_counts_volcano[[i]]$logCPM, normalized_counts_volcano[[i]]$logFC_deg, col = normalized_counts_volcano[[i]]$color, pch = 16, cex = 0.8)
text(
  normalized_counts_volcano[[i]]$logCPM, normalized_counts_volcano[[i]]$logFC_deg,
  labels = normalized_counts_volcano[[i]]$label_volcano,
  pos = 4,  
  cex = 0.7, col = "black")}

venn_tumor_up <- list()
venn_tumor_down <- list()
venn_blood_up <- list()
venn_blood_down <- list()

venn_tumor_up[1] <- deg_FT1_tumor_0.05less %>% filter(logFC > 0 & FDR < 0.05) %>% select(gene_name) 
venn_tumor_up[2] <- deg_EO771_tumor_0.05less %>% filter(logFC > 0 & FDR < 0.05) %>% select(gene_name)
names(venn_tumor_up) <- c("FT1", "EO771")
common_gene_tumor_up <- intersect(venn_tumor_up[1]$FT1, venn_tumor_up[2]$EO771)

venn_tumor_down[1] <- deg_FT1_tumor_0.05less %>% filter(logFC < 0 & FDR < 0.05) %>% select(gene_name)
venn_tumor_down[2] <- deg_EO771_tumor_0.05less %>% filter(logFC < 0 & FDR < 0.05) %>% select(gene_name)
names(venn_tumor_down) <- c("FT1", "EO771")
common_gene_tumor_down <- intersect(venn_tumor_down[1]$FT1, venn_tumor_down[2]$EO771)


venn_blood_up[1] <- deg_FT1_blood_0.05less %>% filter(logFC > 0 & FDR < 0.05) %>% select(gene_name)
venn_blood_up[2] <- deg_EO771_blood_0.05less %>% filter(logFC > 0 & FDR < 0.05) %>% select(gene_name)
names(venn_blood_up) <- c("FT1", "EO771")
common_gene_blood_up <- intersect(venn_blood_up[1]$FT1, venn_blood_up[2]$EO771)
 
venn_blood_down[1] <- deg_FT1_blood_0.05less %>% filter(logFC < 0 & FDR < 0.05) %>% select(gene_name)
venn_blood_down[2] <- deg_EO771_blood_0.05less %>% filter(logFC < 0 & FDR < 0.05) %>% select(gene_name)
names(venn_blood_down) <- c("FT1", "EO771")
common_gene_blood_down <- intersect(venn_blood_down[1]$FT1, venn_blood_down[2]$EO771)

grid.newpage()
grid.draw(venn.diagram(venn_tumor_up, filename = NULL, main = "Upregulated genes in tumor", fill=c(2,4), scaled=TRUE))
grid.newpage()
grid.draw(venn.diagram(venn_tumor_down, filename = NULL, main = "Downregulated genes in tumor", fill=c(2,4), scaled=FALSE))
grid.newpage()
grid.draw(venn.diagram(venn_blood_up, filename = NULL, main = "Upregulated genes in blood", fill=c(2,4), scaled=FALSE))
grid.newpage()
grid.draw(venn.diagram(venn_blood_down, filename = NULL, main = "Downregulated genes in blood", fill=c(2,4), scaled=FALSE))

# To output to image file
venn.diagram(venn_tumor_up, filename = "~/output/venn_tumor_up.png", main = "Upregulated genes in tumor", fill=c(2,4), scaled=FALSE)
venn.diagram(venn_tumor_down, filename = "~/output/venn_tumor_down.png", main = "Downregulated genes in tumor", fill=c(2,4), scaled=FALSE)
venn.diagram(venn_blood_up, filename = "~/output/venn_blood_up.png", main = "Upregulated genes in blood", fill=c(2,4), scaled=FALSE)
venn.diagram(venn_blood_down, filename = "~/output/venn_blood_down.png", main = "Downregulated genes in blood", fill=c(2,4), scaled=FALSE)

venn <- list()
venn[[1]] <- venn.diagram(venn_tumor_up, filename = NULL, main = "Upregulated genes in tumor", fill=c(2,4), scaled=TRUE)
venn[[2]] <- venn.diagram(venn_tumor_down, filename = NULL, main = "Downregulated genes in tumor", fill=c(2,4), scaled=FALSE)
venn[[3]] <- venn.diagram(venn_blood_up, filename = NULL, main = "Upregulated genes in blood", fill=c(2,4), scaled=FALSE)
venn[[4]] <- venn.diagram(venn_blood_down, filename = NULL, main = "Downregulated genes in blood", fill=c(2,4), scaled=FALSE)

grid.arrange(grobs = venn, ncol = 2, nrow = 2)

grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))  

vp1 <- viewport(layout.pos.row = 1, layout.pos.col = 1)
vp2 <- viewport(layout.pos.row = 1, layout.pos.col = 2)
vp3 <- viewport(layout.pos.row = 2, layout.pos.col = 1)
vp4 <- viewport(layout.pos.row = 2, layout.pos.col = 2)

venn1 <- venn.diagram(venn_tumor_up, filename = NULL, main = "Upregulated genes in tumor", fill = c(2,4), scaled = TRUE)
venn2 <- venn.diagram(venn_tumor_down, filename = NULL, main = "Downregulated genes in tumor", fill = c(2,4), scaled = FALSE)
venn3 <- venn.diagram(venn_blood_up, filename = NULL, main = "Upregulated genes in blood", fill = c(2,4), scaled = FALSE)
venn4 <- venn.diagram(venn_blood_down, filename = NULL, main = "Downregulated genes in blood", fill = c(2,4), scaled = FALSE)

pushViewport(vp1)
grid.draw(venn1)
popViewport()

pushViewport(vp2)
grid.draw(venn2)
popViewport()

pushViewport(vp3)
grid.draw(venn3)
popViewport()

pushViewport(vp4)
grid.draw(venn4)
popViewport()

# These are dataset 

deg_graph <- list()

deg_graph[[1]] <- deg_FT1_tumor_0.05less %>% dplyr::filter(logFC > log2(1.5))
deg_graph[[2]] <- deg_FT1_tumor_0.05less %>% dplyr::filter(logFC < log2(0.75))

deg_graph[[3]] <- deg_EO771_tumor_0.05less %>% dplyr::filter(logFC > log2(1.5))
deg_graph[[4]] <- deg_EO771_tumor_0.05less %>% dplyr::filter(logFC < log2(0.75)) 

deg_graph[[5]] <- deg_FT1_blood_0.05less %>% dplyr::filter(logFC > log2(1.5))
deg_graph[[6]] <- deg_FT1_blood_0.05less %>% dplyr::filter(logFC < log2(0.75)) 

deg_graph[[7]] <- deg_EO771_blood_0.05less %>% dplyr::filter(logFC > log2(1.5))
deg_graph[[8]] <- deg_EO771_blood_0.05less %>% dplyr::filter(logFC < log2(0.75)) 


names(deg_graph) <- c("FT1_tumor_up", "FT1_tumor_down", "EO771_tumor_up", "EO771_tumor_down", "FT1_blood_up", "FT1_blood_down", "EO771_blood_up", "EO771_blood_down")

deg_graph_output <- list()

for(i in 1:length(deg_graph)){
deg_graph_output[[i]] <- enrichGO(gene = deg_graph[[i]]$gene_name,
                          OrgDb = org.Mm.eg.db, 
                          keyType = "SYMBOL", 
                          ont = "BP", 
                          pAdjustMethod = "BH", 
                          pvalueCutoff = 0.05, 
                          qvalueCutoff = 0.05, 
                          readable = TRUE)}

names(deg_graph_output) <- c("FT1_tumor_up", "FT1_tumor_down", "EO771_tumor_up", "EO771_tumor_down", "FT1_blood_up", "FT1_blood_down", "EO771_blood_up", "EO771_blood_down")

for(i in 1:length(deg_graph_output)){
  cnetplot(deg_graph_output[[i]], node_label="all", 
         cex_category = 0.5,
         cex_gene = 0.3,
         cex_label_category = 0.5,
         cex_label_gene = 0.3) %>% plot()
}

# Extract DEG genes from normalized count matrix.
gsea_FT1_tumor <- normalized_raw_counts_df[normalized_raw_counts_df$gene %in% deg_FT1_tumor_0.05less$gene_id == TRUE,1:11]
gsea_EO771_tumor <- normalized_raw_counts_df[normalized_raw_counts_df$gene %in% deg_EO771_tumor_0.05less$gene_id == TRUE,c(1, 12:21)]
gsea_FT1_blood <- normalized_raw_counts_df[normalized_raw_counts_df$gene %in% deg_FT1_blood_0.05less$gene_id == TRUE,c(1, 22:31)]
gsea_EO771_blood <- normalized_raw_counts_df[normalized_raw_counts_df$gene %in% deg_EO771_blood_0.05less$gene_id == TRUE,c(1, 32:41)]

# Add gene name into count matirx
gsea_FT1_tumor <-inner_join(gsea_FT1_tumor, gene_id_gene_name, by = "gene_id")
gsea_EO771_tumor <-inner_join(gsea_EO771_tumor, gene_id_gene_name, by = "gene_id")
gsea_FT1_blood <-inner_join(gsea_FT1_blood, gene_id_gene_name, by = "gene_id")
gsea_EO771_blood <-inner_join(gsea_EO771_blood, gene_id_gene_name, by = "gene_id")

# select required colomns.
gsea_FT1_tumor <- gsea_FT1_tumor %>%  select(all_of(c("gene_name", "gene_id", "T1", "T2", "T3", "T4", "T5", "T6", "T7", "T8", "T9", "T10")))
gsea_EO771_tumor <- gsea_EO771_tumor %>%  select(all_of(c("gene_name", "gene_id","T11", "T12", "T13", "T14", "T15", "T16", "T17", "T18", "T19", "T20")))
gsea_FT1_blood <- gsea_FT1_blood %>% select(all_of(c("gene_name", "gene_id","B1", "B2", "B3", "B4", "B5", "B6", "B7", "B8", "B9", "B10")))
gsea_EO771_blood <- gsea_EO771_blood %>% select(all_of(c("gene_name", "gene_id","B11", "B12", "B13", "B14", "B15", "B16", "B17", "B18", "B19", "B20")))

file_path <- c("~/output/gsea/gsea_FT1_tumor.gct",
               "~/output/gsea/gsea_EO771_tumor.gct",
               "~/output/gsea/gsea_FT1_blood.gct",
               "~/output/gsea/gsea_EO771_blood.gct")

gsea_data <- list()
gsea_data[[1]] <- gsea_FT1_tumor
gsea_data[[2]] <- gsea_EO771_tumor
gsea_data[[3]] <- gsea_FT1_blood
gsea_data[[4]] <- gsea_EO771_blood

for(i in 1:length(file_path)){
  file.remove(file_path[i])
  cat("#1.2", file =file_path[i], sep = "\t", append = FALSE)
  cat("\n", file =file_path[i], sep = "\t", append = TRUE)
  cat(paste0(nrow(gsea_data[[i]]),"\t", ncol(gsea_data[[i]][,3:ncol(gsea_data[[i]])])), file =file_path[i], sep = "\t", append = TRUE)
  cat("\n", file =file_path[i], sep = "\t", append = TRUE)
  write_tsv(gsea_data[[i]], file =file_path[i], col_names = TRUE, append = TRUE)
}

file.remove("~/output/gsea/DMSO_vs_CPD.cls")
#
cat("10 2 1", file ="~/output/gsea/DMSO_vs_CPD.cls", sep = "\t", append = FALSE)
cat("\n", file ="~/output/gsea/DMSO_vs_CPD.cls",sep = "\t", append = TRUE)
cat("# DMSO FIM001", file ="~/output/gsea/DMSO_vs_CPD.cls", sep = "\t", append = TRUE)
cat("\n", file ="~/output/gsea/DMSO_vs_CPD.cls",sep = "\t", append = TRUE)
cat("DMSO DMSO DMSO DMSO DMSO CPD CPD CPD CPD CPD", file ="~/output/gsea/DMSO_vs_CPD.cls", sep = "\t", append = TRUE)

# FT1_tumor
gsea-cli.sh GSEA -res /home/kats/output/gsea/gsea_FT1_tumor.gct -cls /home/kats/output/gsea/DMSO_vs_CPD.cls#CPD_versus_DMSO -gmx /home/kats/output/gsea/msigdb.v2024.1.Mm.symbols.gmt -collapse No_Collapse -mode Max_probe -norm meandiv -nperm 1000 -permute gene_set -rnd_seed 500000 -rnd_type no_balance -scoring_scheme weighted -rpt_label FT1_tumor -metric Signal2Noise -sort real -order descending -create_gcts true -create_svgs true -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 20 -save_rnd_lists false -set_max 500 -set_min 3 -zip_report true -out /home/kats/output/gsea/result

# EO771_tumor
gsea-cli.sh GSEA -res /home/kats/output/gsea/gsea_EO771_tumor.gct -cls /home/kats/output/gsea/DMSO_vs_CPD.cls#CPD_versus_DMSO -gmx /home/kats/output/gsea/msigdb.v2024.1.Mm.symbols.gmt -collapse No_Collapse -mode Max_probe -norm meandiv -nperm 1000 -permute gene_set -rnd_seed 500000 -rnd_type no_balance -scoring_scheme weighted -rpt_label EO771_tumor -metric Signal2Noise -sort real -order descending -create_gcts true -create_svgs true -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 20 -save_rnd_lists false -set_max 500 -set_min 3 -zip_report true -out /home/kats/output/gsea/result

# FT1_blood
gsea-cli.sh GSEA -res /home/kats/output/gsea/gsea_FT1_blood.gct -cls /home/kats/output/gsea/DMSO_vs_CPD.cls#CPD_versus_DMSO -gmx /home/kats/output/gsea/msigdb.v2024.1.Mm.symbols.gmt -collapse No_Collapse -mode Max_probe -norm meandiv -nperm 1000 -permute gene_set -rnd_seed 500000 -rnd_type no_balance -scoring_scheme weighted -rpt_label FT1_blood -metric Signal2Noise -sort real -order descending -create_gcts true -create_svgs true -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 20 -save_rnd_lists false -set_max 500 -set_min 3 -zip_report true -out /home/kats/output/gsea/result

# EO771_blood
gsea-cli.sh GSEA -res /home/kats/output/gsea/gsea_EO771_blood.gct -cls /home/kats/output/gsea/DMSO_vs_CPD.cls#CPD_versus_DMSO -gmx /home/kats/output/gsea/msigdb.v2024.1.Mm.symbols.gmt -collapse No_Collapse -mode Max_probe -norm meandiv -nperm 1000 -permute gene_set -rnd_seed 500000 -rnd_type no_balance -scoring_scheme weighted -rpt_label EO771_blood -metric Signal2Noise -sort real -order descending -create_gcts true -create_svgs true -include_only_symbols true -make_sets true -median false -num 100 -plot_top_x 20 -save_rnd_lists false -set_max 500 -set_min 3 -zip_report true -out /home/kats/output/gsea/result

sudo apt update
sudo apt install openjdk-17-jre
sudo sh '/home/kats/Cytoscape_3_10_3_unix.sh'

# FT1_blood

## CPD 
em_FT1_blood <- read_csv("~/output/output_cytoscape/FT1_blood.csv")

colnames(em_FT1_blood)

em_FT1_blood_CPD <- em_FT1_blood %>% dplyr::filter(
  em_FT1_blood$`EnrichmentMap::fdr_qvalue (Dataset 1)` < 0.05 &
  em_FT1_blood$`EnrichmentMap::Colouring (Dataset 1)` > 0 )


em_FT1_blood_CPD_2 <- em_FT1_blood_CPD %>% dplyr::select(`EnrichmentMap::Genes`) %>% dplyr::rename("gene" =`EnrichmentMap::Genes`)

gene_FT1_blood_CPD <- 1 # initnalize dataframe

for(i in 1:nrow(em_FT1_blood_CPD_2))
{
  gene_FT1_blood_CPD[i] <- em_FT1_blood_CPD_2$gene[i] %>% strsplit("\\|") 
}

gene_FT1_blood_CPD <- unlist(gene_FT1_blood_CPD)

leadingedge_FT1_blood_CPD <- table(gene_FT1_blood_CPD) %>% data.frame() %>% rename("leading_edge_gene" = "gene_FT1_blood_CPD")

leadingedge_FT1_blood_CPD <- leadingedge_FT1_blood_CPD %>% arrange(desc(Freq))

leadingedge_FT1_blood_CPD$FIM001_or_DMSO <- "CPD"


## DMSO

em_FT1_blood_DMSO <- em_FT1_blood %>% dplyr::filter(
  em_FT1_blood$`EnrichmentMap::fdr_qvalue (Dataset 1)` < 0.05 &
  em_FT1_blood$`EnrichmentMap::Colouring (Dataset 1)` < 0 )


em_FT1_blood_DMSO_2 <- em_FT1_blood_DMSO %>% dplyr::select(`EnrichmentMap::Genes`) %>% dplyr::rename("gene" =`EnrichmentMap::Genes`)

gene_FT1_blood_DMSO <- 1 # initnalize dataframe

for(i in 1:nrow(em_FT1_blood_DMSO_2))
{
  gene_FT1_blood_DMSO[i] <- em_FT1_blood_DMSO_2$gene[i] %>% strsplit("\\|") 
}

gene_FT1_blood_DMSO <- unlist(gene_FT1_blood_DMSO)

leadingedge_FT1_blood_DMSO <- table(gene_FT1_blood_DMSO) %>% data.frame() %>% rename("leading_edge_gene" = "gene_FT1_blood_DMSO")

leadingedge_FT1_blood_DMSO <- leadingedge_FT1_blood_DMSO %>% arrange(desc(Freq))

leadingedge_FT1_blood_DMSO$FIM001_or_DMSO <- "DMSO"

# EO771_blood

## CPD 
em_EO771_blood <- read_csv("~/output/output_cytoscape/EO771_blood.csv")

colnames(em_EO771_blood)

em_EO771_blood_CPD <- em_EO771_blood %>% dplyr::filter(
  em_EO771_blood$`EnrichmentMap::fdr_qvalue (Dataset 1)` < 0.05 &
  em_EO771_blood$`EnrichmentMap::Colouring (Dataset 1)` > 0 )


em_EO771_blood_CPD_2 <- em_EO771_blood_CPD %>% dplyr::select(`EnrichmentMap::Genes`) %>% dplyr::rename("gene" =`EnrichmentMap::Genes`)

gene_EO771_blood_CPD <- 1 # initnalize dataframe

for(i in 1:nrow(em_EO771_blood_CPD_2))
{
  gene_EO771_blood_CPD[i] <- em_EO771_blood_CPD_2$gene[i] %>% strsplit("\\|") 
}

gene_EO771_blood_CPD <- unlist(gene_EO771_blood_CPD)

leadingedge_EO771_blood_CPD <- table(gene_EO771_blood_CPD) %>% data.frame() %>% rename("leading_edge_gene" = "gene_EO771_blood_CPD")

leadingedge_EO771_blood_CPD <- leadingedge_EO771_blood_CPD %>% arrange(desc(Freq))

leadingedge_EO771_blood_CPD$FIM001_or_DMSO <- "CPD"


## DMSO

em_EO771_blood_DMSO <- em_EO771_blood %>% dplyr::filter(
  em_EO771_blood$`EnrichmentMap::fdr_qvalue (Dataset 1)` < 0.05 &
  em_EO771_blood$`EnrichmentMap::Colouring (Dataset 1)` < 0 )


em_EO771_blood_DMSO_2 <- em_EO771_blood_DMSO %>% dplyr::select(`EnrichmentMap::Genes`) %>% dplyr::rename("gene" =`EnrichmentMap::Genes`)

gene_EO771_blood_DMSO <- 1 # initnalize dataframe

for(i in 1:nrow(em_EO771_blood_DMSO_2))
{
  gene_EO771_blood_DMSO[i] <- em_EO771_blood_DMSO_2$gene[i] %>% strsplit("\\|") 
}

gene_EO771_blood_DMSO <- unlist(gene_EO771_blood_DMSO)

leadingedge_EO771_blood_DMSO <- table(gene_EO771_blood_DMSO) %>% data.frame() %>% rename("leading_edge_gene" = "gene_EO771_blood_DMSO")

leadingedge_EO771_blood_DMSO <- leadingedge_EO771_blood_DMSO %>% arrange(desc(Freq))

leadingedge_EO771_blood_DMSO$FIM001_or_DMSO <- "DMSO"

leadingedge_FT1_blood_CPD$leading_edge_gene <- str_to_title(leadingedge_FT1_blood_CPD$leading_edge_gene)
leadingedge_EO771_blood_CPD$leading_edge_gene <- str_to_title(leadingedge_EO771_blood_CPD$leading_edge_gene)

leadingedge_FT1_blood_DMSO$leading_edge_gene <- str_to_title(leadingedge_FT1_blood_DMSO$leading_edge_gene)
leadingedge_EO771_blood_DMSO$leading_edge_gene <- str_to_title(leadingedge_EO771_blood_DMSO$leading_edge_gene)

venn_leadingedge_blood_CPD<- list()
venn_leadingedge_blood_DMSO<- list()

venn_leadingedge_blood_CPD[1] <- leadingedge_FT1_blood_CPD %>% select(leading_edge_gene) 
venn_leadingedge_blood_CPD[2] <- leadingedge_EO771_blood_CPD %>% select(leading_edge_gene) 
names(venn_leadingedge_blood_CPD) <- c("FT1", "EO771")
common_gene_leadingedge_blood_CPD <- intersect(venn_leadingedge_blood_CPD[1]$FT1, venn_leadingedge_blood_CPD[2]$EO771)

grid.newpage()
grid.draw(venn.diagram(venn_leadingedge_blood_CPD, filename = NULL, main = "Leading edge gene in blood with CPD", fill=c(2,4), scaled=FALSE))

venn.diagram(venn_leadingedge_blood_CPD, filename = "~/output/output_venndiagram/venn_leading_edge_blood_CPD.png", main = "Leading edge gene in blood with CPD", fill=c(2,4), scaled=FALSE)

venn_leadingedge_blood_DMSO[1] <- leadingedge_FT1_blood_DMSO %>% select(leading_edge_gene) 
venn_leadingedge_blood_DMSO[2] <- leadingedge_EO771_blood_DMSO %>% select(leading_edge_gene) 
names(venn_leadingedge_blood_DMSO) <- c("FT1", "EO771")
common_gene_leadingedge_blood_DMSO <- intersect(venn_leadingedge_blood_DMSO[1]$FT1, venn_leadingedge_blood_DMSO[2]$EO771)


grid.newpage()
grid.draw(venn.diagram(venn_leadingedge_blood_DMSO, filename = NULL, main = "Leading edge gene in blood with DMSO", fill=c(2,4), scaled=FALSE))
venn.diagram(venn_leadingedge_blood_DMSO, filename = "~/output/output_venndiagram/venn_leading_edge_blood_DMSO.png", main = "Leading edge gene in blood with DMSO", fill=c(2,4), scaled=FALSE)

common_gene_leadingedge_blood_CPD %>% data.frame() %>% rename("gene" = ".") %>% write_tsv("~/output/output_leadingedge/common_gene_leadingedge_blood_CPD")
common_gene_leadingedge_blood_DMSO %>% data.frame() %>% rename("gene" = ".") %>% write_tsv("~/output/output_leadingedge/common_gene_leadingedge_blood_DMSO")

common_gene_leadingedge_blood_CPD
heatmap_common_gene_leadingedge_blood_CPD  <- normalized_counts_gene_name[normalized_counts_gene_name$gene_name %in% common_gene_leadingedge_blood_CPD == TRUE, ]
heatmap_common_gene_leadingedge_blood_CPD <- heatmap_common_gene_leadingedge_blood_CPD[,-c(1:21)]

data_pheatmap <- as.matrix(heatmap_common_gene_leadingedge_blood_CPD[,1:20])
rownames(data_pheatmap) <- heatmap_common_gene_leadingedge_blood_CPD$gene_name

pheatmap(data_pheatmap,
         cluster_rows = TRUE,
         cluster_cols = FALSE,
         scale = "row",
         show_rownames = FALSE,
         show_colnames = TRUE,
         clustering_method = "ward.D2",
         cutree_rows = 2,
         legend = FALSE)


pheatmap_value <- pheatmap(data_pheatmap,
                           cluster_rows = TRUE,
                           cluster_cols = FALSE,
                           scale = "row",
                           show_rownames = TRUE,
                           show_colnames = TRUE,
                           clustering_method = "ward.D2",
                           cutree_rows = 2,
                           legend = FALSE)

pheatmap_annotation <- data.frame("gene_name" = pheatmap_value@row_names_param$labels)
pheatmap_annotation$order <- c(1:nrow(pheatmap_annotation))
pheatmap_annotation$gene_name <- factor(x = pheatmap_annotation$gene_name,
                                        levels = pheatmap_annotation$gene_name,
                                        labels = pheatmap_annotation$gene_name,
                                        exclude = NA, ordered = is.ordered(pheatmap_annotation$gene_name), nmax = NA)

pheatmap_annotation$order <- c(1:nrow(pheatmap_annotation))

heatmap_common_gene_leadingedge_blood_CPD_check_gene_name <- inner_join(
  heatmap_common_gene_leadingedge_blood_CPD, pheatmap_annotation, by = "gene_name")
heatmap_common_gene_leadingedge_blood_CPD_check_gene_name <- heatmap_common_gene_leadingedge_blood_CPD_check_gene_name %>% mutate(
  mean_FT1_DMSO = apply(heatmap_common_gene_leadingedge_blood_CPD_check_gene_name[,1:5], MARGIN = 1, FUN = mean),
  mean_FT1_CPD = apply(heatmap_common_gene_leadingedge_blood_CPD_check_gene_name[,6:10], MARGIN = 1, FUN = mean),
  mean_EO771_DMSO = apply(heatmap_common_gene_leadingedge_blood_CPD_check_gene_name[,11:15], MARGIN = 1, FUN = mean),
  mean_EO771_CPD = apply(heatmap_common_gene_leadingedge_blood_CPD_check_gene_name[,16:20], MARGIN = 1, FUN = mean),
  FC_FT1 = apply(heatmap_common_gene_leadingedge_blood_CPD_check_gene_name[,6:10], MARGIN = 1, FUN = mean)/apply(heatmap_common_gene_leadingedge_blood_CPD_check_gene_name[,1:5], MARGIN = 1, FUN = mean),
  FC_EO771 = apply(heatmap_common_gene_leadingedge_blood_CPD_check_gene_name[,16:20], MARGIN = 1, FUN = mean)/apply(heatmap_common_gene_leadingedge_blood_CPD_check_gene_name[,11:15], MARGIN = 1, FUN = mean),
)

heatmap_Up <- heatmap_common_gene_leadingedge_blood_CPD_check_gene_name %>% dplyr::filter(FC_FT1 > 1 & FC_EO771 > 1)
heatmap_Down <- heatmap_common_gene_leadingedge_blood_CPD_check_gene_name %>% dplyr::filter(FC_FT1 < 1 & FC_EO771 < 1)


tiff(filename = "~/output/data_pheatmap.tiff",
     width = 800, height = 800, units = "px", pointsize = 12,
     compression = "none",
     bg = "white", res = 150)

pheatmap(data_pheatmap,
         cluster_rows = TRUE,
         cluster_cols = FALSE,
         scale = "row",
         show_rownames = FALSE,
         show_colnames = TRUE,
         clustering_method = "ward.D2",
         cutree_rows = 2,
         legend = FALSE)

dev.off()