Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Association tests yield weird results in single lineage case #17

Open
SONGDONGYUAN1994 opened this issue Oct 7, 2019 · 4 comments
Open
Labels
bug Something isn't working

Comments

@SONGDONGYUAN1994
Copy link

Hi,
Thanks for this nice package! I am trying to apply tradeSeq to single lineage simulation data, but the results seem to be incorrect.
Here is my code:

suppressMessages({library(devtools)
library(tidyverse, quietly = TRUE)
library(TSCAN)
library(plyr)
library(slingshot, quietly = TRUE)
library(RColorBrewer)
library(SingleCellExperiment)
library(scales)
library(gridExtra)
library(splatter)
library(dyntoy)
library(dynwrap)
library(mgcv)
library(mclust, quietly = TRUE)
library(tradeSeq, quietly = TRUE)
library(dynplot)
library(destiny, quietly = TRUE)})

set.seed(12)

dataset <- generate_dataset(
  model = model_linear(num_milestones = 2),
  num_cells = 500,
  num_features = 5000,
 differentially_expressed_rate = 0.01,
  normalise = T
)

counts <- dataset$counts %>% as.matrix() %>% t()
colnames(counts) <- dataset$cell_ids
norms <- dataset$expression %>% t()

simu_dat <- SingleCellExperiment(assays = List(counts = counts, norms = norms))

pca <- prcomp(t(assays(simu_dat)$norms), scale. = FALSE)
rd <- pca$x[,1:3]

reducedDims(simu_dat) <- SimpleList(PCA = rd)
cl <- Mclust(rd, G = 4)$classification
colData(simu_dat)$GMM <- cl

simu_slingshot <- slingshot(simu_dat, clusterLabels = 'GMM', reducedDim = 'PCA')

lin <- getLineages(simu_dat, clusterLabels = colData(simu_slingshot)$GMM)

crv <- getCurves(lin) %>% SlingshotDataSet()

norms <- assays(simu_slingshot)$norms %>% as.matrix()

sce <- fitGAM(counts = counts, pseudotime = slingPseudotime(crv, na = FALSE),
                  sds = crv, nknots = 3, sce = FALSE)

assoRes <- associationTest(sce)
assoRes %>% as_tibble(rownames = "gene") %>% dplyr::arrange(desc(waldStat))

The results show that all genes are significant, and the Wald statistics are extremely high, which does not make sence. Would you mind check this? Thanks!

@koenvandenberge koenvandenberge added the bug Something isn't working label Oct 7, 2019
@koenvandenberge
Copy link
Member

Thank you for letting us know.
It seems like mgcv is not respecting the number of knots that are provided as input. Even though we're telling it to fit a smoother with 3 knots, it fits a smoother with only 2 knots, which sends the associationTest astray since it thinks that the model has 3 knots. Interestingly, this only ever happens when only a single lineage is fitted.
We will look into the best way to solve this.

@koenvandenberge
Copy link
Member

Hi, this should now be fixed.
Thanks again for letting us know about this issue, and please let us know if you encounter any further issues.

Best,
Koen

zouter added a commit to dynverse/tradeSeq that referenced this issue Oct 7, 2019
fix single lineage issues, assocTest and startVsEnd, statOmics#17
@SONGDONGYUAN1994
Copy link
Author

Thanks guys!

@cstrlln
Copy link

cstrlln commented Mar 14, 2024

I'm actually having this issue while running a single lineage, any suggestion with the number of knots?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
bug Something isn't working
Projects
None yet
Development

No branches or pull requests

3 participants