# Predicting cell cycle phase using peco

#### Joyce Hsiao

#### 2020-04-17

## Installation

To install and load the package, run:

    install.packages("devtools")
    library(devtools)
    install_github("jhsiao999/peco")

`peco` uses `SingleCellExperiment` class objects.

    library(peco)
    library(SingleCellExperiment)
    #> Loading required package: SummarizedExperiment
    #> Loading required package: GenomicRanges
    #> Loading required package: stats4
    #> Loading required package: BiocGenerics
    #> Loading required package: parallel
    #>
    #> Attaching package: 'BiocGenerics'
    #> The following objects are masked from 'package:parallel':
    #>
    #>     clusterApply, clusterApplyLB, clusterCall, clusterEvalQ,
    #>     clusterExport, clusterMap, parApply, parCapply, parLapply,
    #>     parLapplyLB, parRapply, parSapply, parSapplyLB
    #> The following objects are masked from 'package:stats':
    #>
    #>     IQR, mad, sd, var, xtabs
    #> The following objects are masked from 'package:base':
    #>
    #>     Filter, Find, Map, Position, Reduce, anyDuplicated, append,
    #>     as.data.frame, basename, cbind, colnames, dirname, do.call,
    #>     duplicated, eval, evalq, get, grep, grepl, intersect, is.unsorted,
    #>     lapply, mapply, match, mget, order, paste, pmax, pmax.int, pmin,
    #>     pmin.int, rank, rbind, rownames, sapply, setdiff, sort, table,
    #>     tapply, union, unique, unsplit, which, which.max, which.min
    #> Loading required package: S4Vectors
    #>
    #> Attaching package: 'S4Vectors'
    #> The following object is masked from 'package:base':
    #>
    #>     expand.grid
    #> Loading required package: IRanges
    #>
    #> Attaching package: 'IRanges'
    #> The following object is masked from 'package:grDevices':
    #>
    #>     windows
    #> Loading required package: GenomeInfoDb
    #> Loading required package: Biobase
    #> Welcome to Bioconductor
    #>
    #>     Vignettes contain introductory material; view with
    #>     'browseVignettes()'. To cite Bioconductor, see
    #>     'citation("Biobase")', and for packages 'citation("pkgname")'.
    #> Loading required package: DelayedArray
    #> Loading required package: matrixStats
    #>
    #> Attaching package: 'matrixStats'
    #> The following objects are masked from 'package:Biobase':
    #>
    #>     anyMissing, rowMedians
    #>
    #> Attaching package: 'DelayedArray'
    #> The following objects are masked from 'package:matrixStats':
    #>
    #>     colMaxs, colMins, colRanges, rowMaxs, rowMins, rowRanges
    #> The following objects are masked from 'package:base':
    #>
    #>     aperm, apply, rowsum
    library(doParallel)
    #> Loading required package: foreach
    #> Loading required package: iterators
    library(foreach)

## Overview

`peco` is a supervised approach for PrEdicting cell cycle phase in a COntinuum using single-cell RNA sequencing data. The R package provides functions to build training dataset and also functions to use existing training data to predict cell cycle on a continuum.

Our work demonstrated that peco is able to predict continuous cell cylce phase using a small set of cylcic genes: CDK1, UBE2C, TOP2A, HISTH1E, and HISTH1C (identified as cell cycle marker genes in studies of yeast (Spellman et al., 1998) and HeLa cells (Whitfield et al., 2002)).

Below we provide two use cases. Vignette 1 shows how to use the built-training dataset to predict continuous cell cycle. Vignette 2 shows how to make a training datast and build a predictor using training data.

Users can also view the vigenettes via `browseVignettes("peco")`.

## About the training dataset

`training_human` stores built-in training data of 101 significant cyclic genes. Below are the slots contained in `training_human`:

- `predict.yy`: a gene by sample matrix (101 by 888) that stores predict cyclic expression values.
- `cellcycle_peco_reordered`: cell cycle phase in a unit circle (angle), ordered from 0 to 2\(pi\)
- `cellcycle_function`: lists of 101 function corresponding to the top 101 cyclic genes identified in our dataset
- `sigma`: standard error associated with cyclic trends of gene expression
- `pve`: proportion of variance explained by the cyclic trend

    data("training_human")

## Predict cell cycle phase using gene expression data

`peco` is integrated with `SingleCellExperiment` object in Bioconductor. Below shows an example of inputting `SingleCellExperiment` object to perform cell cycle phase prediction.

`sce_top101genes` includes 101 genes and 888 single-cell samples and one assay slot of `counts`.

    data("sce_top101genes")
    assays(sce_top101genes)
    #> List of length 1
    #> names(1): counts

Transform the expression values to quantile-normalizesd counts-per-million values. `peco` uses the `cpm_quantNormed` slot as input data for predictions.

    sce_top101genes <- data_transform_quantile(sce_top101genes)
    #> computing on 2 cores
    assays(sce_top101genes)
    #> List of length 3
    #> names(3): counts cpm cpm_quantNormed

Apply the prediction model using function `cycle_npreg_outsample` and generate prediction results contained in a list object `pred_top101genes`.

    pred_top101genes <- cycle_npreg_outsample(
        Y_test=sce_top101genes,
        sigma_est=training_human$sigma[rownames(sce_top101genes),],
        funs_est=training_human$cellcycle_function[rownames(sce_top101genes)],
        method.trend="trendfilter",
        ncores=1,
        get_trend_estimates=FALSE)

The `pred_top101genes$Y` contains a SingleCellExperiment object with the predict cell cycle phase in the `colData` slot.

    head(colData(pred_top101genes$Y)$cellcycle_peco)
    #> 20170905-A01 20170905-A02 20170905-A03 20170905-A06 20170905-A07 20170905-A08
    #>     1.099557     4.680973     2.607522     4.303982     4.052655     1.413717

Visualize results of prediction for one gene. Below we choose CDK1 (“ENSG00000170312”). Because CDK1 is a known cell cycle gene, this visualization serves as a sanity check for the results of fitting. The fitted function `training_human$cellcycle_function[[1]]` was obtained from our training data.

    plot(y=assay(pred_top101genes$Y,"cpm_quantNormed")["ENSG00000170312",],
         x=colData(pred_top101genes$Y)$theta_shifted, main = "CDK1",
         ylab = "quantile normalized expression")
    points(y=training_human$cellcycle_function[["ENSG00000170312"]](seq(0,2*pi, length.out=100)),
           x=seq(0,2*pi, length.out=100), col = "blue", pch =16)

## Visualize cyclic expression trend based on predicted phase

Visualize results of prediction for the top 10 genesone genes. Use `fit_cyclical_many` to estimate cyclic function based on the input data.

    # predicted cell time in the input data
    theta_predict = colData(pred_top101genes$Y)$cellcycle_peco
    names(theta_predict) = rownames(colData(pred_top101genes$Y))

    # expression values of 10 genes in the input data
    yy_input = assay(pred_top101genes$Y,"cpm_quantNormed")[1:6,]

    # apply trendfilter to estimate cyclic gene expression trend
    fit_cyclic <- fit_cyclical_many(Y=yy_input,
                                    theta=theta_predict)
    #> computing on 2 cores

    gene_symbols = rowData(pred_top101genes$Y)$hgnc[rownames(yy_input)]

    par(mfrow=c(2,3))
    for (i in 1:6) {
    plot(y=yy_input[i,],
         x=fit_cyclic$cellcycle_peco_ordered,
         main = gene_symbols[i],
         ylab = "quantile normalized expression")
    points(y=fit_cyclic$cellcycle_function[[i]](seq(0,2*pi, length.out=100)),
           x=seq(0,2*pi, length.out=100), col = "blue", pch =16)
    }

## Session information

    sessionInfo()
    #> R version 4.0.0 beta (2020-04-14 r78227)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server 2012 R2 x64 (build 9600)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.1252
    #> [3] LC_MONETARY=English_United States.1252
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.1252
    #>
    #> attached base packages:
    #> [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    #> [8] methods   base
    #>
    #> other attached packages:
    #>  [1] doParallel_1.0.15           iterators_1.0.12
    #>  [3] foreach_1.5.0               SingleCellExperiment_1.9.3
    #>  [5] SummarizedExperiment_1.17.5 DelayedArray_0.13.12
    #>  [7] matrixStats_0.56.0          Biobase_2.47.3
    #>  [9] GenomicRanges_1.39.3        GenomeInfoDb_1.23.17
    #> [11] IRanges_2.21.8              S4Vectors_0.25.15
    #> [13] BiocGenerics_0.33.3         peco_0.99.16
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] viridis_0.5.1            BiocSingular_1.3.2       genlasso_1.4
    #>  [4] viridisLite_0.3.0        DelayedMatrixStats_1.9.1 assertthat_0.2.1
    #>  [7] GenomeInfoDbData_1.2.3   vipor_0.4.5              yaml_2.2.1
    #> [10] pillar_1.4.3             lattice_0.20-41          glue_1.4.0
    #> [13] digest_0.6.25            XVector_0.27.2           colorspace_1.4-1
    #> [16] htmltools_0.4.0          Matrix_1.2-18            pkgconfig_2.0.3
    #> [19] zlibbioc_1.33.1          purrr_0.3.3              mvtnorm_1.1-0
    #> [22] scales_1.1.0             BiocParallel_1.21.2      pracma_2.2.9
    #> [25] tibble_3.0.0             ggplot2_3.3.0            conicfit_1.0.4
    #> [28] ellipsis_0.3.0           cli_2.0.2                magrittr_1.5
    #> [31] crayon_1.3.4             evaluate_0.14            fansi_0.4.1
    #> [34] MASS_7.3-51.5            beeswarm_0.2.3           geigen_2.3
    #> [37] tools_4.0.0              scater_1.15.30           lifecycle_0.2.0
    #> [40] stringr_1.4.0            munsell_0.5.0            irlba_2.3.3
    #> [43] compiler_4.0.0           rsvd_1.0.3               rlang_0.4.5
    #> [46] grid_4.0.0               RCurl_1.98-1.1           BiocNeighbors_1.5.4
    #> [49] circular_0.4-93          igraph_1.2.5             bitops_1.0-6
    #> [52] rmarkdown_2.1            boot_1.3-24              gtable_0.3.0
    #> [55] codetools_0.2-16         R6_2.4.1                 gridExtra_2.3
    #> [58] knitr_1.28               dplyr_0.8.5              stringi_1.4.6
    #> [61] ggbeeswarm_0.6.0         Rcpp_1.0.4.6             vctrs_0.2.4
    #> [64] tidyselect_1.0.0         xfun_0.13