added script to calculate diffusion map and pseudo diffusion time with the `destiny` package

sc_workflow/calc_dm_dpt.R

+#+ include=FALSE
+
+#***********************************************************************************************************************
+#' # Diffusion Map and Pseudo Diffusion Time
+#' Calculation of both diffusion map and pseudo time is done using the `destiny` package [https://bioconductor.org/packages/release/bioc/html/destiny.html](https://bioconductor.org/packages/release/bioc/html/destiny.html)
+#'
+#' Created by: `r system("whoami", intern=T)`
+#'
+#' Created at: `r format(Sys.Date(), format="%B %d %Y")`
+#***********************************************************************************************************************
+
+#+ include=FALSE
+suppressMessages(require(docopt))
+
+doc = '
+Diffusion map calculation. Please specify count data.
+Usage: calc_dm_dpt.R [options] <count_data>
+
+Options:
+--results_prefix <results_prefix>   result prefix which will be added to each output file [default: ]
+'
+
+
+# commandArgs <- function(x) c("../segerstolpe_counts.txt")
+opts = docopt(doc, commandArgs(TRUE))
+
+
+# loading packages -----------------------------------------------------------------------------------------------------
+devtools::source_url("https://raw.githubusercontent.com/holgerbrandl/datautils/v1.45/R/core_commons.R")
+load_pack(destiny)
+load_pack(knitr)
+load_pack(Biobase)
+load_pack(data.table)
+load_pack(plotly)
+
+
+# handling input data --------------------------------------------------------------------------------------------------
+results_prefix <- opts$results_prefix
+countData <- fread(opts$count_data)
+
+countMatrix <- countData %>%
+    column2rownames(colnames(countData)[str_detect(colnames(countData), "gene_id")]) %>%
+    as.matrix()
+countMatrix <- countMatrix[rowSums(countMatrix)>0, ]
+
+
+# Seurat-based normalization (https://rdrr.io/cran/Seurat/src/R/preprocessing.R) by using its LogNormalize function:
+normCountMatrix <- Seurat::LogNormalize(countMatrix, scale.factor = 10000, display.progress = T) %>% as.matrix()
+
+rm(countData)
+rm(countMatrix)
+
+# calculate and save diffusion map and pseudo diffusion time plot
+es <- ExpressionSet(normCountMatrix)
+
+dm <- DiffusionMap(es)
+# saveRDS(dm, "dm.rds")
+# dm <- readRDS("dm.rds")
+
+dpt <- DPT(dm)
+# saveRDS(dpt, "dpt.rds")
+# dpt <- readRDS("dpt.rds")
+
+
+saveRDS(dm, "dm.rds")
+# readRDS("dm.rds")
+saveRDS(dpt, "dpt.rds")
+# readRDS("dpt.rds")
+
+
+#' ## Diffusion Map
+dm_data <- as.data.frame(dm)
+plot_ly(dm_data, x = ~DC1, y = ~DC2, z = ~DC3, size = I(5), type = "scatter3d")
+
+
+#' ## Pseudo Diffusion Time
+plot(dpt)
+
+
+#-----------------------------------------------------------------------------------------------------------------------
+# save copy of the used sc_quality_check.R script
+system("cp ${NGS_TOOLS}/sc_workflow/calc_dm_dpt.R ./.ngs_tools_calc_dm_dpt_$(date +%Y%m%d_%H%M).R")
+
+# get R version and package infos
+writeLines(capture.output(devtools::session_info()), ".sessionInfo.txt")
+
+session::save.session(".calc_dm_dpt.R.dat")
+# session::restore.session(".calc_dm_dpt.R.dat")
+
+
+
+