fit_trendfilter {peco} | R Documentation

## Using trendfiltering to estimate cyclic trend of gene expression

### Description

We applied quadratic (second order) trend filtering using the trendfilter function in the genlasso package (Tibshirani, 2014). The trendfilter function implements a nonparametric smoothing method which chooses the smoothing parameter by cross-validation and fits a piecewise polynomial regression. In more specifics: The trendfilter method determines the folds in cross-validation in a nonrandom manner. Every k-th data point in the ordered sample is placed in the k-th fold, so the folds contain ordered subsamples. We applied five-fold cross-validation and chose the smoothing penalty using the option lambda.1se: among all possible values of the penalty term, the largest value such that the cross-validation standard error is within one standard error of the minimum. Furthermore, we desired that the estimated expression trend be cyclical. To encourage this, we concatenated the ordered gene expression data three times, with one added after another. The quadratic trend filtering was applied to the concatenated data series of each gene.

### Usage

    fit_trendfilter(yy, polyorder = 2)

### Arguments

`yy` | A vector of gene expression values for one gene that are ordered by cell cycle phase. Also, the expression values are normalized and transformed to standard normal distribution.
`polyorder` | We estimate cyclic trends of gene expression levels using nonparamtric trend filtering. The default fits second degree polynomials.

### Value

A list with two elements:

`trend.yy` | The estimated cyclic trend.
`pve` | Proportion of variance explained by the cyclic trend in the gene expression levels.

### Author(s)

Joyce Hsiao

### Examples

    library(SingleCellExperiment)
    data(sce_top101genes)

    # select top 10 cyclic genes
    sce_top10 <- sce_top101genes[order(rowData(sce_top101genes)$pve_fucci,
                                      decreasing=TRUE)[1:10],]
    coldata <- colData(sce_top10)

    # cell cycle phase based on FUCCI scores
    theta <- coldata$theta
    names(theta) <- rownames(coldata)

    # normalize expression counts
    sce_top10 <- data_transform_quantile(sce_top10, ncores=2)
    exprs_quant <- assay(sce_top10, "cpm_quantNormed")

    # order FUCCI phase and expression
    theta_ordered <- theta[order(theta)]
    yy_ordered <- exprs_quant[1, names(theta_ordered)]

    fit <- fit_trendfilter(yy_ordered)

    plot(x=theta_ordered, y=yy_ordered, pch=16, cex=.7, axes=FALSE,
      ylab="quantile-normalized expression values", xlab="FUCCI phase",
      main = "trendfilter fit")
    points(x=theta_ordered, y=fit$trend.yy, col="blue", pch=16, cex=.7)
    axis(2)
    axis(1,at=c(0,pi/2, pi, 3*pi/2, 2*pi),
      labels=c(0,expression(pi/2), expression(pi), expression(3*pi/2),
      expression(2*pi)))
    abline(h=0, lty=1, col="black", lwd=.7)

[Package peco version 1.2.0 Index]