In statOmics/tradeSeq: trajectory-based differential expression analysis for sequencing data

library(reticulate)

p2 <- "/Library/Frameworks/Python.framework/Versions/3.8/bin/python3.8"
reticulate::use_python(p2)

import cellrank as cr
import pandas as pd

adata = cr.datasets.pancreas_preprocessed()
adata = adata[adata.obs['dpt_pseudotime'].argsort()].copy()

cr.tl.transition_matrix(adata)
cr.tl.terminal_states(adata, n_states=3, cluster_key="clusters")
cr.tl.lineages(adata)

adata.obs['dpt_pseudotime']

cr.pl.lineages(adata, same_plot=True)

Convert to SCE and extract relevant info

suppressPackageStartupMessages({
  library(SingleCellExperiment)
  library(zellkonverter)
})

adata <- py$adata
sce <- AnnData2SCE(adata)
sce

table(colData(sce)$clusters)

# cell-level weights
cw <- reducedDims(sce)$terminal_states_memberships 
#apply(cw, 2, function(x) print(hist(x)))
boxplot(cw[,1] ~ colData(sce)$clusters)
boxplot(cw[,2] ~ colData(sce)$clusters)
boxplot(cw[,3] ~ colData(sce)$clusters)
colnames(cw) <- c("Epsilon", "Alpha", "Beta")

# pseudotime
pt <- colData(sce)$dpt_pseudotime 
pseudotime <- cbind(pt, pt, pt)
colnames(pseudotime) <- colnames(cw)

Fit GAM using tradeSeq

library(tradeSeq)

## note that counts in object are not integers!
## we will proceed anyway using rounded coutns, as proof-of-concept.
assays(sce)$X[1:5, 1:5]

sceGAM <- fitGAM(round(assays(sce)$X), 
                 pseudotime=pseudotime,
                 cellWeights=cw,
                 nknots=6,
                 verbose=FALSE)

Compare start vs end points for each lineage

startEndRes <- startVsEndTest(sceGAM, 
                              l2fc = log2(1.5))
topGene <- rownames(sce)[order(startEndRes$waldStat, decreasing=TRUE)[3]]
plotSmoothers(sceGAM, assays(sce)$X, topGene) +
  ggplot2::ggtitle(topGene)

Compare end points between lineages

diffEndRes <- diffEndTest(sceGAM, 
                          l2fc = log2(1.5))
topGene <- rownames(sce)[order(diffEndRes$waldStat, decreasing=TRUE)[1]]
plotSmoothers(sceGAM, assays(sce)$X, topGene) +
  ggplot2::ggtitle(topGene)

earlyDE test at later stages

edRes <- earlyDETest(sceGAM,
                     l2fc = log2(1.5),
                     knots = c(4,6))
topGene <- rownames(sce)[order(edRes$waldStat, decreasing=TRUE)[2]]
plotSmoothers(sceGAM, assays(sce)$X, topGene) +
  ggplot2::ggtitle(topGene)

statOmics/tradeSeq documentation built on Jan. 19, 2024, 8:26 p.m.