---

title: "NGS_ID_to_PIS" author: "zagor-mju" date: "18 02 2022" output: html_document editor_options: chunk_output_type: console

knitr::opts_chunk$set(echo = TRUE)

library(plyr)
library(dplyr)
library(data.table)

setwd("O:/DEJAVNOSTI/OMIKE/pISA-Projects/_p_SUSPHIRE/_I_T21_SXPsysbio/_S_P4_DiNAR/_A_CKN_NbL35-DiNAR/input/Prioritization_newScript")
pval.cutoff = 0.01

NGS.ID.to.PIS = read.table(file = "input/NODES_long.txt",
                           header = TRUE, 
                           sep = "\t", 
                           quote = NULL,
                           dec = ".", 
                           stringsAsFactors = FALSE,
                           na.strings = "NA",
                           fill = TRUE,
                           comment.char = "#")

NGS.ID.to.PIS = NGS.ID.to.PIS[!duplicated(NGS.ID.to.PIS), ]
colnames(NGS.ID.to.PIS) = c("clade", "probeID")
NGS.ID.to.PIS <- data.table(NGS.ID.to.PIS, key = "probeID") 
dim(NGS.ID.to.PIS)


organism.NGS = list.files("input/NGS_DE/", pattern = ".txt")
organism.NGS.files = NGS.ID.to.PIS

for (i in organism.NGS) {
  tmp = read.table(file = paste0("input/NGS_DE/", i),
                   header = TRUE, 
                   sep = "\t", 
                   quote = NULL,
                   dec = ".", 
                   stringsAsFactors = FALSE,
                   na.strings = "NA",
                   fill = TRUE,
                   comment.char = "#")  
  ind = c(grep("probeID", colnames(tmp)), grep("pAdj", colnames(tmp)), grep("LogFC", colnames(tmp)))
  tmp = tmp[,ind]
  ind = which(tmp[,2] < pval.cutoff)
  tmp = tmp[ind,c(1,3)]
  organism.NGS.files = merge(organism.NGS.files, tmp, by = "probeID", all.x = TRUE, sort = FALSE)
}
dim(organism.NGS.files)

logFC.starts.from.col = 3

NGS.ID.to.PIS[which(!is.na(stringr::str_match(NGS.ID.to.PIS$clade, 'AT3G21510'))),]
NGS.ID.to.PIS[which(!is.na(stringr::str_match(NGS.ID.to.PIS$clade, 'AT5G23660'))),]
NGS.ID.to.PIS[which(!is.na(stringr::str_match(NGS.ID.to.PIS$clade, 'AT3G54700'))),]


organism.NGS.files[which(organism.NGS.files$clade == 'AT3G21510'), ]
organism.NGS.files[which(organism.NGS.files$clade == 'AT5G23660'), ]
organism.NGS.files[which(organism.NGS.files$clade == 'AT3G54700'), ]

lapply(organism.NGS.files, class)

clades = sort(unique(organism.NGS.files$clade))
tmp = NULL

for (i in clades) {
  ind = which(organism.NGS.files$clade == i)
  if (length(ind) > 1) { # (length(ind > 1)) for found, not length(ind) > 1 for multiple (1 or more, doesnt matter here)
    mysubset = organism.NGS.files[ind, ]
    e = sapply(1:nrow(mysubset), function(x) sum(!is.na(mysubset[x, logFC.starts.from.col:ncol(mysubset)])))
    f1 = sapply(1:nrow(mysubset), 
                function(x) ifelse(all(is.na(mysubset[x, logFC.starts.from.col:ncol(mysubset)])), 
                                   0, 
                                   mean(colMeans(abs(mysubset[x, logFC.starts.from.col:ncol(mysubset)]), na.rm = TRUE), na.rm = TRUE)))
    f2 = sapply(1:nrow(mysubset), 
                function(x) ifelse(all(is.na(mysubset[x, logFC.starts.from.col:ncol(mysubset)])), 
                                   0, 
                                   max(abs(mysubset[x, logFC.starts.from.col:ncol(mysubset)]), na.rm = TRUE)))
    g = which(e == max(e))
    h1 = which(f1[g] == max(f1[g]))
    h2 = which(f2[g] == max(f2[g]))
    k2 = ifelse(length(h1) > 1, h2, h1)
    k = ifelse(length(g) > 1, g[k2], g)
    
    tmp = rbind(tmp, mysubset[k,])
  } else { # no match
    tmp = rbind(tmp, organism.NGS.files[ind,]) # take first
  }
}
nrow(tmp) == length(clades)
clades.organism.NGS =  tmp

write.table(x = clades.organism.NGS, file = "output/clades.organism.NGS.priority_MJU_p001.txt", 
            append = FALSE, quote = FALSE, sep = "\t",
            eol = "\n", na = "NA", dec = ".", row.names = FALSE,
            col.names = TRUE)

dir.create("output/organism_NGS_to_PIS_MJU_p001")
getGenes = clades.organism.NGS$probeID
j = !sapply(1:nrow(clades.organism.NGS), 
            function(x) all(is.na(clades.organism.NGS[x, logFC.starts.from.col:ncol(clades.organism.NGS)])))
getGenes = getGenes[j]

organism.NGS = list.files("input/NGS_DE/", pattern = ".txt")
for (i in organism.NGS) {
  tmp = read.table(file = paste0("input/NGS_DE/", i),
                   header = TRUE, 
                   sep = "\t", 
                   quote = NULL,
                   dec = ".", 
                   stringsAsFactors = FALSE,
                   na.strings = "NA",
                   fill = TRUE,
                   comment.char = "#")
  ind = match(getGenes, tmp$probeID)
  temp = merge(NGS.ID.to.PIS, tmp[ind,], by = "probeID", all.y = TRUE, sort = FALSE)
  temp = arrange(temp, clade, probeID)
  write.table(x = temp, file = paste0("output/organism_NGS_to_PIS_MJU_p001/", "PIS_", i), 
              append = FALSE, quote = FALSE, sep = "\t",
              eol = "\n", na = "NA", dec = ".", row.names = FALSE,
              col.names = TRUE)
}