Evaluate the withheld sample on top X cyclical

Run on supervised model
Load previously computed results
Before removing moisy labels, on top 10 cyclical genes
Before removing moisy labels, on top 101 cyclical genes
Session information

Last updated: 2018-07-14

Code version: c115f90

Run on supervised model

test cluster code

ngenes=101

dir <-"/project2/gilad/joycehsiao/fucci-seq"
source(file.path(dir,"code/working/run_methods.R"))

data_training <- readRDS(file=file.path(dir, "data/results/data_training.rds"))
data_withheld <-readRDS(file=file.path(dir, "data/results/data_withheld.rds"))

sig.genes <- readRDS(file=file.path(dir,
    "output/npreg-trendfilter-quantile.Rmd/out.stats.ordered.sig.476.rds"))

# make prediction parameters
Y_train_topX <- data_training$log2cpm.quant.nonvalid[
  rownames(data_training$log2cpm.quant.nonvalid) %in% rownames(sig.genes)[1:ngenes], ]

training_topX <- cycle.npreg.insample(Y = Y_train_topX,
                                        theta = data_training$theta.nonvalid,
                                        polyorder=2,
                                        ncores=15,
                                        method.trend="trendfilter")


seurat.genes <- readLines(
  con = file.path(dir,                                    "data/cellcycle-genes-previous-studies/seurat_cellcycle/regev_lab_cell_cycle_genes.txt"))
seurat.genes <- list(s.genes=seurat.genes[1:43],
                     g2m.genes=seurat.genes[44:97])

training_model <- training_topX
Y_test_normed <- data_withheld$log2cpm.quant.valid
  
cycle.genes <- rownames(training_model$Y)
Y_test_normed.cycle <- Y_test_normed[which(rownames(Y_test_normed) %in% cycle.genes),]
#Y_test.cycle <- Y_test[which(rownames(Y_test) %in% cycle.genes),]

### supervised methods
message("Begin supervised method...")
fit.supervised <- cycle.npreg.outsample(Y_test=Y_test_normed.cycle,
                                        sigma_est=training_model$sigma_est,
                                        funs_est=training_model$funs_est,
                                        method.grid = "uniform",
                                        method.trend="trendfilter",
                                        polyorder=2,
                                        ncores=ncores)
  
pred_time <- fit.supervised$cell_times_est
pred_time_shift <- with(fit.supervised, rotation(data_withheld$theta.valid,
                                                 cell_times_est)$y2shift)
#plot(pred_time, pred_time_shift)
source("../peco/R/utility.R")
pv <- get.pve(data_withheld$pdata.valid$gfp.median.log10sum.adjust[order(pred_time_shift)])


bb <- readRDS("../data/results/results_train_top101.rds")

bb1.pve <- get.pve(bb$fold.1$fit.supervised$dapi[order(bb$fold.1$fit.supervised$pred_time_shift)])
bb2.pve <- get.pve(bb$fold.2$fit.supervised$dapi[order(bb$fold.2$fit.supervised$pred_time_shift)])
bb3.pve <- get.pve(bb$fold.3$fit.supervised$dapi[order(bb$fold.3$fit.supervised$pred_time_shift)])
bb4.pve <- get.pve(bb$fold.4$fit.supervised$dapi[order(bb$fold.4$fit.supervised$pred_time_shift)])
bb5.pve <- get.pve(bb$fold.5$fit.supervised$dapi[order(bb$fold.5$fit.supervised$pred_time_shift)])


mean(c(bb1.pve$pve, bb2.pve$pve, bb3.pve$pve, bb4.pve$pve, bb5.pve$pve))

Load previously computed results

results_eval_top101 <- readRDS("../data/results/results_eval_top101.rds")
results_eval_top10 <- readRDS("../data/results/results_eval_top10.rds")

Before removing moisy 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_eval_top10$fit.supervised$diff_time)/2/pi
mean(results_eval_top10$fit.trend2.unsup$diff_time)/2/pi
mean(results_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_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_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.Rmd/pve_eval_top10.rds")

pve_eval_top10 <- readRDS("../output/method-eval-withheld-explore.Rmd/pve_eval_top10.rds")
print(pve_eval_top10)

                  dapi          gfp       rfp genes_used    methods
supervised 0.125585308 0.0148850520 0.3092444      top10 supervised
trend2     0.113966791 0.0884220119 0.2200606      top10     trend2
bspline    0.101424270 0.1342596991 0.2563613      top10    bspline
loess      0.133801337 0.1079941543 0.2980371      top10      loess
seurat     0.007779095 0.0005450732 0.1126729      top10     seurat

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

get.aov(yy=results_eval_top10$fit.seurat$dapi,
        xx=results_eval_top10$fit.seurat$assignments)

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

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

Seurat scores vs cell time predictions

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

seurat.S.owntime <- with(results_eval_top10,
                     get.pve(fit.seurat$S[order(fit.seurat$cell_times_est)]))
seurat.G2M.owntime <- with(results_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.Rmd/seurat.time.top10.rda")

load(file="../output/method-eval-withheld-explore.Rmd/seurat.time.top10.rda")

c(seurat.S.sup$pve, seurat.G2M.sup$pve)

[1] 0.3896577 0.5287918

c(seurat.S.unsup$pve, seurat.G2M.unsup$pve)

[1] 0.2623815 0.5239003

c(seurat.S.owntime$pve, seurat.G2M.owntime$pve)

[1] 0.8272036 0.8565118

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

Seurat cell class vs unsupervsied predicted time

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

Seurat clases vs supervsied cell times

cols <- c("orange", "red", "brown")
par(mfrow=c(1,3))
with(results_eval_top10, 
     hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="G1"],
     nclass=5, col=cols[1], main = "G1, 10 cells", xlim=c(0,2*pi), ylim=c(0,30),
     xlab="Predicted cell time"))
with(results_eval_top10, hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="S"],
     nclass=10, col=cols[2],
     main = "S, 66 cells", xlim=c(0,2*pi), ylim=c(0,30),
     xlab="Predicted cell time"))
with(results_eval_top10, hist(fit.supervised$pred_time_shift[fit.seurat$assignments=="G2M"],
     nclass=10, col=cols[3],
     main = "G2M, 57 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)

Before 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_eval_top101$fit.supervised$diff_time)/2/pi
mean(results_eval_top101$fit.trend2.unsup$diff_time)/2/pi
mean(results_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_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_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.Rmd/pve_eval_top101.rds")

pve_eval_top101 <- readRDS("../output/method-eval-withheld-explore.Rmd/pve_eval_top101.rds")
print(pve_eval_top101)

                  dapi          gfp       rfp genes_used    methods
supervised 0.099929197 0.1208146060 0.2861203     top101 supervised
trend2     0.147145939 0.0948834158 0.1423536     top101     trend2
bspline    0.183418205 0.0020355877 0.3354409     top101    bspline
loess      0.004744912 0.0649387956 0.1642361     top101      loess
seurat     0.007779095 0.0005450732 0.1126729     top101     seurat

get.aov(yy=results_eval_top101$fit.seurat$dapi,
        xx=results_eval_top101$fit.seurat$assignments)

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

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

Seurat scores vs cell time predictions

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

seurat.S.owntime <- with(results_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.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_eval_top101,
#   get.aov(yy=fit.seurat$G2M,xx=fit.seurat$assignments))
# with(results_eval_top101,
#   get.aov(yy=fit.seurat$S,xx=fit.seurat$assignments))

Seurat cell class vs unsupervsied predicted time

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

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, 10 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, 66 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, 57 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

Evaluate the withheld sample on top X cyclical

Joyce Hsiao

Run on supervised model

Load previously computed results

Before removing moisy labels, on top 10 cyclical genes

Seurat scores vs cell time predictions

Seurat cell class vs unsupervsied predicted time

Seurat clases vs supervsied cell times

Before removing moisy labels, on top 101 cyclical genes

Seurat scores vs cell time predictions

Seurat cell class vs unsupervsied predicted time

Seurat clases vs supervsied cell times

Session information