Hi all,

I am trying to create a predictor matrix for a range of pseudotime to plot a heatmap like your fig 6c using my dataset. After following ./fletcher/20180830_caseStudyFletcher_v5_tradeR.Rmd, I could run fitGAM on my dataset (~ 14k genes and 1k cells) with only one lineage specified, as well as run the associationTest. However, I got error for line 365 in 20180830_caseStudyFletcher_v5_tradeR.Rmd, i.e., df <- .getPredictRangeDf(gamList[[1]], lineageId = 1, nPoints = 100).

The error message I got is that Error in vars[, grep(colnames(vars), pattern = "t[1-9]")] <- 0 : incorrect number of subscripts on matrix. I think it is because the structure & data type of my dm = gamList[[1]] is different from the data frame that you had in your original analysis.

Could you specify the list structure of gamList in your analysis? How would you recommend to extract the data frame (which is dm = gamList[[1]] in your analysis) from a gamList when there is only one lineage in my analysis? Or could you suggest how to modify the utility function .getPredictRangeDf? I had a long and hard time trying to interpret what each line in .getPredictRangeDf is trying to do, without knowing your data structure from the very beginning. I could not reproduce your analysis and see what your data structure was like because of the concern of time, for instance, zinbFit and fitGAM would take a very long time to run.

My gamList is a large list of length 14k, equal to the number of genes. Thus, gamList[[1]] is a gam glm lm object for one gene in my output, rather than a data frame. I suspect that it is because I only specified one lineage for my dataset, for the sake of biological interpretation.

Thank you!

Update: After setting sce = T, my gamList object is a SingleCellExperiment object now. Then I have two questions in order to create the predictor matrix:

(1) Does df <- .getPredictRangeDf(gam_list[[1]], lineageId = 1, nPoints = 100) intend to create evenly-distributed grids across the pseudotime, and also set the starting point to be zero? If this is what .getPredictRangeDf intended to do, I have done this step correctly.
(2) When I tried to use predict function to create the predictor matrix following line 366 in 20180830_caseStudyFletcher_v5_tradeR.Rmd, i.e., y <- do.call(rbind, lapply(gam_list[top200Assoc], predict, newdata = df, type = "link")), I got an error because predict cannot be applied to SingleCellExperiment object, which is gam_list[top200Assoc] in the original code. I assume that predict is stats::predict and there is some way to extract the regression model from gam_list and feed the model to predict, is this correct? If so, any suggestion how I could find such model in the gam_list object? Or should I actually use tradeSeq::predictSmooth instead?

> y <- do.call(rbind, lapply(gam_list[top200Assoc], stats::predict, newdata = df, type = "link"))
Error in h(simpleError(msg, call)) : 
  error in evaluating the argument 'args' in selecting a method for function 'do.call': no applicable method for 'predict' applied to an object of class "c('SingleCellExperiment', 'RangedSummarizedExperiment', 'SummarizedExperiment', 'Vector', 'Annotated', 'vector_OR_Vector')"

Hi there,

I'm running the code for bulk RNA seq time course data. But my gamList is not a large list composed by each gene. So I cannot not proceed to look at fit for each gene.

fit GAM

gamList <- fitGAM(counts, pseudotime=time, cellWeights=weights, nknots=6)

May I ask do you know how to fix it?

Hello,

Could you post where you got the clustering results for the Fletcher paper? Or is the answer to just re-run their analysis? The data they posted has no cell types (or its mostly NA).

Thanks!

Hello,

Thank you for producing this tool. I'm running through the bulk RNA-seq tutorial and have gotten up to line 53 without issue. I would like to carry out the following steps (line 55-59) but receive some errors:

look at fit for six random genes

set.seed(81)
id <- sample(1:length(gamList), size=6)
par(mfrow=c(2,3))
for(ii in 1:6) plotSmoothers(gamList[[id[ii]]], main=resPat[id[ii],"pvalue"], ylim=c(3,9))

This told me that there was no object called "resPat", so I ran the command below (line 63):

test for different expression pattern

resPat <- patternTest(gamList, nPoints=8)

When I try to re-run lines 55-59 it gives me this error:
Error in .local(models, ...) :
unused arguments (main = 3.14194401196399e-05, ylim = c(3, 9))

Are you able to help me with this?

Many thanks in advance,

Rosanna

I'm trying to run your robustK.R script, but keep getting an error for line 10, where it tries to use the function updateCDS. Where is that function defined? I haven't found it in monocle, monocle3, nor in rafalib.

Thank you for all of your work on TradeSeq!