# Predicting cell cycle phase using peco

#### Joyce Hsiao

#### 2022-04-26

## 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
    #>
    #> Attaching package: 'BiocGenerics'
    #> 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
    #> Loading required package: parallel
    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.481858     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.2.0 RC (2022-04-19 r82224 ucrt)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server x64 (build 20348)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.utf8
    #> [3] LC_MONETARY=English_United States.utf8
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.utf8
    #>
    #> attached base packages:
    #> [1] parallel  stats4    stats     graphics  grDevices utils     datasets
    #> [8] methods   base
    #>
    #> other attached packages:
    #>  [1] doParallel_1.0.17           iterators_1.0.14
    #>  [3] foreach_1.5.2               SingleCellExperiment_1.18.0
    #>  [5] SummarizedExperiment_1.26.0 Biobase_2.56.0
    #>  [7] GenomicRanges_1.48.0        GenomeInfoDb_1.32.0
    #>  [9] IRanges_2.30.0              S4Vectors_0.34.0
    #> [11] BiocGenerics_0.42.0         MatrixGenerics_1.8.0
    #> [13] matrixStats_0.62.0          peco_1.8.0
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] bitops_1.0-7              tools_4.2.0
    #>  [3] bslib_0.3.1               utf8_1.2.2
    #>  [5] R6_2.5.1                  irlba_2.3.5
    #>  [7] vipor_0.4.5               DBI_1.1.2
    #>  [9] colorspace_2.0-3          gridExtra_2.3
    #> [11] tidyselect_1.1.2          compiler_4.2.0
    #> [13] cli_3.3.0                 BiocNeighbors_1.14.0
    #> [15] DelayedArray_0.22.0       sass_0.4.1
    #> [17] scales_1.2.0              mvtnorm_1.1-3
    #> [19] stringr_1.4.0             digest_0.6.29
    #> [21] rmarkdown_2.14            geigen_2.3
    #> [23] XVector_0.36.0            scater_1.24.0
    #> [25] pkgconfig_2.0.3           htmltools_0.5.2
    #> [27] sparseMatrixStats_1.8.0   highr_0.9
    #> [29] fastmap_1.1.0             conicfit_1.0.4
    #> [31] rlang_1.0.2               circular_0.4-94.1
    #> [33] DelayedMatrixStats_1.18.0 jquerylib_0.1.4
    #> [35] generics_0.1.2            jsonlite_1.8.0
    #> [37] BiocParallel_1.30.0       dplyr_1.0.8
    #> [39] RCurl_1.98-1.6            magrittr_2.0.3
    #> [41] BiocSingular_1.12.0       GenomeInfoDbData_1.2.8
    #> [43] scuttle_1.6.0             Matrix_1.4-1
    #> [45] Rcpp_1.0.8.3              ggbeeswarm_0.6.0
    #> [47] munsell_0.5.0             fansi_1.0.3
    #> [49] viridis_0.6.2             lifecycle_1.0.1
    #> [51] stringi_1.7.6             yaml_2.3.5
    #> [53] zlibbioc_1.42.0           grid_4.2.0
    #> [55] ggrepel_0.9.1             crayon_1.5.1
    #> [57] lattice_0.20-45           beachmat_2.12.0
    #> [59] knitr_1.38                pillar_1.7.0
    #> [61] igraph_1.3.1              boot_1.3-28
    #> [63] codetools_0.2-18          ScaledMatrix_1.4.0
    #> [65] glue_1.6.2                evaluate_0.15
    #> [67] genlasso_1.5              vctrs_0.4.1
    #> [69] gtable_0.3.0              purrr_0.3.4
    #> [71] assertthat_0.2.1          ggplot2_3.3.5
    #> [73] xfun_0.30                 rsvd_1.0.5
    #> [75] pracma_2.3.8              viridisLite_0.4.0
    #> [77] tibble_3.1.6              beeswarm_0.4.0
    #> [79] ellipsis_0.3.2