Hi everyone!

I am currently dealing with pseudotime trajectory analysis and trying to find out genes that vary along pseudotime axis. I have two objects, one of the counts ~20k cells in state_1 and ~2k cells in state_2, another counts ~2k cells in state_1 and ~1k cells in state_2, the first one possibly being biased by the number of cells.

My code is the following:

test_df <- graph_test(cds, neighbor_graph = "principal_graph", k=500, alternative="two.sided", cores=18)
test_df$gene_name <- row.names(test_df)
test_df <- test_df[test_df$q_value<0.05,]
gene_module_df <- find_gene_modules(cds[test_df$gene_name,])
ggplot(gene_module_df, aes(x=dim_1, y=dim_2, color=module))+
    geom_point(size=0.5)+
    ggtitle("UMAP for genes")
agg <- aggregate_gene_expression(cds, gene_module_df)
pheatmap(agg,
     annotation_col = cds@colData[,c("pseudotime", "sample")],
     cluster_cols = FALSE,
     breaks=c(-2,0,2),
     scale="row",
     cluster_rows = T,
     show_colnames = FALSE,
     show_rownames = T,
     use_raster=FALSE)

The number of genes that pass threshold of q_value<0.05 depends on the number of nearest neighbors (k) that I choose. My idea is that k should depend on the cells number to explain differences driven by pseudotime, not by fluctuations. However, it is still difficult to understand which k would be the most reasonable.

One another question is whether I am right to use alternative="two.sided", as I want to see not only genes that are upregulated along the axis?

Any advice would be appreciated! Thank you in advance!