hello，大家好其障，今天我們來分享一個(gè)新的分析點(diǎn)锋八，跨樣本的基因表達(dá)動(dòng)態(tài)磺浙。文章在A statistical framework for differential pseudotime analysis with multiple single-cell RNA-seq samples，很有意義的一個(gè)分析點(diǎn)跑慕。

看看簡(jiǎn)介

基于單細(xì)胞 RNA-seq (scRNA-seq) 數(shù)據(jù)的偽時(shí)間分析已被廣泛用于研究沿連續(xù)生物過程的動(dòng)態(tài)基因調(diào)控程序万皿，如細(xì)胞分化、免疫反應(yīng)和疾病發(fā)展核行。 現(xiàn)有的偽時(shí)間分析方法主要解決重建細(xì)胞偽時(shí)間軌跡和推斷一個(gè)生物樣本中沿重建軌跡的基因表達(dá)變化的問題牢硅。 隨著越來越多地對(duì)多個(gè)患者樣本進(jìn)行 scRNA-seq 研究，比較不同樣本的基因表達(dá)動(dòng)態(tài)已成為一種新的需求芝雪，但缺乏系統(tǒng)的分析解決方案减余。 （個(gè)人翻譯，總的來說惩系，其實(shí)我們更加關(guān)注相對(duì)于normal位岔，疾病樣本在軌跡分化的過程中的基因動(dòng)態(tài)）。

其實(shí)軟件Lamian就是對(duì)多樣本軌跡分析的補(bǔ)充堡牡，考慮到單細(xì)胞數(shù)據(jù)樣本包含多種生物學(xué)條件（例如抒抬，年齡、性別晤柄、樣本類型擦剑、疾病狀態(tài)等），這個(gè)分析軟件框架有以下功能芥颈。

• (1) 構(gòu)建細(xì)胞偽時(shí)間軌跡惠勒，評(píng)估軌跡分支結(jié)構(gòu)的不確定性。
• (2) 評(píng)估與樣本生物學(xué)因素相關(guān)的分支結(jié)構(gòu)的變化浇借。
• (3) 確定沿偽時(shí)間的細(xì)胞豐度和基因表達(dá)的變化捉撮。
• (4) 表征樣本生物學(xué)差異如何修改細(xì)胞豐度和基因表達(dá)的偽時(shí)間動(dòng)態(tài)。

重要的是妇垢，在識(shí)別與偽時(shí)間和樣本生物學(xué)差異相關(guān)的細(xì)胞豐度或基因表達(dá)變化時(shí)巾遭，Lamian 考慮了其他現(xiàn)有偽時(shí)間方法未考慮的生物樣本的變異性。 因此闯估，在分析多樣本數(shù)據(jù)時(shí)灼舍，Lamian 能夠更恰當(dāng)?shù)乜刂棋e(cuò)誤發(fā)現(xiàn)率 (FDR)，這是任何其他方法都無法提供的特性涨薪。 （這個(gè)可能確實(shí)需要關(guān)注）骑素。

接下來我們?cè)敿?xì)看看分析過程(使用示例數(shù)據(jù))

options(warn=-1)
suppressMessages(library(Lamian))

Module 1: tree variability

Lamian 的模塊 1 設(shè)計(jì)用于檢測(cè)偽時(shí)間樹形結(jié)構(gòu)中分支的穩(wěn)定性。 自動(dòng)隨機(jī)選取偽時(shí)間樹結(jié)構(gòu)中的分支刚夺，然后測(cè)試它們的檢測(cè)率（這算是檢驗(yàn)分支的穩(wěn)定性）献丑。 在每次cell-bootstrapping后末捣，重建樹結(jié)構(gòu)并重新識(shí)別所有分支。 同時(shí)應(yīng)用 Jaccard 統(tǒng)計(jì)量和重疊系數(shù)作為評(píng)估引導(dǎo)設(shè)置中的任何分支是否與原始分支之一匹配的統(tǒng)計(jì)量创橄。 分支的檢測(cè)率定義為分支發(fā)現(xiàn)它匹配的bootstrap settings的百分比箩做。 模塊 1 還報(bào)告和測(cè)試每個(gè)分支中的樣本比例。 （有點(diǎn)費(fèi)勁）妥畏。

數(shù)據(jù)準(zhǔn)備

需要一個(gè)基因的細(xì)胞表達(dá)矩陣邦邦、細(xì)胞的低維表示和細(xì)胞注釋（每個(gè)細(xì)胞屬于哪個(gè)樣本）。這里我們就用示例數(shù)據(jù)分析看看

data(man_tree_data)

man_tree_data 是一個(gè)包含基因按細(xì)胞表達(dá)矩陣醉蚁、細(xì)胞的低維表示（PCA）和細(xì)胞注釋的列表燃辖，其中第一列是細(xì)胞名稱，第二列是樣本名稱网棍，行名稱是細(xì)胞 名稱黔龟。

推斷樹形結(jié)構(gòu)

set.seed(12345)
res = infer_tree_structure(pca = man_tree_data[['pca']], 
                           expression = man_tree_data[['expression']], 
                           cellanno = man_tree_data[['cellanno']], 
                           origin.marker = c('CD34'), 
                           number.cluster = 5,
                           xlab='Principal component 1', 
                           ylab = 'Principal component 2')

這個(gè)軟件看來是按照pca的成分來構(gòu)建樹形結(jié)構(gòu)，我想我們還是可以替換成tSNE或者UMAP确沸，還要指定起始位點(diǎn)的基因特征捌锭，分支數(shù)量，看來這個(gè)軟件罗捎，需要的前期知識(shí)有點(diǎn)多观谦，這里的cluster我們就對(duì)應(yīng)我們聚類的結(jié)果或者細(xì)胞注釋的結(jié)果。

繪圖

plotmclust(res, cell_point_size = 0.5, 
           x.lab = 'Principal component 1', 
           y.lab = 'Principal component 2')

圖片.png

這里我們還是真的需要UMAP的降維結(jié)果最好桨菜。

估計(jì)分支的不穩(wěn)定性

We can call the function evaluate_uncertainty() to evaluate the tree topology uncertainty. We suggested that users set n.permute = 100 or more to ensure enough randomness to construct the null distribution, but here we set n.permute = 5 in order to provide a simplified example.

result <- evaluate_uncertainty(res, n.permute=5)

這里作為簡(jiǎn)單的示例豁状，n.permute 設(shè)置為5，在實(shí)際應(yīng)用中設(shè)置 n.permute = 100 或更多時(shí)倒得，結(jié)果會(huì)更有意義泻红。

Module 2: Evaluate differential topoloty

Lamian 的模塊 2 首先識(shí)別樣本之間樹拓?fù)涞淖兓缓笤u(píng)估是否存在與樣本生物學(xué)差異相關(guān)的差異拓?fù)渥兓?對(duì)于每個(gè)樣本霞掺，Lamian 計(jì)算每個(gè)樹分支中的cell比例谊路，稱為分支cell比例。 由于零或低比例可以反映分支的缺失或耗盡菩彬，因此可以使用分支cell比例變化來描述樹拓?fù)涞淖兓?對(duì)于多個(gè)樣本缠劝，Lamian 通過估計(jì)跨樣本的分支cell比例的方差來表征每個(gè)分支的跨樣本變異。 此外骗灶，可以擬合回歸模型來測(cè)試分支單元格比例是否與樣本生物學(xué)差異相關(guān)惨恭。 這允許識(shí)別不同條件之間的樹拓?fù)渥兓缭诓±鄬?duì)于normal耙旦。

data = result[[2]]
rownames(data) <- c('HSC->lymphocyte', 'HSC->myeloid', 'HSC->erythroid')
design = data.frame(sample= paste0('BM',1,8), sex = c(0, rep(1,4), rep(0,3)))
branchPropTest(data = data, design = design)

Module 3: Trajectory differential tests about gene expression

XDE test(這個(gè)大家自行查看脱羡，就不多做解釋了)。

Res <- lamian.test(expr = expdata$expr, 
                   cellanno = expdata$cellanno, 
                   pseudotime = expdata$pseudotime, 
                   design = expdata$design, 
                   test.type = 'variable', 
                   permuiter = 5)

downstream analysis 1: visualize and cluster XDE genes based on their multiple-sample temporal patterns across sample covariates

stat <- Res$statistics
stat <- stat[order(stat[,1], -stat[,3]), ]
diffgene <- rownames(stat[stat[, grep('^fdr.*overall$', colnames(stat))] < 0.05,])
Res$populationFit <- getPopulationFit(Res, gene = diffgene, type = 'variable')
## clustering
Res$covariateGroupDiff <- getCovariateGroupDiff(testobj = Res, gene = diffgene)
Res$cluster <- clusterGene(Res, gene = diffgene, type = 'variable', k=5)
colnames(Res$populationFit[[1]]) <- colnames(Res$populationFit[[2]]) <- colnames(Res$expr) 
plotXDEHm(Res, cellWidthTotal = 180, cellHeightTotal = 350, subsampleCell = FALSE, sep = ':.*')

圖片.png

圖片還是不錯(cuò)的。

We can plot the cluster mean and difference.

## plotClusterMeanAndDiff
plotClusterMeanAndDiff(Res)

We can also plot the cluster difference seperately by calling plotClusterDiff(testobj=Res, gene = diffgene).

downstream analysis 2: Gene ontology (GO) enrichment analysis

We can further check out whether the DDG in each cluster have any enriched genome functions by applying Gene Ontology (GO) analysis library(topGO); goRes <- GOEnrich(testobj = Res, type = 'variable'). Please note that there are no enriched GO terms in any of the clusters in this simplified example since we only subset 1,000 genes. In practice, if there are significant GO terms in the clusters, we can generate a heatmap of the top GO terms using the command plotGOEnrich(goRes = goRes).

TDE test(這個(gè)大家也自行查看)

Res <- lamian.test(expr = expdata$expr, 
                   cellanno = expdata$cellanno, 
                   pseudotime = expdata$pseudotime, 
                   design = expdata$design, 
                   test.type = 'time', 
                   permuiter = 5)
diffgene <- rownames(Res$statistics)[Res$statistics[,1] < 0.05]
Res$populationFit <- getPopulationFit(Res, gene = diffgene, type = 'time')
Res$cluster <- clusterGene(Res, gene = diffgene, type = 'time', k=3)
plotTDEHm(
  Res,
  subsampleCell  = FALSE,
  showCluster = TRUE,
  type = 'time',
  cellWidthTotal = 200,
  cellHeightTotal = 200
)

圖片.png

plotClusterMean(testobj=Res, cluster = Res$cluster, type = 'time')

圖片.png

downstream analysis 2: GO analysis for TDE genes

Similar to XDE test, we can further proceed to identify the enriched GO terms for each gene cluster using commands library(topGO); goRes <- GOEnrich(testobj = Res, type = 'time', sep = ':.*'). Please note that there are no enriched GO terms in any of the clusters in this simplified example since we only subset 1,000 genes. In practice, if there are enriched GO terms, we can plot them using commands plotGOEnrich(goRes = goRes).

Module 4: trajectory differential test based on cell composition(這部分跟上面的類似锉罐，就不多介紹了)

總之帆竹，軟件有其特點(diǎn)，可以借鑒氓鄙，取其長(zhǎng)處馆揉，為己所用业舍。

生活很好抖拦，有你更好