From f20884645b2ffc1089b39d1b8e2aee9cb02ce400 Mon Sep 17 00:00:00 2001 From: juuussi Date: Tue, 19 Nov 2019 02:23:03 +0200 Subject: [PATCH 1/2] Re-build project structure. Fixed several typos and potential code errors. --- .Rbuildignore | 1 - .gitignore | 5 + DESCRIPTION | 0 NAMESPACE | 0 R/data_description.R | 0 R/script_disease_relevant_tissues.R | 618 +++++++++++------- R/script_integration_efficacy_scores.R | 0 R/script_modulation_score.R | 0 R/script_node_centrality_scores.R | 232 ++++--- R/script_tissue_specific_efficacy_score.R | 0 README.md | 2 +- ThETA.Rproj | 3 + build/partial.rdb | Bin 23769 -> 0 bytes build/vignette.rds | Bin 272 -> 0 bytes data/centrality_score.rda | Bin data/dis_vrnts.rda | Bin data/disease_tissue_zscores.rda | Bin data/geo_gene_sets.rda | Bin data/gtexv7_zscore.rda | Bin data/ppi_strdb_700.rda | Bin inst/REFERENCES.bib | 0 ...15_v7_RNASeQCv1.1.8_gene_median_tpm.gct.gz | Bin inst/{ => extdata}/conversion_enrichr_efo.csv | 0 inst/{ => extdata}/efo_feb2019.obo | 0 man/Borda.Rd | 0 man/centrality_score.Rd | 0 man/dis_vrnts.Rd | 0 man/disease_tissue_zscores.Rd | 0 man/geo.mean.Rd | 0 man/geo_gene_sets.Rd | 0 man/gtexv7_zscore.Rd | 0 man/integrate.scores.Rd | 0 man/l2norm.Rd | 0 man/ppi_strdb_700.Rd | 0 vignettes/DataProcessing.Rmd | 4 +- .../DataProcessing_cache/html/__packages | 1 - vignettes/Introduction.Rmd | 8 +- vignettes/Introduction_cache/html/__packages | 20 - 38 files changed, 553 insertions(+), 341 deletions(-) mode change 100755 => 100644 DESCRIPTION mode change 100755 => 100644 NAMESPACE mode change 100755 => 100644 R/data_description.R mode change 100755 => 100644 R/script_disease_relevant_tissues.R mode change 100755 => 100644 R/script_integration_efficacy_scores.R mode change 100755 => 100644 R/script_modulation_score.R mode change 100755 => 100644 R/script_node_centrality_scores.R mode change 100755 => 100644 R/script_tissue_specific_efficacy_score.R delete mode 100644 build/partial.rdb delete mode 100644 build/vignette.rds mode change 100755 => 100644 data/centrality_score.rda mode change 100755 => 100644 data/dis_vrnts.rda mode change 100755 => 100644 data/disease_tissue_zscores.rda mode change 100755 => 100644 data/geo_gene_sets.rda mode change 100755 => 100644 data/gtexv7_zscore.rda mode change 100755 => 100644 data/ppi_strdb_700.rda mode change 100755 => 100644 inst/REFERENCES.bib rename inst/{ => extdata}/GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_median_tpm.gct.gz (100%) mode change 100755 => 100644 rename inst/{ => extdata}/conversion_enrichr_efo.csv (100%) mode change 100755 => 100644 rename inst/{ => extdata}/efo_feb2019.obo (100%) mode change 100755 => 100644 mode change 100755 => 100644 man/Borda.Rd mode change 100755 => 100644 man/centrality_score.Rd mode change 100755 => 100644 man/dis_vrnts.Rd mode change 100755 => 100644 man/disease_tissue_zscores.Rd mode change 100755 => 100644 man/geo.mean.Rd mode change 100755 => 100644 man/geo_gene_sets.Rd mode change 100755 => 100644 man/gtexv7_zscore.Rd mode change 100755 => 100644 man/integrate.scores.Rd mode change 100755 => 100644 man/l2norm.Rd mode change 100755 => 100644 man/ppi_strdb_700.Rd mode change 100755 => 100644 vignettes/DataProcessing.Rmd mode change 100755 => 100644 vignettes/Introduction.Rmd delete mode 100644 vignettes/Introduction_cache/html/__packages diff --git a/.Rbuildignore b/.Rbuildignore index 272bf8d..91114bf 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -1,3 +1,2 @@ -..Rcheck ^.*\.Rproj$ ^\.Rproj\.user$ diff --git a/.gitignore b/.gitignore index 158593b..d3b7c26 100644 --- a/.gitignore +++ b/.gitignore @@ -20,6 +20,11 @@ # produced vignettes vignettes/*.html vignettes/*.pdf +vignettes/*_cache/* +vignettes/*.R + +.build.timestamp +build/ # OAuth2 token, see https://github.com/hadley/httr/releases/tag/v0.3 .httr-oauth diff --git a/DESCRIPTION b/DESCRIPTION old mode 100755 new mode 100644 diff --git a/NAMESPACE b/NAMESPACE old mode 100755 new mode 100644 diff --git a/R/data_description.R b/R/data_description.R old mode 100755 new mode 100644 diff --git a/R/script_disease_relevant_tissues.R b/R/script_disease_relevant_tissues.R old mode 100755 new mode 100644 index 47a2e79..338226d --- a/R/script_disease_relevant_tissues.R +++ b/R/script_disease_relevant_tissues.R @@ -5,12 +5,16 @@ #'@param tissue_expr_data a data matrix containing gene expression values. #'@return a numeric matrix of gene expression values in the form of Z-scores. Columns are tissues and rows are diseases. #'@export -tissue.expr <- function(tissue_expr_data){ - if(is.null(rownames(tissue_expr_data))|is.null(colnames(tissue_expr_data))){ +tissue.expr <- function(tissue_expr_data) { + if (is.null(rownames(tissue_expr_data)) | + is.null(colnames(tissue_expr_data))) { stop('Both colnames and rownames for tissue_expr_data must be provided!') } - tissue_expr_data <- tissue_expr_data[!rowSums(tissue_expr_data)==0,] - tissue_expr_zscore <- t(apply(tissue_expr_data, 1, function(x) (x-mean(x))/sd(x))) + tissue_expr_data <- + tissue_expr_data[!rowSums(tissue_expr_data) == 0, ] + tissue_expr_zscore <- + t(apply(tissue_expr_data, 1, function(x) + (x - mean(x)) / sd(x))) return(tissue_expr_zscore) } @@ -31,9 +35,12 @@ tissue.expr <- function(tissue_expr_data){ #'@import RCurl #'@import XML #'@import SPARQL -disease.vrnts <- function(x, id_type='efo', min_score=0, curated=F){ - endpoint<-"http://rdf.disgenet.org/sparql/" - partI<-'PREFIX rdf: +disease.vrnts <- function(x, + id_type = 'efo', + min_score = 0, + curated = F) { + endpoint <- "http://rdf.disgenet.org/sparql/" + partI <- 'PREFIX rdf: PREFIX rdfs: PREFIX owl: PREFIX xsd: @@ -57,11 +64,14 @@ disease.vrnts <- function(x, id_type='efo', min_score=0, curated=F){ ?score WHERE {?disease skos:exactMatch ' - if(id_type=='mesh') url<-paste(' .',sep='') - else if(id_type=='omim') url<-paste(' .',sep='') - else if (id_type=='efo') url<-paste(' .',sep='') + if (id_type == 'mesh') + url <- paste(' .', sep = '') + else if (id_type == 'omim') + url <- paste(' .', sep = '') + else if (id_type == 'efo') + url <- paste(' .', sep = '') - partII<-'?gene sio:SIO_000205 ?GeneSymbol . + partII <- '?gene sio:SIO_000205 ?GeneSymbol . ?scoreIRI sio:SIO_000300 ?score . ?variant so:associated_with ?gene . ?vda sio:SIO_000628 ?variant,?disease ; @@ -73,29 +83,46 @@ disease.vrnts <- function(x, id_type='efo', min_score=0, curated=F){ ?scoreIRI rdf:type ncit:C25338 . ?variant a so:0001060 .' - filter<-paste('FILTER(?score >= ',min_score,')',sep='') + filter <- paste('FILTER(?score >= ', min_score, ')', sep = '') - if(curated==T){ - filterII<-'FILTER (!regex(?source,"BEFREE"))' - filter<-paste(filter,filterII,sep='\n ') + if (curated == T) { + filterII <- 'FILTER (!regex(?source,"BEFREE"))' + filter <- paste(filter, filterII, sep = '\n ') } - partIII<-'} ORDER BY ?GeneSymbol' + partIII <- '} ORDER BY ?GeneSymbol' - query<-paste(partI,url,partII,filter,partIII,sep='\n ') + query <- paste(partI, url, partII, filter, partIII, sep = '\n ') - res<-SPARQL::SPARQL(endpoint,query)$results + res <- SPARQL::SPARQL(endpoint, query)$results - if(nrow(res !=0)){ - etz <- sapply(strsplit(res$gene,'/'),function(x)gsub('>','',tail(x,1))) - symb <- sapply(strsplit(res$GeneSymbol,'/'),function(x)gsub('>','',tail(x,1))) + if (nrow(res != 0)) { + etz <- + sapply(strsplit(res$gene, '/'), function(x) + gsub('>', '', tail(x, 1))) + symb <- + sapply(strsplit(res$GeneSymbol, '/'), function(x) + gsub('>', '', tail(x, 1))) scr <- res$score - df <- data.frame(entrez=etz,symbol=symb,score=scr,stringsAsFactors = F) - df <- df[!duplicated(df),] - df <- df[order(df$score,decreasing=T),] + df <- + data.frame( + entrez = etz, + symbol = symb, + score = scr, + stringsAsFactors = F + ) + df <- df[!duplicated(df), ] + df <- df[order(df$score, decreasing = T), ] } - else df <- data.frame(entrez=character(),symbol=character(),score=numeric(),stringsAsFactors = F) + else + df <- + data.frame( + entrez = character(), + symbol = character(), + score = numeric(), + stringsAsFactors = F + ) return(df) } @@ -142,104 +169,170 @@ disease.vrnts <- function(x, id_type='efo', min_score=0, curated=F){ #'@importFrom Rdpack reprompt #'@importFrom igraph V E graph_from_data_frame distances induced.subgraph clusters #'@importFrom scales rescale -dis.rel.tissues <- function(disease_genes, ppi_network, weighted = FALSE, tissue_expr_data, thr = 1, top = 20, rand = 100, verbose=FALSE){ - if(is.null(rownames(tissue_expr_data))|is.null(colnames(tissue_expr_data))){ - stop('Both colnames and rownames for tissue_expr_data must be provided!') - } - if(!is.character(ppi_network[,1])) ppi_network[,1] <- as.character(ppi_network[,1]) - if(!is.character(ppi_network[,2])) ppi_network[,2] <- as.character(ppi_network[,2]) - ppi_network <- ppi_network[!duplicated(ppi_network[,1:2]),] - ppi_network_node<-unique(unlist(ppi_network[,1:2])) - universe <- intersect(ppi_network_node,rownames(tissue_expr_data)) - if(length(universe)==0) stop('No corresponding IDs between ppi_network and tissue_expr_data!') - else if(length(universe)!=length(ppi_network_node)| - length(universe)!=nrow(tissue_expr_data)){ - if (verbose) print(paste(length(universe),'IDs in common between ppi_network and tissue_expr_data will be considered.', sep=' ')) - tissue_expr_data<-tissue_expr_data[universe,] - } - disease_genes_size<-length(disease_genes) - disease_genes<-intersect(disease_genes,universe) - if(length(disease_genes)==0) stop('No disease-associated ID match with ppi_network and tissue_expr_data!') - else if( disease_genes_size!=length(disease_genes)){ - if(verbose) print(paste(length(disease_genes),'disease-associated IDs match with ppi_network and tissue_expr_data. Computing z-scores for disease-tissue pairs...',sep=' ')) - } - gene_sets <- list() - max_comps <- c() - sd_rand <- c() - mean_rand <- c() - zscores <- c() - tissue_genes <- sapply(colnames(tissue_expr_data), function(i) rownames(tissue_expr_data)[tissue_expr_data[,i]>=thr]) - for(tgs in tissue_genes) { - # intersect the ppi network with tissue-specific genes - df <- ppi_network[ppi_network[,1] %in% tgs & ppi_network[,2] %in% tgs,] - g <- igraph::graph_from_data_frame(df[,1:2], directed = FALSE) - genes <- intersect(disease_genes, igraph::V(g)$name) - if(length(genes) > 0){ - if(top > 0) { - # calculate the radial distance between each gene and all disease-relevant genes - if(weighted) igraph::E(g)$weight <- scales::rescale(as.numeric(df[,3]),c(1,.Machine$double.eps)) - rad_dist <- igraph::distances(g, to=intersect(disease_genes,igraph::V(g)$name), mode = "out", weights = NULL, algorithm = "dijkstra") - rad_dist <- rowMeans(rad_dist) - genes <- unique(c(genes, names(sort(rad_dist))[1:top])) # extended set of disease-relevant genes - } +dis.rel.tissues <- + function(disease_genes, + ppi_network, + weighted = FALSE, + tissue_expr_data, + thr = 1, + top = 20, + rand = 100, + verbose = FALSE) { + if (is.null(rownames(tissue_expr_data)) | + is.null(colnames(tissue_expr_data))) { + stop('Both colnames and rownames for tissue_expr_data must be provided!') } - gene_sets[[length(gene_sets)+1]] <- genes - # calculate the maximal components by considering disease-relevant genes - if(length(genes) > 1){ - gs <- igraph::induced.subgraph(graph=g, vids=genes) - max_comps <- c(max_comps, max(igraph::clusters(gs)$csize)) - # random - rand_max_comps <- c() - for(i in 1:rand) { - gr <- igraph::induced.subgraph(graph=g, vids=sample(igraph::V(g)$name, length(genes), replace = T)) - rand_max_comps <- c(rand_max_comps, max(igraph::clusters(gr)$csize)) - } - #hist(rand_max_comps) - mean_rand <- c(mean_rand, mean(rand_max_comps)) - sd_rand <- c(sd_rand, sd(rand_max_comps)) - zscores <- c(zscores, (max_comps[length(max_comps)] - mean(rand_max_comps)) / sd(rand_max_comps)) + if (!is.character(ppi_network[, 1])) + ppi_network[, 1] <- as.character(ppi_network[, 1]) + if (!is.character(ppi_network[, 2])) + ppi_network[, 2] <- as.character(ppi_network[, 2]) + ppi_network <- ppi_network[!duplicated(ppi_network[, 1:2]), ] + ppi_network_node <- unique(unlist(ppi_network[, 1:2])) + universe <- intersect(ppi_network_node, rownames(tissue_expr_data)) + if (length(universe) == 0) + stop('No corresponding IDs between ppi_network and tissue_expr_data!') + else if (length(universe) != length(ppi_network_node) | + length(universe) != nrow(tissue_expr_data)) { + if (verbose) + print( + paste( + length(universe), + 'IDs in common between ppi_network and tissue_expr_data will be considered.', + sep = ' ' + ) + ) + tissue_expr_data <- tissue_expr_data[universe, ] } - else{ - max_comps <- c(max_comps, length(genes)) - mean_rand <- c(mean_rand, length(genes)) - sd_rand <- c(sd_rand, 0) - zscores <- c(zscores, 0) + disease_genes_size <- length(disease_genes) + disease_genes <- intersect(disease_genes, universe) + if (length(disease_genes) == 0) + stop('No disease-associated ID match with ppi_network and tissue_expr_data!') + else if (disease_genes_size != length(disease_genes)) { + if (verbose) + print( + paste( + length(disease_genes), + 'disease-associated IDs match with ppi_network and tissue_expr_data. Computing z-scores for disease-tissue pairs...', + sep = ' ' + ) + ) } + gene_sets <- list() + max_comps <- c() + sd_rand <- c() + mean_rand <- c() + zscores <- c() + tissue_genes <- + sapply(colnames(tissue_expr_data), function(i) + rownames(tissue_expr_data)[tissue_expr_data[, i] >= thr]) + for (tgs in tissue_genes) { + # intersect the ppi network with tissue-specific genes + df <- + ppi_network[ppi_network[, 1] %in% tgs & ppi_network[, 2] %in% tgs, ] + g <- igraph::graph_from_data_frame(df[, 1:2], directed = FALSE) + genes <- intersect(disease_genes, igraph::V(g)$name) + if (length(genes) > 0) { + if (top > 0) { + # calculate the radial distance between each gene and all disease-relevant genes + if (weighted) + igraph::E(g)$weight <- + scales::rescale(as.numeric(df[, 3]), c(1, .Machine$double.eps)) + rad_dist <- + igraph::distances( + g, + to = intersect(disease_genes, igraph::V(g)$name), + mode = "out", + weights = NULL, + algorithm = "dijkstra" + ) + rad_dist <- rowMeans(rad_dist) + genes <- + unique(c(genes, names(sort(rad_dist))[1:top])) # extended set of disease-relevant genes + } + } + gene_sets[[length(gene_sets) + 1]] <- genes + # calculate the maximal components by considering disease-relevant genes + if (length(genes) > 1) { + gs <- igraph::induced.subgraph(graph = g, vids = genes) + max_comps <- c(max_comps, max(igraph::clusters(gs)$csize)) + # random + rand_max_comps <- c() + for (i in 1:rand) { + gr <- + igraph::induced.subgraph(graph = g, + vids = sample(igraph::V(g)$name, length(genes), replace = T)) + rand_max_comps <- + c(rand_max_comps, max(igraph::clusters(gr)$csize)) + } + #hist(rand_max_comps) + mean_rand <- c(mean_rand, mean(rand_max_comps)) + sd_rand <- c(sd_rand, sd(rand_max_comps)) + zscores <- + c(zscores, (max_comps[length(max_comps)] - mean(rand_max_comps)) / sd(rand_max_comps)) + } + else{ + max_comps <- c(max_comps, length(genes)) + mean_rand <- c(mean_rand, length(genes)) + sd_rand <- c(sd_rand, 0) + zscores <- c(zscores, 0) + } + } + final_df <- + data.frame( + maxc = max_comps, + mer = mean_rand, + sdr = sd_rand, + z = zscores, + genes = sapply(gene_sets, length) + ) + rownames(final_df) <- colnames(tissue_expr_data) + return(final_df) } - final_df <- data.frame(maxc = max_comps, mer = mean_rand, sdr = sd_rand, z = zscores, genes = sapply(gene_sets, length)) - rownames(final_df) <- colnames(tissue_expr_data) - return(final_df) -} #'@rdname dis.rel.tissues #'@export #'@importFrom snow makeCluster stopCluster #'@importFrom doParallel registerDoParallel #'@importFrom foreach foreach %dopar% -list.dis.rel.tissues <- function(disease_gene_list, ppi_network, weighted = FALSE, tissue_expr_data, thr = 1, top = 20, rand = 100, parallel = NULL) { - if(!is.list(disease_gene_list))stop('Argument disease_gene_list is not a list!') - if(is.null(names(disease_gene_list)))stop('Names for disease_gene_list must be provided!') - mat_set_tiss_zscore <- NULL - if(!is.null(parallel)) { - cl <- snow::makeCluster(parallel) - doParallel::registerDoParallel(cl) - `%dopar%` <- foreach::`%dopar%` - dis_rlvnt_tiss <- foreach::foreach(i=disease_gene_list, .export = 'get.disease.relevant.tissues') %dopar% - get.disease.relevant.tissues(i, ppi_network, weighted, tissue_expr_data, thr, top, rand) - snow::stopCluster(cl) - } - else { - warning("A parallel computation is highly recommended",immediate. = T) - dis_rlvnt_tiss <- list() - for(d in disease_gene_list) { - dis_rlvnt_tiss[[length(dis_rlvnt_tiss)+1]] <- get.disease.relevant.tissues(d, ppi_network, weighted, tissue_expr_data, thr, top, rand) +list.dis.rel.tissues <- + function(disease_gene_list, + ppi_network, + weighted = FALSE, + tissue_expr_data, + thr = 1, + top = 20, + rand = 100, + parallel = NULL) { + if (!is.list(disease_gene_list)) + stop('Argument disease_gene_list is not a list!') + if (is.null(names(disease_gene_list))) + stop('Names for disease_gene_list must be provided!') + mat_set_tiss_zscore <- NULL + if (!is.null(parallel)) { + cl <- snow::makeCluster(parallel) + doParallel::registerDoParallel(cl) + `%dopar%` <- foreach::`%dopar%` + dis_rlvnt_tiss <- + foreach::foreach(i = disease_gene_list, .export = 'get.disease.relevant.tissues') %dopar% + get.disease.relevant.tissues(i, ppi_network, weighted, tissue_expr_data, thr, top, rand) + snow::stopCluster(cl) + } + else { + warning("A parallel computation is highly recommended", immediate. = T) + dis_rlvnt_tiss <- list() + for (d in disease_gene_list) { + dis_rlvnt_tiss[[length(dis_rlvnt_tiss) + 1]] <- + get.disease.relevant.tissues(d, ppi_network, weighted, tissue_expr_data, thr, top, rand) + } } + names(dis_rlvnt_tiss) <- names(disease_gene_list) + df <- + do.call(function(...) + mapply(rbind, ..., SIMPLIFY = F), dis_rlvnt_tiss) + for (i in 1:length(df)) + colnames(df[[i]]) <- colnames(tissue_expr_data) + return(df) } - names(dis_rlvnt_tiss) <- names(disease_gene_list) - df <- do.call(function(...)mapply(rbind,...,SIMPLIFY=F),dis_rlvnt_tiss) - for (i in 1:length(df)) colnames(df[[i]]) <- colnames(tissue_expr_data) - return(df) -} #'Interactive visualization of tissue-specific networks. @@ -263,127 +356,196 @@ list.dis.rel.tissues <- function(disease_gene_list, ppi_network, weighted = FALS #'@export #'@importFrom clusterProfiler enrichKEGG #'@importFrom clusterProfiler enrichGO -visualize.graph <- function(tissue_scores, disease_genes, ppi_network, directed_network = FALSE, - tissue_expr_data, top_targets = NULL, db='kegg',verbose = FALSE){ - - if(is.null(top_targets)) stop('Please specifiy a set of targets (ENTREZ ids)!') - if(is.null(rownames(tissue_expr_data))|is.null(colnames(tissue_expr_data))){ - stop('Both colnames and rownames for tissue_expr_data must be provided.')} - for(i in 1:2) ppi_network[,i] <- as.character(ppi_network[,i]) - ppi_network <- ppi_network[!duplicated(ppi_network[,1:2]),] - ppi_network_size <- nrow(ppi_network) - if(directed_network) idx <- ppi_network[,1]%in%rownames(tissue_expr_data) - else idx <- ppi_network[,1]%in%rownames(tissue_expr_data) & ppi_network[,2]%in%rownames(tissue_expr_data) - ppi_network <- ppi_network[idx,] - if(nrow(ppi_network)==0) stop('No corresponding IDs between ppi_network and tissue_expr_data.') - else if(ppi_network_size!=nrow(ppi_network)){ - if(verbose) print(paste(nrow(ppi_network),'out of',ppi_network_size,'network edges selected.', sep=' ')) - } - if(!directed_network){ - if(verbose) print('Undirect network. Converting to direct network...') - ppi_network_rev <- ppi_network[,c(2:1)] - colnames(ppi_network_rev) <- colnames(ppi_network)[1:2] - ppi_network <- rbind(ppi_network[,1:2],ppi_network_rev) - } - disease_genes_size <- length(disease_genes) - disease_genes <- intersect(disease_genes,ppi_network[,2]) - if(length(disease_genes)==0) stop('No disease-associated ID match with ppi_network and tissue_expr_data!') - else if( disease_genes_size!=length(disease_genes)){ - if(verbose) print(paste(length(disease_genes),'disease-associated IDs are reachable from the network.',sep=' ')) - } - if(!(db %in% c('kegg','BP','MF','CC'))) stop('Unused value for argument db') - idx <- order(tissue_scores$avg_tissue_score,decreasing=T) - tissue_scores <- tissue_scores[idx,] - g <- igraph::graph_from_edgelist(as.matrix(ppi_network[,1:2]), directed=T) - tissue_expr_data <- scales::rescale(tissue_expr_data,c(1,.Machine$double.eps)) - sign_tiss <- colnames(tissue_scores)[-ncol(tissue_scores)] - sh.path <- lapply(sign_tiss,function(i){ - igraph::E(g)$weight <- tissue_expr_data[ppi_network[,1],i] - lapply(top_targets,function(x) - igraph::shortest_paths(g, from=x, to=disease_genes, mode = "out", weights = NULL,output='epath')$epath) - }) - names(sh.path) <- sign_tiss - for(i in names(sh.path)) names(sh.path[[i]]) <- top_targets - all_target_path <-lapply(sh.path,function(x)lapply(x,function(y)do.call(c,y))) - all_tissue_path <- lapply(all_target_path,function(x)do.call(c,x)) - ids <- lapply(all_tissue_path,igraph::as_ids) - shinyApp( - ui = fluidPage( - fluidRow( +visualize.graph <- + function(tissue_scores, + disease_genes, + ppi_network, + directed_network = FALSE, + tissue_expr_data, + top_targets = NULL, + db = 'kegg', + verbose = FALSE) { + if (is.null(top_targets)) + stop('Please specifiy a set of targets (ENTREZ ids)!') + if (is.null(rownames(tissue_expr_data)) | + is.null(colnames(tissue_expr_data))) { + stop('Both colnames and rownames for tissue_expr_data must be provided.') + } + for (i in 1:2) + ppi_network[, i] <- as.character(ppi_network[, i]) + ppi_network <- ppi_network[!duplicated(ppi_network[, 1:2]), ] + ppi_network_size <- nrow(ppi_network) + if (directed_network) + idx <- ppi_network[, 1] %in% rownames(tissue_expr_data) + else + idx <- + ppi_network[, 1] %in% rownames(tissue_expr_data) & + ppi_network[, 2] %in% rownames(tissue_expr_data) + ppi_network <- ppi_network[idx, ] + if (nrow(ppi_network) == 0) + stop('No corresponding IDs between ppi_network and tissue_expr_data.') + else if (ppi_network_size != nrow(ppi_network)) { + if (verbose) + print(paste( + nrow(ppi_network), + 'out of', + ppi_network_size, + 'network edges selected.', + sep = ' ' + )) + } + if (!directed_network) { + if (verbose) + print('Undirect network. Converting to direct network...') + ppi_network_rev <- ppi_network[, c(2:1)] + colnames(ppi_network_rev) <- colnames(ppi_network)[1:2] + ppi_network <- rbind(ppi_network[, 1:2], ppi_network_rev) + } + disease_genes_size <- length(disease_genes) + disease_genes <- intersect(disease_genes, ppi_network[, 2]) + if (length(disease_genes) == 0) + stop('No disease-associated ID match with ppi_network and tissue_expr_data!') + else if (disease_genes_size != length(disease_genes)) { + if (verbose) + print(paste( + length(disease_genes), + 'disease-associated IDs are reachable from the network.', + sep = ' ' + )) + } + if (!(db %in% c('kegg', 'BP', 'MF', 'CC'))) + stop('Unused value for argument db') + idx <- order(tissue_scores$avg_tissue_score, decreasing = T) + tissue_scores <- tissue_scores[idx, ] + g <- + igraph::graph_from_edgelist(as.matrix(ppi_network[, 1:2]), directed = T) + tissue_expr_data <- + scales::rescale(tissue_expr_data, c(1, .Machine$double.eps)) + sign_tiss <- colnames(tissue_scores)[-ncol(tissue_scores)] + sh.path <- lapply(sign_tiss, function(i) { + igraph::E(g)$weight <- tissue_expr_data[ppi_network[, 1], i] + lapply(top_targets, function(x) + igraph::shortest_paths( + g, + from = x, + to = disease_genes, + mode = "out", + weights = NULL, + output = 'epath' + )$epath) + }) + names(sh.path) <- sign_tiss + for (i in names(sh.path)) + names(sh.path[[i]]) <- top_targets + all_target_path <- + lapply(sh.path, function(x) + lapply(x, function(y) + do.call(c, y))) + all_tissue_path <- lapply(all_target_path, function(x) + do.call(c, x)) + ids <- lapply(all_tissue_path, igraph::as_ids) + shinyApp( + ui = fluidPage(fluidRow( column( width = 2, - checkboxGroupInput("tissue","Select tissue/s",sign_tiss,selected = sign_tiss[1]) + checkboxGroupInput("tissue", "Select tissue/s", sign_tiss, selected = sign_tiss[1]) ), - column( - width = 12, - visNetworkOutput("network") - ) + column(width = 12, + visNetworkOutput("network")) ), - fluidRow( - column( - width = 12, - tableOutput("table") - ) - ) - ), - server = function(input, output) { - output$network <- renderVisNetwork({ - validate(need(input$tissue!='','Please select one or more tissues')) - if(length(input$tissue)>1){ - edge <- strsplit(Reduce(intersect,ids[input$tissue]),split='|',fixed=T) - } - else if (length(input$tissue) ==1){ - edge <- strsplit(ids[[input$tissue]],split='|',fixed=T) - } - edge <- as.data.frame(do.call(rbind,edge)) - edge <- edge[!duplicated(edge),] - colnames(edge) <- c('from','to') - edge$width <- scales::rescale(rowMeans(1-tissue_expr_data[edge$from,input$tissue,drop=F]),c(1,5)) - nb <- unique(unlist(edge[,1:2])) - node <- data.frame(id=nb,label=nb,stringsAsFactors = F) - node$value <- rowMeans(tissue_scores[nb,input$tissue,drop=F]) - gp <- rep('bridge gene',nrow(node)) - gp[node$label%in%top_targets] <- 'target gene' - gp[node$label%in%disease_genes] <- 'disease gene' - node$group <- gp - node$shape <- c("circle","star","diamond")[as.numeric(as.factor(gp))] - visNetwork::visNetwork(node, edge, height = "1500px", width = "500%") %>% - visGroups(groupname = "target gene", color = list(background = "skyblue", border = "deepskyblue")) %>% - visGroups(groupname = "disease gene", color = list(background = "lightcoral", border = "red")) %>% - visGroups(groupname = "bridge gene", shape="circle") %>% - visLegend(addNodes=list( - list(label = "disease gene", shape = "star", - color = list(background = "lightcoral", border = "red")), - list(label = "target gene", shape = "diamond", - color = list(background = "skyblue", border = "deepskyblue")), - list(label = "bridge gene", shape = "circle")), - useGroups = FALSE,width=0.1) %>% - visLayout(hierarchical = TRUE) %>% # visLayout(randomSeed = 123) - visEvents(select = "function(nodes) { + fluidRow(column( + width = 12, + tableOutput("table") + ))), + server = function(input, output) { + output$network <- renderVisNetwork({ + validate(need(input$tissue != '', 'Please select one or more tissues')) + if (length(input$tissue) > 1) { + edge <- + strsplit(Reduce(intersect, ids[input$tissue]), + split = '|', + fixed = T) + } + else if (length(input$tissue) == 1) { + edge <- strsplit(ids[[input$tissue]], split = '|', fixed = T) + } + edge <- as.data.frame(do.call(rbind, edge)) + edge <- edge[!duplicated(edge), ] + colnames(edge) <- c('from', 'to') + edge$width <- + scales::rescale(rowMeans(1 - tissue_expr_data[edge$from, input$tissue, drop = + F]), c(1, 5)) + nb <- unique(unlist(edge[, 1:2])) + node <- data.frame(id = nb, + label = nb, + stringsAsFactors = F) + node$value <- + rowMeans(tissue_scores[nb, input$tissue, drop = F]) + gp <- rep('bridge gene', nrow(node)) + gp[node$label %in% top_targets] <- 'target gene' + gp[node$label %in% disease_genes] <- 'disease gene' + node$group <- gp + node$shape <- + c("circle", "star", "diamond")[as.numeric(as.factor(gp))] + visNetwork::visNetwork(node, edge, height = "1500px", width = "500%") %>% + visGroups( + groupname = "target gene", + color = list(background = "skyblue", border = "deepskyblue") + ) %>% + visGroups( + groupname = "disease gene", + color = list(background = "lightcoral", border = "red") + ) %>% + visGroups(groupname = "bridge gene", shape = "circle") %>% + visLegend( + addNodes = list( + list( + label = "disease gene", + shape = "star", + color = list(background = "lightcoral", border = "red") + ), + list( + label = "target gene", + shape = "diamond", + color = list(background = "skyblue", border = "deepskyblue") + ), + list(label = "bridge gene", shape = "circle") + ), + useGroups = FALSE, + width = 0.1 + ) %>% + visLayout(hierarchical = TRUE) %>% # visLayout(randomSeed = 123) + visEvents(select = "function(nodes) { Shiny.onInputChange('current_node_id', nodes.nodes); ;}") - }) - output$table <- renderTable({ - validate(need(input$tissue!='',NULL), - need(input$current_node_id!=''&input$current_node_id%in%top_targets,'Please select a target gene')) - current_target_path <- sapply(all_target_path[input$tissue],function(x) x[input$current_node_id]) - current_target_path <- Reduce(igraph::intersection,current_target_path) - target_interactors <- unique(c(igraph::ends(g,current_target_path))) - if(db=='kegg') res <- clusterProfiler::enrichKEGG(target_interactors,organism = 'hsa')@result - else res <- clusterProfiler::enrichGO(target_interactors, 'org.Hs.eg.db', ont = db)@result - res <- res[res$p.adjust<0.05,] - if (nrow(res) > 25) res <- res[1:25,] - res - }) - } - ) + }) + output$table <- renderTable({ + validate( + need(input$tissue != '', NULL), + need( + input$current_node_id != '' & + input$current_node_id %in% top_targets, + 'Please select a target gene' + ) + ) + current_target_path <- + sapply(all_target_path[input$tissue], function(x) + x[input$current_node_id]) + current_target_path <- + Reduce(igraph::intersection, current_target_path) + target_interactors <- + unique(c(igraph::ends(g, current_target_path))) + if (db == 'kegg') + res <- + clusterProfiler::enrichKEGG(target_interactors, organism = 'hsa')@result + else + res <- + clusterProfiler::enrichGO(target_interactors, 'org.Hs.eg.db', ont = db)@result + res <- res[res$p.adjust < 0.05, ] + if (nrow(res) > 25) + res <- res[1:25, ] + res + }) + } + ) } - - - - - - - - - diff --git a/R/script_integration_efficacy_scores.R b/R/script_integration_efficacy_scores.R old mode 100755 new mode 100644 diff --git a/R/script_modulation_score.R b/R/script_modulation_score.R old mode 100755 new mode 100644 diff --git a/R/script_node_centrality_scores.R b/R/script_node_centrality_scores.R old mode 100755 new mode 100644 index e2a1b6d..e314407 --- a/R/script_node_centrality_scores.R +++ b/R/script_node_centrality_scores.R @@ -7,8 +7,8 @@ #'@param x numeric, a vector of values. #'@param na.rm logical, whether or not to remove NA values from the calculation. #'@return The L2 norm of \code{x} -l2norm <- function (x, na.rm = TRUE) { - return(sqrt(mean(abs(x)^2, na.rm = TRUE))) +l2norm <- function(x, na.rm = TRUE) { + return(sqrt(mean(abs(x) ^ 2, na.rm = TRUE))) } #'Calculate The Geometric Mean #' @@ -19,7 +19,7 @@ l2norm <- function (x, na.rm = TRUE) { #'@param x numeric, a vector of values. #'@param na.rm logical, whether or not to remove NA values from the calculation. #'@return the geometric mean of \code{x} -geo.mean <- function (x, na.rm = TRUE) { +geo.mean <- function(x, na.rm = TRUE) { return(exp(mean(log(x), na.rm = TRUE))) } @@ -36,16 +36,24 @@ geo.mean <- function (x, na.rm = TRUE) { #'@return a list with two components:\cr #' -\strong{TopK} a matrix with 4 columns corresponding to rankings by each of the 4 aggregation functions;\cr #' -\strong{Scores} a matrix with 4 columns corresponding to the Borda scores from each of the 4 aggregation functions. -Borda <- function(input){ +Borda <- function(input) { aggreg.function = c("mean", "median", "geo.mean", "l2norm") - allranks = data.frame(matrix(0, nrow = nrow(input), ncol = length(aggreg.function))) + allranks = data.frame(matrix( + 0, + nrow = nrow(input), + ncol = length(aggreg.function) + )) names(allranks) = aggreg.function - allfun = data.frame(matrix(0, nrow = nrow(input), ncol = length(aggreg.function))) + allfun = data.frame(matrix( + 0, + nrow = nrow(input), + ncol = length(aggreg.function) + )) names(allfun) = aggreg.function for (i in 1:length(aggreg.function)) { allfun[, i] = apply(input, 1, aggreg.function[i], na.rm = TRUE) - allranks[, i] = rank(allfun[,i],ties.method ='random') + allranks[, i] = rank(allfun[, i], ties.method = 'random') } results = list(allranks, allfun) names(results) = c("TopK", "Scores") @@ -83,81 +91,137 @@ Borda <- function(input){ #'@importFrom doParallel registerDoParallel #'@importFrom foreach foreach %dopar% #'@importFrom mltools bin_data -node.centrality <- function(ppi_network, tissue_expr_data, agg_function = 4, directed_network=F, parallel = NULL, verbose = FALSE) { - if(length(rownames(tissue_expr_data))==0|length(colnames(tissue_expr_data))==0){ - stop('Both colnames and rownames for tissue_expr_data must be provided') +node.centrality <- + function(ppi_network, + tissue_expr_data, + agg_function = 4, + directed_network = F, + parallel = NULL, + verbose = FALSE) { + if (length(rownames(tissue_expr_data)) == 0 | + length(colnames(tissue_expr_data)) == 0) { + stop('Both colnames and rownames for tissue_expr_data must be provided') + } + for (i in 1:2) + ppi_network[, i] <- as.character(ppi_network[, i]) + ppi_network <- ppi_network[!duplicated(ppi_network[, 1:2]), ] + ppi_network_size <- nrow(ppi_network) + if (directed_network) + idx <- ppi_network[, 1] %in% rownames(tissue_expr_data) + else + idx <- + ppi_network[, 1] %in% rownames(tissue_expr_data) & + ppi_network[, 2] %in% rownames(tissue_expr_data) + ppi_network <- ppi_network[idx, ] + if (nrow(ppi_network) == 0) + stop('No corresponding IDs between ppi_network and tissue_expr_data.') + else if (ppi_network_size != nrow(ppi_network)) { + if (verbose) + print(paste( + nrow(ppi_network), + 'out of', + ppi_network_size, + 'network edges selected.', + sep = ' ' + )) + } + if (!directed_network) { + if (verbose) + print('Undirect network. Duplicating network edges...') + ppi_network_rev <- ppi_network[, c(2:1)] + colnames(ppi_network_rev) <- colnames(ppi_network)[1:2] + ppi_network <- rbind(ppi_network[, 1:2], ppi_network_rev) + } + # rescaling z-scores into positive values + tissue_expr_scaled <- + scales::rescale(tissue_expr_data, c(.Machine$double.eps, 1)) + tissue_expr_scaled_btw <- + scales::rescale(tissue_expr_data, c(1, .Machine$double.eps)) + # building a graph from an edge list + g = igraph::graph_from_edgelist(as.matrix(ppi_network[, 1:2]), directed = + T) + # computing node degree, clustering coefficient and betweenness per tissue + node_degree <- c() + node_cc <- c() + node_btw <- c() + if (verbose) + print("Compiling node centrality scores...") + if (is.null(parallel)) { + node_degree <- lapply(colnames(tissue_expr_data), function(i) { + igraph::E(g)$weight <- tissue_expr_scaled[ppi_network[, 1], i] + igraph::strength(g, mode = 'in') + }) + node_cc <- lapply(colnames(tissue_expr_data), function(i) { + igraph::E(g)$weight <- tissue_expr_scaled[ppi_network[, 1], i] + igraph::transitivity(igraph::as.undirected(g), + type = 'barrat', + isolates = 'zero') + }) + node_btw <- lapply(colnames(tissue_expr_data), function(i) { + igraph::E(g)$weight <- tissue_expr_scaled_btw[ppi_network[, 1], i] + igraph::betweenness(g) + }) + } + else { + cl <- snow::makeCluster(parallel) + doParallel::registerDoParallel(cl) + `%dopar%` <- foreach::`%dopar%` + node_degree <- + foreach::foreach(i = colnames(tissue_expr_data)) %dopar% { + igraph::E(g)$weight <- tissue_expr_scaled[ppi_network[, 1], i] + igraph::strength(g, mode = 'in') + } + names(node_degree) <- colnames(tissue_expr_data) + node_cc <- foreach(i = colnames(tissue_expr_data)) %dopar% { + igraph::E(g)$weight <- tissue_expr_scaled[ppi_network[, 1], i] + igraph::transitivity(igraph::as.undirected(g), + type = 'barrat', + isolates = 'zero') + } + names(node_cc) <- colnames(tissue_expr_data) + node_btw <- foreach(i = colnames(tissue_expr_data)) %dopar% { + igraph::E(g)$weight <- tissue_expr_scaled_btw[ppi_network[, 1], i] + igraph::betweenness(g) + } + names(node_btw) <- colnames(tissue_expr_data) + snow::stopCluster(cl) + } + if (verbose) + print("Merging different node centrality scores...") + # aggregating individual centrality rankings with Borda (l2norm aggregation function) + node_cent_scores <- + mapply(function(x, y, z) + cbind(degree = x, cc = y, btw = z), + node_degree, + node_cc, + node_btw, + SIMPLIFY = F) + node_cent_ranks <- + lapply(node_cent_scores, function(x) + apply(x, 2, rank, ties.method = 'random')) + borda_var <- lapply(node_cent_ranks, Borda) + rank_agg_function <- + sapply(borda_var, function(x) + x$TopK[, agg_function]) + rownames(rank_agg_function) <- igraph::V(g)$name + if (verbose) + print("Discretizing the node centrality score...") + # assigning each node a discrete value between 0.25 and 1 depending on its centrality ranking quartile + chunk <- as.data.frame(apply(rank_agg_function, 2, function(x) + mltools::bin_data(x, bins = quantile(x), boundaryType = "[lorc"))) + W <- + lapply(chunk, function(x) + as.numeric(factor(x, levels = rev(levels( + x + )))) / 4) + for (i in 1:length(W)) + names(W[[i]]) <- rownames(node_cent_scores[[1]]) + return( + list( + scores = node_cent_scores, + ranks = node_cent_ranks, + borda = rank_agg_function, + borda.disc = W + ) + ) } - for(i in 1:2) ppi_network[,i]<-as.character(ppi_network[,i]) - ppi_network<-ppi_network[!duplicated(ppi_network[,1:2]),] - ppi_network_size<-nrow(ppi_network) - if(directed_network) idx<-ppi_network[,1]%in%rownames(tissue_expr_data) - else idx<-ppi_network[,1]%in%rownames(tissue_expr_data) & ppi_network[,2]%in%rownames(tissue_expr_data) - ppi_network<-ppi_network[idx,] - if(nrow(ppi_network)==0) stop('No corresponding IDs between ppi_network and tissue_expr_data.') - else if(ppi_network_size!=nrow(ppi_network)){ - if(verbose) print(paste(nrow(ppi_network),'out of',ppi_network_size,'network edges selected.', sep=' ')) - } - if(!directed_network){ - if(verbose) print('Undirect network. Duplicating network edges...') - ppi_network_rev <- ppi_network[,c(2:1)] - colnames(ppi_network_rev) <- colnames(ppi_network)[1:2] - ppi_network <- rbind(ppi_network[,1:2],ppi_network_rev) - } - # rescaling z-scores into positive values - tissue_expr_scaled <- scales::rescale(tissue_expr_data, c(.Machine$double.eps,1)) - tissue_expr_scaled_btw <- scales::rescale(tissue_expr_data, c(1,.Machine$double.eps)) - # building a graph from an edge list - g = igraph::graph_from_edgelist(as.matrix(ppi_network[,1:2]), directed=T) - # computing node degree, clustering coefficient and betweenness per tissue - node_degree <- c() - node_cc <- c() - node_btw <- c() - if(verbose) print("Compiling node centrality scores...") - if(is.null(parallel)) { - node_degree <- lapply(colnames(tissue_expr_data),function(i){ - igraph::E(g)$weight <- tissue_expr_scaled[ppi_network[,1], i] - igraph::strength(g, mode='in')}) - node_cc <- lapply(colnames(tissue_expr_data),function(i){ - igraph::E(g)$weight <- tissue_expr_scaled[ppi_network[,1], i] - igraph::transitivity(igraph::as.undirected(g), type='barrat',isolates='zero')}) - node_btw <- lapply(colnames(tissue_expr_data),function(i){ - igraph::E(g)$weight <- tissue_expr_scaled_btw[ppi_network[,1],i] - igraph::betweenness(g)}) - } - else { - cl <- snow::makeCluster(parallel) - doParallel::registerDoParallel(cl) - `%dopar%` <- foreach::`%dopar%` - node_degree <- foreach::foreach(i=colnames(tissue_expr_data)) %dopar% { - igraph::E(g)$weight <- tissue_expr_scaled[ppi_network[,1], i] - igraph::strength(g, mode='in')} - names(node_degree) <- colnames(tissue_expr_data) - node_cc <- foreach(i=colnames(tissue_expr_data)) %dopar% { - igraph::E(g)$weight <- tissue_expr_scaled[ppi_network[,1], i] - igraph::transitivity(igraph::as.undirected(g), type='barrat',isolates='zero')} - names(node_cc) <- colnames(tissue_expr_data) - node_btw <- foreach(i=colnames(tissue_expr_data)) %dopar% { - igraph::E(g)$weight <- tissue_expr_scaled_btw[ppi_network[,1],i] - igraph::betweenness(g)} - names(node_btw)<-colnames(tissue_expr_data) - snow::stopCluster(cl) - } - if(verbose) print("Merging different node centrality scores...") - # aggregating individual centrality rankings with Borda (l2norm aggregation function) - node_cent_scores <- mapply(function(x,y,z) cbind(degree=x,cc=y,btw=z), node_degree, node_cc, node_btw, SIMPLIFY=F) - node_cent_ranks <- lapply(node_cent_scores, function(x) apply(x, 2, rank, ties.method ='random')) - borda_var <- lapply(node_cent_ranks, Borda) - rank_agg_function <- sapply(borda_var, function(x)x$TopK[,agg_function]) - rownames(rank_agg_function) <- igraph::V(g)$name - if(verbose) print("Discretizing the node centrality score...") - # assigning each node a discrete value between 0.25 and 1 depending on its centrality ranking quartile - chunk <- as.data.frame(apply(rank_agg_function, 2, function(x) - mltools::bin_data(x,bins=quantile(x), boundaryType = "[lorc"))) - W <- lapply(chunk,function(x) as.numeric(factor(x, levels=rev(levels(x))))/4) - for (i in 1:length(W)) names(W[[i]]) <- rownames(node_cent_scores[[1]]) - return(list(scores = node_cent_scores, - ranks = node_cent_ranks, - borda = rank_agg_function, - borda.disc = W)) -} - diff --git a/R/script_tissue_specific_efficacy_score.R b/R/script_tissue_specific_efficacy_score.R old mode 100755 new mode 100644 diff --git a/README.md b/README.md index 7f79662..e2c21a5 100644 --- a/README.md +++ b/README.md @@ -1,7 +1,7 @@ # ThETA ThETA: Transcriptome-driven Efficacy estimates for gene-based TArget discovery -This R package was built to provide the implementation of novel trasncriptome-based efficacy scores of target(gene)-disease associations, which have been recentely described in [Failli et al. 2019](https://www.nature.com/articles/s41598-019-46293-7). It provides utility functions to recompile the scores based on the selection of different disease-relevant gene sets and tissue-specific gene expresstion profiles. Morevoer, it provides the users with an easy access to the disease-gene association scores compiled by the [Open Target platform](https://www.targetvalidation.org/) and functions to merge the OT-based scores with our novel efficacy scores in order to provide a final prioritization of target(gene)-disease associations. Finally, the users can run basic visualization functions in order to visualize the tissue-specific gene networks and biological annotations associated to top drug targets (or genes) and closely related genes slected with a random walk algorithm. +This R package was built to provide the implementation of novel trasncriptome-based efficacy scores of target(gene)-disease associations, which have been recentely described in [Failli et al. 2019](https://www.nature.com/articles/s41598-019-46293-7). It provides utility functions to recompile the scores based on the selection of different disease-relevant gene sets and tissue-specific gene expresstion profiles. Morevoer, it provides the users with an easy access to the disease-gene association scores compiled by the [Open Targets platform](https://www.targetvalidation.org/) and functions to merge the OT-based scores with our novel efficacy scores in order to provide a final prioritization of target(gene)-disease associations. Finally, the users can run basic visualization functions in order to visualize the tissue-specific gene networks and biological annotations associated to top drug targets (or genes) and closely related genes slected with a random walk algorithm. ## Installation diff --git a/ThETA.Rproj b/ThETA.Rproj index 21a4da0..497f8bf 100644 --- a/ThETA.Rproj +++ b/ThETA.Rproj @@ -12,6 +12,9 @@ Encoding: UTF-8 RnwWeave: Sweave LaTeX: pdfLaTeX +AutoAppendNewline: Yes +StripTrailingWhitespace: Yes + BuildType: Package PackageUseDevtools: Yes PackageInstallArgs: --no-multiarch --with-keep.source diff --git a/build/partial.rdb b/build/partial.rdb deleted file mode 100644 index bd62c6db57d7d74f663d93467a68b9918386e4b1..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 23769 zcmZs?bx@oA6ZVTs(ctdx#oe`daWC#pad(H}PVwUIRv@^$ySuxc^!c4Lb7tN*f83KH z6LJ%F_q+SKt|g9w1N-j*dEP1GOfZ;1k|xKb+2Ri+qeLOW+9Knh1O^vK-y*#QfriEf z9-otv^jo-?oDDbzM?;22@@K}YTHmd)v+TmEoanUn;`l_?VBm#iv6d%Bi*UupSD1Li za>5xev#HFh^O||*R&H`ATs*X0%0x79OmuT_O!(#SgR#qJG&Lw)jw*Ax(_VwB9e=9e z&JUI~eH~13+K;%{LBmYZ%)fWd^c?-r@{o(d~)Cw`p<=FbSVOtKKKO)$&{0Y+cwo_ObXY?%i!d8bPC zi?dwF8~zCq@2$@DyG{W=EHdMlx`#4Rm$0V{&4w?HbfzH}cqyg0#Z4?)b9z3l(piDF z6+I>B2!LQNYcuq7twQ`Q_MO?8hCw1xiw`{5YXHgapUns=r-#4N{e^IdM_i=@-HJ~ zcE5z6CmN%PcWm`v7D&jJL}3KB^K~4}a_ox^S$Mc^3(`HB9AEIXiar++&N|b4mcz~< zpKz`Rzi1uf*-^?(E;W2{K)!Ch(}JHwN8Tw2Pd2XC?sw3w{4s60&T{;cz|_%~#s^5W zVF`tVD?!cuKG4|3pWIk>A<5o>u&oBHE}rDH@XZK%f~#X(rRDL$2F2hVcW}$dL&q;QsMiz>Hl9XW5aST#&E)@svk5j_DJIixL zt$cChg8M4J!~IJXdwvvh$liaG{G-!_nK2D2LZ7C2ERd>bzRp0=Le?b`==by22eh!g0g?ZIU3-*-g4 ztUIQbk<7n=xC+8sSb8iNaYR6@Hc||Be{9t5-*QKiX;VBzKD+&xfONIRU`flFhF&kq*C;i-JNxP|CLgnjA%;G7Dk2W` z19QiZay?yWDu+8bcT^oorCOD9cQnUZ zG3bxg&Wz$1>#s6Aq2rKyvAP(zr#`lmv}9agg~z~yrBH?P<%%8-5a!Ut*+;h=#1bf$ zZ#38}+;mb;T$0-)7jq;kyXHC(MrR3_v@fl#+YGG>H$9RuhH#}z0M4GM8@^%X+j zal^)$oTO+Z(95!Y%)yb19I8;~{PozSX0#tv-UmYi%DC;>LYdl7;Cko;sA=nZ;EUkm z%V=y!a{p1-@a5Vaj^sjBHL@x`WmI1N&O{Q`#vEe5GW_HHa6rRF^Ukxp$Lv?uoqUVK>Q9UE%0}yeVB#8i%0uZa2QNi5o#pvsm$FuPDa2A=f={2HHd8Y(9xf;bW@Eo#! zI0}KvMC^#&o2y_q+#+C8skKT?_b^9qjHgPU?2$jN)!t8=rs_Ckk&P|%7Bx%| zZP4y`|8VytKdq(Mws|!Wkn(FnC6Xg^vdPZ$M66nV83{c5TGRN*;pp&FrKyP(rtto}8&; zzC#Vc*;R()d73SU*ddHwRx zuvcshHZ&lR6Xi`wV01Vh zO^%aO2{8$OPYHVZvvTuyzgEmBm&{=3X`>1(j-+bhXPwYzhK-%3>(cZ-G=Z*vUpUGr zGuciZbcWD~h&K?1M7KnH#pJVAc0Zbq*0#xgV;q-)&`?G-e@n;=FY)j(z+gGkt>ccW zB3%zQqkU~@V-JDzaV0Na<6$o%wK1gAnf*$@Z&26I(tfK5@;TD5OTS})wXEjxA&m3M zvuH`qYsn+tXjplq=>fxV{lZb_$X1PU5 zWM&2tETXY$<^8+-0=*f9lb4&1WSiTX#wyJleg7|y3rb5m*b!Rn>aER8p%&|&Ef8)L zC|;#|??bmNG#>$Ji8Mr^Z!DPg$xSD%Rkxg`ext%C9K7+a)|mkTz(RnnL>^#y<6-dO zz!4WWUtPp_L7D)wXx_(Tu(7gOgKUnSg+>4GQVJ(7JN=(LAHuhfF7G0X*0*p%tFLHn zdN)pJ*Y-$_kqV=e4eYM*Ovk)4qIMk919D9Ir_$|5LpmJoQyUh>6O}ah+|pWdgJ3g% z9PJ71Nb&Q1qniWVv}xNcDl#N`%Wp!kBHsv&-a~ft0$~S6Xm#+@sf;m6i+;>A)YsvD zIE7opm)TwQV%yCao*UY!uUcH{+v~Pqhl){lcAg(U@>(xGxoMwQVzSEvI9XV@u5@CW zB=dwe={hUd?=P=bV|Dr%i({spn><6RIHQ}G9dQ4wXr029Bb_o|prC{zhhqy@OW~9@ zA>AE<>aV(G-CO18mlp-975?*EAD9*{7*Q#&yeFs1VPT}ESd81)jqebDepf&j$Njd= zxR~Os4f^$BP0mON9K(C#)(|#*&9mZnyu(zo1#v&pic$H~R=I+BD#}pCw zgejiF@srhvZL!0~x@owfJJj-b*R9{(A7>#nipa_8LEUfplLTN@WL+yu~M*PRm(yG<8St#e+8@U0es`MJUBL6+k0%r{R`*g!#tT2uKIS- z!yyrB4{^iLJBmpYg(?fStgT8}cU&on6UAak?gEs1A6ts{&EoN-Ix89oqBEL^T{j|8 zXN@n6o!U5^+?aHR!Xn6{<(SUV?pRC^)Iviv`4h>e$rg_ zhKpdJ<|@a63zBWBFJu`zr&N8q4}atPEw$@DWVb*qn5B&GZ@w>=4`Dosv#}pgD#44F z@l3bVr;`A`s=`qh+C2}nQoy>3A^$t$ej(OEdEyefzXobI8F@A7D;dEy<)?!Tw8rIN zE@Wpx<@BmwA6$e0Omr!+J9md^8AJNwiNOL(0_0l9w0<#COtCBFH`1%g!hA13 za4$rbICuRB0dadn%=Z+BHoNm{iU+hUHejyItykLhYpb0*4|@bzj!^XEtyOmLk#bA4 z0Ia5X$Dg(DJ9hLA3j1NQ!Y8X0y9hLKPNvQQO=y7$LHL1Y>T_DF4eMzk0bEcCdvZ+_TyGf(p49hsG^)U27P;8 zE8b~xQ94!1+x&CTu%|wMB7k$KOHC$nh3?L1dG!+CM~IFWoA6=}3^e?eT3L6D;p7n+ zyPvM*3jmYlw(_!6=jJ#5c-PcBnmQBP-4 z$Ru$XtN-NS$0p$OvJP=6ffn}0b=+ImjZCQ)&wc7xS4v0U>;s`ds2A`~V-qnb8U_o0 zyVV@-S1vxjjhDeB>3}*H-P65ngqY@yipovL5^9_7z%>M~w|BMI5UDmrztY)um5|k> z1dNh>xrYAG7Vx(Gi6e!&A*^tO};KOUKiW&it`vwxMs!x_;@Sz+jFIP;b(-t{akOSllr78B;I87jM9=)Zw^ zN-9cnIJi);O^ak@RP7u<_26lTY3*WMF`l`SU=U#TOYL|1AR39fT3O&K!H)A6k7XI? zX6kSuncUPX{(Ta`z>c*=TAp9Hc(3u$h%ny*;E_k){k!7o=0j88xJxQ4aYr@pX-aHS z3sgfn-UP;1yAFa+($7Wz8_}}e3k_X$I~3;h1YkX`XtPx?Y8Pmn?!O^_=YTv{(V#(| zhQyjDJd4!5u9pfI`lcX9e*8mXiYzh@e&(9OhSdTs*3@0B)xQRf{56@bJKTRtzDmCo zm0QRA`nlfx8alXIm5P!C@W58aOT?#cTmrO$B5QhwO8a-$54RNzJYn(Zh&$xc$Szbi zk%CJmg2oBl!DINn#FwM&5{a$~MpQNfS~N*NHQ`Y+R|MNgYMxK>;yZjlerdGADGA*Q z%zIw~2O8N4=XyB4Y-v>+5{E=`$B;0P>L*~0;X5?!0Z^c`{`k&L7&6kWGbX)@m)POj zyl2=};!C_xuImpLzO&Gw{lyR8RBEbQoz$9cF@w=>3%i2f5yfINYRCJ9PORS+aYegh zG+I!6(816E(UC6oL~wmL;_l4N?HtZx0dEL?Ey|1N8Bi1MuACuK%ck*{FBtScJ?$ih z!LqVW>k4&;FXn*n<)cL2N#w+pa)83YFhFa_eB-S#I=A|BPs$(o1lTTRP8)M-(_ajk z@N9$mR^ef%_ObW79ghImhF6r8>;RyJa6{ZC*b&X#r|Rvl>pYvV2`%0&=lDLsy+Shw zLz@!L#+fG&ibsms7fi<9B%j4G%E0?Er^jydg6+*QwWPvz8TyI2V<6*Ze~>4%q&=)j zeIEkbjS3z+cU)ZHpDb^DNPZ{)DFz{)P$$}_TknS2o!IpvcI(CyuDD`yg}pkgmFL?8 zA0tCpAqmO^(YN@z`N!5sHP*i>{kM;c0)`tj=ybbH)JG6NA`$M0Clf^msepWZI3dZm z-8~JysWr^t`tZqbQBSl@Qa57^ME2aombc|(h+h6z_`q)Eff|7=WAffo%wQpdY;=p$6$PHksc>EMLqu$Fk@MF(pS# zjKg@2_?yYgEiAiVWYANgJPg+ryn=DI9bds8D8YNuy<*OuV+PHKv%D&gVcBSi=DLH0 z`85VuLK0*sCJSUYKTNL>p9HTiOJSJ{{7`$gzfu@j*#>dOAx<%s^4qHLlilxlh5Dq@ z_$b2{{@K@MnYGOCeyY|kLbWs8PLZ8E2MD#e;rKJ_@`p9i)iVn^VWz^<0>QLIDI(5X zfHFRYC?71y@kwXY+`0oS_Zfkiuu&4)sN!j{O z7pT^icQ0)@h{=4GtZ{ZNB?jREU!^zCj2DQA=?N5={l?8nLoKht2a0Ah@V3lfRXCiHWdveF?j-{;V}%wP;3Oy(oZG)KKp zh@m_`0|T+?=JbT8Sw|K4RrpFv^$sV>s*}HI1maJ*sGo|%KhAQ!VLw8A}D2n_JM!W|B+wZ<_c%dXwq;v?Vr! zTl*h&4z8|7@?rpNRZrbvm*8)pE-C^0^D%xHEEQm&U~uvK`SQQS>xd)&`-64Iudv^d ziz>Z<<_%1dgb4!-y-iwoj^lSk^W)f*iq>PQL{|``9-T-$fiue0$H?z;X2zp6oK%#l zS`|3@aE-$LZDTF$(X>t+F?filuRFxl)iVXU_4VPS23-RSlolW0u-7t~(I}w8BsVl3 z@LV!jcp)~Y?;plmXNT+sp~@J{3^%$Esy;t1Ki;lWI%&>Hs|@|{;OTIq@z}zwIp3H; z43F=X3j59*=m337r}TWkbAY3Y4(9IN)b*eZ%`6xH zjoNPiQN{#Kcw$u1p^z0&`6|&Xicc)^guG*Jy*^Ia zkH%(CLSXVPbKxL1P51}ED6}v1}S_C7K(zQ|F>h4@g-#9)6J|It4EYZ9Dc7xhR z2lgAq@8;}CRjmaiVZ#ZD=^nLFZj2CScttvZUqf3W$5NRf28>7q+XMVEQ@#dE`iVZV z?Fg0FBj!DCI-%y;;r5kxLAbMk8y{ip;)#)%;24$-pyCmT|CzwqI32_3op)Alj?(^B zv3s8GuLm_Fg~~+2WQLTqEJ3nQAM61P$!gnKmawn`Y=c_L)p}RJpEu-f1JVx-OVCg6 z3`pO|cl_iwnM7~>#znl_OtDLRps#*=wzr&Cg7Neb0I5C7;gFEY_#<1CMxu&8x-X{a zZnoicxbZCjju4JN0lc1w+x7tL%v@*}0Ny)cZ@2gx^X#I8#Rm!jutcZL9qIm3o2;!x zg93ewV@|}MJRBig6F5I;{U_*OYUCa{G_v_HOkqYEdM<}~*xEeAS-UP%I?ys4l_B_e z!){N8BUCPl=~QC>u(x!BFg@e#Wy-b_yVU~mN16vMEhjsZH)lVX>t6)P&tWd_CaM{> z0@S+Nm@3qP9*-*2ybM3{dx{jI8FjgS)q))Yay6)QVQ(b5_H2)Sh(=ZXFi90Z!-~WV zaw2n8P03ot-TeJEJZ|q#G-EPSi6fGwaLkuWBLtt8-zN2!cle3AOp4H#OQg~lcEx~` z?+>p++)t(3yf3Ha8pU=sv?#OgVGpibxBY_%ljg)s9nnWOuyBW{xlya8?z_TQF~u{m zR||H!cAVEY*4@l2w^>b^3+H83;|C3gKn5iR--OwmIk~!0fnw?9-%(;>T#fx%4`WVW zJohW1?fIu{efI4PUq1AFPP&XO>Q`tG$24tT!O5?J@gu0o$N{Z1&r^!%<+bd!+`EF6 zYL=-P^&Cq&RmoG?QCEwKF!4ptJcTac;bmR566Xx zv$kTb`C<&w?m(m(ZSg=%!1VE{Mz#FU*`=Q5ioc=anBPBPnMB`D(uaq4GT~pcHJUt6 zMu2n4$2y6z@(UdkkrPfc>qy!AHKEo%zL4GZQ}yaJoEwZ(A>Z!+l-`ooxcattw;|?i zc=QxR!=WHbSiJ6Vc^?_*mJIyytwVriL%@_T?54QU2>rcLBl*1~w;(*~4}W|G#5!lCZlGZPV-&%yjU6%~NhSKX6z40UpLhVukpx1{X$4h$7nyj}Y{F01p!v zeJn1FgF(GZ1bb%e$SE)C9c`~YAoUT*2W(yMSM>APG;e}MA!Tw+t`x5ih~5QdyybZ4 zK`3gh$;rGQ1&=d+wEeHCWa)X9uL056J2IlF9Gn4jvTfyq`mx~4)TmwJZ$<3S0~Nci z-ETjZ9aL;|PB-)%nb=k4{;6sr)Y#1Dn1yNc^vN`fIK1rKIJNgjaWr^nq&YKAtkOWdGBM#lfgR}xNPnNz@ZS)LkoYk;Lo zm9`KC_YE?xX6afkd;6{~T##@I?^zulufH z!o5zdjO9NWd1qWKXAQeDFg0BGalKj2)({N~aWpGz9VqHuIYhP&#$)d%vto=w%I+5f z(n8_^Z6N>}Ik42iAUG5=?Syz21?4#cJIoZ$#a_;@x>zYSKRvV=0uI8s!{f+O5a%kT zF|+_0P(O-aQISnglBF@PP3?(J$i&&(5w;D1C?%5^79Vh+V~$hhcf7YyZ(ley2msR9{NF*DrK0Hu-+Y*@G|aLSrJ{7>=eV%HnlAZMH79K^5pOYd-9k{G`l8%>A8z zmcj@i-4U*JzO)|xQ;AEK)`snM?MfyET+4e%1S4CxW~6mzY539kd|1ueAGWlBlXS+) z21GmkcQyDQf_n&mox^iwbET#_21LA}Jd)0nvW9NQ|I>t#gMyRaqm2I!UD)c)4rq0;o)Td(}cnQ(}bNkE?j=c3>3bCIiYU- z%SyKJFfL9`sz1??<4o=`KwWDvwPsX7k+%p|+|9)5<>YYn5`dt@TPZuU+n}2$48#Gz zH}a`uq^(;GxZsNc>x!{JDFK_d|05NX5W_1zC?K4%N1622&idD0O=EDba~i#M*s(WZ zwF$tmF~3nV+03Yo)1w@>Uw5qvnRFH$4H)5~aXVr5v%kz88#V*Gc0;;nvxVM~J`>OU z?k#`mo&cEUk^W#aTh1MmE2WUHEWn>tI+>|o6|0wCtVAi8?=`bQlEf%48lH6bcJyWqs&D@d7c8yA=3 zHcP37_lCJ(eCn_!wv^cHiK_9{LFY=~lJ zC*QH|ypgi%s;^&tqqV8iO+iVOPZ9bx!8Ytb?b~8gceebqeXs=pV9G?~dXciQ7TK2*w(45vBq<1;p1ao{c)A91Kiw-97ax6D&HOv!A&vZC}&n= zMM%!sx&0r~l*<2Anu-h`SLgdDP0^ygcGpWQ&=1v&=0fum{QqlN%aUt5fu;B~$ z4rs}6f|VLi*T+aHmpmLE^8u;iANdQf0cX2>a3NJGy^?aiy8eZr_F1{HAmCN_j|y>O zSt9Nx2ly5M7}8-XrgJ{oO_hnCLg>osfd-4c#j?|PV;$UL6L1g*YMN`o!*K&$zqDs7 zQx;thqw#PI6PMxY#33aRg#*4PJ^(B<;(sv8%6Pa z_?aqRG)5R)T8aFh3d+C+MuLbGEZL@T@kEnh=@yxYR?GdqnGK0IJyW?88isSk9rKw( z5H3pz3h)tDbO?`A<_l_>QVMBW(TZNt0`?C`MoLV*@%%H5NF0g*!0+Cb){9pNSFmIA5qCYxO36=mx{Yw&?*KZmu z5W)EHWUOev=p)D(TQ&BilDfzdIiOG*QsyaS%ls(ou!?4 z%LR2G()KqEtKsA4ZyMj)@OeimkL#+VWjd9PzxkLa@g0>@Du)F5)-|cjhzM!P;_b)$ zFI)>cuH@_LD34NeE8HVq_U+wTce*jFTtN9m+54zenPwKGt_58EDYY#&lw9S~xV}$s zRE>9R)Ygd~%GGwfno%R8u5$d1W*%@lPq$M(L-pAsxlEJ9vHs@;H=7`gTj9N`;rXfL zjPMg`V%*sN0D0FP+}r;N^4Ql~;Bx<~l}TZnez(Y-v$C)PdBGno-wpEB6^D*!jDO}A zT!XxKLu??jSCLX}8OLn;ui?Nujva$h#qo^K%@(Z0EJF(LEwCNP5p{Vl&^#xaq!rU< zhd+P~P9*X6hxI+WVgPa+iHtw=>8s^H^n`fomK8>KWHpuWzGJW(pnFx|t5HN^EESn# zdpR(cULtGn3N5(-W4lp!-Ox)()jiV(&s?*IaghG6n0%s|8oiQ~@)`I9soUn+%M69HlYF$fx>({>_##zbIWwJgi7PJTP;Xe*(Ocb=G};) zg{o#p%2V#HAre6u*#*ZWDt}UN&Oluc-dR=Vv2AbK0mt>FBzQzN`M^z#j+c9G&#j2P zG2_G?jc~tt0bG2MRy-N;fd|iF$=#7+;R`vl28c%e&{7SD zPk?eaQq;nWAcp=tmjR-&>=NALTQN8$k^%Lztd_bcseN`0VdM|{#hQtQ-JqT)nm zN??TZ2ug|B2I@C>MKVaRGUp&61ck>~Rx&LplE3xgxye*Fx*)R#93FFLjN-m!A>tMf zaFTX*xT)1hi{!GM^&4%q|7HHn*R%u($zPB+-tqsLFq(!iM=irzc>C}LQs^HT9`C82rcxCla<5yu<2EmeKh%Y-ZNN9k^vf z@$*}a=D*kbIQ7#VyB&$A7Vf6A+TEK~Y2L4^TtpV9h+|e+3~Q^&Xf&AE+9Zv}uaFK) z4k;)=d(`Q($4mi)vH=MXL8INZ*sI=-$3faU%xAnU*dM7 z);L?f0#)IEf2kc{g7hr!RiA`*RK#2L+yctE-wJoL66|*py{tv~xl_LtcMDh;fH=i~ z9B7Kb?w1r4pj&(W|H=8;gLdYsGPK@Jwod)70udIgz}ibFC~Im*2I~0U2xC9FtV2p8 z#-P_3R-C?U?^>i|yh9Q}xwa1yGixt0g9?)KvfKUv67nMNS`g{~`V-+WNGY=C8ZJC1mS9V=BC!z#x zCxmawLFYQa{6c_@VDVN5_QB#W5ibgcHMHNb>Z3Ai)CYn!(AD}O^0E@c?Y4t?c8K5= z?po>#&$7(H(qOodLT=*Km?J8&L>50U{us_H6_ekScWntks4`lJ@F2>Rsra&ZlEmTz z5yW^*l$unW#Bck*Adh~X>~0?>8!d`R2GwCzEf#qj?SwjsjHR@MI)juf#u?z3*pa30 z=}M~x*Jg0ci6!!Xv`9}~>@*uI%N~{r7m7PTy&FZUA}n9%QyQz7`!7Aqqxr@N4rzHF z<)2w6P_zqFK9T3cX|kv!E@<24;eN$(E025Du9W%nC?FuxMKV6$XS`u^v1*ixA>s7Q zM5Zg^Es;^?>~d?l>>`Ddj{PMC9@m(|t#1Y~Rdtt|=m)F~*d)lFmm$3k#$z|wucIIV zu?g$y-G##^eu3Hr=K$rHmuw1*r6vf&ER4-XV=rAhF+{YEnj|?idEOa z2aPnCl}OHOR7KkTd`Rpkt$YMeO+PK^6Xc%>gUiFKT@|-Nm4-sm04?uvfzM5R*lncV zfGt{}A|xy(A~{n-gDqb8rSQ=gm&j=ZCJsH{*fe>#DUwhf7-MxRJda(pdg%~k z1u$XGZLy}>x;Hm|@VmO9Cn?8^1MRmo`J>f#I4;1ZDrZBsnU+x2@yQ6{#;}u;ohem! zs8FC--wzVW`b_I5!#ytKejsIo;fl!n^^#(N2xdSiukz z4O44986q!1k&vtJNE=5i`7ABV8OEIfX}yXJ;rS&)ffBGK5V&Gy&MO49@@-~R_Kny) z=gRV2jgw9*_AXWj-4m9%+NuMK60F=7y44)N^qXn^j_d>%YecP31HI*6k#Qe|re)Qq z+4Ze+4;sup^K|FsXGiA@+{Mqz_S>db4_$iX5m#r3JR@GWjbhSP&|gA6y1Hau>9aBo z^o07b$R-_y4(0pQe~wd`8u1SBmgbsS0lgrZT92Bfs%>~-&P+9q+#Yzdf9Gkq#7K`@ zo}Imqv|%}FW8uoBN1p+k)Q9U__arKnNF_xw8F3_}bVyt}O|VA<`U{5qM@+aegnqMX zVVPFzoqlx{7gwd~iv{q)u5~%3>UR1)E+oIY*(+{-j-O}hvBzIk);*2L0{FuqOOkXO}Ndaw4YaCUMp++WauT!XUtk=!FDF{LI1 zrgbSCbQa?0@ip0+4i-!{m0;_5jBy*CJ-3>L)g^nf?7mLz2w0!zc62O!a>rkb&#bF^ z^t5Fp#n1n65G-$r)1px~83&0uU%M!5GOPBIDPWTI)i=5IT}g87_M$Zntw$j2@)dCq zNo22(F-5OSN0VT*k-_1AeRc%@Z|I)S7iM9KZkI)gLP*?GZ*O5qDZ~8NQ_dhj=G|3D z@!*6bL~((A-R;#390OUFLMEH!bv2uff?a!i_~+k3LV{ZN#$V(S-&f&bt)(?G9cuh& z&5XRr2;|~j=eK`+9g?OcsBn%XGqc%tsAui^! z4M?zw2ZcZWy&>2<_Fn^OapSsp zOD+=#h((r7l@wx&-jBkRi*OL1rUwar5p>x=V+_xcTq0KIyVf)WRKH0 z7$b6W3&jWGyH!?1KZG|r6|%Loi6nk+HSyrO6Lqt08Wh~Y3u)mp9%iMULEy(-06sz} zrf@A=FM3520`#|rYj&rOVrsKo-xrdvPn)aUruTZ9M>@JEr8E}5%A#x7>&D&Kt~sS7 zT-Y2jBp04GvtosOcPS|?L`J^A`ss323BrBit+UPAXlmR()cn=MJU-Oc(gkb6^T;xs zqP`Qhr0?r^WlF%!fxW9ypn&XFTt zLq^nrFdQ?TZvRbAv(?;&#%qM#MqMhGod06Ex7r=&bx6S~Bv!8P|2Huaq2AHMn1_i3 zllu10=zp8`?f14q9ROL^I0EAF-;|MI-XeFAgyp(-0B5d}$w|HDe7^lU+?h@H^hC%* zv(@v}Ej6yEpM+2L;5--)^seaUjfQQ7P`Vj#7?{hBjCP0GZ`Htxvt35LsQ1@r(*dlP zXI3pejxU~u&gqUcU2P_{)L>#+?NDCidp7)oAh0n6mpOt0Yszx~1CMyUg|s)XlII`zS%gXn5-EBCN&6)_)@2nlUa$GxGV}j~YWE;~6(FjxDgR z%|kH}TQ~*2ZcCFpqN4BmCaPB7Chu$3-DF8dgvPC9Y1osB9|!mzStqxcmgDFYEN`R#g$);y)-OPD7#9?O zGkht7x&nZ{`MK8~{>F}t0bUy!A^&wjipqiSVh%%#wY+CH6hcYYi};9TTI5XlRO^}+ zYvOH&EvAY^L->mer|N-p=#4PeDOz`DVMSUk0Qe<@80{>Glyc2w44I%w?-`J;m|_68 zNHB?dWZ!|us-!c>8WVUhSGHm`e#ZnPu$?0nkCI__{sBwVh^(++>gAxoMTfb2l3cnv z_5sIdE2beN{seVR{U4;^wT%Nonzzkrm-RV$BdgrK1|Iu+b4vqvyS3)sc1JP49H)0r zx|!-B!nf3Fn;-?Ls_O-n+@p|^YCo#^sxNwlRcLOriCSv)8UNXn2ex(AuWmNKgm_a7 z?r--8icoVfLwWS9#sl@ps77<8b1zL=&^HhHwEpKz4tID~2*rOdXtvkuUXkyE+oW@T zn2_M};sdSaY<6)dlmuK`b5(8R5Mh6Nt}NWlL)T02-5hDb{5EYa@4H#K+Yy4`eo7qu zv(sy&vu8fxPcJ^dRyo}%V_l7OHVcma2G44u2jiWl+3hUl84@N{eKBm+K@DKsmt+deZ zF5ND+v}_nGitT48M(t;ZNMQ=HB(*T6ph)lnLAwfvC-y>aF_GJL3(WmnIn`UOLO@%_G4M<#N{CyK}}KZ7#4r1QIbKhP6?gME|(AWX&qBnl@an zZ6;0*6O3yE(k>4cYN>%O&SY+dAKixmxboRHOHP>@R5jd?M7Ihl3CR!JUiS|gGYfNG zq+QDzQ`}`V1F1S`!N0gb^TKu#jCN&E;}N4_Y^7!PFdDh>mM|G5Qew3yV%ax^ zuOGAiyZ+jXaX3@Z9$O}Jl}g}Y&RI+ks^!HKzf~gFnIl4V3qqaE>}eeEish7DoQ|8~ zJHo78&*&sONS%Ov@s52u{Gzog^Lr9I-0zp5HAv8OKqVGgaP6=~N4f(Lq#%gA6NkLO zMgi5(x17PR5x(q8vizlHY+WY%OI^GMb$O)*`vk<^okE}YrQ-Nnx*Oz;OzJD0)`kvV z*0tCex!n@q4VE^W zS-ZWM$CSG^KyT{{VVaZePfRnympVCVso;~3P_Yj;~)0$zUPvx#}35ak`AN2QPX7inny7(eL|4TvWl;HSE$v{T%@Xe-lD zMk9n)Hb>Bm-Xg@+c-R(zx(Xkacp9xPHLkQIJ$%9K5W!V^cNdNtkSn)Y1KH6F3aMe7 zskhI8#GQBE)dVQ=B(nR=8xi!L0~$<+lF?r6}6lI+zOsD_yNn#aQ#mNhh2-uvz5 zU>|XwFxtOhWw&(UjXZ1uf7fOQr&$JGOdYR21YJ;X4X|~^2d4orke5CFJ&h&`=}7KV z6m;MRs=4dm+elR#^2C$uJ`1+`MGqsnf|()sL!L_aFzw-zk+v*>in6zdd`caXa&$0u zy%dUwyw|UeOM5Z0V`OUlvpO&pp2SY498V=B!7+@nFv+M!ii=atZU8&7F530M-VViD zA+AMg>+gEbT+UZ9$3C4_YIF-QF4^)sQ`Cb7<#sq9ys!+){WDP&oa@H3oL}~#(^1Pxx=?vw9B3HVcYljRdnKazBvIK39)PJe{Y0FtM4p{BbN$c{{$Zb=kZf;&N zPX9f(>8SUfEv-S0S@oprMqr$6ac|;!T3GuyC*N6j{MW|pgkFg+`Pu0xd7#0r}uo0k7^R};i-gOYz)0#R5~n}mN^MxNdG^O>gsvsJc)z zoK|+**XKSU3jF|-92&ZBi3GLev0a!#w2V{*ouN971T#LM+^}zr=?BNzJE-qt`sJ?n zlcFqALw=%+E*x|PHX%$XbN$=+z&+JOY2&c4J;9M*KcfAb2(Pg|antYTRcPFLJmXm` zZxPGuI1W?UV8OB!3m#7LYGT9cmv@6zdu)sUtI_8lsU-Z8^1D=BWna2jrQ4d{DEkWn zFj6Uok#!A%r(5dnkt{G@kP1=yMfM+Le$U+~G7k-=6k|6S;cw=Uw+=i>!d}#*;O$t><@E4it8J?J#&;rWgpeIH{%(0aBOL zf<~HZ%Z?ep>rctk37-O8PX)bU*1n_NnFj~4?{Ht8WL(VTqJ>lJiz2YfhBMP9)>o3r zXw3}|Uwm7YNn(N7FSM;@;MrPd9s0NZW@=b!!+uf9UeKmax`7x8KY8931y3DVz>x1T z);yRoIf=Q6&tTTAKpqDE4}h#y+x|kuX=@9g5?w>nioDYl{c$vjw+NcReI*DK9&bXvN8M`TWffvW zau>?sD!K)NcXy#b{|E2j1OCCg;F5pvjsj}qKX~^IuH%33?sD(nDZMsa{=ZZD8vK8! z^dDKcj&(|=hGotPoYi(i_kv7)5_o)bN_`ToQn4Ijz^AfS^~-Hq7I)#k%2J;j;33`` z@;^Fi6ofX>mcqI0XI?{&tJVjedtPX~&MR`Lz^krk+^jHY^_dHK$P+tq7k%F(QWft*km;M_TANZ9qP$=ahTuEW! zpG8Gj@<_y5IMSGD@Hwes5Fo@PstaUo4oa#RgSMvdS*VGwLEM?euh?d*`+vR2t|`Ab zDD>>KI}G6saS;ngllR8u$GXcpfAh}}EqDND>L;|CILigTE<}`edU8~El3awEcE`3e zZ>dyoNWFE561gUL-jf4`zU78Zy$ymQFFuU09XzU7ACn+u**_<)3S=xUW@h0T=O8GN zyz0_+yU!|4{WkjZv%zbnF-z5K0VK>6V|=XQ%G2X@6GrOqz6{CN?fOOxR& zsYLhc&sX7jcYRpi*?lTC)w0s?Cq14ix;HAHtAk9R@-fx++&;hK##s^@n2Q1}YyQ2y z0fU7sH8oCv%!Zeb94M8LgLo{cfaFgoJ+z2a{4J$v;9p~@SAJ4u+Yw-5Tq*e8R<71j zK^5;(|3l4c{t4bV;L7-)k0vptz1q9kD`h((wt9j4+NnW#h44EAz}!ak_IsrIU+ByU zdzV0+`mub-A|IDI59UvVJ?DV&YB{30YShnfZC-9fFF}>bB-HW(2m};%;wbveNOT36 zvAr3#tw2-<6muyuuuBqOsW!-4=J63Kp_!*T5nW%-a(zs7fC7eM z_{qNrwGiuXo{_Oi5AA}D8UEWB{f}vNInOX9e=uMFC)h8on1H+HT|dus7Y&SX)LdWP`fPgC(NCjS?%Q?@FM4QawZ0G^{B#}n)fd^x2BIE#@zJ}lF7i-(D=v2+W#(A6{fL)F`154uHDo#GDqFRV zgrpV7QKplc>L!_C1GVA#EdE#r%xg7EktMGEWaAha1Y*L`rmEMZ%H0Rm2%`S?p+@|8 zC)9uD1Zt3kJ5jvv>Oa0wz;j+kV-N}enTHEmyUv1=qK&eqQS~L?oOrW%LaqdWn<_o& zmdXF?(3&UxfBHDjXt>_BjSohR=%RO`_e75}dMAPq zqPIk^(L19vT9g=~lMp2m(M4yB7DVqYh;EV$-Yx%g&b!|8;r+7Lnl;aQX79b$-g95~ zb^Y!jf@Hs%UUO{BhZfB2(DV)Ca6TLC(I!R#xgP}IOa(kh;EjVTHA>`W*cT06SFfkK zN-ZtQ(l}W4erS4Te=5wQeFpXn{8Il>;WBuuY-<@=H%atm1?ay^AAr1`bX2Y&rIqR} z=@kI9BEFGiA=jdJp4{!U|1h*QIm|#0!Lqv1Us8bd_5d;!gSxaZ}2oKO#TeE~JG=uBz~NzO>E+

IQTJ^+6$8wX8 z2Po-_IA%0cGpU*l9B3trI|Sj(u*kw}_bIv6eD9fnQcsukKZ>5DyvBBviaO$2Q_ZAF zZ*N_n%c3zd-R8qhvD=ddtGAFN0drgshNYR6fWYXQ{mx^6-27MsK9A(XMtdW z0YB9n&#FqQ+(Ylp=#_8ZB$XR`#q^X`w(G<`*(n45j?Ln;nk2^QP}d-3KWJ?>#>e~~ zy2IJM$|4O+v>)`QrE!X!_6kmN9oTad{Cf`{af1XL9=bh$M#Gb!GKZ;CnW~Lnjo(gX zaKBY1PR)H33vfQ zH~?a`klMj+uo-GS{C9#vSz6rpWqkj#H}gb>afefu^CEk}2J9TajdTE^`X!kmfUIt!d-)(oKBNA(9{z&awEAtw*&pf~u}}crkNem@Ara zNWH|nY3owt8lNdL+!+y3>`6Yc^rs$5XMf}s`3?7@n~$I4az{VXRl0~DHmjVltSW}t z8i&%hTKlA|VZXD-AEgQHZ~qn}9>B*TdbaIVV*yV%oQf zFKG?FYY!YU5CM#JWnk6&DqG;_;{4=V{i7l*DHXmk!B9n6f7@mNcFGhT&ls6E@67F> zR3FUevsf$7zvA#vV&*s{l7`vD>S=5-#IS!>Ovz&w-+g4P6k+x+_mMltZ<_0%iFNO> zh?`&B9hW9;`;CJ|q1JIk#>i~}!=00|*~;5cW#KWQ8N=tZZ=hmkJ{!`)j`(Q>S0yD_ zOpi_~L>~+PrAv}miZiZIn_Z0kKvToWFP(w=b_FUyrn$aF zB_lIHgbd)Ts4y{`(^s4WNmdnlh~xDr9;Q5tE=7F&u*nEd{6(@BJK=!;Kn`ekP}FT$ zrgOwKD`I*r#*&DtorDC#TH{7wZMZoYN7QTm6ApwALi2&?B$zgxs~MZdvj13=v~Sq* zyu3Hl&q@$9jv=>9&f3WSB<-6BlgEhThjbwLDD2mzlb8pPH&aAsxEOv^rZTf&d<1Q= zo6mMCaRB032GI+NLs^KeCt;^hWI0m%mF((m-pt?I=quVk0Q(VY{RDSp3Di9ag&Xz3 z+huHfCch@j{0f>Gu|O@J(SQ+fV0{GPkgtZ$Jl`Bk%Ekx!$lhVcbi-u$N+76t+Hq)}5Y7#fLEOV;qHbe1T9U&b6de(*o> z9pS^^e|Kyqc%{%1+Huo#_m@bQRkOWtN*}_5&~azWlJBKB#04Ad^~$ zO(*mYk}tonYG;@HWY(aM^qkY!_m=Tzo0k{jdaW1kutP(N~f!~ww*nN%Ci*UN{ z*j)q;dTe@+-b>iEe>I6n5uXSwagM#r1}IZt*pT`7LH9gjX0AA2Y`eG2D1Odb=TD=N z7|Vv`c!g=%M^}m!raSlEc<;a;iMEnNZF}2SrF~+p^cVCfc2~EhFwI!osm9l?y~onC z`SF6M{QWa~Hxl39&2rIlBTjIVA$Ehnt!QVh{!1+*Rm<7Ft-RjS=wDR2jpy^^2BHhD z-{R+Lq$ERWtY4qav>k0=+LD&VZ$!Jl8r*g+T6j4mMO&v3pE#$c0*WGKW)p|m>0it}NvUIF^=6~-r5|xhaW!v#c)c?{buM74 zjxsVq_iaK^N%NEGWv;qe&!s6k?tgigdg#7W&pOz7*Up)sbjOO#z1Tl&UfKCUhvf7* z5jXKSJwWK_CD|e>Y%5gi#U2g`8vYeJ+t#?WWA9f_LRhFBmt(;f8DVjKj?R9AYaNT) z9QLu?kNPs&`Ja5&tRBFdvHY0DMwic{oImdO&Yi1E)rJIhC$kPJ8m%UH^!KfxeR&#{ zR8My4ReAbG^M9jZ?o)IVD4~nA#b`Wpji{v8=+7W^Ky?(qFb7eY6UDtL2GghGh$iKl zMCO#`noKZ8*A8yQ*u}o7aTn0HYVMtsN|q_8bhvlPENhKLB1&Kh{&z4e4s$2V=G8oP zSX!9Om8lJgAhiB^xBLfSdn9AL%9E%TE-GzMAxN!r!W=F=1Y>!V5gxhqIv-5(M`8rY zy|#NtJ@UF?gshN1&)c+QF60GgyOVp)KB}p>$GF1 z#95dUs`RHy4M-*um2)NuJcmGr8& z;Xa8XlRXXx1H)}nG2tlrhy4>%|V1C<_inT4cW;y?wXofPc-D)?PSfqsp{KjwV8q? zkJO&47{lo^IrAuXFyX;(7bmwoAq5vLB%^EsP1FGXQCy_5@=bbi3#)6iUYp0(T>Sf9 z{Kc-g-nPt8g<2_Vl&C(nbl@0IV%j|D~5*&2P|nz)_3$Hy#Rg zT(2d2A@_E6wf;phuP5m%0$&pl-E=K?S}*^Y_wZn?4s#eYVBA9;@aX&NSeiU()AHluMP!rc}l#o=e+V}7m{s~|{nk5q6-7t7<5 z&6o`E1|JGw(vQmNj|oSk47&Q&06V*OMQe@gG~n4ipEpBO?r!4_7J` z+w|V}wSqL>y`Kxnv?Qsj?!x#~tv_U*c~f%2Sm1s_TL%Xzvbj2n!~C z<)pSMhK;{9mA;kV2VjnXu^$(OQ*^Qi!uMDUFSlb@;aHS~TSR9CPv~_br4;=_QD4c) zn?$#er7nyg!w>8p$2rAEOlcdvgJt`D6-uB8fVAn~e=_Z&yJFFnJ40CBKAgX(`aYI&uP`xj+w9F&Ly2sMc3((0E7@k0h&^W_->{$n#&MWdc0?goP5 z$DV7%49v+`GlR+cBCAhE5{{Fk_siV}FhU|x;u8zvgdIlwwk6>vyFWaSxS>82c%g})8N zGg#18yu~c%wfmReDV(0tQ-~EFm!%lwO9VxrI8WMAs`bav{z-Bp;Li(21V_2k9NRfz&u$-lS$}~scJv4m@p-ISUo-_Xfy~21DmK?&VDf_n zq@IaKb6Y-Z?N6f;0fnAF<3DVhynbEP2Rt4yp19YAR{aeo~% z85B9lXX$%&|Gvq}lbKs9*_gLWbZdEscJ`&7JE8Yi^UvRtl$bfX*$sOR=q&`|D264Quc)2yp|T0Khqs-JTP5MR0qhzSf+8 zL(nkz2($(X;VF91R`Ij+Od|UwyyKv+uf;a1Kre1$D|ey?ZP@*|VxT=|1q`ZIw3`WBN!7rm@{AInT7NHVxNa;pP~ z*8vQuMuu_*j{*C&swvvv1|}KPO?(7i0lJYH(U4xmUS2zQk8W%_xQ zxhlVO7$$CICM(~s!}5+2$fUy`Zf8e=tp(TTByKR;xLN#`>$wnbvNfV~tDE z*rV4|n|gs5jJGH@_krl5uJbciF}do@_jlh{<)wNoZ*7Tw!7nAya^^xp92&nQNIOfs ztkHY>Y-gG@N@@+_;BfpP>`pQA9dkRZeMjLGpk|9JE+hRH1^maStuP{E>*m)dtOXx} z*DFK;I=?2^b3@A(J*mP~{S>;=El~gLCv^J2)7vtcm;6Jp zzZ(;^aN^i><{`lhB?T8w3!j+*&2o(BEe(f0gX-m(VJ~YVw}qu1Oj~>{yi?H;A7|R( zUqI{K(p~fFS;)n=>R+g}SCq;aE4t8X>Grv*iXn<7iAD3Gd=&*-JsmDwKv{0xs4)Oj z!+0obOOSwxO}<5^$J~<_2eHf8{5=L;r#0k>zOK*HvFKq?C3R^Bol3T`ZpeBw_!H}0 zU+~nO&_#MlnfocMCe1rbvYEZY`eLe#nt|Tnwy;`s@4uSurpX1o|EZcOihF%&EAY#Pw9D!tR7{Jc3GybS8}(5SK;~YJj-p zRbU@3F#i;oNJkpY3|m#X**SPNQ|mh^1lT4o;1cLcx8*TGt!aW!ju z^eDu@`)&H+ZSvKY8_#ZmKyJe+BLa_M+n16Ba%1`K$KuO=P*FF@3w6Q(KPW|v+H{r1 zJ>b}C}5;$mP)iC{+@6DZ8|6f=_S=Qi%EA`TjJ!xiCQUi>Y=eDIDnuO}e61)D2FcyCp z>nObYCHx9kA`$hWzS7?|lDbsH_LKWi)&W`no2QeD6?<)I;b8d%^(bxb@RRZ^>)DLCL>{Vi(B54_dajOqbaX?lHUle-PaJ4 zuN>Yl!$6S=B)87Wex2gYRg@7-E6ZbFb1*LaP#9<>+q1~rGga2ASG4=ynB`S2FSZN%Q<(gPduibjN1sk9lQWp%u! zzYY7mb?C6lt5$Am4v5a$bmAEe?hmX?+BFSVydhBG*6x%jIB?LQv5&Mx_@#Mc}0 z4jD7Z*_8SYIdW6EVlT45#t@vdmtZTkD?)`XEzWgLixjA|xzU5lVOumf{*Z-m5XIfb zm@J_Z15|mp9|p9Fe$&Wz%XDw&Wc{XrwWig;F>Lu+2rJ=f?2*rz`Rq;i9H2T_9)j(K z{(uFZaW2c+o|P#ZX@Bcihmu8YY0?u;A=MHHeVc{~Q7~p;H+PNL@J*|8QEQ!G^9ndd z6qBMvwo@gXWnibE5-uDC!8+pY>3|(R4>%VkVncww`}As^810bEIICSHHcl}uH4ooj zTjPL4#4Xuxq}CMD=iVnS2T~Q`_rbQ#eaq~FKZ69l)wQSXozk(N62M{G2%>Jr|B%o> zGlrJI?MS1f#jO>SWY8hs73H?IzV;(pai0GN`zXLmC2i8e|u2?$b+CXXN5#Sr}=w zh3UR{kAc|5Tede987%)e9kvy^X33l$^cW{zP?`#+F3bO)njuoqHakC@Y?z&RQkAJY zRqF(XE4+wcCl02CPy-9OSGVmqK4_gdsWUvd0gQEE;z7ghxv-%layF1fJy*0tt98jL~;=bv)AiWJ-#;0=2`FX`3^3?TS9O0^H|Lo&<6qMFa5w=f9^ zfesH$xbW8vFkWj+Q2gA}eE;*&p7+lg9 z(Bada--%`7xan`jzI{r*oha(J+Vd!ux8aW;xNl#k8%av#_cN)|dTfUZ51{CW&kPB3 zZ-Bul-kQAvMW*H&O1&|8hR!Ayg7{3$mdxj+rvKs%9U!vfUhhh8YHu9Nx}mj+v;CyeHd1gbsy zwC|G_P-;l>Am_H?F0=Qm$9?XWoOFhUwPK?Pj~%)n&E4}~DXK=2(n9HOb7*Hz8DH%n ziOo7fNOp}kG)j$6I(dKC^#pK=5Vr!bh={_Vr->|aP(3KB2ps2TD@s1S zh-If^ksI86YRGGg)kUXU2YfM<s!~6*YQ7n+6>iQ4IFgZS)fAh;GGZXj!q--uBDrkvN&YT=x<-7jHk43)v zosab5rp^^7;gMy@Hx0Og<6}gv(ucsO?<+#RK;4(qE|Z(G(u1%RD_hjO%3H^q=&u~s z{yBqL3~OS+hcNLU;>X0aaHi-fI*tL#N2KV=V#DgE)vvO3$Q9nQ@Y%l|=dNME8}Tt~ h`CK1#Iz*X#v_m+9g5({KqEBy#O=bB#FYB>D{{!d7Xk`EZ diff --git a/build/vignette.rds b/build/vignette.rds deleted file mode 100644 index 266149b58de2afec293a1a3f87eeeabdda71f496..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 272 zcmV+r0q_1FiwFP!0000019ecbZi6roHGz^QQY+dhyMLhU9V5|5RToq-ban8hSOOav zr%|_ldFVmxC`QU)-`zVuzjybwB!tXJ7zbn)Vj8aSit!{P9+3sc=ZqPaYgIzm3HhCV zHr#fScWJcZhtdcoN8D$=<{Dnw3zJF!4xAcAOVzZZ25Ksx?1zb(?XXHUXLOHvXFph9 z3k)mDR*wWt!Tyo;-tN!{;W`*)0C(Pd2EGtTUTj#T$4uV1fPZPO=7PIyYkr zMq3=zJt%~^6k(oiZ5lW~^9)*$-0kK51AdQ6W9|t>BhgaYCr>Uc|W}S27 Wz$M3Jq Date: Tue, 19 Nov 2019 02:35:21 +0200 Subject: [PATCH 2/2] Several code fixes to fix problems identified by project diagnostics. --- R/script_disease_relevant_tissues.R | 17 +- R/script_integration_efficacy_scores.R | 12 +- R/script_modulation_score.R | 136 ++- R/script_tissue_specific_efficacy_score.R | 1018 +++++++++++++-------- 4 files changed, 770 insertions(+), 413 deletions(-) diff --git a/R/script_disease_relevant_tissues.R b/R/script_disease_relevant_tissues.R index 338226d..7b215ac 100644 --- a/R/script_disease_relevant_tissues.R +++ b/R/script_disease_relevant_tissues.R @@ -314,7 +314,7 @@ list.dis.rel.tissues <- `%dopar%` <- foreach::`%dopar%` dis_rlvnt_tiss <- foreach::foreach(i = disease_gene_list, .export = 'get.disease.relevant.tissues') %dopar% - get.disease.relevant.tissues(i, ppi_network, weighted, tissue_expr_data, thr, top, rand) + ThETA::get.disease.relevant.tissues(i, ppi_network, weighted, tissue_expr_data, thr, top, rand) snow::stopCluster(cl) } else { @@ -322,7 +322,7 @@ list.dis.rel.tissues <- dis_rlvnt_tiss <- list() for (d in disease_gene_list) { dis_rlvnt_tiss[[length(dis_rlvnt_tiss) + 1]] <- - get.disease.relevant.tissues(d, ppi_network, weighted, tissue_expr_data, thr, top, rand) + ThETA::get.disease.relevant.tissues(d, ppi_network, weighted, tissue_expr_data, thr, top, rand) } } names(dis_rlvnt_tiss) <- names(disease_gene_list) @@ -365,6 +365,7 @@ visualize.graph <- top_targets = NULL, db = 'kegg', verbose = FALSE) { + library(shiny) if (is.null(top_targets)) stop('Please specifiy a set of targets (ENTREZ ids)!') if (is.null(rownames(tissue_expr_data)) | @@ -488,16 +489,16 @@ visualize.graph <- node$shape <- c("circle", "star", "diamond")[as.numeric(as.factor(gp))] visNetwork::visNetwork(node, edge, height = "1500px", width = "500%") %>% - visGroups( + visNetwork::visGroups( groupname = "target gene", color = list(background = "skyblue", border = "deepskyblue") ) %>% - visGroups( + visNetwork::visGroups( groupname = "disease gene", color = list(background = "lightcoral", border = "red") ) %>% - visGroups(groupname = "bridge gene", shape = "circle") %>% - visLegend( + visNetwork::visGroups(groupname = "bridge gene", shape = "circle") %>% + visNetwork::visLegend( addNodes = list( list( label = "disease gene", @@ -514,8 +515,8 @@ visualize.graph <- useGroups = FALSE, width = 0.1 ) %>% - visLayout(hierarchical = TRUE) %>% # visLayout(randomSeed = 123) - visEvents(select = "function(nodes) { + visNetwork::visLayout(hierarchical = TRUE) %>% # visLayout(randomSeed = 123) + visNetwork::visEvents(select = "function(nodes) { Shiny.onInputChange('current_node_id', nodes.nodes); ;}") }) diff --git a/R/script_integration_efficacy_scores.R b/R/script_integration_efficacy_scores.R index f80ca52..79686ad 100644 --- a/R/script_integration_efficacy_scores.R +++ b/R/script_integration_efficacy_scores.R @@ -9,16 +9,18 @@ #'using max and harmonic sum functions (see \insertRef{Failli2019}{ThETA} for details). #'@export #'@importFrom Rdpack reprompt -integrate.scores <- function(data, col.scores = NULL){ - if(is.null(col.scores) | length(col.scores) < 2) - stop("Please indicate, at least, two different column-scores.") +integrate.scores <- function(data, col.scores = NULL) { + if (is.null(col.scores) | length(col.scores) < 2) + stop("Please indicate, at least, two different column-scores.") data[is.na(data)] <- 0 print(dim(data)) data$MAX <- rep(0, nrow(data)) data$HS <- rep(0, nrow(data)) - for(i in 1:nrow(data)){ + for (i in 1:nrow(data)) { data$MAX[i] <- max(data[i, col.scores]) - data$HS[i] <- sum(sort(data[i,col.scores],decreasing=T)/(1:length(data[i,col.scores]))^2) + data$HS[i] <- + sum(sort(data[i, col.scores], decreasing = T) / (1:length(data[i, col.scores])) ^ + 2) } return(data) } diff --git a/R/script_modulation_score.R b/R/script_modulation_score.R index 87ef250..a294f57 100644 --- a/R/script_modulation_score.R +++ b/R/script_modulation_score.R @@ -20,71 +20,121 @@ #'@importFrom reshape2 melt #'@importFrom data.table as.data.table #'@importFrom scales rescale -modulation.score <- function(geneSets = NULL){ - if(is.null(geneSets)){ - diff_exp_gene_sets <- MIGSA::downloadEnrichrGeneSets(c('Disease_Perturbations_from_GEO_up', - 'Disease_Perturbations_from_GEO_down', - 'Single_Gene_Perturbations_from_GEO_up', - 'Single_Gene_Perturbations_from_GEO_down')) +modulation.score <- function(geneSets = NULL) { + if (is.null(geneSets)) { + diff_exp_gene_sets <- + MIGSA::downloadEnrichrGeneSets( + c( + 'Disease_Perturbations_from_GEO_up', + 'Disease_Perturbations_from_GEO_down', + 'Single_Gene_Perturbations_from_GEO_up', + 'Single_Gene_Perturbations_from_GEO_down' + ) + ) geneSets <- list() - for(i in names(diff_exp_gene_sets)){ + for (i in names(diff_exp_gene_sets)) { x <- diff_exp_gene_sets[[i]] geneSets[[i]] <- list() - is_disease<-grepl('Disease_Perturbations',i) - for(j in 1:length(x)) { - geneSets[[i]][[j]] <- list(geneIds=x[[j]]@geneIds,Description=x[[j]]@setName) - temp <- strsplit(geneSets[[i]][[j]]$Description,' ')[[1]] - idx <- length(temp)-4 - if(is_disease) { - geneSets[[i]][[j]]$Description <- tolower(paste(temp[1:(idx-1)],collapse=' ')) + is_disease <- grepl('Disease_Perturbations', i) + for (j in 1:length(x)) { + geneSets[[i]][[j]] <- + list(geneIds = x[[j]]@geneIds, + Description = x[[j]]@setName) + temp <- strsplit(geneSets[[i]][[j]]$Description, ' ')[[1]] + idx <- length(temp) - 4 + if (is_disease) { + geneSets[[i]][[j]]$Description <- + tolower(paste(temp[1:(idx - 1)], collapse = ' ')) names(geneSets[[i]])[j] <- temp[idx] } else{ - geneSets[[i]][[j]]$Description <- paste(temp[2:idx],collapse=' ') + geneSets[[i]][[j]]$Description <- paste(temp[2:idx], collapse = ' ') names(geneSets[[i]])[j] <- toupper(temp[1]) } - if(is_disease){ - names(geneSets[[i]])[grepl("^[0-9]+$",names(geneSets[[i]]))] <- - paste('DOID:',names(geneSets[[i]])[grepl("^[0-9]+$",names(geneSets[[i]]))],sep='') - names(geneSets[[i]])[grepl("^C[0-9]+$",names(geneSets[[i]]))] <- - paste('CUI:',names(geneSets[[i]])[grepl("^C[0-9]+$",names(geneSets[[i]]))],sep='') - names(geneSets[[i]])<-gsub('-',':',names(geneSets[[i]])) + if (is_disease) { + names(geneSets[[i]])[grepl("^[0-9]+$", names(geneSets[[i]]))] <- + paste('DOID:', names(geneSets[[i]])[grepl("^[0-9]+$", names(geneSets[[i]]))], sep = + '') + names(geneSets[[i]])[grepl("^C[0-9]+$", names(geneSets[[i]]))] <- + paste('CUI:', names(geneSets[[i]])[grepl("^C[0-9]+$", names(geneSets[[i]]))], sep = + '') + names(geneSets[[i]]) <- gsub('-', ':', names(geneSets[[i]])) } } } } - prtrb_specular <- lapply(seq(1,4,by=2),function(i) geneSets[c(i,i+1)]) + prtrb_specular <- + lapply(seq(1, 4, by = 2), function(i) + geneSets[c(i, i + 1)]) prtrb_specular[[2]] <- rev(prtrb_specular[[2]]) - occ.<-mapply(function(x,y) sapply(x,function(i)sapply(y,function(j) length(intersect(i$geneIds,j$geneIds)))), - prtrb_specular[[1]], prtrb_specular[[2]],SIMPLIFY=F) - names(occ.)<-c('up-down','down-up') + occ. <- + mapply( + function(x, y) + sapply(x, function(i) + sapply(y, function(j) + length( + intersect(i$geneIds, j$geneIds) + ))), + prtrb_specular[[1]], + prtrb_specular[[2]], + SIMPLIFY = F + ) + names(occ.) <- c('up-down', 'down-up') ### Calculating the composite z-scores for each disease-gene perturbation interaction - z1 <-lapply(occ.,function(x)t(apply(x,1,function(y)(y-mean(y))/sd(y)))) - z2 <-lapply(occ.,function(x)apply(x,2,function(y)(y-mean(y))/sd(y))) - Z <- mapply('+',z1,z2,SIMPLIFY = F) + z1 <- + lapply(occ., function(x) + t(apply(x, 1, function(y) + (y - mean(y))/sd(y)))) + z2 <- + lapply(occ., function(x) + apply(x, 2, function(y) + (y - mean(y))/sd(y))) + Z <- mapply('+', z1, z2, SIMPLIFY = F) ### Composite score - Z[['both']] <- (Z[[1]]+Z[[2]])/2 + Z[['both']] <- (Z[[1]] + Z[[2]]) / 2 ### Removing indexes of perturbations whose gene signatures correspond to few human orthologs - len <- lapply(geneSets,function(x)sapply(x,function(y)length(y$geneIds))) - out <- sapply(len,function(x){q<-quantile(x);q[2]-1.5*(q[4]-q[2])}) - idx <- mapply(function(x,i)x>i,len,out,SIMPLIFY=F) - idx <- lapply(seq(1,4,by=2),function(i)do.call('&',idx[c(i,i+1)])) - mat <- Z$both[idx[[2]],idx[[1]]] + len <- + lapply(geneSets, function(x) + sapply(x, function(y) + length(y$geneIds))) + out <- + sapply(len, function(x) { + q <- quantile(x) + q[2] - 1.5 * (q[4] - q[2]) + }) + idx <- mapply(function(x, i) + x > i, len, out, SIMPLIFY = F) + idx <- lapply(seq(1, 4, by = 2), function(i) + do.call('&', idx[c(i, i + 1)])) + mat <- Z$both[idx[[2]], idx[[1]]] ### Aggregating target-disease pairs by max score df <- reshape2::melt(mat, value.name = "modscore") - colnames(df)[1:2] <- c('target.id','disease.id') - annotation <- expand.grid(target.modulationType=unlist(lapply(geneSets[[3]][idx[[2]]],'[',2),use.names = F), - disease.name=unlist(lapply(geneSets[[1]][idx[[1]]],'[',2),use.names = F),KEEP.OUT.ATTRS = F) - df <- cbind(df,annotation) + colnames(df)[1:2] <- c('target.id', 'disease.id') + annotation <- + expand.grid( + target.modulationType = unlist(lapply(geneSets[[3]][idx[[2]]], '[', 2), use.names = F), + disease.name = unlist(lapply(geneSets[[1]][idx[[1]]], '[', 2), use.names = F), + KEEP.OUT.ATTRS = F + ) + df <- cbind(df, annotation) dt <- data.table::as.data.table(df) - dt <- dt[, .SD[which.max(modscore)], by=list(target.id,disease.id)] + dt <- + dt[, data.table::.SD[which.max(modscore)], by = list(target.id, disease.id)] ### Rescaling scores to [0-1] | Setting outliers above Q3 + 1.5 IQR equal to 1 q <- quantile(dt$modscore) - outliers <- q[4]+(1.5*(q[4]-q[2])) - dt$modscore <- scales::rescale(dt$modscore,from=c(min(dt$modscore),outliers),to=c(0,1)) - dt$modscore[dt$modscore>1] <- 1 + outliers <- q[4] + (1.5 * (q[4] - q[2])) + dt$modscore <- + scales::rescale(dt$modscore, + from = c(min(dt$modscore), outliers), + to = c(0, 1)) + dt$modscore[dt$modscore > 1] <- 1 ### Generating final data frame df <- as.data.frame(dt) - df[] <- lapply(df, function(x) if(is.factor(x)) as.character(x) else x) + df[] <- + lapply(df, function(x) + if (is.factor(x)) + as.character(x) + else + x) return(df) } diff --git a/R/script_tissue_specific_efficacy_score.R b/R/script_tissue_specific_efficacy_score.R index 00698b5..e0ea7cb 100644 --- a/R/script_tissue_specific_efficacy_score.R +++ b/R/script_tissue_specific_efficacy_score.R @@ -25,63 +25,109 @@ #'@importFrom Rdpack reprompt #'@importFrom igraph E graph_from_edgelist distances #'@importFrom scales rescale -weighted.shortest.path <- function(disease_genes, ppi_network, directed_network = F, - tissue_expr_data, dis_relevant_tissues, W, cutoff = 1.6, - verbose=F) { - # - if(is.null(rownames(tissue_expr_data))|is.null(colnames(tissue_expr_data))){ - stop('Both colnames and rownames for tissue_expr_data must be provided.') - } - for(i in 1:2) ppi_network[,i] <- as.character(ppi_network[,i]) - ppi_network <- ppi_network[!duplicated(ppi_network[,1:2]),] - ppi_network_size <- nrow(ppi_network) - if(directed_network) idx <- ppi_network[,1]%in%rownames(tissue_expr_data) - else idx <- ppi_network[,1]%in%rownames(tissue_expr_data) & ppi_network[,2]%in%rownames(tissue_expr_data) - ppi_network <- ppi_network[idx,] - if(nrow(ppi_network)==0) stop('No corresponding IDs between ppi_network and tissue_expr_data.') - else if(ppi_network_size!=nrow(ppi_network)){ - if(verbose) print(paste(nrow(ppi_network),'out of',ppi_network_size,'network edges selected.', sep=' ')) - } - if(!directed_network){ - #if(verbose) print('Undirect network. Converting to direct network...') - ppi_network_rev <- ppi_network[,c(2:1)] - colnames(ppi_network_rev) <- colnames(ppi_network)[1:2] - ppi_network <- rbind(ppi_network[,1:2],ppi_network_rev) - } - disease_genes_size <- length(disease_genes) - disease_genes <- intersect(disease_genes,ppi_network[,2]) - if(length(disease_genes)==0) stop('No disease-associated ID match with ppi_network and tissue_expr_data!') - else if( disease_genes_size!=length(disease_genes)){ - if(verbose) print(paste(length(disease_genes),'disease-associated IDs are reachable from the network.',sep=' ')) - } - tissue_expr_data <- scales::rescale(tissue_expr_data,c(1,.Machine$double.eps)) - g <- igraph::graph_from_edgelist(as.matrix(ppi_network[,1:2]), directed=T) - if(is.vector(dis_relevant_tissues) == FALSE) stop('Argument dis_relevant_tissues is not a vector!') - if(is.null(names(dis_relevant_tissues))) stop('Names for dis_relevant_tissues must be provided!') - # selecting relevenat tissues - sign_tiss <- names(which(dis_relevant_tissues >= cutoff)) - if(length(sign_tiss) != 0){ - if(verbose) print(paste(length(sign_tiss),' tissue/s significant for given disease.',sep='')) - target_path_mean <- NULL - sh.path <- lapply(sign_tiss,function(i){ - igraph::E(g)$weight <- tissue_expr_data[ppi_network[,1],i] - igraph::distances(g, to=disease_genes, mode = "out", weights = NULL, algorithm = "dijkstra") - }) - names(sh.path) <- sign_tiss - sh.path <- lapply(sh.path,function(x){ - x[is.infinite(x)]<-(3*max(x[which(is.finite(x))],na.rm = TRUE)) +weighted.shortest.path <- + function(disease_genes, + ppi_network, + directed_network = F, + tissue_expr_data, + dis_relevant_tissues, + W, + cutoff = 1.6, + verbose = F) { + # + if (is.null(rownames(tissue_expr_data)) | + is.null(colnames(tissue_expr_data))) { + stop('Both colnames and rownames for tissue_expr_data must be provided.') + } + for (i in 1:2) + ppi_network[, i] <- as.character(ppi_network[, i]) + ppi_network <- ppi_network[!duplicated(ppi_network[, 1:2]), ] + ppi_network_size <- nrow(ppi_network) + if (directed_network) + idx <- ppi_network[, 1] %in% rownames(tissue_expr_data) + else + idx <- + ppi_network[, 1] %in% rownames(tissue_expr_data) & + ppi_network[, 2] %in% rownames(tissue_expr_data) + ppi_network <- ppi_network[idx, ] + if (nrow(ppi_network) == 0) + stop('No corresponding IDs between ppi_network and tissue_expr_data.') + else if (ppi_network_size != nrow(ppi_network)) { + if (verbose) + print(paste( + nrow(ppi_network), + 'out of', + ppi_network_size, + 'network edges selected.', + sep = ' ' + )) + } + if (!directed_network) { + #if(verbose) print('Undirect network. Converting to direct network...') + ppi_network_rev <- ppi_network[, c(2:1)] + colnames(ppi_network_rev) <- colnames(ppi_network)[1:2] + ppi_network <- rbind(ppi_network[, 1:2], ppi_network_rev) + } + disease_genes_size <- length(disease_genes) + disease_genes <- intersect(disease_genes, ppi_network[, 2]) + if (length(disease_genes) == 0) + stop('No disease-associated ID match with ppi_network and tissue_expr_data!') + else if (disease_genes_size != length(disease_genes)) { + if (verbose) + print(paste( + length(disease_genes), + 'disease-associated IDs are reachable from the network.', + sep = ' ' + )) + } + tissue_expr_data <- + scales::rescale(tissue_expr_data, c(1, .Machine$double.eps)) + g <- + igraph::graph_from_edgelist(as.matrix(ppi_network[, 1:2]), directed = T) + if (is.vector(dis_relevant_tissues) == FALSE) + stop('Argument dis_relevant_tissues is not a vector!') + if (is.null(names(dis_relevant_tissues))) + stop('Names for dis_relevant_tissues must be provided!') + # selecting relevenat tissues + sign_tiss <- names(which(dis_relevant_tissues >= cutoff)) + if (length(sign_tiss) != 0) { + if (verbose) + print(paste( + length(sign_tiss), + ' tissue/s significant for given disease.', + sep = '' + )) + target_path_mean <- NULL + sh.path <- lapply(sign_tiss, function(i) { + igraph::E(g)$weight <- tissue_expr_data[ppi_network[, 1], i] + igraph::distances( + g, + to = disease_genes, + mode = "out", + weights = NULL, + algorithm = "dijkstra" + ) + }) + names(sh.path) <- sign_tiss + sh.path <- lapply(sh.path, function(x) { + x[is.infinite(x)] <- (3 * max(x[which(is.finite(x))], na.rm = TRUE)) x }) - # weighting shortest paths according to disease-associeted genes discrete values - new_sh.path <- mapply(function(x,y) t(x)*y[colnames(x)],sh.path, W[sign_tiss], SIMPLIFY=F) - # averaging shortest paths across the disease-associated genes and the disease-relevant tissues - target_path <- sapply(new_sh.path,colMeans) - target_path_mean <- rowMeans(target_path) - return(list(shortest_paths = target_path, - shortest_paths_avg = target_path_mean)) + # weighting shortest paths according to disease-associeted genes discrete values + new_sh.path <- + mapply(function(x, y) + t(x) * y[colnames(x)], sh.path, W[sign_tiss], SIMPLIFY = F) + # averaging shortest paths across the disease-associated genes and the disease-relevant tissues + target_path <- sapply(new_sh.path, colMeans) + target_path_mean <- rowMeans(target_path) + return(list( + shortest_paths = target_path, + shortest_paths_avg = target_path_mean + )) + } + else + stop('No tissue significant for given disease!') } - else stop('No tissue significant for given disease!') -} #'Compile tissue-specific scores for a given disease. #' @@ -115,71 +161,137 @@ weighted.shortest.path <- function(disease_genes, ppi_network, directed_network #'@export #'@importFrom Rdpack reprompt #'@importFrom scales rescale -tissue.specific.scores <- function(disease_genes, ppi_network, directed_network = F, - tissue_expr_data, dis_relevant_tissues, W, cutoff = 1.6, - verbose = FALSE) { - if(is.null(disease_genes) | is.null(tissue_expr_data) | is.null(ppi_network) | is.null(dis_relevant_tissues)) - stop('Incomplete data input.') - # compile the wsp - wsp_list = weighted.shortest.path(disease_genes,ppi_network, - directed_network,tissue_expr_data, - dis_relevant_tissues,W, cutoff, verbose) - tissue_scores <- apply(wsp_list$shortest_paths, 2, function(x) { - q <- quantile(x) - outliers <- q[4]+(1.5*(q[4]-q[2])) - ts <- scales::rescale(x, from=c(min(x), outliers), to=c(1,0)) - ts[ts< 0] <- 0 - ts - }) - dt <- data.frame(tissue_scores, avg_tissue_score = rowMeans(tissue_scores[,,drop=F])) - return(dt) -} +tissue.specific.scores <- + function(disease_genes, + ppi_network, + directed_network = F, + tissue_expr_data, + dis_relevant_tissues, + W, + cutoff = 1.6, + verbose = FALSE) { + if (is.null(disease_genes) | + is.null(tissue_expr_data) | + is.null(ppi_network) | is.null(dis_relevant_tissues)) + stop('Incomplete data input.') + # compile the wsp + wsp_list = weighted.shortest.path( + disease_genes, + ppi_network, + directed_network, + tissue_expr_data, + dis_relevant_tissues, + W, + cutoff, + verbose + ) + tissue_scores <- apply(wsp_list$shortest_paths, 2, function(x) { + q <- quantile(x) + outliers <- q[4] + (1.5 * (q[4] - q[2])) + ts <- scales::rescale(x, from = c(min(x), outliers), to = c(1, 0)) + ts[ts < 0] <- 0 + ts + }) + dt <- + data.frame(tissue_scores, avg_tissue_score = rowMeans(tissue_scores[, , drop = + F])) + return(dt) + } #'@rdname tissue.specific.scores #'@export #'@importFrom scales rescale #'@importFrom snow makeCluster stopCluster #'@importFrom doParallel registerDoParallel #'@importFrom foreach foreach %dopar% -list.tissue.specific.scores <- function(disease_gene_list, ppi_network, directed_network = F, - tissue_expr_data, dis_relevant_tissues, W, cutoff = 1.6, - verbose = FALSE, parallel = NULL) { - require(foreach) - require(doParallel) - if(is.list(disease_gene_list) == FALSE) stop('Argument disease_gene_list is not a list!') - if(is.null(names(disease_gene_list))) stop('Names for disease_gene_list must be provided!') - if(is.matrix(dis_relevant_tissues) == FALSE) stop('Argument dis_relevant_tissues is not a matrix!') - if(is.null(rownames(dis_relevant_tissues))|is.null(colnames(dis_relevant_tissues))){ - stop('Both colnames and rownames for dis_relevant_tissues must be provided!') - } - common_diseases <- intersect(names(disease_gene_list), rownames(dis_relevant_tissues)) - if(length(common_diseases)==0) stop('No disease in common between disease_gene_list and dis_relevant_tissues!') - else if(length(common_diseases) != length(names(disease_gene_list))){ - if(verbose) print(paste(length(common_diseases),'diseases in common between disease_gene_list and dis_relevant_tissues.',sep=' ')) - } - tissue_scores <- NULL - if(!is.null(parallel)) { - cl <- snow::makeCluster(parallel) - doParallel::registerDoParallel(cl) - `%dopar%` <- foreach::`%dopar%` - wsp <- foreach::foreach(i=common_diseases,.export = 'weighted.shortest.path',.final = function(i) setNames(i, common_diseases)) %dopar% - weighted.shortest.path(disease_gene_list[[i]], ppi_network, directed_network, tissue_expr_data, dis_relevant_tissues[i,], W, cutoff = 1.6) - snow::stopCluster(cl) - } - else { - warning("A parallel computation is highly recommended",immediate. = T) - wsp <- sapply(common_diseases,function(i) weighted.shortest.path(disease_gene_list[[i]], ppi_network, directed_network, tissue_expr_data, dis_relevant_tissues[i,], W, cutoff = 1.6),simplify=F) +list.tissue.specific.scores <- + function(disease_gene_list, + ppi_network, + directed_network = F, + tissue_expr_data, + dis_relevant_tissues, + W, + cutoff = 1.6, + verbose = FALSE, + parallel = NULL) { + require(foreach) + require(doParallel) + if (is.list(disease_gene_list) == FALSE) + stop('Argument disease_gene_list is not a list!') + if (is.null(names(disease_gene_list))) + stop('Names for disease_gene_list must be provided!') + if (is.matrix(dis_relevant_tissues) == FALSE) + stop('Argument dis_relevant_tissues is not a matrix!') + if (is.null(rownames(dis_relevant_tissues)) | + is.null(colnames(dis_relevant_tissues))) { + stop('Both colnames and rownames for dis_relevant_tissues must be provided!') + } + common_diseases <- + intersect(names(disease_gene_list), rownames(dis_relevant_tissues)) + if (length(common_diseases) == 0) + stop('No disease in common between disease_gene_list and dis_relevant_tissues!') + else if (length(common_diseases) != length(names(disease_gene_list))) { + if (verbose) + print( + paste( + length(common_diseases), + 'diseases in common between disease_gene_list and dis_relevant_tissues.', + sep = ' ' + ) + ) + } + tissue_scores <- NULL + if (!is.null(parallel)) { + cl <- snow::makeCluster(parallel) + doParallel::registerDoParallel(cl) + `%dopar%` <- foreach::`%dopar%` + wsp <- + foreach::foreach( + i = common_diseases, + .export = 'weighted.shortest.path', + .final = function(i) + setNames(i, common_diseases) + ) %dopar% + weighted.shortest.path( + disease_gene_list[[i]], + ppi_network, + directed_network, + tissue_expr_data, + dis_relevant_tissues[i, ], + W, + cutoff = 1.6 + ) + snow::stopCluster(cl) + } + else { + warning("A parallel computation is highly recommended", immediate. = T) + wsp <- + sapply(common_diseases, function(i) + weighted.shortest.path( + disease_gene_list[[i]], + ppi_network, + directed_network, + tissue_expr_data, + dis_relevant_tissues[i, ], + W, + cutoff = 1.6 + ), simplify = F) + } + wsp <- do.call(function(...) + mapply(list, ..., SIMPLIFY = F), wsp) + tissue_scores <- + lapply(wsp$shortest_paths, function(x) + apply(x, 2, function(y) { + q <- quantile(y) + outliers <- q[4] + (1.5 * (q[4] - q[2])) + ts <- scales::rescale(y, from = c(min(x), outliers), to = c(1, 0)) + ts[ts < 0] <- 0 + ts + })) + tissue_scores <- + lapply(tissue_scores, function(x) + data.frame(x, avg_tissue_score = rowMeans(x[, , drop = F]))) + return(tissue_scores) } - wsp <- do.call(function(...)mapply(list,...,SIMPLIFY = F),wsp) - tissue_scores <- lapply(wsp$shortest_paths, function(x) apply(x, 2, function(y) { - q <- quantile(y) - outliers <- q[4]+(1.5*(q[4]-q[2])) - ts <- scales::rescale(y, from=c(min(x), outliers), to=c(1,0)) - ts[ts< 0] <- 0 - ts - })) - tissue_scores <- lapply(tissue_scores,function(x)data.frame(x, avg_tissue_score = rowMeans(x[,,drop=F]))) - return(tissue_scores) -} #'Compile tissue-specific networks and, within these, set of genes that are relevant for the selected targets. @@ -215,60 +327,125 @@ list.tissue.specific.scores <- function(disease_gene_list, ppi_network, directed #'@importFrom igraph V #'@importFrom igraph E #'@importFrom igraph shortest_paths -build.tissue.specific.networks <- function(tissue_scores, disease_genes, ppi_network, - directed_network = F, tissue_expr_data, - top_targets = NULL, rwr_restart = 0.75, - rwr_norm = "quantile", rwr_cutoff = 0.001, - verbose = FALSE){ - # check - if(is.null(top_targets)) stop('Please specifiy a set of targets (ENTREZ ids)!') - if(!is.character(ppi_network[,1])) ppi_network[,1] <- as.character(ppi_network[,1]) - if(!is.character(ppi_network[,2])) ppi_network[,2] <- as.character(ppi_network[,2]) - ppi_network <- ppi_network[!duplicated(ppi_network[,1:2]),] - ppi_network_node<-unique(unlist(ppi_network[,1:2])) - universe <- intersect(ppi_network_node,rownames(tissue_expr_data)) - if(length(universe)==0) stop('No corresponding IDs between ppi_network and tissue_expr_data!') - else if(length(universe)!=length(ppi_network_node)| - length(universe)!=nrow(tissue_expr_data)){ - if (verbose) print(paste(length(universe),'IDs in common between ppi_network and tissue_expr_data will be considered.', sep=' ')) - tissue_expr_data<-tissue_expr_data[universe,] - } - ## scale tissue expression scores - tissue_expr_data <- scales::rescale(tissue_expr_data,c(1,.Machine$double.eps)) - ## select the sign tissues - sign_tiss <- colnames(tissue_scores)[-ncol(tissue_scores)] - tissue_spec_network <- list() - shp_top_genes <- list() - rwr_top_genes <- list() - for(i in 1: length(sign_tiss)) { - g <- igraph::graph_from_edgelist(as.matrix(ppi_network[,1:2]), directed=T) - igraph::E(g)$weight <- tissue_expr_data[ppi_network[,1], sign_tiss[i]] - # select disease genes active in this tissue and search for shortest paths - tissue_disease_genes <- intersect(disease_genes, igraph::V(g)$name) - shp_top_genes[[length(shp_top_genes) + 1]] <- lapply(top_targets, function(x) { - ss = igraph::distances(g, v=x, to=tissue_disease_genes, mode = "out", weights = NULL, algorithm = "dijkstra") - if(verbose) print(paste("Targets ", x, " cannot reach ", length(which(is.infinite(ss)))," nodes", sep="")) - suppressWarnings(selected_nodes <- unique(unlist(igraph::shortest_paths(g, from=x, to=tissue_disease_genes, mode = "out", weights = NULL, output='vpath')$vpath))) - igraph::V(g)$name[selected_nodes] - }) - names(shp_top_genes[[length(shp_top_genes)]]) <- top_targets - # apply random walk on the top targets - aSeeds <- rep(0, length(igraph::V(g)$name)) - aSeeds[which((igraph::V(g)$name%in%top_targets) == TRUE)] <- 1 - setSeeds <- data.frame(aSeeds) - rownames(setSeeds) <- igraph::V(g)$name - rwr_top_genes[[length(rwr_top_genes)+1]] <- as.numeric(dnet::dRWR(g, normalise = "none", setSeeds = setSeeds, restart = rwr_restart, - normalise.affinity.matrix = rwr_norm, parallel = FALSE, multicores = NULL, verbose = F)) - names(rwr_top_genes[[length(rwr_top_genes)]]) = igraph::V(g)$name - rwr_top_genes[[length(rwr_top_genes)]] = names(rwr_top_genes[[length(rwr_top_genes)]])[which(rwr_top_genes[[length(rwr_top_genes)]] >= rwr_cutoff)] - # - tissue_spec_network[[length(tissue_spec_network) + 1]] = g +build.tissue.specific.networks <- + function(tissue_scores, + disease_genes, + ppi_network, + directed_network = F, + tissue_expr_data, + top_targets = NULL, + rwr_restart = 0.75, + rwr_norm = "quantile", + rwr_cutoff = 0.001, + verbose = FALSE) { + # check + if (is.null(top_targets)) + stop('Please specifiy a set of targets (ENTREZ ids)!') + if (!is.character(ppi_network[, 1])) + ppi_network[, 1] <- as.character(ppi_network[, 1]) + if (!is.character(ppi_network[, 2])) + ppi_network[, 2] <- as.character(ppi_network[, 2]) + ppi_network <- ppi_network[!duplicated(ppi_network[, 1:2]), ] + ppi_network_node <- unique(unlist(ppi_network[, 1:2])) + universe <- intersect(ppi_network_node, rownames(tissue_expr_data)) + if (length(universe) == 0) + stop('No corresponding IDs between ppi_network and tissue_expr_data!') + else if (length(universe) != length(ppi_network_node) | + length(universe) != nrow(tissue_expr_data)) { + if (verbose) + print( + paste( + length(universe), + 'IDs in common between ppi_network and tissue_expr_data will be considered.', + sep = ' ' + ) + ) + tissue_expr_data <- tissue_expr_data[universe, ] + } + ## scale tissue expression scores + tissue_expr_data <- + scales::rescale(tissue_expr_data, c(1, .Machine$double.eps)) + ## select the sign tissues + sign_tiss <- colnames(tissue_scores)[-ncol(tissue_scores)] + tissue_spec_network <- list() + shp_top_genes <- list() + rwr_top_genes <- list() + for (i in 1:length(sign_tiss)) { + g <- + igraph::graph_from_edgelist(as.matrix(ppi_network[, 1:2]), directed = T) + igraph::E(g)$weight <- + tissue_expr_data[ppi_network[, 1], sign_tiss[i]] + # select disease genes active in this tissue and search for shortest paths + tissue_disease_genes <- + intersect(disease_genes, igraph::V(g)$name) + shp_top_genes[[length(shp_top_genes) + 1]] <- + lapply(top_targets, function(x) { + ss = igraph::distances( + g, + v = x, + to = tissue_disease_genes, + mode = "out", + weights = NULL, + algorithm = "dijkstra" + ) + if (verbose) + print(paste( + "Targets ", + x, + " cannot reach ", + length(which(is.infinite(ss))), + " nodes", + sep = "" + )) + suppressWarnings(selected_nodes <- + unique(unlist( + igraph::shortest_paths( + g, + from = x, + to = tissue_disease_genes, + mode = "out", + weights = NULL, + output = 'vpath' + )$vpath + ))) + igraph::V(g)$name[selected_nodes] + }) + names(shp_top_genes[[length(shp_top_genes)]]) <- top_targets + # apply random walk on the top targets + aSeeds <- rep(0, length(igraph::V(g)$name)) + aSeeds[which((igraph::V(g)$name %in% top_targets) == TRUE)] <- 1 + setSeeds <- data.frame(aSeeds) + rownames(setSeeds) <- igraph::V(g)$name + rwr_top_genes[[length(rwr_top_genes) + 1]] <- + as.numeric( + dnet::dRWR( + g, + normalise = "none", + setSeeds = setSeeds, + restart = rwr_restart, + normalise.affinity.matrix = rwr_norm, + parallel = FALSE, + multicores = NULL, + verbose = F + ) + ) + names(rwr_top_genes[[length(rwr_top_genes)]]) = igraph::V(g)$name + rwr_top_genes[[length(rwr_top_genes)]] = names(rwr_top_genes[[length(rwr_top_genes)]])[which(rwr_top_genes[[length(rwr_top_genes)]] >= rwr_cutoff)] + # + tissue_spec_network[[length(tissue_spec_network) + 1]] = g + } + names(rwr_top_genes) <- sign_tiss + names(shp_top_genes) <- sign_tiss + names(tissue_spec_network) <- sign_tiss + return( + list( + "tsn" = tissue_spec_network, + "shp" = shp_top_genes, + "rwr" = rwr_top_genes, + "universe" = universe + ) + ) } - names(rwr_top_genes) <- sign_tiss - names(shp_top_genes) <- sign_tiss - names(tissue_spec_network) <- sign_tiss - return(list("tsn"=tissue_spec_network, "shp" = shp_top_genes, "rwr" = rwr_top_genes, "universe" = universe)) -} #'Over-Representation Analysis. #' @@ -287,38 +464,53 @@ build.tissue.specific.networks <- function(tissue_scores, disease_genes, ppi_net #'@importFrom clusterProfiler enrichKEGG #'@importFrom clusterProfiler enrichGO #'@importFrom ReactomePA enrichPathway -generate.ora.data <- function(input, databases = c("GO", "KEGG"), - orgdb_go = 'org.Hs.eg.db', orgdb_kegg = 'hsa', - apval = 0.01, verbose = TRUE) { +generate.ora.data <- function(input, + databases = c("GO", "KEGG"), + orgdb_go = 'org.Hs.eg.db', + orgdb_kegg = 'hsa', + apval = 0.01, + verbose = TRUE) { list_annots <- list() annots <- c() - for(i in 1:length(input)) { - if("GO" %in% databases) { - list_annots[[length(list_annots) + 1]] <- clusterProfiler::enrichGO(unique(unlist(input[[i]])), - orgdb_go, ont = "BP", - pvalueCutoff = apval, - universe = input$universe) - annots <- c(annots,"GO") - if(verbose) print(paste("ORA completed for GO-", names(input)[i], ".", sep="")) + for (i in 1:length(input)) { + if ("GO" %in% databases) { + list_annots[[length(list_annots) + 1]] <- + clusterProfiler::enrichGO( + unique(unlist(input[[i]])), + orgdb_go, + ont = "BP", + pvalueCutoff = apval, + universe = input$universe + ) + annots <- c(annots, "GO") + if (verbose) + print(paste("ORA completed for GO-", names(input)[i], ".", sep = "")) } - if("KEGG" %in% databases) { - list_annots[[length(list_annots) + 1]] <- clusterProfiler::enrichKEGG(unique(unlist(input[[i]])), - organism = orgdb_kegg, - keyType = "kegg", - pvalueCutoff = apval, - universe = input$universe) - annots <- c(annots,"KEGG") - if(verbose) print(paste("ORA completed for KEGG-", names(input)[i], ".", sep="")) + if ("KEGG" %in% databases) { + list_annots[[length(list_annots) + 1]] <- + clusterProfiler::enrichKEGG( + unique(unlist(input[[i]])), + organism = orgdb_kegg, + keyType = "kegg", + pvalueCutoff = apval, + universe = input$universe + ) + annots <- c(annots, "KEGG") + if (verbose) + print(paste("ORA completed for KEGG-", names(input)[i], ".", sep = "")) } - if("REACTOME" %in% databases) { - list_annots[[length(list_annots) + 1]] <- ReactomePA::enrichPathway(unique(unlist(input[[i]])), - pvalueCutoff = apval, - readable = T) - annots <- c(annots,"REACTOME") - if(verbose) print(paste("ORA completed for REACTOME-", names(input)[i], ".", sep="")) + if ("REACTOME" %in% databases) { + list_annots[[length(list_annots) + 1]] <- + ReactomePA::enrichPathway(unique(unlist(input[[i]])), + pvalueCutoff = apval, + readable = T) + annots <- c(annots, "REACTOME") + if (verbose) + print(paste("ORA completed for REACTOME-", names(input)[i], ".", sep = + "")) } } - names(list_annots) <- paste(annots, names(input), sep="-") + names(list_annots) <- paste(annots, names(input), sep = "-") return(list_annots) } @@ -340,40 +532,60 @@ generate.ora.data <- function(input, databases = c("GO", "KEGG"), #'@importFrom enrichplot upsetplot #'@importFrom ggplot2 theme #'@importFrom ggplot2 element_text -generate.ora.plots <- function(ora_data, set_plots = c("dotplot", "emapplot", "cnetplot", "upsetplot"), - showCategory = 20, font_size = 8){ - list_plots <- list() - plots <- c() - for(i in 1:length(ora_data)) { - if("dotplot" %in% set_plots) { - list_plots[[length(list_plots)+1]] <- enrichplot::dotplot(ora_data[[i]], showCategory = showCategory, font.size = font_size) + - ggplot2::theme(legend.title=ggplot2::element_text(size=font_size), - legend.text=ggplot2::element_text(size=font_size)) - plots <- c(plots, paste("dotplot",names(ora_data)[i],sep="-")) - } - if("emapplot" %in% set_plots) { - list_plots[[length(list_plots)+1]] <- enrichplot::emapplot(ora_data[[i]], showCategory = showCategory) + - ggplot2::theme(legend.title=ggplot2::element_text(size=font_size), - legend.text=ggplot2::element_text(size=font_size)) - plots <- c(plots, paste("emapplot",names(ora_data)[i],sep="-")) - } - if("cnetplot" %in% set_plots) { - list_plots[[length(list_plots)+1]] <- enrichplot::cnetplot(ora_data[[i]], showCategory = showCategory) + - ggplot2::theme(legend.title=ggplot2::element_text(size=font_size), - legend.text=ggplot2::element_text(size=font_size)) - plots <- c(plots, paste("cnetplot",names(ora_data)[i],sep="-")) - } - if("upsetplot" %in% set_plots) { - list_plots[[length(list_plots)+1]] <- enrichplot::upsetplot(ora_data[[i]]) + - ggplot2::theme(legend.title=ggplot2::element_text(size=font_size), - legend.text=ggplot2::element_text(size=font_size)) - plots <- c(plots, paste("upsetplot",names(ora_data)[i],sep="-")) +generate.ora.plots <- + function(ora_data, + set_plots = c("dotplot", "emapplot", "cnetplot", "upsetplot"), + showCategory = 20, + font_size = 8) { + list_plots <- list() + plots <- c() + for (i in 1:length(ora_data)) { + if ("dotplot" %in% set_plots) { + list_plots[[length(list_plots) + 1]] <- + enrichplot::dotplot(ora_data[[i]], + showCategory = showCategory, + font.size = font_size) + + ggplot2::theme( + legend.title = ggplot2::element_text(size = font_size), + legend.text = ggplot2::element_text(size = font_size) + ) + plots <- c(plots, paste("dotplot", names(ora_data)[i], sep = "-")) + } + if ("emapplot" %in% set_plots) { + list_plots[[length(list_plots) + 1]] <- + enrichplot::emapplot(ora_data[[i]], showCategory = showCategory) + + ggplot2::theme( + legend.title = ggplot2::element_text(size = font_size), + legend.text = ggplot2::element_text(size = font_size) + ) + plots <- + c(plots, paste("emapplot", names(ora_data)[i], sep = "-")) + } + if ("cnetplot" %in% set_plots) { + list_plots[[length(list_plots) + 1]] <- + enrichplot::cnetplot(ora_data[[i]], showCategory = showCategory) + + ggplot2::theme( + legend.title = ggplot2::element_text(size = font_size), + legend.text = ggplot2::element_text(size = font_size) + ) + plots <- + c(plots, paste("cnetplot", names(ora_data)[i], sep = "-")) + } + if ("upsetplot" %in% set_plots) { + list_plots[[length(list_plots) + 1]] <- + enrichplot::upsetplot(ora_data[[i]]) + + ggplot2::theme( + legend.title = ggplot2::element_text(size = font_size), + legend.text = ggplot2::element_text(size = font_size) + ) + plots <- + c(plots, paste("upsetplot", names(ora_data)[i], sep = "-")) + } } + print(plots) + names(list_plots) <- plots + return(list_plots) } - print(plots) - names(list_plots) <- plots - return(list_plots) -} #'Generate PubMed Central Trend plots. @@ -392,16 +604,27 @@ generate.ora.plots <- function(ora_data, set_plots = c("dotplot", "emapplot", "c #'@importFrom ggplot2 theme #'@importFrom ggplot2 theme_minimal #'@importFrom ggplot2 element_text -novelty.plots <- function(list_genes, gene_id = 'ENTREZID', orgdb = NULL, pubmed = c(2010,2019), font_size = 8){ - pmc.genes = as.character(AnnotationDbi::mapIds(orgdb, list_genes, 'SYMBOL', gene_id)) - print(pmc.genes) - plot <- enrichplot::pmcplot(pmc.genes, pubmed[1]:pubmed[2]) + ggplot2::theme_minimal() + - ggplot2::theme(legend.title=ggplot2::element_text(size=font_size), - legend.text=element_text(size=font_size), - axis.text.x = element_text(size=font_size), - axis.text.y = element_text(size=font_size)) - return(plot) -} +novelty.plots <- + function(list_genes, + gene_id = 'ENTREZID', + orgdb = NULL, + pubmed = c(2010, 2019), + font_size = 8) { + pmc.genes = as.character(AnnotationDbi::mapIds(orgdb, list_genes, 'SYMBOL', gene_id)) + print(pmc.genes) + plot <- + enrichplot::pmcplot(pmc.genes, pubmed[1]:pubmed[2]) + ggplot2::theme_minimal() + + ggplot2::theme( + legend.title = ggplot2::element_text(size = font_size), + legend.text = element_text(size = + font_size), + axis.text.x = element_text(size = + font_size), + axis.text.y = element_text(size = + font_size) + ) + return(plot) + } #'Interactive visualization of tissue-specific networks. #' @@ -422,132 +645,213 @@ novelty.plots <- function(list_genes, gene_id = 'ENTREZID', orgdb = NULL, pubmed #'Possible values are: \code{kegg}, \code{BP}, \code{MF} and \code{CC}. Defaults to \code{kegg}. #'@param verbose logical indicating whether the messages will be displayed or not in the screen. #'@importFrom AnnotationDbi mapIds -visualize.graph <- function(tissue_scores, disease_genes, ppi_network, directed_network = FALSE, - tissue_expr_data, top_target = NULL, db='kegg', verbose = FALSE){ - if(is.null(top_target)) stop('Please specifiy a set of targets (ENTREZ ids)!') - if(is.null(rownames(tissue_expr_data))|is.null(colnames(tissue_expr_data))){ - stop('Both colnames and rownames for tissue_expr_data must be provided.') - } - for(i in 1:2) ppi_network[,i] <- as.character(ppi_network[,i]) - ppi_network <- ppi_network[!duplicated(ppi_network[,1:2]),] - ppi_network_size <- nrow(ppi_network) - if(directed_network) idx <- ppi_network[,1]%in%rownames(tissue_expr_data) - else idx <- ppi_network[,1]%in%rownames(tissue_expr_data) & ppi_network[,2]%in%rownames(tissue_expr_data) - ppi_network <- ppi_network[idx,] - if(nrow(ppi_network)==0) stop('No corresponding IDs between ppi_network and tissue_expr_data.') - else if(ppi_network_size!=nrow(ppi_network)){ - if(verbose) print(paste(nrow(ppi_network),'out of',ppi_network_size,'network edges selected.', sep=' ')) - } - if(!directed_network){ - if(verbose) print('Undirect network. Converting to direct network...') - ppi_network_rev <- ppi_network[,c(2:1)] - colnames(ppi_network_rev) <- colnames(ppi_network)[1:2] - ppi_network <- rbind(ppi_network[,1:2],ppi_network_rev) - } - disease_genes_size <- length(disease_genes) - disease_genes <- intersect(disease_genes,ppi_network[,2]) - if(length(disease_genes)==0) stop('No disease-associated ID match with ppi_network and tissue_expr_data!') - else if( disease_genes_size!=length(disease_genes)){ - if(verbose) print(paste(length(disease_genes),'disease-associated IDs are reachable from the network.',sep=' ')) - } - if(!(db %in% c('kegg','BP','MF','CC'))) stop('Unused value for argument db') - idx <- order(tissue_scores$avg_tissue_score,decreasing=T) - tissue_scores <- tissue_scores[idx,] - g <- igraph::graph_from_edgelist(as.matrix(ppi_network[,1:2]), directed=T) - tissue_expr_data <- scales::rescale(tissue_expr_data,c(1,.Machine$double.eps)) - sign_tiss <- colnames(tissue_scores)[-ncol(tissue_scores)] - sh.path <- lapply(sign_tiss,function(i){ - igraph::E(g)$weight <- tissue_expr_data[ppi_network[,1],i] - lapply(top_target,function(x) - igraph::shortest_paths(g, from=x, to=disease_genes, mode = "out", weights = NULL,output='epath')$epath) - }) - names(sh.path) <- sign_tiss - for(i in names(sh.path)) names(sh.path[[i]]) <- top_target - all_target_path <-lapply(sh.path,function(x)lapply(x,function(y)do.call(c,y))) - all_tissue_path <- lapply(all_target_path,function(x)do.call(c,x)) - ids <- lapply(all_tissue_path,igraph::as_ids) - shinyApp( - ui = fluidPage( - fluidRow( - column( - width = 2, - checkboxGroupInput("tissue","Select tissue/s",sign_tiss,selected = sign_tiss[1]) +visualize.graph <- + function(tissue_scores, + disease_genes, + ppi_network, + directed_network = FALSE, + tissue_expr_data, + top_target = NULL, + db = 'kegg', + verbose = FALSE) { + if (is.null(top_target)) + stop('Please specifiy a set of targets (ENTREZ ids)!') + if (is.null(rownames(tissue_expr_data)) | + is.null(colnames(tissue_expr_data))) { + stop('Both colnames and rownames for tissue_expr_data must be provided.') + } + for (i in 1:2) + ppi_network[, i] <- as.character(ppi_network[, i]) + ppi_network <- ppi_network[!duplicated(ppi_network[, 1:2]), ] + ppi_network_size <- nrow(ppi_network) + if (directed_network) + idx <- ppi_network[, 1] %in% rownames(tissue_expr_data) + else + idx <- + ppi_network[, 1] %in% rownames(tissue_expr_data) & + ppi_network[, 2] %in% rownames(tissue_expr_data) + ppi_network <- ppi_network[idx, ] + if (nrow(ppi_network) == 0) + stop('No corresponding IDs between ppi_network and tissue_expr_data.') + else if (ppi_network_size != nrow(ppi_network)) { + if (verbose) + print(paste( + nrow(ppi_network), + 'out of', + ppi_network_size, + 'network edges selected.', + sep = ' ' + )) + } + if (!directed_network) { + if (verbose) + print('Undirect network. Converting to direct network...') + ppi_network_rev <- ppi_network[, c(2:1)] + colnames(ppi_network_rev) <- colnames(ppi_network)[1:2] + ppi_network <- rbind(ppi_network[, 1:2], ppi_network_rev) + } + disease_genes_size <- length(disease_genes) + disease_genes <- intersect(disease_genes, ppi_network[, 2]) + if (length(disease_genes) == 0) + stop('No disease-associated ID match with ppi_network and tissue_expr_data!') + else if (disease_genes_size != length(disease_genes)) { + if (verbose) + print(paste( + length(disease_genes), + 'disease-associated IDs are reachable from the network.', + sep = ' ' + )) + } + if (!(db %in% c('kegg', 'BP', 'MF', 'CC'))) + stop('Unused value for argument db') + idx <- order(tissue_scores$avg_tissue_score, decreasing = T) + tissue_scores <- tissue_scores[idx, ] + g <- + igraph::graph_from_edgelist(as.matrix(ppi_network[, 1:2]), directed = T) + tissue_expr_data <- + scales::rescale(tissue_expr_data, c(1, .Machine$double.eps)) + sign_tiss <- colnames(tissue_scores)[-ncol(tissue_scores)] + sh.path <- lapply(sign_tiss, function(i) { + igraph::E(g)$weight <- tissue_expr_data[ppi_network[, 1], i] + lapply(top_target, function(x) + igraph::shortest_paths( + g, + from = x, + to = disease_genes, + mode = "out", + weights = NULL, + output = 'epath' + )$epath) + }) + names(sh.path) <- sign_tiss + for (i in names(sh.path)) + names(sh.path[[i]]) <- top_target + all_target_path <- + lapply(sh.path, function(x) + lapply(x, function(y) + do.call(c, y))) + all_tissue_path <- lapply(all_target_path, function(x) + do.call(c, x)) + ids <- lapply(all_tissue_path, igraph::as_ids) + + require(shiny) + require(visNetwork) + shinyApp( + ui = fluidPage( + fluidRow( + column( + width = 2, + checkboxGroupInput("tissue", "Select tissue/s", sign_tiss, selected = sign_tiss[1]) + ), + column(width = 10, + visNetworkOutput("network")) ), - column( - width = 10, - visNetworkOutput("network") - ) - ), - fluidRow( - column( - width = 12, - tableOutput("table") - ) + fluidRow(column(width = 12, + tableOutput("table"))), + textOutput("text") ), - textOutput("text") - ), - server = function(input, output,session) { - last_id <- NULL - output$network <- renderVisNetwork({ - validate(need(input$tissue!='','Please select one or more tissues')) - if(length(input$tissue)>1){ - edge <- strsplit(Reduce(intersect,ids[input$tissue]),split='|',fixed=T) - } - else if (length(input$tissue) ==1){ - edge <- strsplit(ids[[input$tissue]],split='|',fixed=T) - } - edge <- as.data.frame(do.call(rbind,edge)) - edge <- edge[!duplicated(edge),] - colnames(edge) <- c('from','to') - edge$width <- scales::rescale(rowMeans(1-tissue_expr_data[edge$from,input$tissue,drop=F]),c(1,5)) - edge$arrows <- 'to' - nb <- unique(unlist(edge[,1:2])) - node <- data.frame(id=nb,label=nb,stringsAsFactors = F) - node$title <- AnnotationDbi::mapIds(org.Hs.eg.db,as.character(nb),'SYMBOL','ENTREZID', verbose = FALSE) - node$size <- rowMeans(tissue_scores[nb,input$tissue,drop=F])*20 - gp <- rep('bridge gene',nrow(node)) - gp[node$label%in%top_target] <- 'target gene' - gp[node$label%in%disease_genes] <- 'disease gene' - node$group <- gp - node$shape <- rep('dot',length(gp)) - node$shape[gp=='target gene'] <- 'triangle' - node$shape[gp=='disease gene'] <- 'star' - visNetwork::visNetwork(node, edge, height = "1500px", width = "500%") %>% - visGroups(groupname = "target gene", color = list(background = "skyblue", border = "deepskyblue")) %>% - visGroups(groupname = "disease gene", color = list(background = "lightcoral", border = "red")) %>% - visGroups(groupname = "bridge gene", shape="dot") %>% - visLegend(addNodes=list( - list(label = "disease gene", shape = "star", - color = list(background = "lightcoral", border = "red")), - list(label = "target gene", shape = "triangle", - color = list(background = "skyblue", border = "deepskyblue")), - list(label = "bridge gene", shape = "dot")), - useGroups = FALSE,width=0.1) %>% - visIgraphLayout(layout = "layout_with_sugiyama") %>% - visEvents(select = "function(nodes) { + server = function(input, output, session) { + last_id <- NULL + output$network <- renderVisNetwork({ + validate(need(input$tissue != '', 'Please select one or more tissues')) + if (length(input$tissue) > 1) { + edge <- + strsplit(Reduce(intersect, ids[input$tissue]), + split = '|', + fixed = T) + } + else if (length(input$tissue) == 1) { + edge <- strsplit(ids[[input$tissue]], split = '|', fixed = T) + } + edge <- as.data.frame(do.call(rbind, edge)) + edge <- edge[!duplicated(edge), ] + colnames(edge) <- c('from', 'to') + edge$width <- + scales::rescale(rowMeans(1 - tissue_expr_data[edge$from, input$tissue, drop = + F]), c(1, 5)) + edge$arrows <- 'to' + nb <- unique(unlist(edge[, 1:2])) + node <- data.frame(id = nb, + label = nb, + stringsAsFactors = F) + node$title <- + AnnotationDbi::mapIds(org.Hs.eg.db, + as.character(nb), + 'SYMBOL', + 'ENTREZID', + verbose = FALSE) + node$size <- + rowMeans(tissue_scores[nb, input$tissue, drop = F]) * 20 + gp <- rep('bridge gene', nrow(node)) + gp[node$label %in% top_target] <- 'target gene' + gp[node$label %in% disease_genes] <- 'disease gene' + node$group <- gp + node$shape <- rep('dot', length(gp)) + node$shape[gp == 'target gene'] <- 'triangle' + node$shape[gp == 'disease gene'] <- 'star' + visNetwork::visNetwork(node, edge, height = "1500px", width = "500%") %>% + visGroups( + groupname = "target gene", + color = list(background = "skyblue", border = "deepskyblue") + ) %>% + visGroups( + groupname = "disease gene", + color = list(background = "lightcoral", border = "red") + ) %>% + visGroups(groupname = "bridge gene", shape = "dot") %>% + visLegend( + addNodes = list( + list( + label = "disease gene", + shape = "star", + color = list(background = "lightcoral", border = "red") + ), + list( + label = "target gene", + shape = "triangle", + color = list(background = "skyblue", border = "deepskyblue") + ), + list(label = "bridge gene", shape = "dot") + ), + useGroups = FALSE, + width = 0.1 + ) %>% + visIgraphLayout(layout = "layout_with_sugiyama") %>% + visEvents(select = "function(nodes) { Shiny.onInputChange('current_node_id', nodes.nodes); ;}") - }) - output$table <- renderTable({ - res = NULL - if (is.null(last_id) | - !identical(last_id,input$current_node_id) ) { - validate(need(input$current_node_id%in%top_target,'Please select a target gene')) - last_id <<- input$current_node_id - current_target_path <- sapply(all_target_path[input$tissue],function(x) x[input$current_node_id]) - current_target_path <- Reduce(igraph::intersection,current_target_path) - target_interactors <- unique(c(igraph::ends(g,current_target_path))) - if(db=='kegg') res <- clusterProfiler::enrichKEGG(target_interactors,organism = 'hsa')@result - else res <- clusterProfiler::enrichGO(target_interactors, 'org.Hs.eg.db', ont = db)@result - res <- res[res$p.adjust<0.05,] - if (nrow(res) > 25) res <- res[1:25,] - res - } - else{ - validate(need(!is.null(res),'Please select a target gene')) - } - }) - } - ) + }) + output$table <- renderTable({ + res = NULL + if (is.null(last_id) | + !identical(last_id, input$current_node_id)) { + validate(need( + input$current_node_id %in% top_target, + 'Please select a target gene' + )) + last_id <<- input$current_node_id + current_target_path <- + sapply(all_target_path[input$tissue], function(x) + x[input$current_node_id]) + current_target_path <- + Reduce(igraph::intersection, current_target_path) + target_interactors <- + unique(c(igraph::ends(g, current_target_path))) + if (db == 'kegg') + res <- + clusterProfiler::enrichKEGG(target_interactors, organism = 'hsa')@result + else + res <- + clusterProfiler::enrichGO(target_interactors, 'org.Hs.eg.db', ont = db)@result + res <- res[res$p.adjust < 0.05, ] + if (nrow(res) > 25) + res <- res[1:25, ] + res + } + else{ + validate(need(!is.null(res), 'Please select a target gene')) + } + }) + } + ) } -