Pseudotemporal Analysis

Setup knitr and load utility functions

knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir="E:/DISC/reproducibility")

source_path = "./source/trajectory_analysis.r"
source(source_path)

Registered S3 method overwritten by 'dplyr':
  method           from
  print.rowwise_df     
Registered S3 methods overwritten by 'htmltools':
  method               from         
  print.html           tools:rstudio
  print.shiny.tag      tools:rstudio
  print.shiny.tag.list tools:rstudio
Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio
Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     
reldist: Relative Distribution Methods
Version 1.6-6 created on 2016-10-07.
copyright (c) 2003, Mark S. Handcock, University of California-Los Angeles
 For citation information, type citation("reldist").
 Type help(package="reldist") to get started.

搼㸴搼㸸挼㸸攼㹢搼㸰攼㸸Ҫ戼㸵ij̼愼㹤戼㸰昼㹣愼㸳戼㹡Biobase
搼㸴搼㸸挼㸸攼㹢搼㸰攼㸸Ҫ戼㸵ij̼愼㹤戼㸰昼㹣愼㸳戼㹡BiocGenerics

搼㸴搼㸸挼㸸攼㹢戼㸳̼愼㹤戼㸰昼㹣愼㸳戼㹡愼㸱愼㹥BiocGenerics愼㸱愼㹦

The following objects are masked from 愼㸱愼㹥package:parallel愼㸱愼㹦:

    clusterApply, clusterApplyLB, clusterCall, clusterEvalQ, clusterExport, clusterMap, parApply, parCapply, parLapply,
    parLapplyLB, parRapply, parSapply, parSapplyLB

The following object is masked from 愼㸱愼㹥package:Matrix愼㸱愼㹦:

    which

The following objects are masked from 愼㸱愼㹥package:stats愼㸱愼㹦:

    IQR, mad, sd, var, xtabs

The following objects are masked from 愼㸱愼㹥package:base愼㸱愼㹦:

    anyDuplicated, 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, 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")'.

搼㸴搼㸸挼㸸攼㹢搼㸰攼㸸Ҫ戼㸵ij̼愼㹤戼㸰昼㹣愼㸳戼㹡VGAM
戼㸳̼愼㹤戼㸰昼㹣愼㸱愼㹥VGAM愼㸱愼㹦挼㹡挼㸷搼㸳挼㸳R戼㸰汾3.6.3 挼㸰戼㸴戼㹤愼㸸搼㸴攼㹣戼㸵挼㸴搼㸴搼㸸挼㸸攼㹢搼㸰攼㸸Ҫ戼㸵ij̼愼㹤戼㸰昼㹣愼㸳戼㹡stats4
搼㸴搼㸸挼㸸攼㹢搼㸰攼㸸Ҫ戼㸵ij̼愼㹤戼㸰昼㹣愼㸳戼㹡splines
搼㸴搼㸸挼㸸攼㹢搼㸰攼㸸Ҫ戼㸵ij̼愼㹤戼㸰昼㹣愼㸳戼㹡DDRTree
戼㸳̼愼㹤戼㸰昼㹣愼㸱愼㹥DDRTree愼㸱愼㹦挼㹡挼㸷搼㸳挼㸳R戼㸰汾3.6.3 挼㸰戼㸴戼㹤愼㸸搼㸴攼㹣戼㸵挼㸴搼㸴搼㸸挼㸸攼㹢搼㸰攼㸸Ҫ戼㸵ij̼愼㹤戼㸰昼㹣愼㸳戼㹡irlba

Settings.

method_names = c("Raw", "DISC", "scVI", "MAGIC", "DCA", "scScope", "DeepImpute", "VIPER", "scImpute")
method_color = c("#A5A5A5", "#E83828", "#278BC4", "#EADE36", "#198B41", "#920783", "#F8B62D", "#8E5E32", "#1E2B68")
names(method_color) = method_names
text_color = rep("black", length(method_names))
names(text_color) = method_names
cell_type = readRDS("./data/BONE_MARROW/cell_type.rds")
raw_mat = readh5_loom("./data/BONE_MARROW/raw.loom")

[1] "./data/BONE_MARROW/raw.loom"
[1] 32738  6939

raw_filt = raw_mat[gene_selection(raw_mat, 10), ]
dim(raw_filt) #  13813, 6939

[1] 13813  6939

used_genes = rownames(raw_filt)
result_list = list()

The percentage of correct order is calculated for every possible root state.

for(ii in method_names){
  switch(ii, 
         "Raw" = {
           gene_bc_filt = raw_filt
         },
         "DISC" = {
           gene_bc_mat = readh5_loom(paste0("./data/BONE_MARROW/", ii, ".loom"))
           gene_bc_filt = gene_bc_mat[used_genes, ]
         }, {
           gene_bc_mat = readh5_imputation(paste0("./data/BONE_MARROW/", ii, ".hdf5"))
           gene_bc_filt = gene_bc_mat[used_genes, ]
         }
  )
  used_cells = names(cell_type)[!is.na(cell_type)]
  cds = get_cds_monocle2(gene_bc_filt[, used_cells])
  result_mat = get_score_monocle2(cds, cell_type, correct_order = correct_order_all, wrong_order = wrong_order_all, type_level = type_level, show_interactively = T)
  print(ii)
  print(result_mat[which.max(result_mat[, 2]), ])
  result_list[[ii]] = result_mat
}

Removing 49 outliers

[1] "Reducing dimension..."
[1] "Looking for the root state..."
[1] "Raw"
 pair_number          acc 
3.301487e+06 7.783153e-01 
[1] "./data/BONE_MARROW/DISC.loom"
[1] 32738  6939
[1] "Reducing dimension..."
[1] "Looking for the root state..."
[1] "DISC"
 pair_number          acc 
1.815709e+06 8.165014e-01 
[1] "./data/BONE_MARROW/scVI.hdf5"
[1] 13813  6939
[1] "Reducing dimension..."
[1] "Looking for the root state..."
[1] "scVI"
 pair_number          acc 
4.316787e+06 8.392429e-01 
[1] "./data/BONE_MARROW/MAGIC.hdf5"
[1] 13813  6939
[1] "Reducing dimension..."
[1] "Looking for the root state..."

错误: 无法分配大小为170.3 Mb的矢量

We show results of DISC and no imputation here since our pc cannot not handle this analysis for the imputation results of some methods, the complete results are shown in our paper which were using a 128GB machine for computation. The highest percentage of correct orders is used for the comparison among methods.

percentage_value = sapply(result_list, function(x){
  return(x[which.max(x[, 2]), ])
})["acc", c(1, 2)]


barplot_usage(percentage_value, main = "", cex.main = 1.5, bar_color = method_color[c(1, 2)], text_color = text_color[c(1, 2)], use_data_order = T, ylab = "Percentage", cex.lab = 1.5, font.main = 1, ylim = c(-0.1, 1), use_border = F)