rm(list = ls())
library(Seurat)
library(Matrix)Read in Seurat Object file.
obj = readRDS("Output/SeuratObject.RDS")First, extract the cluster-specific expression data. As an example, we will use data from cells in cluster 5.
clusterid = 5
#extract data and filter lowly expressed genes
expr_dat = GetAssayData(obj[, Idents(obj) == clusterid], slot = 'data')
expr_dat = expr_dat[rowSums(expr_dat > 0) >= 0.2 * ncol(expr_dat), ]
expr_dat = expr_dat[apply(expr_dat, 1, function(x) length(unique(x))) > 3, ]Prepare output directory.
outdir = paste0('BN/Cluster_', clusterid, '/')
if(!file.exists(outdir)) dir.create(outdir, recursive = TRUE)
setwd(outdir)Define a discretization function. Here we will use k-means clustering (k = 3) to discretize the expression data.
discretizeK = function(x, k = 3){
y = kmeans(x, k)
order(y$centers)[y$cluster] - 1
}
expr_dat_3 <- t(apply(expr_dat, 1, discretizeK))Write discretized data into file data.discretized.txt.
write.table(expr_dat_3, file = paste0(outdir, 'data.discretized.txt'), col.names = FALSE, row.names = TRUE, quote = FALSE, sep = "\t")Now run the RIMBANet analysis with data.discretized.txt by following a pipeline at https://github.com/mw201608/BayesianNetwork
sessionInfo()## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: CentOS Linux 7 (Core)
##
## Matrix products: default
## BLAS/LAPACK: /hpc/packages/minerva-common/intel/parallel_studio_xe_2018/compilers_and_libraries_2018.1.163/linux/mkl/lib/intel64_lin/libmkl_gf_lp64.so
##
## locale:
## [1] C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] Matrix_1.3-2 Seurat_3.1.4
##
## loaded via a namespace (and not attached):
## [1] TH.data_1.0-10 Rtsne_0.15 colorspace_2.0-1
## [4] ellipsis_0.3.2 ggridges_0.5.1 leiden_0.3.1
## [7] listenv_0.8.0 npsurv_0.4-0 ggrepel_0.9.1
## [10] fansi_0.5.0 mvtnorm_1.0-11 codetools_0.2-16
## [13] splines_3.6.0 mnormt_1.5-5 lsei_1.2-0
## [16] knitr_1.26 TFisher_0.2.0 jsonlite_1.7.2
## [19] ica_1.0-2 cluster_2.0.8 png_0.1-7
## [22] uwot_0.1.5 sctransform_0.2.1 compiler_3.6.0
## [25] httr_1.4.1 assertthat_0.2.1 lazyeval_0.2.2
## [28] htmltools_0.5.0 tools_3.6.0 rsvd_1.0.2
## [31] igraph_1.2.6 gtable_0.3.0 glue_1.4.2
## [34] RANN_2.6.1 reshape2_1.4.4 dplyr_1.0.6
## [37] rappdirs_0.3.3 Rcpp_1.0.6 Biobase_2.46.0
## [40] vctrs_0.3.8 multtest_2.42.0 gdata_2.18.0
## [43] ape_5.3 nlme_3.1-139 gbRd_0.4-11
## [46] lmtest_0.9-37 xfun_0.12 stringr_1.4.0
## [49] globals_0.14.0 lifecycle_1.0.0 irlba_2.3.3
## [52] gtools_3.8.1 future_1.21.0 MASS_7.3-51.4
## [55] zoo_1.8-6 scales_1.1.1 parallel_3.6.0
## [58] sandwich_2.5-1 RColorBrewer_1.1-2 yaml_2.2.1
## [61] reticulate_1.20 pbapply_1.4-3 gridExtra_2.3
## [64] ggplot2_3.3.3 stringi_1.6.2 mutoss_0.1-12
## [67] plotrix_3.7-7 caTools_1.18.0 BiocGenerics_0.32.0
## [70] bibtex_0.4.2.1 Rdpack_0.11-1 rlang_0.4.11
## [73] pkgconfig_2.0.3 bitops_1.0-6 evaluate_0.14
## [76] lattice_0.20-38 ROCR_1.0-7 purrr_0.3.4
## [79] patchwork_1.1.1 htmlwidgets_1.5.1 cowplot_1.1.1
## [82] tidyselect_1.1.1 parallelly_1.25.0 RcppAnnoy_0.0.14
## [85] plyr_1.8.6 magrittr_2.0.1 R6_2.5.0
## [88] gplots_3.0.3 generics_0.1.0 multcomp_1.4-11
## [91] DBI_1.1.0 pillar_1.6.1 sn_1.5-4
## [94] fitdistrplus_1.0-14 survival_2.44-1.1 tibble_3.1.2
## [97] future.apply_1.7.0 tsne_0.1-3 crayon_1.4.1
## [100] KernSmooth_2.23-15 utf8_1.2.1 plotly_4.9.2.1
## [103] rmarkdown_1.12 grid_3.6.0 data.table_1.12.8
## [106] blob_1.2.0 metap_1.2 digest_0.6.27
## [109] tidyr_1.1.3 numDeriv_2016.8-1.1 RcppParallel_4.4.4
## [112] stats4_3.6.0 munsell_0.5.0 viridisLite_0.4.0