x <- 1:100
y <- cos(x/15)+rnorm(length(x), sd=.3)

library(splines)
b <- bs(x, knots = c(25, 50, 75))

f <- lm(y~b)

head(b)

plot(x,y)
lines(x,predict(f))


# for gene expression (e.g., limma/edgeR), could make a design matrix with
# mm <- model.matrix(~b) .. where 'x' above corresponds to the inferred pseudotime
# and you'd have to match the knots to the scale of the pseudotime