scRNA-seq擬時(shí)間分析數(shù)據(jù)處理

擬時(shí)間分析簡介

擬時(shí)序（pseudotime）分析，又稱細(xì)胞軌跡（cell trajectory）分析碴巾，通過擬時(shí)分析可以推斷出發(fā)育過程細(xì)胞的分化軌跡或細(xì)胞亞型的演化過程，在發(fā)育相關(guān)研究中使用頻率較高。

主要基于關(guān)鍵基因的表達(dá)模式法梯，通過學(xué)習(xí)每個(gè)細(xì)胞必須經(jīng)歷的基因表達(dá)變化的序列，在擬時(shí)間中對單個(gè)細(xì)胞進(jìn)行排序，模擬出時(shí)間發(fā)育過程的動(dòng)態(tài)變化立哑。而這個(gè)排序技術(shù)表現(xiàn)是一種在低維空間排布高維數(shù)據(jù)的降維技術(shù)夜惭。（具體描敘請參考：https://www.meiwen.com.cn/subject/zlgsvqtx.html）

如何利用單細(xì)胞數(shù)據(jù)來分析細(xì)胞的分化軌跡或細(xì)胞亞型的演化過程？這里給大家介紹monocle軟件铛绰，該軟件能從數(shù)據(jù)中用無監(jiān)督或半監(jiān)督學(xué)習(xí)獲得這個(gè)軌跡诈茧。無監(jiān)督方式是利用seurat對細(xì)胞進(jìn)行聚類，通過cluster間的差異基因?qū)?xì)胞進(jìn)行排序捂掰。

在做擬時(shí)序分析之前敢会，先要明白幾個(gè)研究目的：

1.?對什么類型的細(xì)胞群進(jìn)行擬時(shí)許分析？這類細(xì)胞群在不同的分化軌跡或細(xì)胞亞型的區(qū)別是否明顯这嚣？

2.?monocle做擬時(shí)序分析時(shí)鸥昏，可與seurat聚類結(jié)果進(jìn)行對接，因此可對seurat的分群結(jié)果進(jìn)行定性的解釋姐帚。

3.?可分析比較組間的差異或者多組樣本間擬時(shí)間分析差別吏垮。

Seurat對象轉(zhuǎn)換成monocle對象

如果你已經(jīng)完成了seurat分類(數(shù)據(jù)已經(jīng)經(jīng)過質(zhì)控)，而且你希望通過seurat的分類后cluster的差異基因來對細(xì)胞進(jìn)行排序罐旗。以下方法可供選擇膳汪。

Seurat2的版本∮容海可直接使用importRDS進(jìn)行seurat對象和monocle對象間的轉(zhuǎn)換旅敷。

????spleen <-readRDS("seurat.analysis.rds")?

????monocle_cds?<-importCDS(spleen,import_all=TRUE)??

Seurat3版本。Seurat3版本沒有importRDS函數(shù)了颤霎，因此需要手動(dòng)構(gòu)建clustered_spleen_monocle對象媳谁。??

????spleen <-readRDS("seurat.analysis.rds")

????data <- as(as.matrix(spleen@assays$RNA@counts), 'sparseMatrix')??##此處需要提供絕對的count值。??

????pd <- new('AnnotatedDataFrame', data = spleen@meta.data)??

????fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data)??

????fd <- new('AnnotatedDataFrame', data = fData)???????

????monocle_cds?<- newCellDataSet(data,?

????????????????????????phenoData = pd,??

????????????????????????featureData = fd,??

????????????????????????lowerDetectionLimit = 0.5,??

????????????????????????expressionFamily = negbinomial.size()) ##newCellDataSet構(gòu)建monocle對象?

使用差異基因進(jìn)行細(xì)胞排序

對細(xì)胞進(jìn)行排序之前友酱，首先要進(jìn)行降維晴音。降維后形成細(xì)胞cluster，通過cluster形成的細(xì)胞集缔杉，對細(xì)胞軌跡進(jìn)行訓(xùn)練锤躁。

這里根據(jù)需要，提供兩種思路進(jìn)行分析或详。

由于做擬時(shí)間分析時(shí)系羞，構(gòu)建monocle_cds對象采用的是基因的counts值，因此要對基因的表達(dá)量進(jìn)行標(biāo)準(zhǔn)化處理霸琴。而作者也建議即使seurat做了標(biāo)準(zhǔn)化椒振，在這里還要進(jìn)行標(biāo)準(zhǔn)化。（見 http://www.reibang.com/p/aaa0c768c861）

monocle_cds? <- estimateSizeFactors(monocle_cds?) #計(jì)算SizeFactor

monocle_cds? <- estimateDispersions(monocle_cds?) #計(jì)算Dispersions

1 重新進(jìn)行聚類分析來對細(xì)胞軌跡進(jìn)行無監(jiān)督學(xué)習(xí)分析梧乘。

??????#先繪制主成分對應(yīng)的解釋方差的值可視化圖形澎迎。??????

????plot_pc_variance_explained(monocle_cds?, return_all = F, max_components = 25)

????# 這里的num_dim選擇多少個(gè)主成分來進(jìn)行降維庐杨，需要根據(jù)PCA分布圖進(jìn)行確定，一般是選擇斜率第一次出現(xiàn)最小值時(shí)對應(yīng)的主成分點(diǎn)夹供。

????HSMM_myo <- reduceDimension(my_cds,max_components = 2,norm_method = 'log',num_dim = 3,reduction_method = 'tSNE',verbose = T)

????# number指要將細(xì)胞分成多少個(gè)cluster灵份。

????HSMM_myo <- clusterCells(HSMM_myo, verbose = F,num_clusters = number)

????graph <- paste(prefix,'cell_clusters.pdf',sep='_')

????pdf(graph, w=12, h=8)

????plot_cell_clusters(HSMM_myo, color_by = 'as.factor(Cluster)')

????dev.off()

2 根據(jù)seurat已經(jīng)聚好的類對細(xì)胞軌跡進(jìn)行無監(jiān)督學(xué)習(xí)分析。

????# 如果導(dǎo)入的是Seurat3版本的rds文件哮洽，需在pData(HSMM_myo)矩陣中添加一列Cluster填渠，即細(xì)胞群編號。

????HSMM_myo <- moocle_cds

????pData(HSMM_myo)$Cluster <- pData(HSMM_myo)$seurat_clusters

????clustering_DEG_genes <-row.names(subset(fData(HSMM_myo),num_cells_expressed >= 10))?

3 細(xì)胞排序

? ??這一步很慢袁铐，可調(diào)整所使用的的線程cores加快運(yùn)行速度揭蜒。

????diff_test_res <- differentialGeneTest(HSMM_myo,fullModelFormulaStr = "~Cluster"，,cores = 4)

????# 選擇用來排序的基因剔桨。官網(wǎng)給的是q-value值在top1000的基因，此處也可以進(jìn)行調(diào)整徙融。

????HSMM_ordering_genes <-row.names(clustering_DEG_genes)[order(clustering_DEG_genes$qval)][1:1000]

????# 調(diào)整用來排序的基因洒缀，選擇qval< 1e-4。

????HSMM_ordering_genes <- row.names (subset(clustering_DEG_genes, qval < 1e-4))

????# 對細(xì)胞進(jìn)行排序

????HSMM_myo <-setOrderingFilter(HSMM_myo,ordering_genes = HSMM_ordering_genes)

????HSMM_myo <-reduceDimension(HSMM_myo, method = 'DDRTree')

????HSMM_myo <-orderCells(HSMM_myo)

4 對排序好的結(jié)果進(jìn)行可視化輸出

????Cluster_num<-ceiling(length(unique(pData(HSMM_myo)$Cluster))/4)

????pData(HSMM_myo)$seurat_clusters <- as.factor(pData(HSMM_myo)$seurat_clusters)

????graph <- paste(prefix,'cell_trajectory_cluster.pdf',sep='_')

????pdf(graph, w=12, h=8)

????#按state對細(xì)胞進(jìn)行著色

????plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75)????

????plot_cell_trajectory(HSMM_myo, color_by = "State",cell_size = 0.75)+facet_wrap(~State, nrow = Cluster_num)

????#按擬時(shí)間值對細(xì)胞進(jìn)行著色

????plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime",cell_size = 0.75)

????#按cluster id進(jìn)行著色

????plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters",cell_size = 0.75)

????plot_cell_trajectory(HSMM_myo, color_by = "seurat_clusters",cell_size = 0.75)+facet_wrap(~seurat_clusters, nrow =Cluster_num)

????#按樣本進(jìn)行著色

????plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75)

????plot_cell_trajectory(HSMM_myo, color_by = "stim",cell_size = 0.75)+facet_wrap(~stim, ncol =opt$sample_number)

????dev.off()

參考文章見：

1?https://www.jieandze1314.com/post/cnposts/scrna-notes-18/

2?http://www.reibang.com/p/aaa0c768c861

3?http://www.reibang.com/p/2b9982cb7d89