#!/usr/bin/env Rscript

# Copyright Saez Lab 2020
#
# Author: Denes Turei
# Contact: turei.denes@gmail.com
#
# License: MIT License https://opensource.org/licenses/MIT
#
# Combines the OmniPath PPI and transcriptional regulatory network
# with the human-SARS CoV-2 interactions from Gordon et al. 2020

require(dplyr)
require(tibble)
require(purrr)
require(openxlsx)
require(igraph)
require(qgraph)

if(!'OmnipathR' %in% installed.packages()[,"Package"]){
    require(devtools)
    install_github('saezlab/OmnipathR')
}

require(OmnipathR)

# URLs for supplementary tables S2 and S3 from Gordon et al. 2020
# https://www.biorxiv.org/content/10.1101/2020.03.22.002386v2
url_s2 <- paste0(
    'https://www.biorxiv.org/content/biorxiv/early/',
    '2020/03/23/2020.03.22.002386/DC4/embed/media-4.xlsx?download=true'
)
url_s3 <- paste0(
    'https://www.biorxiv.org/content/biorxiv/early/',
    '2020/03/23/2020.03.22.002386/DC5/embed/media-5.xlsx?download=true'
)

# the OmniPath PPI interaction network
ia_omnipath <- import_Omnipath_Interactions() %>% as_tibble()

# in the example we don't use these but if you need higher coverage on PPI
# add these with `bind_rows` below
ia_ligrec <- import_LigrecExtra_Interactions() %>% as_tibble()
ia_pwextra <- import_PathwayExtra_Interactions() %>% as_tibble()
ia_kinaseextra <- import_KinaseExtra_Interactions() %>% as_tibble()

# transcriptional regulation
# if you want more interactions, pass the argument
# `confidence_level = c('A', 'B', 'C', 'D')`
ia_transcriptional <- import_TFregulons_Interactions() %>% as_tibble()

# post-transcriptional regulation
# here we don't use it but if you are interested
# add these by `bind_rows` below
ia_post_transcriptional <- import_miRNAtarget_Interactions() %>% as_tibble()

# downloading the host-pathogen interactions from the paper
human_sarscov2_raw <- read.xlsx(url_s2, startRow = 2) %>% as_tibble()
# downloading the compound-host-pathogen interactions from the paper
# here we don't use it but if you are interested combine with the other
# interactions a similar way as the host-pathogen ones
human_sarscov2_compounds_raw <- read.xlsx(url_s3) %>% as_tibble()

# combining all the interactions to one data frame
ia_human_sarscov2 <- human_sarscov2_raw %>%
    select(
        source = Bait,
        target = Preys,
        source_genesymbol = Bait,
        target_genesymbol = PreyGene,
        MIST,
        Saint_BFDR,
        AvgSpec,
        FoldChange
    ) %>%
    # manually adding the known S--ACE2 and S--BSG interactions
    add_row(
        source = 'SARS-CoV2 Spike',
        target = 'Q9BYF1',
        source_genesymbol = 'Spike',
        target_genesymbol = 'ACE2'
    ) %>%
    add_row(
        source = 'SARS-CoV2 Spike',
        target = 'P35613',
        source_genesymbol = 'Spike',
        target_genesymbol = 'BSG'
    )

interactions <- as_tibble(
    bind_rows(
        ia_omnipath %>% mutate(type = 'ppi'),
        ia_pwextra %>% mutate(type = 'ppi'),
        ia_kinaseextra %>% mutate(type = 'ppi'),
        ia_ligrec %>% mutate(type = 'ppi'),
        ia_transcriptional %>% mutate(type = 'transcriptional'),
        ia_human_sarscov2 %>% mutate(
            source_genesymbol = sub('SARS-CoV2 ', '', source_genesymbol),
            type = 'host_pathogen',
            # the direction or effect of these interactions is unknown
            is_directed = 0,
            is_stimulation = 0,
            is_inhibition = 0,
            consensus_direction = 0,
            consensus_stimulation = 0,
            consensus_inhibition = 0,
        )
    )
)

# creating an igraph network from the interactions data frame
net_human_sarscov2 <- interaction_graph(interactions = interactions) %>%
    simplify(remove.multiple = FALSE, remove.loops = TRUE)

# labeling the proteins by organism
sarscov2_proteins <- ia_human_sarscov2 %>% pull(source) %>% unique()

V(net_human_sarscov2)$organism <- ifelse(
    V(net_human_sarscov2)$up_ids %in% sarscov2_proteins,
    'SARS-CoV2',
    'human'
)

# labeling direct viral targets and transcription factors
V(net_human_sarscov2)$viral_target <- (
    V(net_human_sarscov2)$up_ids %in% (
        ia_human_sarscov2 %>% pull(target) %>% unique()
    )
)

V(net_human_sarscov2)$tf <- (
    V(net_human_sarscov2)$up_ids %in% (
        ia_transcriptional %>% pull(source) %>% unique()
    )
)

# node indices of the SARS-CoV2 proteins
nodes_sarscov2 <- which(V(net_human_sarscov2)$organism == 'SARS-CoV2')


# creating subgraphs of the neighborhood of each virus protein
neighborhoods_human_sarscov2_1 <- make_ego_graph(
    net_human_sarscov2,
    order = 1,
    nodes = nodes_sarscov2,
    mode = 'all'
)

nodes_in_neighborhoods <- neighborhoods_human_sarscov2_1 %>%
    map(function(g){V(g)$name}) %>%
    unlist() %>%
    unique()

# reducing the network to the neighborhoods
net_human_sarscov2_1 <- net_human_sarscov2 %>%
    delete_vertices(
        which(
            !V(net_human_sarscov2)$name %in% nodes_in_neighborhoods
        )
    )

net_human_sarscov2_1 <- net_human_sarscov2_1 %>%
    delete_vertices(
        which(!(
            V(net_human_sarscov2_1)$viral_target |
            V(net_human_sarscov2_1)$organism == 'SARS-CoV2'
        ))
    )

nodes_having_ppi <- net_human_sarscov2_1 %>%
    get.edges(
        E(net_human_sarscov2_1)[
            E(net_human_sarscov2_1)$type %in% c('ppi', 'transcriptional')
        ]
    ) %>%
    c() %>%
    unique()

sarscov2_nodes <- which(V(net_human_sarscov2_1)$organism == 'SARS-CoV2')

nodes_to_keep <- c(nodes_having_ppi, sarscov2_nodes) %>% unique()

net_human_sarscov2_1 <- net_human_sarscov2_1 %>%
    induced_subgraph(nodes_to_keep)

net_human_sarscov2_1 <- net_human_sarscov2_1 %>%
    delete_vertices(degree(.) == 0)

# creating a figure of the direct neighborhood network
fr_layout <- qgraph.layout.fruchtermanreingold(
    get.edgelist(net_human_sarscov2_1, names = FALSE),
    vcount = vcount(net_human_sarscov2_1),
    area = vcount(net_human_sarscov2_1) ** 2.3,
    repulse.rad = vcount(net_human_sarscov2_1) ** 2.1,
    niter = 3000
)

png(
    filename = 'SARS-CoV2_OmniPath_neighborhood.png',
    width = 7,
    height = 7,
    units = 'in',
    res = 600
)

plot(
    net_human_sarscov2_1,
    layout = fr_layout,
    vertex.size = 5,
    vertex.label.cex = .33,
    vertex.color = ifelse(
        V(net_human_sarscov2_1)$organism == 'SARS-CoV2',
        '#FCCC06',
        ifelse(
            V(net_human_sarscov2_1)$viral_target,
            '#97BE73',
            '#B6B7B9'
        )
    ),
    vertex.label.family = 'DINPro',
    vertex.label.color = '#454447',
    vertex.frame.width = 0,
    vertex.frame.color = NA,
    edge.color = ifelse(
        E(net_human_sarscov2_1)$type == 'host_pathogen',
        '#FCCC06',
        ifelse(
            E(net_human_sarscov2_1)$is_inhibition,
            '#E25C49AA',
            ifelse(
                E(net_human_sarscov2_1)$is_stimulation,
                '#49969AAA',
                '#646567AA'
            )
        )
    ),
    edge.arrow.size = ifelse(
        E(net_human_sarscov2_1)$is_directed,
        .2,
        .0
    ),
#     edge.lty = ifelse(
#         E(net_human_sarscov2_1)$type == 'transcriptional',
#         'dashed',
#         ifelse(
#             E(net_human_sarscov2_1)$type == 'host_pathogen',
#             'dotted',
#             'solid'
#         )
#     ),
    edge.width = .2
)

dev.off()