diff --git a/DESCRIPTION b/DESCRIPTION
index 14170830..87a3c193 100644
--- a/DESCRIPTION
+++ b/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
diff --git a/NAMESPACE b/NAMESPACE
index c73e7eca..0c113c17 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -39,6 +39,7 @@ export(extractSDRF)
 export(getProcessed)
 export(getSamplesInfo)
 export(getSelectedProteins)
+export(getSubnetworkFromIndra)
 export(groupComparison)
 export(groupComparisonPlots)
 export(makePeptidesDictionary)
@@ -46,6 +47,7 @@ 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)
diff --git a/R/getSubnetworkFromIndra.R b/R/getSubnetworkFromIndra.R
new file mode 100644
index 00000000..eb3ac7c0
--- /dev/null
+++ b/R/getSubnetworkFromIndra.R
@@ -0,0 +1,40 @@
+#' 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
+#' @param stmt_types list of statement types to filter, default is c("Complex")
+#' @importFrom r2r hashmap query keys
+#' 
+#' @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)
+#' # head(subnetwork$nodes)
+#' # head(subnetwork$edges)
+#'
+#'   
+getSubnetworkFromIndra = function(input, pvalue_cutoff = 0.05, stmt_types = c("Complex")) {
+    input = .filterAndProcessInput(input, pvalue_cutoff)
+    res = .callIndraCogexApi(input$HgncId)
+    res = .filterIndraResults(res, stmt_types)
+    res = .collapseIndraResultsIntoEdgeToEvidenceMapping(res, input)
+    
+    nodes = data.frame(id=input$HgncId,
+                       uniprot_id=input$Protein,
+                       logFC=input$log2FC, 
+                       pvalue=input$adj.pvalue,
+                       stringsAsFactors=FALSE)
+    edges = data.frame(source=sapply(keys(res), function(x) query(res, x)$source_id),
+                       target=sapply(keys(res), function(x) query(res, x)$target_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)$data$evidence_list),
+                       stringsAsFactors=FALSE)
+    
+    return(list(nodes=nodes, edges=edges))
+}
\ No newline at end of file
diff --git a/R/utils_getSubnetworkFromIndra.R b/R/utils_getSubnetworkFromIndra.R
new file mode 100644
index 00000000..cab2511f
--- /dev/null
+++ b/R/utils_getSubnetworkFromIndra.R
@@ -0,0 +1,78 @@
+#' 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 and process 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
+.filterAndProcessInput = function(input, pvalue_cutoff) {
+    input = input[input$adj.pvalue < pvalue_cutoff,]
+    input = input[is.na(input$issue),]
+    input$HgncId = as.character(input$HgncId)
+    return(input)
+}
+
+#' Filter INDRA results based on user preferences
+#' @param res INDRA response
+#' @param stmt_types statement types filter
+#' @return filtered INDRA response
+#' @keywords internal
+#' @noRd
+.filterIndraResults = function(res, stmt_types) {
+    res = Filter(function(x) x$data$stmt_type %in% stmt_types, res)
+    return(res)
+}
+
+#' Collapse INDRA statements into a mapping of edge to evidence data
+#' @param res INDRA response
+#' @param input filtered groupComparison result
+#' @return processed INDRA statements
+#' @keywords internal
+#' @noRd
+.collapseIndraResultsIntoEdgeToEvidenceMapping = function(res, input) {
+    gene_id_map = hashmap()
+    gene_id_map[input$HgncId] = input$HgncName
+    
+    edge_to_evidence_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_evidence_mapping)) {
+            # Collapse duplicate edges and keep track of evidence count
+            edge_to_evidence_mapping[[key]]$data$evidence_count = 
+                edge_to_evidence_mapping[[key]]$data$evidence_count + edge$data$evidence_count
+        } else {
+            # Add INDRA DB link for new edges
+            edge$data$evidence_list = paste(
+                "https://db.indra.bio/statements/from_agents?subject=", 
+                query(gene_id_map, edge$source_id), "&object=", 
+                query(gene_id_map, edge$target_id), "&type=",
+                edge$data$stmt_type, "&format=html", sep="")
+            edge_to_evidence_mapping[[key]] = edge
+        }
+    }
+    
+    return(edge_to_evidence_mapping)
+}
\ No newline at end of file
diff --git a/R/visualizeSubnetwork.R b/R/visualizeSubnetwork.R
new file mode 100644
index 00000000..c1288503
--- /dev/null
+++ b/R/visualizeSubnetwork.R
@@ -0,0 +1,33 @@
+#' Create visualization of subnetwork in cytoscape
+#' 
+#' @param nodes dataframe of nodes
+#' @param edges dataframe of edges
+#' @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) {
+    DEFAULT_VISUAL_STYLE = list(
+        NODE_FILL_COLOR="lightblue",
+        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"
+    
+    createNetworkFromDataFrames(nodes, edges)
+    
+    VISUAL_STYLE_MAPPINGS = list(
+        mapVisualProperty('Node Label','uniprot_id','p')
+    )
+    createVisualStyle(VISUAL_STYLE_NAME, DEFAULT_VISUAL_STYLE, VISUAL_STYLE_MAPPINGS)
+    setVisualStyle(VISUAL_STYLE_NAME)
+}
\ No newline at end of file
diff --git a/inst/tinytest/processed_data/groupComparisonModel.csv b/inst/tinytest/processed_data/groupComparisonModel.csv
new file mode 100644
index 00000000..cafef55c
--- /dev/null
+++ b/inst/tinytest/processed_data/groupComparisonModel.csv
@@ -0,0 +1,4 @@
+,Protein,Label,log2FC,SE,Tvalue,DF,pvalue,adj.pvalue,issue,MissingPercentage,ImputationPercentage,HgncId,HgncName
+1,BRD2_HUMAN,DMSO-DbET6,2.046185244,0.114338622,17.89583616,260,0,0,NA,0.310067114,0,1103,BRD2
+2,BRD3_HUMAN,DMSO-DbET6,3.333427936,0.126570545,26.33652185,257,0,0,NA,0.252348993,0,1104,BRD3
+3,BRD4_HUMAN,DMSO-DbET6,2.668934662,0.101282782,26.35131665,257,0,0,NA,0.118120805,0,13575,BRD4
\ No newline at end of file
diff --git a/inst/tinytest/processed_data/indraResponse.rds b/inst/tinytest/processed_data/indraResponse.rds
new file mode 100644
index 00000000..8ba5b2d2
Binary files /dev/null and b/inst/tinytest/processed_data/indraResponse.rds differ
diff --git a/man/getSubnetworkFromIndra.Rd b/man/getSubnetworkFromIndra.Rd
new file mode 100644
index 00000000..a74546eb
--- /dev/null
+++ b/man/getSubnetworkFromIndra.Rd
@@ -0,0 +1,31 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/getSubnetworkFromIndra.R
+\name{getSubnetworkFromIndra}
+\alias{getSubnetworkFromIndra}
+\title{Get subnetwork from INDRA database with differential analysis results}
+\usage{
+getSubnetworkFromIndra(input, pvalue_cutoff = 0.05, stmt_types = c("Complex"))
+}
+\arguments{
+\item{input}{groupComparison comparisionResult table with additional HGNC ID 
+and HGNC name columns}
+
+\item{pvalue_cutoff}{p-value cutoff for filtering}
+
+\item{stmt_types}{list of statement types to filter, default is c("Complex")}
+}
+\value{
+list of 2 data.frames, nodes and edges
+}
+\description{
+Get subnetwork from INDRA database with differential analysis results
+}
+\examples{
+input = data.table::fread(system.file("tinytest/processed_data/groupComparisonModel.csv", 
+                             package = "MSstats"))
+# subnetwork = getSubnetworkFromIndra(input)
+# head(subnetwork$nodes)
+# head(subnetwork$edges)
+
+  
+}
diff --git a/man/visualizeSubnetwork.Rd b/man/visualizeSubnetwork.Rd
new file mode 100644
index 00000000..a065569f
--- /dev/null
+++ b/man/visualizeSubnetwork.Rd
@@ -0,0 +1,24 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/visualizeSubnetwork.R
+\name{visualizeSubnetwork}
+\alias{visualizeSubnetwork}
+\title{Create visualization of subnetwork in cytoscape}
+\usage{
+visualizeSubnetwork(nodes, edges)
+}
+\arguments{
+\item{nodes}{dataframe of nodes}
+
+\item{edges}{dataframe of edges}
+}
+\description{
+Create visualization of subnetwork in cytoscape
+}
+\examples{
+input = data.table::fread(system.file("tinytest/processed_data/groupComparisonModel.csv", 
+                             package = "MSstats"))
+# subnetwork = getSubnetworkFromIndra(input)
+# visualizeSubnetwork(subnetwork$nodes, subnetwork$edges)
+
+              
+}
diff --git a/tests/testthat.R b/tests/testthat.R
new file mode 100644
index 00000000..73f1cfff
--- /dev/null
+++ b/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")
diff --git a/tests/testthat/test-getSubnetworkFromIndra.R b/tests/testthat/test-getSubnetworkFromIndra.R
new file mode 100644
index 00000000..a17ab528
--- /dev/null
+++ b/tests/testthat/test-getSubnetworkFromIndra.R
@@ -0,0 +1,11 @@
+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), 3)
+    expect_equal(nrow(subnetwork$edges), 6)
+})
\ No newline at end of file