Merge pull request #106 from jr-leary7/dev

Dev

R/testDynamic.R

@@ -115,6 +115,9 @@ testDynamic <- function(expr.mat = NULL,
   if (inherits(pt, "SlingshotDataSet")) {
     pt <- as.data.frame(slingshot::slingPseudotime(pt))
   }
+  if (inherits(pt, "tbl") || inherits(pt, "tbl_df")) {
+    pt <- as.data.frame(pt)
+  }
   if (!inherits(pt, "data.frame")) { stop("pt must be of class data.frame.") }
   # set pseudotime lineage column names automatically to prevent user error (uses e.g., "Lineage_A", "Lineage_B", etc.)
   n_lineages <- ncol(pt)
@@ -180,7 +183,7 @@ testDynamic <- function(expr.mat = NULL,
     if (!is.null(log.file) && i %in% print_nums) {
       withr::with_output_sink(log.file,
                               code = {
-                                cat(paste(Sys.time(), "- Starting iteration:", 20, "\n"))
+                                cat(paste(Sys.time(), "- Starting iteration:", i, "\n"))
                               },
                               append = TRUE)
     }

README.Rmd

@@ -1,6 +1,11 @@
 ---
 output: 
-  github_document
+  github_document:
+    toc: true
+    toc_depth: 2
+    fig_width: 9
+    fig_height: 6
+    dev: png
 always_allow_html: true
 ---
 
@@ -23,7 +28,6 @@ knitr::opts_chunk$set(message = FALSE,
 [![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)
 [![License: MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
 ![release](https://img.shields.io/github/v/release/jr-leary7/scLANE?color=purple)
-![downloads](https://img.shields.io/github/downloads/jr-leary7/scLANE/total?color=blue)
 [![CodeFactor](https://www.codefactor.io/repository/github/jr-leary7/sclane/badge)](https://www.codefactor.io/repository/github/jr-leary7/sclane)
 ![last commit](https://img.shields.io/github/last-commit/jr-leary7/scLANE/main?color=darkgreen)
 <!-- badges: end -->
@@ -89,7 +93,7 @@ plotPCA(sim_data, colour_by = "subject")
 plotPCA(sim_data, colour_by = "cell_time_normed")
 ```
 
-## Differential Expression Testing
+## Trajectory DE Testing
 
 Since we have multi-subject data, we can use any of the three model backends to run our DE testing. We'll start with the simplest model, the GLM, then work our way through the other options in order of increasing complexity. We first prepare our inputs - a dense matrix with cells as rows & genes as columns (i.e., transposed from the way it's stored in `SingleCellExperiment` & `Seurat` objects), a dataframe containing our cell ordering, a set of genes to build models for, and a vector of per-cell size factors to be used as offsets during estimation. In reality, it's usually unnecessary to fit a model for every single gene in a dataset, as trajectories are usually estimated using a subset of the entire set of genes (usually a few thousand most highly variable genes). For the purpose of demonstration, we'll select 50 genes each from the dynamic and non-dynamic populations. **Note**: in this case we're working with a single pseudotime lineage, though in real datasets several lineages often exist; in order to fit models for a subset of lineages simply remove the corresponding columns from the cell ordering dataframe passed as input to `testDynamic()`. 
 

README.md

@@ -1,4 +1,21 @@
 
+- <a href="#sclane" id="toc-sclane"><code>scLANE</code></a>
+  - <a href="#installation" id="toc-installation">Installation</a>
+  - <a href="#model-structure" id="toc-model-structure">Model Structure</a>
+- <a href="#usage" id="toc-usage">Usage</a>
+  - <a href="#input-data" id="toc-input-data">Input Data</a>
+  - <a href="#trajectory-de-testing"
+    id="toc-trajectory-de-testing">Trajectory DE Testing</a>
+  - <a href="#downstream-analysis--visualization"
+    id="toc-downstream-analysis--visualization">Downstream Analysis &amp;
+    Visualization</a>
+- <a href="#conclusions--best-practices"
+  id="toc-conclusions--best-practices">Conclusions &amp; Best
+  Practices</a>
+- <a href="#contact-information" id="toc-contact-information">Contact
+  Information</a>
+- <a href="#references" id="toc-references">References</a>
+
 <!-- README.md is generated from README.Rmd. Please edit that file -->
 
 # `scLANE`
@@ -10,7 +27,6 @@
 [![License:
 MIT](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT)
 ![release](https://img.shields.io/github/v/release/jr-leary7/scLANE?color=purple)
-![downloads](https://img.shields.io/github/downloads/jr-leary7/scLANE/total?color=blue)
 [![CodeFactor](https://www.codefactor.io/repository/github/jr-leary7/sclane/badge)](https://www.codefactor.io/repository/github/jr-leary7/sclane)
 ![last
 commit](https://img.shields.io/github/last-commit/jr-leary7/scLANE/main?color=darkgreen)
@@ -139,7 +155,7 @@ plotPCA(sim_data, colour_by = "cell_time_normed")
 
 <img src="man/figures/README-plot-sims-2.png" width="100%" />
 
-## Differential Expression Testing
+## Trajectory DE Testing
 
 Since we have multi-subject data, we can use any of the three model
 backends to run our DE testing. We’ll start with the simplest model, the
@@ -184,7 +200,7 @@ de_test_glm <- testDynamic(expr.mat = counts_mat,
                            n.potential.basis.fns = 5, 
                            n.cores = 2, 
                            track.time = TRUE)
-#> [1] "testDynamic evaluated 100 genes with 1 lineage apiece in 26.322 secs"
+#> [1] "testDynamic evaluated 100 genes with 1 lineage apiece in 26.588 secs"
 ```
 
 After the function finishes running, we use `getResultsDE()` to generate
@@ -280,7 +296,7 @@ de_test_gee <- testDynamic(expr.mat = counts_mat,
                            cor.structure = "ar1", 
                            n.cores = 2, 
                            track.time = TRUE)
-#> [1] "testDynamic evaluated 100 genes with 1 lineage apiece in 3.105 mins"
+#> [1] "testDynamic evaluated 100 genes with 1 lineage apiece in 3.849 mins"
 ```
 
 We again generate the table of DE test results. The variance of the
@@ -374,7 +390,7 @@ de_test_glmm <- testDynamic(expr.mat = counts_mat,
                             id.vec = sim_data$subject, 
                             n.cores = 2, 
                             track.time = TRUE)
-#> [1] "testDynamic evaluated 100 genes with 1 lineage apiece in 3.343 mins"
+#> [1] "testDynamic evaluated 100 genes with 1 lineage apiece in 5.306 mins"
 ```
 
 Like the GLM backend, the GLMM backend use a likelihood ratio test to
@@ -518,13 +534,13 @@ slice_sample(gene_clust_table, n = 5) %>%
                col.names = c("Gene", "Lineage", "Cell", "Fitted (link)", "Fitted (response)", "Pseudotime", "Cluster"))
 ```
 
-| Gene   | Lineage | Cell | Fitted (link) | Fitted (response) | Pseudotime | Cluster |
-|:-------|:--------|:-----|--------------:|------------------:|-----------:|:--------|
-| SFMBT2 | A       | 63   |        -2.677 |             0.043 |      0.158 | 1       |
-| UBIAD1 | A       | 945  |        -2.227 |             0.188 |      0.362 | 1       |
-| TSPAN1 | A       | 106  |         0.279 |             0.807 |      0.265 | 1       |
-| CPA3   | A       | 79   |        -4.372 |             0.120 |      0.198 | 1       |
-| WAPAL  | A       | 1150 |         0.755 |             1.023 |      0.875 | 1       |
+| Gene    | Lineage | Cell | Fitted (link) | Fitted (response) | Pseudotime | Cluster |
+|:--------|:--------|:-----|--------------:|------------------:|-----------:|:--------|
+| CKAP4   | A       | 1130 |         0.933 |             1.458 |      0.825 | 2       |
+| CBX6    | A       | 607  |         1.396 |             2.407 |      0.517 | 2       |
+| NLGN4Y  | A       | 390  |        -2.941 |             0.056 |      0.975 | 1       |
+| MPG     | A       | 1150 |        -0.147 |             0.415 |      0.875 | 1       |
+| TRAPPC1 | A       | 91   |         0.538 |             1.683 |      0.228 | 2       |
 
 The results can then be plotted as desired using `ggplot2` or another
 visualization package. Upon visual inspection, the genes seem to cluster