diff --git a/README.md b/README.md index 587b7fc..877be80 100644 --- a/README.md +++ b/README.md @@ -24,6 +24,7 @@ [![R-CMD-check](https://github.com/jr-leary7/scLANE/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/jr-leary7/scLANE/actions/workflows/R-CMD-check.yaml) +[![Bioc-check](https://github.com/jr-leary7/scLANE/actions/workflows/bioc-check.yaml/badge.svg)](https://github.com/jr-leary7/scLANE/actions/workflows/bioc-check.yaml) ![release](https://img.shields.io/github/v/release/jr-leary7/scLANE?color=purple) ![last commit](https://img.shields.io/github/last-commit/jr-leary7/scLANE/main?color=darkgreen) @@ -94,11 +95,8 @@ conflict. ``` r library(dplyr) -library(scater) library(scLANE) library(ggplot2) -select <- dplyr::select -filter <- dplyr::filter ``` ## Input data @@ -117,8 +115,101 @@ The PCA embeddings show us a pretty simple trajectory that’s strongly correlated with the first principal component. ``` r -plotPCA(sim_data, colour_by = "cell_time_normed") + +data.frame(sim_data@int_colData$reducedDims@listData$PCA[, 1:2]) %>% + mutate(pseudotime = sim_data$cell_time_normed) %>% + ggplot(aes(x = PC1, y = PC2, color = pseudotime)) + + geom_point(size = 2, alpha = 0.75, stroke = 0) + + scale_color_gradientn(colors = viridisLite::plasma(n = 20)) + + labs(x = "PC 1", y = "PC 2", color = "Pseudotime") + theme_scLANE(umap = TRUE) +#> Loading required package: S4Vectors +#> Loading required package: stats4 +#> Loading required package: BiocGenerics +#> +#> Attaching package: 'BiocGenerics' +#> The following objects are masked from 'package:dplyr': +#> +#> combine, intersect, setdiff, union +#> The following objects are masked from 'package:stats': +#> +#> IQR, mad, sd, var, xtabs +#> The following objects are masked from 'package:base': +#> +#> anyDuplicated, aperm, append, as.data.frame, basename, cbind, +#> colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find, +#> get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply, +#> match, mget, order, paste, pmax, pmax.int, pmin, pmin.int, +#> Position, rank, rbind, Reduce, rownames, sapply, setdiff, sort, +#> table, tapply, union, unique, unsplit, which.max, which.min +#> +#> Attaching package: 'S4Vectors' +#> The following objects are masked from 'package:dplyr': +#> +#> first, rename +#> The following object is masked from 'package:utils': +#> +#> findMatches +#> The following objects are masked from 'package:base': +#> +#> expand.grid, I, unname +#> Loading required package: SingleCellExperiment +#> Loading required package: SummarizedExperiment +#> Loading required package: MatrixGenerics +#> Loading required package: matrixStats +#> +#> Attaching package: 'matrixStats' +#> The following object is masked from 'package:dplyr': +#> +#> count +#> +#> 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: IRanges +#> +#> Attaching package: 'IRanges' +#> The following objects are masked from 'package:dplyr': +#> +#> collapse, desc, slice +#> Loading required package: GenomeInfoDb +#> +#> Attaching package: 'GenomicRanges' +#> The following object is masked from 'package:magrittr': +#> +#> subtract +#> 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 +#> Registered S3 method overwritten by 'SparseArray': +#> method from +#> rowsum.dgCMatrix DelayedArray ``` @@ -127,7 +218,11 @@ We also see that the data are not clustered by subject, which indicates that gene dynamics are mostly homogeneous across subjects. ``` r -plotPCA(sim_data, colour_by = "subject") + +data.frame(sim_data@int_colData$reducedDims@listData$PCA[, 1:2]) %>% + mutate(subject = sim_data$subject) %>% + ggplot(aes(x = PC1, y = PC2, color = subject)) + + geom_point(size = 2, alpha = 0.75, stroke = 0) + + labs(x = "PC 1", y = "PC 2", color = "Subject ID") + theme_scLANE(umap = TRUE) ``` @@ -181,7 +276,7 @@ scLANE_models_glm <- testDynamic(sim_data, #> Registered S3 method overwritten by 'bit': #> method from #> print.ri gamlss -#> scLANE testing completed for 100 genes across 1 lineage in 35.114 secs +#> scLANE testing completed for 100 genes across 1 lineage in 29.977 secs ``` After the function finishes running, we use `getResultsDE()` to generate @@ -200,13 +295,13 @@ select(scLANE_res_glm, Gene, Lineage, Test_Stat, P_Val, P_Val_Adj, Gene_Dynamic_ col.names = c("Gene", "Lineage", "LRT stat.", "P-value", "Adj. p-value", "Predicted dynamic status")) ``` -| Gene | Lineage | LRT stat. | P-value | Adj. p-value | Predicted dynamic status | -|:---------|:--------|----------:|--------:|-------------:|-------------------------:| -| MFSD2B | A | 209.755 | 0.000 | 0.000 | 1 | -| GOLGA8EP | A | 4.390 | 0.036 | 0.537 | 0 | -| NOP14 | A | 7.712 | 0.005 | 0.132 | 0 | -| TMCO3 | A | 167.582 | 0.000 | 0.000 | 1 | -| BATF2 | A | 5.216 | 0.074 | 0.537 | 0 | +| Gene | Lineage | LRT stat. | P-value | Adj. p-value | Predicted dynamic status | +|:-------|:--------|----------:|--------:|-------------:|-------------------------:| +| MFSD2B | A | 209.755 | 0.000 | 0.000 | 1 | +| SMG1 | A | 5.062 | 0.080 | 0.419 | 0 | +| UAP1L1 | A | 9.118 | 0.010 | 0.199 | 0 | +| TMCO3 | A | 167.582 | 0.000 | 0.000 | 1 | +| MYOF | A | 4.015 | 0.045 | 0.419 | 0 | ### GEE mode @@ -229,7 +324,7 @@ scLANE_models_gee <- testDynamic(sim_data, id.vec = sim_data$subject, cor.structure = "ar1", n.cores = 4) -#> scLANE testing completed for 100 genes across 1 lineage in 2.201 mins +#> scLANE testing completed for 100 genes across 1 lineage in 3.076 mins ``` We again generate the table of DE test results. The variance of the @@ -248,11 +343,11 @@ select(scLANE_res_gee, Gene, Lineage, Test_Stat, P_Val, P_Val_Adj, Gene_Dynamic_ | Gene | Lineage | Wald stat. | P-value | Adj. p-value | Predicted dynamic status | |:---------|:--------|-----------:|--------:|-------------:|-------------------------:| -| DGUOK | A | 200675.460 | 0.000 | 0.000 | 1 | -| VDAC1 | A | 13.025 | 0.001 | 0.025 | 0 | +| CKAP4 | A | 48922.644 | 0 | 0 | 1 | +| TBCC | A | 32.151 | 0 | 0 | 1 | | GOLGA8EP | A | NA | NA | NA | 0 | -| WAPAL | A | 3254.351 | 0.000 | 0.000 | 1 | -| IDH3G | A | 863.479 | 0.000 | 0.000 | 1 | +| PFDN2 | A | 2131.763 | 0 | 0 | 1 | +| IDH3G | A | 863.479 | 0 | 0 | 1 | ### GLMM mode @@ -276,7 +371,7 @@ scLANE_models_glmm <- testDynamic(sim_data, glmm.adaptive = TRUE, id.vec = sim_data$subject, n.cores = 4) -#> scLANE testing completed for 100 genes across 1 lineage in 3.427 mins +#> scLANE testing completed for 100 genes across 1 lineage in 4.863 mins ``` {% note %} @@ -306,8 +401,8 @@ select(scLANE_res_glmm, Gene, Lineage, Test_Stat, P_Val, P_Val_Adj, Gene_Dynamic | Gene | Lineage | LRT stat. | P-value | Adj. p-value | Predicted dynamic status | |:--------|:--------|----------:|--------:|-------------:|-------------------------:| -| DDX1 | A | 131.707 | 0.000 | 0 | 1 | -| TSPAN1 | A | 78.616 | 0.000 | 0 | 1 | +| PGLS | A | 129.086 | 0.000 | 0 | 1 | +| TSPAN1 | A | 85.987 | 0.000 | 0 | 1 | | WDSUB1 | A | NA | NA | NA | 0 | | FAM135B | A | NA | NA | NA | 0 | | NLGN4Y | A | 9.878 | 0.627 | 1 | 0 | diff --git a/man/figures/README-plot-models-gee-1.png b/man/figures/README-plot-models-gee-1.png index c1d170c..104d0b7 100644 Binary files a/man/figures/README-plot-models-gee-1.png and b/man/figures/README-plot-models-gee-1.png differ diff --git a/man/figures/README-plot-models-glm-1.png b/man/figures/README-plot-models-glm-1.png index 5554ee8..55df417 100644 Binary files a/man/figures/README-plot-models-glm-1.png and b/man/figures/README-plot-models-glm-1.png differ diff --git a/man/figures/README-plot-models-glmm-1.png b/man/figures/README-plot-models-glmm-1.png index 0bc197b..57ac989 100644 Binary files a/man/figures/README-plot-models-glmm-1.png and b/man/figures/README-plot-models-glmm-1.png differ diff --git a/man/figures/README-plot-sims-pt-1.png b/man/figures/README-plot-sims-pt-1.png index 0734055..bdc2fbc 100644 Binary files a/man/figures/README-plot-sims-pt-1.png and b/man/figures/README-plot-sims-pt-1.png differ diff --git a/man/figures/README-plot-sims-subj-1.png b/man/figures/README-plot-sims-subj-1.png index 07248a3..2b8641f 100644 Binary files a/man/figures/README-plot-sims-subj-1.png and b/man/figures/README-plot-sims-subj-1.png differ diff --git a/man/figures/README-unnamed-chunk-2-1.png b/man/figures/README-unnamed-chunk-2-1.png index 3666dea..8eeb103 100644 Binary files a/man/figures/README-unnamed-chunk-2-1.png and b/man/figures/README-unnamed-chunk-2-1.png differ