tradeSeq是一個(gè)可以沿著軌跡分析基因表達(dá)的R包，也可以用來(lái)分析特定pathway沿軌跡的表達(dá)，具體方法見后續(xù)分享。
1.準(zhǔn)備文件

library(BiocParallel)
library(tradeSeq)
#load slingshot sds文件
load(file = "./scsds.RData")
#load sc counts
load(file = "./sccounts.RData")
pseudotime <- slingPseudotime(sds, na = FALSE)
cellWeights <- slingCurveWeights(sds)
crv <- SlingshotDataSet(sds)

2.擬合負(fù)二項(xiàng)式模型確定nknots，這一步很慢，一定要并行運(yùn)算

icMat <- evaluateK(counts = counts, 
                   sds = crv, nGenes = 200,
                   k = 3:10,verbose = T, plot = TRUE)
sce <- fitGAM(counts = counts, 
              pseudotime = pseudotime, 
              cellWeights = cellWeights,
              nknots = 6, verbose = TRUE,
              BPPARAM = MulticoreParam(20),parallel=T)
table(rowData(sce)$tradeSeq$converged)#查看收斂基因個(gè)數(shù)

3.譜系內(nèi)比較

startRes <- startVsEndTest(sce,lineages=T)
oStart <- order(startRes$waldStat, decreasing = TRUE)
sigGeneStart <- names(sce)[oStart[1]]
plotSmoothers(sce, counts, gene = sigGeneStart)
plotGeneCount(crv, counts, gene = sigGeneStart)

4.譜系間比較

endRes <- diffEndTest(sce,pairwise=T)
o <- order(endRes$waldStat, decreasing = TRUE)
sigGene <- names(sce)[o[1]]
plotSmoothers(sce, counts, sigGene)
plotGeneCount(crv1, counts, gene = sigGene)

出的圖和monocle出圖類似，接下來(lái)可以篩選譜系間和譜系內(nèi)都顯著表達(dá)的基因，確定每個(gè)譜系特異表達(dá)的基因。

#Lineage1
tmp1 <- startRes[startRes$pvalue_lineage1<0.01,] %>% rownames()
tmp2 <- endRes[endRes$pvalue_1vs2<0.01,]%>%rownames()
intersect(tmp1,tmp2)

參考：The tradeSeq workflow ? tradeSeq (statomics.github.io)