# from "test_exp_raw" file.
# first - identify those problematic samples and delete

source("quant_assist.R")
## note: higher level functions are depreciated here,
## since being straight is more than being neat.

# load raw file SRRX*SYM
exp_raw <- read_rds('test_exp_raw.rds')
m_raw <- read_rds('test_m_raw.rds')

# load the smaple dictionary
homo <- read_rds("srr_withgtfs.rds")
h <- homo %>% dplyr::select(wecall, sra_acc)
h$source_abbr <- mapply(function(a,b) gsub(paste0("_",a), "", b), a = h$sra_acc, b = h$wecall)

# split into three matrices
m_ip_raw <- m_raw[h$sra_acc[grepl("_ip_", h$wecall)],]
m_input_raw <- m_raw[h$sra_acc[grepl("_input_", h$wecall)],]

# identify the samples to be deleted, e.g. humanBrain
srrs_to_delete <- c('SRR456936','SRR456937')
exp_del <- delete.from.matrix(exp_raw, srrs_to_delete)
m_ip_del <- delete.from.matrix(m_ip_raw, srrs_to_delete)
m_input_del <- delete.from.matrix(m_input_raw, srrs_to_delete)

# for display and check by eyes (only differ by 'BR' ids, presumably bio-reps)
# check_abbr(input_rep_srrs)
# check_abbr(ip_rep_srrs)

exp_merge <- merge_reps(input_rep_srrs, exp_del, "_input")
ip_merge <- merge_reps(ip_rep_srrs, m_ip_del, "_ip")
input_merge <- merge_reps(input_rep_srrs, m_input_del, "_input")

# calculate the methylation level
meth_merge <- log2((ip_merge + 0.01) / (input_merge + 0.01))
Note <- "here the three stored matrices have gone thorough stringtie quantification, reading from gtf files and merging the techreps and bioreps. However these data has not been normalised. Start your journey here!"
save(exp_merge, ip_merge, input_merge, Note, file = "raw_merge.RData")


### ---------------- ######
# remove the entire workspace except three matrices
#rm(ls()[!ls() %in% c("meth_merge", "exp_merge")])

# filter by mean/sd quantile, and normalize by experiments
meth_norm <- norm_merge(meth_merge, 0.5, 0.25) # dim = 38*31705
exp_norm <- norm_merge(exp_merge, 0.75, 0.75) # dim = 38*16356

# merge the sites by expression correlation and physical location


g <- merge_sites(m = meth_norm, mdict0 = refer_18w, distance_cutoff = 200, cor_cutoff = 0.7)
write_rds(g, "try-graph.rds")
#


# find the clusters
g <- read_rds("try-graph.rds")

# make meth_exp
g_list <- cluster_fast_greedy(g) %>% communities()
single_groupings <- unlist(g_list[lengths(g_list) == 1]) %>% unname()
multi_groupings <- g_list[! lengths(g_list) == 1]
cs <- comb_site(multi_groupings, meth_norm)
ss <- meth_norm[,colnames(meth_norm) %in% single_groupings]
meth_exp <- cbind(cs, ss) %>% reorder


# make network
met_exp_adj <- adjmake(x = meth_exp, y = exp_norm, quant = 0.05, pcut = 0.01)



# characterize network


# play with test_exp and test_m
nmf(meth_exp + 18, rank = 3) -> res