#!/usr/bin/env Rscript

# Dénes Türei EMBL 2017
# turei.denes@gmail.com

require(igraph)
require(qgraph)
require(dplyr)
require(stringi)

make_fr_layout <- function(g){
    # layout with qgraph
    # g is an igraph object
    el <- get.edgelist(g, names = FALSE)
    w  <- ifelse('weight' %in% names(edge.attributes(g)), E(g)$weight, NULL)
    w  <- 1 / w
    w  <- w**4
    lo <- qgraph.layout.fruchtermanreingold(el,
                                            weights = w,
                                            vcount = vcount(g),
                                            area = vcount(g)^2.3,
                                            repulse.rad = vcount(g)^0.3,
                                            niter = 30000)
    return(lo)
}

###

infile        <- '170829_charlotte_dataset.csv'
outfile       <- 'lima-network.pdf'
ewidth_min    <- .5
ewidth_max    <- 3.0

###

data <- read.table(infile, sep = ',', dec = '.', header = TRUE)

data15 <- data %>%
    mutate(
        prop = as.numeric(stri_sub(membrane_ID, -2)),
        base = sapply(strsplit(as.character(membrane_ID), '_'), '[', 1),
        spec = gsub('^[^_]*_', '', gsub('_[^_]*$', '', membrane_ID))
    ) %>%
    filter(log10fc > 1.5 & (prop == 5 | is.na(prop))) %>%
    group_by(protein, membrane_ID) %>%
    mutate(mlog10fc = mean(log10fc)) %>%
    summarize_all(first)

# see how many connections we have:
# dim(data15)
# 114 20

proteins  <- unique(as.character(data15$protein))
membranes <- unique(as.character(data15$membrane_ID))

edgelist <- data15 %>%
    select(protein, membrane_ID, mlog10fc)

g <- graph_from_data_frame(edgelist, directed = TRUE)
V(g)$type <- ifelse(V(g)$name %in% proteins, 'protein', 'membrane')
V(g)$base <- ifelse(V(g)$type == 'membrane',
                    sapply(strsplit(as.character(V(g)$name), '_'), '[', 1),
                    '')
V(g)$spec <- ifelse(V(g)$type == 'membrane',
                    gsub('^[^_]*_', '', gsub('_[^_]*$', '', V(g)$name)),
                    '')
V(g)$label <- ifelse(V(g)$base == V(g)$spec,
                     V(g)$name,
                     paste(V(g)$base, V(g)$spec, sep = '_'))
V(g)$color <- ifelse(V(g)$type == 'protein',
                     '#007B7FAA', # embl dark green
                     '#6EA945AA'  # embl light green
                     )
E(g)$width <- (E(g)$mlog10fc - min(E(g)$mlog10fc)) /
              max(E(g)$mlog10fc - min(E(g)$mlog10fc)) *
              (ewidth_max - ewidth_min) +
              ewidth_min
E(g)$weight <- E(g)$mlog10fc
lo <- make_fr_layout(g)


cairo_pdf(filename = outfile, onefile = TRUE, width = 12, height = 12, pointsize = 12, family = 'DINPro')
    
    plot(g,
         layout = lo,
         edge.alpha = .7,
         edge.width = E(g)$width,
         vertex.alpha = .7,
         vertex.frame.color = NA,
         vertex.label.family = 'DINPro',
         vertex.label.color = '#454447'
    )
    
dev.off()