For each gene, test if it is more variable than expected under the null
Arguments
- x
model fit from
lucida()- ...
other args
Details
Under the null model of no expression variation, the deviance residuals from the NB model fit have a chisq distribution. Hypothesis testing is performed for each gene comparing the observed variance to the theoretical null
Examples
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: generics
#>
#> Attaching package: ‘generics’
#> The following objects are masked from ‘package:base’:
#>
#> as.difftime, as.factor, as.ordered, intersect, is.element, setdiff,
#> setequal, union
#>
#> 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, aperm, append,
#> as.data.frame, basename, cbind, colnames, dirname, do.call,
#> duplicated, eval, evalq, get, grep, grepl, is.unsorted, lapply,
#> mapply, match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
#> rank, rbind, rownames, sapply, saveRDS, table, tapply, unique,
#> unsplit, which.max, which.min
#> Loading required package: S4Vectors
#>
#> Attaching package: ‘S4Vectors’
#> The following object is masked from ‘package:utils’:
#>
#> findMatches
#> The following objects are masked from ‘package:base’:
#>
#> I, expand.grid, unname
#> Loading required package: IRanges
#> Loading required package: Seqinfo
#> 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
# Load example data
data(example_sce, package="muscat")
sce <- example_sce
# Compute library size for each cell
sce$libSize <- colSums(counts(sce))
# run NB model on all cells
fit <- lucida(sce, ~ group_id, verbose=FALSE)
#> all
# Find highly variable genes
findHVGs( fit )
#> # A tibble: 1,205 × 5
#> ID var P.Value FDR isHVG
#> <chr> <dbl> <dbl> <dbl> <lgl>
#> 1 HES4 1.08 1.19 e- 2 2.60e- 2 TRUE
#> 2 ISG15 1.59 5.45 e-45 2.63e-43 TRUE
#> 3 AURKAIP1 1.06 4.60 e- 2 8.99e- 2 FALSE
#> 4 MRPL20 1.19 2.04 e- 7 1.02e- 6 TRUE
#> 5 SSU72 1.10 2.72 e- 3 6.64e- 3 TRUE
#> 6 RER1 0.967 8.21 e- 1 1 e+ 0 FALSE
#> 7 RPL22 1.23 6.37 e-10 4.44e- 9 TRUE
#> 8 PARK7 1.19 2.92 e- 7 1.42e- 6 TRUE
#> 9 ENO1 1.15 2.07 e- 5 7.35e- 5 TRUE
#> 10 FBXO6 0.847 1.000e+ 0 1 e+ 0 FALSE
#> # ℹ 1,195 more rows