tradeSeq | Slingshot下游 沿轨迹分析基因表达

tradeSeq是一个可以沿着轨迹分析基因表达的R包，也可以用来分析特定pathway沿轨迹的表达，具体方法见后续分享。
1.准备文件

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.拟合负二项式模型确定nknots，这一步很慢，一定要并行运算

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)#查看收敛基因个数

3.谱系内比较

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出图类似，接下来可以筛选谱系间和谱系内都显著表达的基因，确定每个谱系特异表达的基因。

#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)