# Predicting cell cycle phase using peco

#### Joyce Hsiao

#### 2021-05-19

## 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: MatrixGenerics
    #> Loading required package: matrixStats
    #>
    #> Attaching package: 'MatrixGenerics'
    #> The following objects are masked from 'package:matrixStats':
    #>
    #>     colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
    #>     colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
    #>     colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
    #>     colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
    #>     colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
    #>     colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
    #>     colWeightedMeans, colWeightedMedians, colWeightedSds,
    #>     colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
    #>     rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
    #>     rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
    #>     rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
    #>     rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
    #>     rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
    #>     rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
    #>     rowWeightedSds, rowWeightedVars
    #> 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.max, which.min
    #> Loading required package: S4Vectors
    #>
    #> Attaching package: 'S4Vectors'
    #> The following objects are masked from 'package:base':
    #>
    #>     I, expand.grid, unname
    #> 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")'.
    #>
    #> Attaching package: 'Biobase'
    #> The following object is masked from 'package:MatrixGenerics':
    #>
    #>     rowMedians
    #> The following objects are masked from 'package:matrixStats':
    #>
    #>     anyMissing, rowMedians
    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.1.0 RC (2021-05-10 r80283)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server x64 (build 17763)
    #>
    #> 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.16           iterators_1.0.13
    #>  [3] foreach_1.5.1               SingleCellExperiment_1.14.0
    #>  [5] SummarizedExperiment_1.22.0 Biobase_2.52.0
    #>  [7] GenomicRanges_1.44.0        GenomeInfoDb_1.28.0
    #>  [9] IRanges_2.26.0              S4Vectors_0.30.0
    #> [11] BiocGenerics_0.38.0         MatrixGenerics_1.4.0
    #> [13] matrixStats_0.58.0          peco_1.4.0
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] bitops_1.0-7              tools_4.1.0
    #>  [3] bslib_0.2.5.1             utf8_1.2.1
    #>  [5] R6_2.5.0                  irlba_2.3.3
    #>  [7] vipor_0.4.5               DBI_1.1.1
    #>  [9] colorspace_2.0-1          tidyselect_1.1.1
    #> [11] gridExtra_2.3             compiler_4.1.0
    #> [13] BiocNeighbors_1.10.0      DelayedArray_0.18.0
    #> [15] sass_0.4.0                scales_1.1.1
    #> [17] mvtnorm_1.1-1             stringr_1.4.0
    #> [19] digest_0.6.27             rmarkdown_2.8
    #> [21] geigen_2.3                XVector_0.32.0
    #> [23] scater_1.20.0             pkgconfig_2.0.3
    #> [25] htmltools_0.5.1.1         sparseMatrixStats_1.4.0
    #> [27] highr_0.9                 conicfit_1.0.4
    #> [29] rlang_0.4.11              circular_0.4-93
    #> [31] DelayedMatrixStats_1.14.0 jquerylib_0.1.4
    #> [33] generics_0.1.0            jsonlite_1.7.2
    #> [35] BiocParallel_1.26.0       dplyr_1.0.6
    #> [37] RCurl_1.98-1.3            magrittr_2.0.1
    #> [39] BiocSingular_1.8.0        GenomeInfoDbData_1.2.6
    #> [41] scuttle_1.2.0             Matrix_1.3-3
    #> [43] Rcpp_1.0.6                ggbeeswarm_0.6.0
    #> [45] munsell_0.5.0             fansi_0.4.2
    #> [47] viridis_0.6.1             lifecycle_1.0.0
    #> [49] stringi_1.6.2             yaml_2.2.1
    #> [51] zlibbioc_1.38.0           grid_4.1.0
    #> [53] crayon_1.4.1              lattice_0.20-44
    #> [55] beachmat_2.8.0            knitr_1.33
    #> [57] pillar_1.6.1              igraph_1.2.6
    #> [59] boot_1.3-28               codetools_0.2-18
    #> [61] ScaledMatrix_1.0.0        glue_1.4.2
    #> [63] evaluate_0.14             genlasso_1.5
    #> [65] vctrs_0.3.8               gtable_0.3.0
    #> [67] purrr_0.3.4               assertthat_0.2.1
    #> [69] ggplot2_3.3.3             xfun_0.23
    #> [71] rsvd_1.0.5                pracma_2.3.3
    #> [73] viridisLite_0.4.0         tibble_3.1.2
    #> [75] beeswarm_0.3.1            ellipsis_0.3.2