netZooR application with TB dataset
Tian Wang
2019-10-16
Source:vignettes/ApplicationwithTBdataset.Rmd
ApplicationwithTBdataset.RmdIntroduction
netZooR is an R package which consists of seven main algorithms and is able to construct, analyse and plot gene regulatory networks.
PANDA(Passing Attributes between Networks for Data Assimilation) is a message-passing model to gene regulatory network reconstruction. It integrates multiple sources of biological data, including protein-protein interaction, gene expression, and transcription factor binding motifs information to reconstruct genome-wide, condition-specific regulatory networks.[Glass et al. 2013]
LIONESS(Linear Interpolation to Obtain Network Estimates for Single Samples) is a method to estimate sample-specific regulatory networks by applying linear interpolation to the predictions made by existing aggregate network inference approaches.[[Kuijjer et al. 2019]]](https://www.sciencedirect.com/science/article/pii/S2589004219300872)
CONDOR (COmplex Network Description Of Regulators) implements methods to cluster biapartite networks and estimatiing the contribution of each node to its community’s modularity.[Platig et al. 2016]
ALPACA(ALtered Partitions Across Community Architectures) is a method to compare two genome-scale networks derived from different phenotypic states to identify condition-specific modules. [Padi and Quackenbush 2018]
SAMBAR(Subtyping Agglomerated Mutations By Annotation Relations) is a method to identify subtypes based on somatic mutation data.[Kuijjer et al.].
MONSTER(Modeling Network State Transitions from Expression and Regulatory data)[Schlauch et al.]: infers transcription factor which drivers of cell state conditions at the gene regulatory network level.
OTTER(Optimization to Estimate Regulation) [publication in preparation]: models gene regulation estimation as a graph matrching problem
Installation
Prerequisites
Using this pacakage requires Python (3.X) and some Python libraries, R (>= 3.3.3), and stable Internet access.
Some plotting functions will require the Cytoscape installed.
Required Python libraries
How to install Python libraries depends varies from different platforms. More instructions could be find here.
The following Python libraries (or packages) are required by running PANDA and LIONESS algorithms:
The required Python packages are: pandas, numpy, networkx, matplotlib.pyplot.
Installing
This package could be downloaded via install_github()
function from devtools package.
# install.packages("devtools")
library(devtools)
# install netZooR pkg with vignettes, otherwise remove the "build_vignettes = TRUE" argument.
devtools::install_github("netZoo/netZooR", build_vignettes = TRUE)
library(viridisLite)#To visualize communitiesData Resources
Motif data
Here is some pre-prepared specie-sepcific PANDA-ready transcription factor binding motifs data stored in our AWS bucket https://s3.console.aws.amazon.com/s3/buckets/netzoo/netZooR/example_datasets/PANDA_ready_motif_prior/?region=us-east-2&tab=overview, which are derived from motif scan and motif info files located on https://sites.google.com/a/channing.harvard.edu/kimberlyglass/tools/resourcesby .
PPI
This package includes a function source.PPI may source a
Protein-Protein Interactions (PPI) througt STRING database given a list
of proteins of interest. The STRINGdb
is already loaded while loading netZooR.
# TF is a data frame with single column filled with TFs of Mycobacterium tuberculosis H37Rv.
PPI <- source.PPI(TF, STRING.version="10", species.index=83332, score_threshold=0)Running the sample TB datasets
library(netZooR)
#> Loading required package: igraph
#>
#> Attaching package: 'igraph'
#> The following objects are masked from 'package:stats':
#>
#> decompose, spectrum
#> The following object is masked from 'package:base':
#>
#> union
#> Loading required package: reticulate
#> Loading required package: pandaR
#> Loading required package: Biobase
#> Loading required package: BiocGenerics
#>
#> Attaching package: 'BiocGenerics'
#> The following objects are masked from 'package:igraph':
#>
#> normalize, path, union
#> The following objects are masked from 'package:stats':
#>
#> IQR, mad, sd, var, xtabs
#> The following objects are masked from 'package:base':
#>
#> anyDuplicated, aperm, append, as.data.frame, basename, cbind,
#> colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
#> get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
#> match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#> Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort,
#> table, tapply, union, unique, unsplit, which.max, which.min
#> Welcome to Bioconductor
#>
#> Vignettes contain introductory material; view with
#> 'browseVignettes()'. To cite Bioconductor, see
#> 'citation("Biobase")', and for packages 'citation("pkgname")'.
#> Loading required package: yarn
#> Setting options('download.file.method.GEOquery'='auto')
#> Setting options('GEOquery.inmemory.gpl'=FALSE)
#> Loading required package: matrixcalc
#>
#> Attaching package: 'matrixcalc'
#> The following object is masked from 'package:igraph':
#>
#> %s%
#>
#> Accessing the help pages for the usage of core functions.
?pandaPy
?createCondorObject
?pandaToCondorObject
?lionessPy
?alpaca
?pandaToAlpaca
?sambarThis package will invoke the Python in R environment through reticulate package. Configure which version of Python to use if necessary, here in netZooR, Python 3.X is required. More details can be found here
#check your Python configuration and the specific version of Python in use currently
py_config()
# reset to Python 3.X if necessary, like below:
use_python("/usr/local/bin/python3")The previous command is necessary to bind R to Python since we are calling PANDA from Python because netZooPy has an optimized implementation of PANDA. Check this tutorial for an example using a pure R implementation of PANDA. Use example data sets within package to test this package. Refer to four input datasets files: one TB expression dataset control group , one TB expression dataset treated, one transcription factor binding motifs dataset, and one protein-protein interaction datasets from either inst/extdat or AWS.
retrieve the file path of these files came with the netZooR package.
# retrieve the file path of these files
treated_expression_file_path <- system.file("extdata", "expr4_matched.txt", package = "netZooR", mustWork = TRUE)
control_expression_file_path <- system.file("extdata", "expr10_matched.txt", package = "netZooR", mustWork = TRUE)
motif_file_path <- system.file("extdata", "chip_matched.txt", package = "netZooR", mustWork = TRUE)
ppi_file_path <- system.file("extdata", "ppi_matched.txt", package = "netZooR", mustWork = TRUE)PANDA algorithm
Assign the above file paths to flag e(refers to
“expression dataset”), m(refers to “motif dataset”), and
ppi(refers to “PPI” dataset), respectively. Then set option
rm_missing to TRUE to run
PANDA to generate an aggregate network without
unmatched TF and genes.
Repeat with control group.
treated_all_panda_result <- pandaPy(expr_file = treated_expression_file_path, motif_file = motif_file_path, ppi_file= ppi_file_path,modeProcess="legacy", remove_missing = TRUE )
control_all_panda_result <- pandaPy(expr_file = control_expression_file_path,motif_file = motif_file_path, ppi_file= ppi_file_path,modeProcess="legacy", remove_missing = TRUE )Vector treated_all_panda_result and vector
control_all_panda_result below are large lists with three
elements: the entire PANDA network, indegree (“to” nodes) nodes and
score, outdegree (“from” nodes) nodes and score. Use
$panda,$indegree and $outdegree
to access each list item resepctively.
Use $pandato access the entire PANDA network.
treated_net <- treated_all_panda_result$panda
control_net <- control_all_panda_result$pandaPANDA Cytoscape Plotting
Cytoscape is an interactivity network visualization tool highly
recommanded to explore the PANDA network. Before using this function
plot.panda.in.cytoscape, please install and launch
Cytoscape (3.6.1 or greater) and keep it running whenever using.
LIONESS Algorithm
How to run LIONESS is mostly idential with method how to run PANDA in
this package, unless the return values of lionessPy() is a
data frame where first two columns represent TFs (regulators) and Genes
(targets) while the rest columns represent each sample. each cell filled
with estimated score calculated by LIONESS.
# Run LIONESS algorithm for the first two samples
# removing start_sample and end_sample arguments to generate whole LIONESS network with all samples.
control_lioness_result <- lionessPy(expr_file = control_expression_file_path,motif_file = motif_file_path, ppi_file= ppi_file_path,modeProcess="legacy", remove_missing = TRUE, start_sample=1, end_sample=2)CONDOR Algorithm and plotting
PANDA network can simply be converted into condor.object by
pandaToCondorObject(panda.net, threshold) Defaults option
threshold is the average of [median weight of non-prior
edges] and [median weight of prior edges], all weights mentioned
previous are transformationed with formula w'=ln(e^w+1)
before calculating the median and average. But all the edges selected
will remain the orginal weights calculated by PANDA.
treated_condor_object <- pandaToCondorObject(treated_net, threshold = 0)The communities structure can be plotted by igraph.
library(viridisLite)
treated_condor_object <-condorCluster(treated_condor_object,project = FALSE)
#> [1] "modularity of projected graph 0.308479519056006"
#> [1] "Q = 0.308566500564539"
#> [1] "Q = 0.308658883596709"
#> [1] "Q = 0.308658883596709"
treated_color_num <- max(treated_condor_object$red.memb$com)
treated_color <- viridis(treated_color_num, alpha = 1, begin = 0, end = 1, direction = 1, option = "D")
condorPlotCommunities(treated_condor_object, color_list=treated_color, point.size=0.04, xlab="Genes", ylab="TFs")
ALPACA Algorithm
ALPACA community structure can also be generated from two PANDA
network by pandaToAlpaca
alpaca<- pandaToAlpaca(treated_net, control_net, NULL, verbose=FALSE)
#> [1] "Detecting communities in control network..."
#> Weights detected. Building condor object with weighted edges.
#> [1] "modularity of projected graph 0.482825466382495"
#> [1] "Q = 0.49524047821125"
#> [1] "Q = 0.49524047821125"
#> [1] "Computing differential modularity matrix..."
#> [1] "Computing differential modules..."
#> [1] "Merging 2644 communities"
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] "Merging 138 communities"
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 5
#> [1] 6
#> [1] "Merging 23 communities"
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] "Merging 13 communities"
#> [1] 1
#> [1] "Computing node scores..."
#> [1] 1
#> [1] 2
#> [1] 3
#> [1] 4
#> [1] 5
#> [1] 6
#> [1] 7
#> [1] 8
#> [1] 9
#> [1] 10
#> [1] 11
#> [1] 12
#> [1] 13Information
sessionInfo()
#> R version 4.2.2 (2022-10-31)
#> Platform: x86_64-apple-darwin17.0 (64-bit)
#> Running under: macOS Big Sur ... 10.16
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.2/Resources/lib/libRlapack.dylib
#>
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] viridisLite_0.4.2 netZooR_1.5.0 matrixcalc_1.0-6
#> [4] yarn_1.24.0 pandaR_1.19.6 Biobase_2.58.0
#> [7] BiocGenerics_0.44.0 reticulate_1.32.0 igraph_1.5.1
#>
#> loaded via a namespace (and not attached):
#> [1] rappdirs_0.3.3 rtracklayer_1.58.0
#> [3] pbdZMQ_0.3-10 AnnotationForge_1.40.0
#> [5] ragg_1.2.5 tidyr_1.3.0
#> [7] bumphunter_1.40.0 minfi_1.44.0
#> [9] ggplot2_3.4.3 bit64_4.0.5
#> [11] knitr_1.44 DelayedArray_0.24.0
#> [13] data.table_1.14.8 KEGGREST_1.38.0
#> [15] RCurl_1.98-1.12 GEOquery_2.66.0
#> [17] doParallel_1.0.17 generics_0.1.3
#> [19] GenomicFeatures_1.50.3 preprocessCore_1.60.1
#> [21] RSQLite_2.3.1 chron_2.3-61
#> [23] bit_4.0.5 tzdb_0.4.0
#> [25] base64url_1.4 xml2_1.3.5
#> [27] SummarizedExperiment_1.28.0 assertthat_0.2.1
#> [29] STRINGdb_2.10.0 xfun_0.40
#> [31] hms_1.1.3 jquerylib_0.1.4
#> [33] evaluate_0.21 fansi_1.0.4
#> [35] restfulr_0.0.15 scrime_1.3.5
#> [37] progress_1.2.2 caTools_1.18.2
#> [39] dbplyr_2.3.3 Rgraphviz_2.42.0
#> [41] DBI_1.1.3 reshape_0.8.9
#> [43] tensorA_0.36.2 stats4_4.2.2
#> [45] purrr_1.0.2 hash_2.2.6.3
#> [47] dplyr_1.1.3 backports_1.4.1
#> [49] permute_0.9-7 annotate_1.76.0
#> [51] biomaRt_2.54.0 sparseMatrixStats_1.10.0
#> [53] MatrixGenerics_1.10.0 vctrs_0.6.3
#> [55] cmdstanr_0.6.1.9000 penalized_0.9-52
#> [57] here_1.0.1 abind_1.4-5
#> [59] cachem_1.0.8 withr_2.5.0
#> [61] checkmate_2.2.0 vegan_2.6-4
#> [63] GenomicAlignments_1.34.0 prettyunits_1.1.1
#> [65] mclust_6.0.0 cluster_2.1.4
#> [67] IRdisplay_1.1 crayon_1.5.2
#> [69] uchardet_1.1.1 genefilter_1.80.2
#> [71] edgeR_3.40.1 pkgconfig_2.0.3
#> [73] GenomeInfoDb_1.34.4 nlme_3.1-163
#> [75] nnet_7.3-19 rlang_1.1.1
#> [77] RJSONIO_1.3-1.8 lifecycle_1.0.3
#> [79] downloader_0.4 filelock_1.0.2
#> [81] BiocFileCache_2.6.0 GOstats_2.64.0
#> [83] quantro_1.32.0 distributional_0.3.2
#> [85] rprojroot_2.0.3 matrixStats_1.0.0
#> [87] graph_1.76.0 rngtools_1.5.2
#> [89] base64_2.0.1 Matrix_1.6-1.1
#> [91] IRkernel_1.3.2 Rhdf5lib_1.20.0
#> [93] base64enc_0.1-3 processx_3.8.2
#> [95] png_0.1-8 rjson_0.2.21
#> [97] bitops_1.0-7 KernSmooth_2.23-22
#> [99] rhdf5filters_1.10.0 Biostrings_2.66.0
#> [101] blob_1.2.4 DelayedMatrixStats_1.20.0
#> [103] doRNG_1.8.6 stringr_1.5.0
#> [105] nor1mix_1.3-0 readr_2.1.4
#> [107] S4Vectors_0.36.1 scales_1.2.1
#> [109] memoise_2.0.1 GSEABase_1.60.0
#> [111] magrittr_2.0.3 plyr_1.8.8
#> [113] hexbin_1.28.3 gplots_3.1.3
#> [115] zlibbioc_1.44.0 compiler_4.2.2
#> [117] BiocIO_1.8.0 RColorBrewer_1.1-3
#> [119] illuminaio_0.40.0 plotrix_3.8-2
#> [121] Rsamtools_2.14.0 cli_3.6.1
#> [123] XVector_0.38.0 Category_2.64.0
#> [125] ps_1.7.5 MASS_7.3-60
#> [127] mgcv_1.9-0 tidyselect_1.2.0
#> [129] stringi_1.7.12 textshaping_0.3.6
#> [131] yaml_2.3.7 askpass_1.2.0
#> [133] locfit_1.5-9.8 grid_4.2.2
#> [135] sass_0.4.7 tools_4.2.2
#> [137] parallel_4.2.2 rstudioapi_0.15.0
#> [139] uuid_1.1-1 foreach_1.5.2
#> [141] posterior_1.4.1 farver_2.1.1
#> [143] digest_0.6.33 proto_1.0.0
#> [145] quadprog_1.5-8 Rcpp_1.0.11
#> [147] GenomicRanges_1.50.2 siggenes_1.72.0
#> [149] org.Hs.eg.db_3.16.0 httr_1.4.7
#> [151] ggdendro_0.1.23 AnnotationDbi_1.60.0
#> [153] colorspace_2.1-0 XML_3.99-0.14
#> [155] fs_1.6.3 IRanges_2.32.0
#> [157] splines_4.2.2 RBGL_1.74.0
#> [159] pkgdown_2.0.7 multtest_2.54.0
#> [161] systemfonts_1.0.4 xtable_1.8-4
#> [163] jsonlite_1.8.7 R6_2.5.1
#> [165] RUnit_0.4.32 gsubfn_0.7
#> [167] pillar_1.9.0 htmltools_0.5.6
#> [169] glue_1.6.2 fastmap_1.1.1
#> [171] BiocParallel_1.32.4 beanplot_1.3.1
#> [173] codetools_0.2-19 utf8_1.2.3
#> [175] lattice_0.21-8 bslib_0.5.1
#> [177] tibble_3.2.1 sqldf_0.4-11
#> [179] curl_5.0.2 gtools_3.9.4
#> [181] GO.db_3.16.0 openssl_2.1.0
#> [183] survival_3.5-7 limma_3.54.0
#> [185] rmarkdown_2.25 repr_1.1.6
#> [187] desc_1.4.2 munsell_0.5.0
#> [189] rhdf5_2.42.0 GenomeInfoDbData_1.2.9
#> [191] iterators_1.0.14 RCy3_2.18.0
#> [193] HDF5Array_1.26.0 reshape2_1.4.4
#> [195] gtable_0.3.4Note
If there is an error like
Error in fetch(key) : lazy-load database.rdb' is corrupt
when accessing the help pages of functions in this package after being
loaded. It’s a
limitation of base R and has not been solved yet. Restart R session
and re-load this package will help.