# Interaction analysis

library(DESeq2)
library(ggplot2)
library(openxlsx)
library(pheatmap)
library(RColorBrewer)
library(clusterProfiler)
library(DOSE)          
library(AnnotationHub)
library(org.Hs.eg.db)


out_dir <- "../00-Interaction/"

read.counts <- read.table(file = "../04-counts/featureCounts_results_corrected.txt", header = T)
rownames(read.counts) <- read.counts$Geneid
read.counts <- read.counts[,c(7:ncol(read.counts))]
colnames(read.counts) <-  unlist(strsplit(x = colnames(read.counts),split = "\\.", ))[seq.int(2, 42, by = 3)]
sample_info <- read.table(file = "SampleInfo.txt", header = T)
sample_info$OFP <- substr(x = sample_info$Group, start = 1, stop = 4)
sample_info$Treatment <- substr(x = sample_info$Group, start = 5, stop = length(sample_info$Group))
rownames(sample_info) <- sample_info$SampleID

dds <- DESeqDataSetFromMatrix(countData = read.counts,
                              colData = sample_info,
                              design = as.formula("~ OFP + Treatment + OFP:Treatment"))

keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]

dds$OFP <- relevel(dds$OFP, ref = "OFPn")
dds$Treatment <- relevel(dds$Treatment, ref = "NT")

dds <- DESeq(dds)
sizeFactors(dds)
resultsNames(dds)

rlog <- rlog(dds)

for (contr in c("OFPOFPp.TreatmentT")) {
  res <- results(dds, name = contr)
  res <- res[!is.na(res$padj), ]
  res$padj <- as.numeric(res$padj)
  res <- res[res$padj < 0.05 & res$baseMean > 5 & abs(res$log2FoldChange) > 1, ]
  write.table(x = as.data.frame(res), file = paste(out_dir, paste0("DESeq2_", contr, "_interaction.tsv"), sep = "/"),
              sep = "\t", quote = FALSE, col.names = NA)
  res$Gene <- rownames(res)
  write.xlsx(x = res, file = paste(out_dir, paste0("DESeq2_", contr, "_interaction.xlsx"), sep = "/"))
}


sample_info$SampleID <- NULL

pheatmap(mat = as.matrix(assay(rlog)[res[res$log2FoldChange > 0,]$Gene,]), 
         color = rev(brewer.pal(9, "RdBu")), 
         scale = "row", annotation_col = sample_info, 
         show_rownames = F, cluster_cols = F)

saveRDS(object = list(dds = dds, 
                      read.counts = read.counts, 
                      res = res, 
                      rlog = rlog,
                      sample_info = sample_info), file = "00-Interaction.rds")

### Performing downstream analysis


interaction_genes <- list()

interaction_genes[["up"]] <- subset(res, log2FoldChange > 0)
interaction_genes[["down"]] <- subset(res, log2FoldChange < 0)

gene_list <- list()
gene_list[["up"]] <- interaction_genes$up$log2FoldChange
names(gene_list[["up"]]) <- interaction_genes$up$Gene
gene_list[["down"]] <- interaction_genes$down$log2FoldChange
names(gene_list[["down"]]) <- interaction_genes$down$Gene


contr.enrich_go_up <- enrichGO(gene = names(gene_list[["up"]]),
                            universe = rownames(dds),
                            OrgDb = 'org.Hs.eg.db',
                            ont = "BP",
                            pAdjustMethod = "BH",
                            pvalueCutoff = 0.05,
                            keyType = 'SYMBOL')

contr.enrich_go_down <- enrichGO(gene = names(gene_list[["down"]]),
                               universe = rownames(dds),
                               OrgDb = 'org.Hs.eg.db',
                               ont = "BP",
                               pAdjustMethod = "BH",
                               pvalueCutoff = 0.05,
                               keyType = 'SYMBOL')

contr.enrich_go <- list(UP = contr.enrich_go_up@result, 
                        DOWN = contr.enrich_go_down@result)

write.xlsx(contr.enrich_go,file = "GO_BP_Enrichment_interaction_OFP_Treatment.xlsx", 
           asTable = T, firstRow = T, colWidths = "auto") 

# GSEA 


for (contr in c("OFPOFPp.TreatmentT")) {
  res_all <- results(dds, name = contr)
  res_all <- res_all[!is.na(res_all$padj), ]
  res_all$padj <- as.numeric(res_all$padj)
 # res_all <- res[res$padj < 0.05 & res$baseMean > 5 & abs(res$log2FoldChange) > 1, ]
  write.table(x = as.data.frame(res_all), file = paste(out_dir, paste0("DESeq2_", contr, "_interaction_full.tsv"), sep = "/"),
              sep = "\t", quote = FALSE, col.names = NA)
  res_all$Gene <- rownames(res_all)
  write.xlsx(x = res_all, file = paste(out_dir, paste0("DESeq2_", contr, "_interaction_full.xlsx"), sep = "/"))
}

gene_list[["full"]] <- res_all$log2FoldChange
names(gene_list[["full"]]) <- res_all$Gene
gene_list[["full"]] <- sort(gene_list[["full"]], decreasing = T)

t2gene <- read.gmt(gmtfile = paste0("/opt/common/tools/ric.tiget/GSEA/GSEA_human_v7.4/h.all.v7.4.symbols.gmt"))

gsea_res <- GSEA(gene_list[["full"]], TERM2GENE = t2gene,
                                      verbose = FALSE,minGSSize = 10,maxGSSize = 1000,
                                      pvalueCutoff = 0.05)

gseaplot2(gsea_res, geneSetID = 1:3, pvalue_table = TRUE,
          color = c("#E495A5", "#86B875", "#7DB0DD"), ES_geom = "dot")

a <- ridgeplot(gsea_res, fill = "NES", showCategory = 50) + scale_fill_gradient2() 
a$theme$axis.text.y$size <- 10
a$theme$axis.text.x$size <- 10

png(filename = "GSEA_ridges_enrich_core_distribution.png",width = 9, height = 6, units = "in", res = 200)
a
dev.off()