Last updated: 2018-07-14

Code version: eed7f05


Load previously computed results

results_filt_eval_top10 <- readRDS("../data/results/results_filt_eval_top10.rds")
results_filt_eval_top101 <- readRDS("../data/results/results_filt_eval_top101.rds")

After removing noisy labels, on top 10 cyclical genes

source("../peco/R/utility.R")
source("../peco/R/fit.trendfilter.generic.R")
source("../peco/R/run_seurat.R")

mean(results_filt_eval_top10$fit.supervised$diff_time)/2/pi
mean(results_filt_eval_top10$fit.trend2.unsup$diff_time)/2/pi
mean(results_filt_eval_top10$fit.seurat$diff_time)/2/pi


pve.wrapper <- function(results_list, methods_list) {
  res <- lapply(1:length(results_list), 
       function(i) {
          obj <- results_list[[i]]
          out <- data.frame(
            dapi=with(obj, get.pve(dapi[order(pred_time_shift)])$pve),
            gfp=with(obj, get.pve(gfp[order(pred_time_shift)])$pve),
            rfp=with(obj, get.pve(rfp[order(pred_time_shift)])$pve) )
  })
  names(res) <- methods_list
  return(res)
}


results_list <- results_filt_eval_top10
methods_list <- sapply(names(results_list), function(x) strsplit(x, split=".", fixed=TRUE)[[1]][2])

pve_eval_top10 <- do.call(rbind, pve.wrapper(results_list=results_filt_eval_top10, 
                               methods_list=methods_list))
pve_eval_top10$genes_used <- "top10"
pve_eval_top10$methods <- methods_list
saveRDS(pve_eval_top10,
        "../output/method-eval-withheld-explore-removenoisy.Rmd/pve_eval_top10.rds")
pve_eval_top10 <- readRDS("../output/method-eval-withheld-explore-removenoisy.Rmd/pve_eval_top10.rds")
print(pve_eval_top10)
                 dapi          gfp       rfp genes_used    methods
supervised 0.07345013 0.0489577790 0.3138331      top10 supervised
trend2     0.14252546 0.0006007315 0.2895497      top10     trend2
bspline    0.06255537 0.0062581397 0.3511136      top10    bspline
loess      0.12301345 0.0166849294 0.2876617      top10      loess
seurat     0.09360356 0.0006058249 0.3565288      top10     seurat
get.aov(yy=results_filt_eval_top10$fit.seurat$dapi,
        xx=results_filt_eval_top10$fit.seurat$assignments)

get.aov(yy=results_filt_eval_top10$fit.seurat$gfp,
        xx=results_filt_eval_top10$fit.seurat$assignments)

get.aov(yy=results_filt_eval_top10$fit.seurat$rfp,
        xx=results_filt_eval_top10$fit.seurat$assignments)

Seurat scores vs cell time predictions

seurat.S.sup <- with(results_filt_eval_top10,
                     get.pve(fit.seurat$S[order(fit.supervised$pred_time_shift)]))
seurat.S.unsup <- with(results_filt_eval_top10,
                     get.pve(fit.seurat$S[order(fit.trend2.unsup$pred_time_shift)]))
seurat.G2M.sup <- with(results_filt_eval_top10,
                     get.pve(fit.seurat$G2M[order(fit.supervised$pred_time_shift)]))
seurat.G2M.unsup <- with(results_filt_eval_top10,
                     get.pve(fit.seurat$G2M[order(fit.trend2.unsup$pred_time_shift)]))

seurat.S.owntime <- with(results_filt_eval_top10,
                     get.pve(fit.seurat$S[order(fit.seurat$cell_times_est)]))
seurat.G2M.owntime <- with(results_filt_eval_top10,
                     get.pve(fit.seurat$G2M[order(fit.seurat$cell_times_est)]))

save(seurat.S.sup, seurat.S.unsup,
     seurat.G2M.sup, seurat.G2M.unsup,
     seurat.S.owntime, seurat.G2M.owntime,
     file = "../output/method-eval-withheld-explore-removenoisy.Rmd/seurat.time.top10.rda")
load(file="../output/method-eval-withheld-explore-removenoisy.Rmd/seurat.time.top10.rda")

c(seurat.S.sup$pve, seurat.G2M.sup$pve)
[1] 0.4109253 0.5839934
c(seurat.S.unsup$pve, seurat.G2M.unsup$pve)
[1] 0.3196591 0.6121693
c(seurat.S.owntime$pve, seurat.G2M.owntime$pve)
[1] 0.8383764 0.8031149
# with(results_filt_eval_top10,
#   get.aov(yy=fit.seurat$G2M,xx=fit.seurat$assignments))

# with(results_filt_eval_top10,
#   get.aov(yy=fit.seurat$S,xx=fit.seurat$assignments))

Seurat clases vs supervsied cell times

cols <- c("orange", "red", "brown")
par(mfrow=c(1,3))
with(results_filt_eval_top10, 
     hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="G1"],
     nclass=5, col=cols[1], main = "G1, 9 cells", xlim=c(0,2*pi), ylim=c(0,30),
     xlab="Predicted cell time"))
with(results_filt_eval_top10, hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="S"],
     nclass=10, col=cols[2],
     main = "S, 53 cells", xlim=c(0,2*pi), ylim=c(0,30),
     xlab="Predicted cell time"))
with(results_filt_eval_top10, hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="G2M"],
     nclass=10, col=cols[3],
     main = "G2M, 55 cells", xlim=c(0,2*pi), ylim=c(0,30),
     xlab="Predicted cell time"))
title("Supervised cell time and Seurat classes", outer=TRUE, line=-1)

After removing moisy labels, on top 101 cyclical genes

source("../peco/R/utility.R")
source("../peco/R/fit.trendfilter.generic.R")
source("../peco/R/run_seurat.R")

mean(results_filt_eval_top101$fit.supervised$diff_time)/2/pi
mean(results_filt_eval_top101$fit.trend2.unsup$diff_time)/2/pi
mean(results_filt_eval_top101$fit.seurat$diff_time)/2/pi


pve.wrapper <- function(results_list, methods_list) {
  res <- lapply(1:length(results_list), 
       function(i) {
          obj <- results_list[[i]]
          out <- data.frame(
            dapi=with(obj, get.pve(dapi[order(pred_time_shift)])$pve),
            gfp=with(obj, get.pve(gfp[order(pred_time_shift)])$pve),
            rfp=with(obj, get.pve(rfp[order(pred_time_shift)])$pve) )
  })
  names(res) <- methods_list
  return(res)
}


results_list <- results_filt_eval_top101
methods_list <- sapply(names(results_list), function(x) strsplit(x, split=".", fixed=TRUE)[[1]][2])

pve_eval_top101 <- do.call(rbind, pve.wrapper(results_list=results_filt_eval_top101, 
                               methods_list=methods_list))
pve_eval_top101$genes_used <- "top101"
pve_eval_top101$methods <- methods_list
saveRDS(pve_eval_top101,
        "../output/method-eval-withheld-explore-removenoisy.Rmd/pve_eval_top101.rds")
pve_eval_top101 <- readRDS("../output/method-eval-withheld-explore-removenoisy.Rmd/pve_eval_top101.rds")
print(pve_eval_top101)
                   dapi          gfp       rfp genes_used    methods
supervised  0.012566735 0.0538205897 0.4034993     top101 supervised
trend2     -0.003921578 0.0045257798 0.3864573     top101     trend2
bspline     0.070241404 0.0015846912 0.4212742     top101    bspline
loess       0.025346711 0.0091878022 0.1472481     top101      loess
seurat      0.093603559 0.0006058249 0.3565288     top101     seurat
get.aov(yy=results_filt_eval_top101$fit.seurat$dapi,
        xx=results_filt_eval_top101$fit.seurat$assignments)

get.aov(yy=results_filt_eval_top101$fit.seurat$gfp,
        xx=results_filt_eval_top101$fit.seurat$assignments)

get.aov(yy=results_filt_eval_top101$fit.seurat$rfp,
        xx=results_filt_eval_top101$fit.seurat$assignments)

Seurat scores vs cell time predictions

seurat.S.sup <- with(results_filt_eval_top101,
                     get.pve(fit.seurat$S[order(fit.supervised$pred_time_shift)]))
seurat.S.unsup <- with(results_filt_eval_top101,
                     get.pve(fit.seurat$S[order(fit.trend2.unsup$pred_time_shift)]))
seurat.G2M.sup <- with(results_filt_eval_top101,
                     get.pve(fit.seurat$G2M[order(fit.supervised$pred_time_shift)]))
seurat.G2M.unsup <- with(results_filt_eval_top101,
                     get.pve(fit.seurat$G2M[order(fit.trend2.unsup$pred_time_shift)]))

seurat.S.owntime <- with(results_filt_eval_top101,
                     get.pve(fit.seurat$S[order(fit.seurat$cell_times_est)]))
seurat.G2M.owntime <- with(results_eval_top101,
                     get.pve(fit.seurat$G2M[order(fit.seurat$cell_times_est)]))

save(seurat.S.sup, seurat.S.unsup,
     seurat.G2M.sup, seurat.G2M.unsup,
     seurat.S.owntime, seurat.G2M.owntime,
     file = "../output/method-eval-withheld-explore-removenoisy.Rmd/seurat.time.top101.rda")
load(file="../output/method-eval-withheld-explore.Rmd/seurat.time.top101.rda")

c(seurat.S.sup$pve, seurat.G2M.sup$pve)
[1] 0.4832478 0.6548333
c(seurat.S.unsup$pve, seurat.G2M.unsup$pve)
[1] 0.5051351 0.5050090
c(seurat.S.owntime$pve, seurat.G2M.owntime$pve)
[1] 0.8272036 0.8565118
# with(results_filt_eval_top101,
#   get.aov(yy=fit.seurat$G2M,xx=fit.seurat$assignments))

# with(results_filt_eval_top101,
#   get.aov(yy=fit.seurat$S,xx=fit.seurat$assignments))

Seurat clases vs supervsied cell times

cols <- c("orange", "red", "brown")
par(mfrow=c(1,3))
with(results_eval_top101, 
     hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="G1"],
     nclass=5, col=cols[1], main = "G1, 9 cells", xlim=c(0,2*pi), ylim=c(0,30),
     xlab="Predicted cell time"))
with(results_eval_top101, hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="S"],
     nclass=10, col=cols[2],
     main = "S, 53 cells", xlim=c(0,2*pi), ylim=c(0,30),
     xlab="Predicted cell time"))
with(results_eval_top101, hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="G2M"],
     nclass=10, col=cols[3],
     main = "G2M, 55 cells", xlim=c(0,2*pi), ylim=c(0,30),
     xlab="Predicted cell time"))
title("Supervised cell time and Seurat classes", outer=TRUE, line=-1)

Session information

sessionInfo()
R version 3.4.3 (2017-11-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)

Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
 [9] LC_ADDRESS=C               LC_TELEPHONE=C            
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

loaded via a namespace (and not attached):
 [1] compiler_3.4.3  backports_1.1.2 magrittr_1.5    rprojroot_1.3-2
 [5] tools_3.4.3     htmltools_0.3.6 yaml_2.1.16     Rcpp_0.12.17   
 [9] stringi_1.1.6   rmarkdown_1.10  knitr_1.20      git2r_0.21.0   
[13] stringr_1.2.0   digest_0.6.15   evaluate_0.10.1

This R Markdown site was created with workflowr