library(ggplot2)
library(ggrepel)
library(gprofiler2)
library(stringr)
#install.packages("readr")   #if needed
library(readr)

#########################################
#
#    Required files:
#        1 - Biogrid file is very big, download it from ->
#            https://downloads.thebiogrid.org/File/BioGRID/Release-Archive/BIOGRID-5.0.257/BIOGRID-ORGANISM-5.0.257.tab3.zip
#            After extracting you can remve all except BIOGRID-ORGANISM-Homo_sapiens-5.0.257.tab3.txt
#        2 - DESeq2_results_basal_vs_classical.csv from -> 
#            output of lab 2a
#        3 - ensg_gprofiler_results_foldchange0_ordered.tsv from ->
#            output of lab 2b
#
#########################################




DEG_table = read.table(file.path("outputs","DESeq2_results_basal_vs_classical.csv"), sep=',',  header = TRUE)

DEG_table$entrez_id = str_split_i(DEG_table$gene, '\\|', 2)  # get the second column. | would mean OR, so we search for "slash"|
DEG_table$gene_name = str_split_i(DEG_table$gene, '\\|', 1)  # get the second column. | would mean OR, so we search for "slash"|

#sort by p-value (usually it it better than sorting by FDR, but right now it wouldn't matter)
DEG_table = DEG_table[order(DEG_table$pvalue),]
DEG_table = DEG_table[1:500,]  # take top500 as an example.

entrez_list = DEG_table$entrez_id


result_list = list()
chunk_num = 1
chunk_size = 1000  # read only 1000 lines at a time.

# Open a file connection
con = file(file.path("data","BIOGRID-ORGANISM-Homo_sapiens-5.0.257.tab3.txt"), "r")
print(con)

# Read header if present
chunk = readLines(con, n = 1)  #get onle the header first
header = chunk

while(length(chunk) != 0){
    chunk = read.table(con, nrows = chunk_size, sep = "\t", header = FALSE, quote = "", comment.char = "")
    colnames(chunk) = str_split_1(header, '\\t')
    filtered_chunk = chunk[((chunk$"Entrez Gene Interactor A" %in% entrez_list) & (chunk$"Entrez Gene Interactor B" %in% entrez_list)),]
    if(nrow(filtered_chunk) > 0) {
        result_list[[chunk_num]] = filtered_chunk
        chunk_num = chunk_num + 1
    }
}

# Combine results
final_result = do.call(rbind, result_list)
colnames(final_result) = str_split_1(header, '\\t')
final_result <- final_result[final_result$"Entrez Gene Interactor A" != final_result$"Entrez Gene Interactor B",]

final_result_sorted <- data.frame(
  'Entrez Gene Interactor A' = pmin(final_result[,1], final_result[,2]),
  'Entrez Gene Interactor B' = pmax(final_result[,1], final_result[,2])
)

# Then remove duplicates
final_result_sorted <- unique(final_result_sorted)

write.table(final_result_sorted, file=file.path("outputs","mini_biogrid.tsv"), quote=FALSE, sep='\t', row.names = FALSE)

###################
#  Create example for enrichmentmap

pathway_table = read.table(file.path("outputs","ensg_gprofiler_results_foldchange0_ordered.tsv"), sep='\t',  header = TRUE)

custom_GMT = pathway_table[, c("term_id", "term_name", "intersection")]
custom_GMT$intersection = gsub(',', '\t', custom_GMT$intersection)

write.table(custom_GMT, file=file.path("outputs","mini_gmt.gmt"), quote=FALSE, sep='\t', row.names = FALSE)

