在發(fā)育過程中，細(xì)胞會(huì)對(duì)刺激做出反應(yīng)，在整個(gè)生命過程中，從一種功能性“狀態(tài)”轉(zhuǎn)變?yōu)榱硪环N功能性“狀態(tài)”。擬時(shí)序分析（pseudo-time），它指通過構(gòu)建細(xì)胞間的變化軌跡來重塑細(xì)胞隨著時(shí)間的變化過程。用于擬時(shí)序分析的軟件和算法很多，目前文章中用到 多的是 Monocle 軟件，這也是是 CNS 主流期刊的常用軟件。
對(duì)于單細(xì)胞的數(shù)據(jù)做擬時(shí)序分析基本的流程是挑選一些感興趣并且可能有分化關(guān)系的亞群，然后來分析它們之間的分化關(guān)系。 那么對(duì)于空間轉(zhuǎn)錄組數(shù)據(jù)有沒有什么新花樣呢？ 其實(shí)目前的空間轉(zhuǎn)錄組測(cè)序已發(fā)表的文章還沒太注重細(xì)胞分化這一內(nèi)容，可能大家做空轉(zhuǎn)數(shù)據(jù)挖掘的時(shí)候還沒太把關(guān)注點(diǎn)放到分化關(guān)系這上面吧，那么今天就來給大家寫寫自己的突發(fā)奇想吧！
既然是做的空間轉(zhuǎn)錄組，那么我們就要有效的利用起來空間位置信息。選擇的亞群不能太離散。其次，我們可以分析相同細(xì)胞類型的亞群在空間位置的分化關(guān)系，也可以按照病理狀態(tài)分布的漸進(jìn)變化來選擇區(qū)域做擬時(shí)序分化。
讀取 seurat 對(duì)象
library(monocle)
# 讀取前面保存的 seurat 對(duì)象文件
combin.data <- readRDS("combin.data.RDS")
這里我們選擇大腦皮層區(qū)的亞群作為示例進(jìn)行細(xì)胞分化分析。
# 展示亞群分布
SpatialDimPlot(combin.data, label = TRUE, label.size = 3)
# 展示皮層 marker 的分布
SpatialFeaturePlot(combin.data, features = c("Stx1a"))
根據(jù)亞群和皮層 marker STX1A 的表達(dá)分布，我們選擇 2、5、7、9、10 號(hào)群作為皮層的亞群。
# 選擇要分析的亞群 subdata <- subset(combin.data, idents = c(2,5,7,9,10))
導(dǎo)入 seurat 對(duì)象，構(gòu)建 CellDataSet 對(duì)象
這里我們使用 monocle2 做擬時(shí)序分析。monocle2 做擬時(shí)序分析需要三個(gè)矩陣文件：expression_matrix（表達(dá)矩陣，一般是 raw count）、phenoData（細(xì)胞 meta 文件，可以包括細(xì)胞的樣本、亞群等信息）、featureData（gene 的 meta 文件，注意要包括 gene_short_name 這一列）。
# 單細(xì)胞 count 文件、細(xì)胞類型注釋文件、基因注釋文件
expression_matrix = combin.data@assays$Spatial@counts
cell_metadata <- data.frame(group = combin.data[['orig.ident']],clusters = Idents(combin.data))
gene_annotation <- data.frame(gene_short_name = rownames(expression_matrix), stringsAsFactors = F)
rownames(gene_annotation) <- rownames(expression_matrix)
##### 新建 CellDataSet object
pd <- new("AnnotatedDataFrame", data = cell_metadata)
fd <- new("AnnotatedDataFrame", data = gene_annotation)
HSMM <- newCellDataSet(expression_matrix,
phenoData = pd,
featureData = fd,
expressionFamily=negbinomial.size())
monocle2 做擬時(shí)序分析主要包括三個(gè)步驟：
1、基因篩選。選擇要作為擬時(shí)序分化依據(jù)的基因，軟件官方提供了幾種可選的方法。
A、選擇不同時(shí)間點(diǎn)分化差異基因，這應(yīng)該是比較理想的情況，需要你提供的數(shù)據(jù)是根據(jù)不同分化時(shí)間點(diǎn)取樣的；
B、選擇表達(dá)離散度高的基因，也就是在所有細(xì)胞里變化比較大的基因，如果只是想看某個(gè)亞群里細(xì)胞之間的分化選擇這種方法是比較合適的；
C、選擇亞群間差異大的基因，一般想比較多個(gè)亞群間的分化關(guān)系選擇這種方法效果會(huì)好一點(diǎn)；
D、選擇一些已知跟分化相關(guān)的基因。
2、數(shù)據(jù)降維。Monocle2 使用 DDRTree 的非線性重建算法對(duì)細(xì)胞進(jìn)行排序。
3、細(xì)胞排序。根據(jù)細(xì)胞的降維后的結(jié)果給每個(gè)細(xì)胞計(jì)算一個(gè)分化時(shí)間。不過這里有一點(diǎn)需要注意，做細(xì)胞排序這一步可以自己指定分化起點(diǎn)，要不然算法會(huì)隨機(jī)選擇一端作為起點(diǎn)，也就是說計(jì)算出來的分化時(shí)間有可能是倒過來的，即起點(diǎn)是終點(diǎn)，終點(diǎn)是起點(diǎn)。
# 評(píng)估 SizeFactors
HSMM_myo <- estimateSizeFactors(HSMM)
# 計(jì)算離散度
HSMM_myo <- estimateDispersions(HSMM_myo)
# 計(jì)算差異
diff_test_res1 <- differentialGeneTest(HSMM_myo,fullModelFormulaStr = '~clusters', cores = 4)
# 選擇差異基因
ordering_genes <- subset(diff_test_res1, qval < 0.05)[,'gene_short_name']
# 基因過濾
HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
# 降維
HSMM_myo <- reduceDimension(HSMM_myo, max_components=2, method = 'DDRTree')
# 細(xì)胞排序
HSMM_myo <- orderCells(HSMM_myo)
Monocle 結(jié)果可視化：
# 查看亞群分化關(guān)系
plot_cell_trajectory(HSMM_myo, color_by = "clusters",show_branch_points = F,cell_size =0.5)
# 按亞群分開展示 plot_cell_trajectory(HSMM_myo, color_by = "clusters",cell_size =0.5) + facet_wrap(~clusters, nrow = 4)
# 分樣本展示 plot_cell_trajectory(HSMM_myo, color_by = "orig.ident",show_branch_points = F,cell_size =0.5)
從亞群和樣本的分布來看其實(shí)還是體現(xiàn)了一定的樣本的異質(zhì)性的，同一樣本的亞群更容易分布到同一個(gè)分支上，這里的 2、7、10 號(hào)群有點(diǎn)分不太開，這時(shí)候我們就要思考一下，對(duì)于空轉(zhuǎn)的數(shù)據(jù)有時(shí)候是不是單獨(dú)分析某個(gè)樣本上的分化是不是更合適一些。
選擇單個(gè)樣本的數(shù)據(jù)進(jìn)行分析
這里我們選擇小鼠大腦 Sagittal-Anterior1 的樣本的皮層的亞群?jiǎn)为?dú)進(jìn)行分析。
#### 單獨(dú)繪制某個(gè)樣本的圖片
library("cowplot")
p1 <- SpatialPlot(combin.data,crop=FALSE,images='image')# 這里 Sagittal-Anterior1 樣本的鏡像圖片的名字是 image，具體名稱應(yīng)跟前期讀取鏡像文件時(shí)設(shè)置的名字對(duì)應(yīng)。
p2 <- SpatialFeaturePlot(combin.data, features = c("Stx1a"),images='image')
plot_grid(p1, p2)
選擇皮層 marker 高表達(dá)的 2、7、10 群進(jìn)行擬時(shí)序分析。
# 選擇單個(gè)樣本 Sagittal-Anterior1 的亞群進(jìn)行分析
subdata <- subset(combin.data, orig.ident == 'Sagittal-Anterior1')
subdata <- subset(subdata, idents = c(2,4,10))
# 單細(xì)胞 count 文件、細(xì)胞類型注釋文件、基因注釋文件
expression_matrix = subdata@assays$Spatial@counts
cell_metadata <- data.frame(group = subdata[['orig.ident']],clusters = Idents(subdata))
gene_annotation <- data.frame(gene_short_name = rownames(expression_matrix), stringsAsFactors = F)
rownames(gene_annotation) <- rownames(expression_matrix)
##### 新建 CellDataSet object
pd <- new("AnnotatedDataFrame", data = cell_metadata)
fd <- new("AnnotatedDataFrame", data = gene_annotation)
HSMM <- newCellDataSet(expression_matrix,
phenoData = pd,
featureData = fd,
expressionFamily=negbinomial.size())
# 評(píng)估 SizeFactors
HSMM_myo <- estimateSizeFactors(HSMM)
# 計(jì)算離散度
HSMM_myo <- estimateDispersions(HSMM_myo)
# 計(jì)算差異
diff_test_res1 <- differentialGeneTest(HSMM_myo,fullModelFormulaStr = '~clusters', cores = 4)
# 選擇差異基因
ordering_genes <- subset(diff_test_res1, qval < 0.05)[,'gene_short_name']
# 基因過濾
HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
# 降維
HSMM_myo <- reduceDimension(HSMM_myo, max_components=2, method = 'DDRTree')
# 細(xì)胞排序
HSMM_myo <- orderCells(HSMM_myo)
結(jié)果可視化：
# 查看亞群分化關(guān)系
plot_cell_trajectory(HSMM_myo, color_by = "clusters",show_branch_points = F,cell_size =0.5)
# 按亞群分開展示
plot_cell_trajectory(HSMM_myo, color_by = "clusters",cell_size =0.5) + facet_wrap(~clusters, nrow = 4)
這里我們可以發(fā)現(xiàn)，之前選擇兩個(gè)樣本一起的 5 個(gè)亞群分析的時(shí)候這 3 個(gè)群是有點(diǎn)區(qū)分不開的，在這次重新分析的時(shí)候卻區(qū)分的很明顯。這也說明選擇的亞群范圍不一樣 monocle 得到的結(jié)果差異是很大的，究其原因可能是加入某些差異比較大的亞群或細(xì)胞后會(huì)進(jìn)一步壓縮原本差異比較小的亞群之間的差異分布。因?yàn)槭撬惴ㄍ瑫r(shí)計(jì)算多組差異變化難免會(huì)出現(xiàn)這種情況，所以我們?cè)谝婚_始選擇亞群上就要稍微注意一點(diǎn)。
# 繪制分化時(shí)間
cell_Pseudotime <- data.frame(pData(HSMM_myo)$Pseudotime)
rownames(cell_Pseudotime) <- rownames(cell_metadata)
plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime",cell_size =0.5)
從分化時(shí)間的分布來看三個(gè)亞群的分化大概順序是從 7 號(hào)群一 >10 號(hào)群一 >2 號(hào)群，由于 monocle 的判斷的分化起點(diǎn)是隨機(jī)的，所以也有可能是倒過來的順序。
我們?cè)賮砜匆幌缕?marker 基因 Stx1a 的在分化過程中的表達(dá)分布
# 繪制 Stx1a 基因分化時(shí)間的變化
data_subset <- HSMM_myo['Stx1a',]
plot_genes_in_pseudotime(data_subset, color_by = "clusters")
marker 基因的表達(dá)分布基本是由低到高再降低的趨勢(shì)。
下面重點(diǎn)來了！我們把細(xì)胞的分化時(shí)間對(duì)應(yīng)到組織切片上。
# 把分化時(shí)間對(duì)應(yīng)到到組織切片位置上
combin.data[['Pseudotime']] <- 0
combin.data[['Pseudotime']][rownames(cell_Pseudotime),] <- cell_Pseudotime
SpatialFeaturePlot(combin.data, features = c("Pseudotime"),images='image')
這時(shí)候我們就發(fā)現(xiàn)結(jié)果有點(diǎn)意思了，皮層細(xì)胞的分化在組織切片上是從外向內(nèi)的方向的，當(dāng)然也有可能實(shí)際上是反過來的，也就是從里往外的，總之就是細(xì)胞的分化是存在一定的空間位置規(guī)律的。
接著，我們?cè)賮砜匆幌录?xì)胞分化過程中表達(dá)變化顯著的基因有哪些。
# 分析分化時(shí)間顯著變化的基因
diff_test_res2 <- differentialGeneTest(HSMM_myo[ordering_genes,],fullModelFormulaStr = "~sm.ns(Pseudotime)",cores = 4)
# 選擇 top50 基因繪圖
top50_gene = diff_test_res2[order(diff_test_res2$qval),][1:50,'gene_short_name']
plot_pseudotime_heatmap(HSMM_myo[top50_gene,], num_clusters = 6,cores = 1,show_rownames = T)
從 top50 差異基因的情況來看，血紅蛋白基因和線粒體基因在小鼠大腦皮層切片外層的表達(dá)比較高。我們也可以把所有差異基因都輸出來，將不同表達(dá)趨勢(shì)的基因模塊分別進(jìn)行功能富集，這樣就可以知道隨著皮層細(xì)胞的在空間位置的分化哪些功能發(fā)生了變化。
sig_gene_names <- subset(diff_test_res2, qval < 0.05)[,'gene_short_name']
heatmap = plot_pseudotime_heatmap(HSMM_myo[sig_gene_names,], num_clusters = 6,cores = 1,show_rownames = F,return_heatmap=T)
###get cluster info
row_cluster = cutree(heatmap$tree_row,k=6)
write.table(row_cluster,file=paste("_gene_clusters.xls",sep=""),sep="\t",row.names=T)
細(xì)胞更有意思呢！因?yàn)槭纠龜?shù)據(jù)里沒有病理信息，如果結(jié)合病理特征來做擬時(shí)序分化分析應(yīng)該能發(fā)現(xiàn)更多有價(jià)值的東西！
轉(zhuǎn)載：簡(jiǎn)生信（侵刪）
更多伯豪生物人工服務(wù)：
