feat(indra): fetch networks from indra cogex and visualize with cytoscape

DESCRIPTION

@@ -34,14 +34,19 @@ Imports:
     graphics,
     methods,
     statmod,
-    parallel
+    parallel,
+    RCy3,
+    httr,
+    jsonlite,
+    r2r
 Suggests: 
     BiocStyle,
     knitr,
     rmarkdown,
     tinytest,
     covr,
-    markdown
+    markdown,
+    testthat (>= 3.0.0)
 VignetteBuilder: knitr
 biocViews: ImmunoOncology, MassSpectrometry, Proteomics, Software, Normalization,
         QualityControl, TimeCourse
@@ -55,3 +60,4 @@ Packaged: 2017-10-20 02:13:12 UTC; meenachoi
 LinkingTo: 
     Rcpp,
     RcppArmadillo
+Config/testthat/edition: 3

NAMESPACE

@@ -39,13 +39,15 @@ export(extractSDRF)
 export(getProcessed)
 export(getSamplesInfo)
 export(getSelectedProteins)
+export(getSubnetworkFromIndra)
 export(groupComparison)
 export(groupComparisonPlots)
 export(makePeptidesDictionary)
 export(modelBasedQCPlots)
 export(quantification)
 export(savePlot)
 export(theme_msstats)
+export(visualizeSubnetwork)
 import(data.table)
 import(ggplot2)
 import(limma)
@@ -57,6 +59,10 @@ importFrom(MSstatsConvert,MSstatsImport)
 importFrom(MSstatsConvert,MSstatsLogsSettings)
 importFrom(MSstatsConvert,MSstatsMakeAnnotation)
 importFrom(MSstatsConvert,MSstatsPreprocess)
+importFrom(RCy3,createNetworkFromDataFrames)
+importFrom(RCy3,createVisualStyle)
+importFrom(RCy3,mapVisualProperty)
+importFrom(RCy3,setVisualStyle)
 importFrom(Rcpp,sourceCpp)
 importFrom(data.table,as.data.table)
 importFrom(data.table,melt)
@@ -79,6 +85,10 @@ importFrom(graphics,title)
 importFrom(htmltools,div)
 importFrom(htmltools,save_html)
 importFrom(htmltools,tagList)
+importFrom(httr,POST)
+importFrom(httr,add_headers)
+importFrom(httr,content)
+importFrom(jsonlite,toJSON)
 importFrom(limma,squeezeVar)
 importFrom(lme4,lmer)
 importFrom(marray,maPalette)
@@ -94,6 +104,9 @@ importFrom(plotly,plot_ly)
 importFrom(plotly,style)
 importFrom(plotly,subplot)
 importFrom(preprocessCore,normalize.quantiles)
+importFrom(r2r,hashmap)
+importFrom(r2r,keys)
+importFrom(r2r,query)
 importFrom(stats,dist)
 importFrom(stats,fitted)
 importFrom(stats,formula)

R/getSubnetworkFromIndra.R

@@ -0,0 +1,26 @@
+#' Get subnetwork from INDRA database with differential analysis results
+#' 
+#' @param input groupComparison comparisionResult table with additional HGNC ID 
+#' and HGNC name columns
+#' @param pvalue_cutoff p-value cutoff for filtering. Default is NULL, i.e. no 
+#' filtering
+#' 
+#' @return list of 2 data.frames, nodes and edges
+#' 
+#' @export
+#' 
+#' @examples 
+#' input = data.table::fread(system.file("tinytest/processed_data/groupComparisonModel.csv", 
+#'                              package = "MSstats"))
+#' # subnetwork = getSubnetworkFromIndra(input, pvalue_cutoff = 0.05)
+#' # head(subnetwork$nodes)
+#' # head(subnetwork$edges)
+#'
+#'   
+getSubnetworkFromIndra = function(input, pvalue_cutoff = NULL) {
+    input = .filterGetSubnetworkFromIndraInput(input, pvalue_cutoff)
+    res = .callIndraCogexApi(input$HgncId)
+    nodes = .constructNodesDataFrame(input)
+    edges = .constructEdgesDataFrame(res, input)
+    return(list(nodes=nodes, edges=edges))
+}

R/utils_getSubnetworkFromIndra.R

@@ -0,0 +1,111 @@
+#' Call INDRA Cogex API and return response
+#' @param hgnc_ids list of hgnc_ids
+#' @return list of INDRA statements
+#' @importFrom jsonlite toJSON
+#' @importFrom httr POST add_headers content
+#' @keywords internal
+#' @noRd
+.callIndraCogexApi = function(hgnc_ids) {
+    INDRA_COGEX_URL = "https://discovery.indra.bio/api/indra_subnetwork_relations"
+
+    groundings = lapply(hgnc_ids, function(x) list("HGNC", x))
+    groundings = list(nodes = groundings)
+    groundings = jsonlite::toJSON(groundings, auto_unbox = TRUE)
+
+    res = POST(
+        INDRA_COGEX_URL, 
+        body = groundings, 
+        add_headers("Content-Type" = "application/json"), 
+        encode = "raw"
+    )
+    res = content(res)
+    return(res)
+}
+
+#' Filter groupComparison result input based on user-defined cutoffs
+#' @param input groupComparison result
+#' @param pvalue_cutoff p-value cutoff
+#' @return filtered groupComparison result
+#' @keywords internal
+#' @noRd
+.filterGetSubnetworkFromIndraInput = function(input, pvalue_cutoff) {
+    if (!is.null(pvalue_cutoff)) {
+        input = input[input$adj.pvalue < pvalue_cutoff,]
+    }
+    input = input[is.na(input$issue),]
+    return(input)
+}
+
+#' Add additional metadata to an edge
+#' @param edge object representation of an INDRA statement
+#' @param input filtered groupComparison result
+#' @return edge with additional metadata
+#' @keywords internal
+#' @noRd
+.addAdditionalMetadataToIndraEdge = function(edge, input) {
+    edge$evidence_list = paste(
+        "https://db.indra.bio/statements/from_agents?subject=", 
+        edge$source_id, "@HGNC&object=", 
+        edge$target_id, "@HGNC&type=",
+        edge$data$stmt_type, "&format=html", sep="")
+    edge$source_uniprot_id = input[input$HgncId == edge$source_id,]$Protein
+    edge$target_uniprot_id = input[input$HgncId == edge$target_id,]$Protein
+    return(edge)
+}
+
+
+#' Collapse duplicate INDRA statements into a mapping of edge to metadata
+#' @param res INDRA response
+#' @param input filtered groupComparison result
+#' @importFrom r2r hashmap keys
+#' @return processed edge to metadata mapping
+#' @keywords internal
+#' @noRd
+.collapseDuplicateEdgesIntoEdgeToMetadataMapping = function(res, input) {
+    edge_to_metadata_mapping = hashmap()
+
+    for (edge in res) {
+        key = paste(edge$source_id, edge$target_id, edge$data$stmt_type, sep="_")
+        if (key %in% keys(edge_to_metadata_mapping)) {
+            edge_to_metadata_mapping[[key]]$data$evidence_count = 
+                edge_to_metadata_mapping[[key]]$data$evidence_count + 
+                    edge$data$evidence_count
+        } else {
+            edge = .addAdditionalMetadataToIndraEdge(edge, input)
+            edge_to_metadata_mapping[[key]] = edge
+        }
+    }
+
+    return(edge_to_metadata_mapping)
+}
+
+#' Construct edges data.frame from INDRA response
+#' @param res INDRA response
+#' @param input filtered groupComparison result
+#' @importFrom r2r query keys
+#' @return edge data.frame
+#' @keywords internal
+#' @noRd
+.constructEdgesDataFrame = function(res, input) {
+    res = .collapseDuplicateEdgesIntoEdgeToMetadataMapping(res, input)
+    edges = data.frame(source=sapply(keys(res), function(x) query(res, x)$source_uniprot_id),
+                       target=sapply(keys(res), function(x) query(res, x)$target_uniprot_id),
+                       interaction=sapply(keys(res), function(x) query(res, x)$data$stmt_type), 
+                       evidenceCount=sapply(keys(res), function(x) query(res, x)$data$evidence_count),
+                       evidenceLink=sapply(keys(res), function(x) query(res, x)$evidence_list),
+                       stringsAsFactors=FALSE)
+    return(edges)
+}
+
+#' Construct nodes data.frame from groupComparison output
+#' @param input filtered groupComparison result
+#' @return nodes data.frame
+#' @keywords internal
+#' @noRd
+.constructNodesDataFrame = function(input) {
+    nodes = data.frame(id=input$Protein,
+                       logFC=input$log2FC, 
+                       pvalue=input$adj.pvalue,
+                       stringsAsFactors=FALSE)
+    return(nodes)
+}

R/visualizeSubnetwork.R

@@ -0,0 +1,47 @@
+#' Create visualization of subnetwork in cytoscape
+#' 
+#' @param nodes dataframe of nodes
+#' @param edges dataframe of edges
+#' @param pvalue_cutoff p-value cutoff for coloring significant proteins. Default is 0.05
+#' @param logfc_cutoff log fold change cutoff for coloring significant proteins. Default is 0.5
+#' @importFrom RCy3 createNetworkFromDataFrames mapVisualProperty createVisualStyle setVisualStyle
+#' 
+#' @export
+#' 
+#' @examples 
+#' input = data.table::fread(system.file("tinytest/processed_data/groupComparisonModel.csv", 
+#'                              package = "MSstats"))
+#' # subnetwork = getSubnetworkFromIndra(input)
+#' # visualizeSubnetwork(subnetwork$nodes, subnetwork$edges)
+#'
+#'               
+visualizeSubnetwork = function(nodes, edges, 
+                               pvalue_cutoff = 0.05, logfc_cutoff = 0.5) {
+
+    # Add additional columns for visualization
+    nodes$logFC_color = nodes$logFC
+    nodes$logFC_color[nodes$pvalue > pvalue_cutoff | 
+                          abs(nodes$logFC) < logfc_cutoff] = 0
+
+    # Create network
+    createNetworkFromDataFrames(nodes, edges)
+
+    # Apply visual style
+    DEFAULT_VISUAL_STYLE = list(
+        NODE_SHAPE="ROUNDRECT",
+        NODE_SIZE=50,
+        NODE_LABEL_FONT_SIZE=6,
+        NODE_LABEL_POSITION="center",
+        EDGE_TARGET_ARROW_SHAPE="Arrow")
+    VISUAL_STYLE_NAME = "MSstats-Indra Visual Style"
+
+    VISUAL_STYLE_MAPPINGS = list(
+        mapVisualProperty('Node Label','id','p'),
+        mapVisualProperty('Node Fill Color','logFC_color','c',
+                          c(-logfc_cutoff, 0.0, logfc_cutoff),
+                          c('#5588DD','#5588DD','#D3D3D3','#DD8855','#DD8855'))
+
+    )
+    createVisualStyle(VISUAL_STYLE_NAME, DEFAULT_VISUAL_STYLE, VISUAL_STYLE_MAPPINGS)
+    setVisualStyle(VISUAL_STYLE_NAME)
+}

inst/tinytest/processed_data/groupComparisonModel.csv

@@ -0,0 +1,5 @@
+,Protein,Label,log2FC,SE,Tvalue,DF,pvalue,adj.pvalue,issue,MissingPercentage,ImputationPercentage,HgncId
+1,BRD2_HUMAN,DMSO-DbET6,2.046185244,0.114338622,17.89583616,260,0,0,NA,0.310067114,0,1103
+2,BRD3_HUMAN,DMSO-DbET6,3.333427936,0.126570545,26.33652185,257,0,0,NA,0.252348993,0,1104
+3,BRD4_HUMAN,DMSO-DbET6,2.668934662,0.101282782,26.35131665,257,0,0,NA,0.118120805,0,13575
+4,CLH1_HUMAN,DMSO-DbET7,0.668934662,0.101282782,2.35131665,257,0.99,0.99,NA,0.118120805,0,2092

tests/testthat.R

@@ -0,0 +1,12 @@
+# This file is part of the standard setup for testthat.
+# It is recommended that you do not modify it.
+#
+# Where should you do additional test configuration?
+# Learn more about the roles of various files in:
+# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview
+# * https://testthat.r-lib.org/articles/special-files.html
+
+library(testthat)
+library(MSstats)
+
+test_check("MSstats")

tests/testthat/test-getSubnetworkFromIndra.R

@@ -0,0 +1,23 @@
+test_that("getSubnetworkFromIndra works correctly", {
+    input = data.table::fread(
+        system.file("tinytest/processed_data/groupComparisonModel.csv", package = "MSstats")
+    )
+    local_mocked_bindings(.callIndraCogexApi = function(x) {
+        return(readRDS(system.file("tinytest/processed_data/indraResponse.rds", package = "MSstats")))
+    })
+    subnetwork = getSubnetworkFromIndra(input)
+    expect_equal(nrow(subnetwork$nodes), 4)
+    expect_equal(nrow(subnetwork$edges), 16)
+})
+
+test_that("getSubnetworkFromIndra with pvalue filter works correctly", {
+    input = data.table::fread(
+        system.file("tinytest/processed_data/groupComparisonModel.csv", package = "MSstats")
+    )
+    local_mocked_bindings(.callIndraCogexApi = function(x) {
+        return(readRDS(system.file("tinytest/processed_data/indraResponse.rds", package = "MSstats")))
+    })
+    subnetwork = getSubnetworkFromIndra(input, pvalue_cutoff = 0.05)
+    expect_equal(nrow(subnetwork$nodes), 3)
+    expect_equal(nrow(subnetwork$edges), 16)
+})