1.定義

做擬時序分析是為了探索自己感興趣的幾種細(xì)胞之間的發(fā)育關(guān)系钧唐，一般不是用全部類型的細(xì)胞來做的溪烤。以下是進(jìn)行擬時序分析的幾個主要目的：細(xì)胞狀態(tài)推斷县昂、細(xì)胞軌跡追蹤邢疙、細(xì)胞命運(yùn)預(yù)測棍弄、動態(tài)過程建模、細(xì)胞亞群識別疟游、基因調(diào)控網(wǎng)絡(luò)推斷呼畸、疾病機(jī)理探索、藥物作用機(jī)制研究颁虐。
輸入數(shù)據(jù)是seurat做完降維聚類分群注釋的數(shù)據(jù)蛮原，從調(diào)用之前注釋好的數(shù)據(jù)開始

rm(list = ls())
library(Seurat)
library(monocle)
library(dplyr)
load("sce.Rdata")
scRNA = sce
table(Idents(scRNA))
head(scRNA@meta.data)
DimPlot(scRNA,label = T)
scRNA$celltype = Idents(scRNA)  #celltype是細(xì)胞類型注釋，添加到數(shù)據(jù)中

做擬時序分析通常不是拿全部的細(xì)胞另绩，而是拿感興趣的一部分儒陨。用subset提取子集即可。因為要使用差異基因來排序板熊，所以要兩類及以上細(xì)胞框全。基于背景知識選擇有進(jìn)化關(guān)系的細(xì)胞類型干签。

levels(Idents(scRNA)) #打出來細(xì)胞類型供復(fù)制
scRNA = subset(scRNA,idents = c("NK","CD8 T"))
table(Idents(scRNA))
set.seed(1234)      #抽樣津辩，實戰(zhàn)時不能抽樣！！４亍Ｕ⒍取！
scRNA = subset(scRNA,downsample = 100)

因為monocle和seurat是兩個不同的體系蚜印，所以需要將seurat對象轉(zhuǎn)換為monocle可以接受的CellDataSet對象莺禁。雖然monocle3已經(jīng)出來很久了，但大家都不約而同的選擇monocle2窄赋，大概就是習(xí)慣了吧哟冬。。

2.創(chuàng)建CellDataSet對象

以下代碼無需改動

# count矩陣忆绰，官方建議用count
ct <- scRNA@assays$RNA$counts
# 基因注釋
gene_ann <- data.frame(
  gene_short_name = row.names(ct), 
  row.names = row.names(ct)
)
fd <- new("AnnotatedDataFrame",
          data=gene_ann)
# 臨床信息
pd <- new("AnnotatedDataFrame",
          data=scRNA@meta.data)
#新建CellDataSet對象
sc_cds <- newCellDataSet(
  ct, 
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=1)
sc_cds

3.構(gòu)建細(xì)胞發(fā)育軌跡

1.使用seurat給出的高變化基因
2.按照平均表達(dá)量大于某個數(shù)字(比如0.1浩峡，官網(wǎng)用的是這個)的基因。
3.使用不同細(xì)胞類型之間的差異基因错敢，differentialGeneTest計算翰灾。
根據(jù)以往經(jīng)驗，選擇使用最后一個策略效果可能會更好稚茅。
以下代碼無需改動

sc_cds <- estimateSizeFactors(sc_cds)
sc_cds <- estimateDispersions(sc_cds)
fdif = "diff_test_res.Rdata"
if(!file.exists(fdif)){
  diff_test_res <- differentialGeneTest(sc_cds,
                                        fullModelFormulaStr = "~celltype",
                                        cores = 4)
  save(diff_test_res,file = fdif)
}
load(fdif)
ordering_genes <- row.names(subset(diff_test_res, qval < 0.01))
#查看基因纸淮，篩選適合用于排序的，設(shè)置為排序要使用的基因
head(ordering_genes)
length(ordering_genes)
sc_cds <- setOrderingFilter(sc_cds, ordering_genes)
#畫出選擇的基因
plot_ordering_genes(sc_cds)

plot_ordering_genes畫的圖縱坐標(biāo)是基因表達(dá)量的變異性（離散性）,橫坐標(biāo)是每個基因在所有細(xì)胞種的平均表達(dá)量亚享。
紅線表示Monocle2基于橫縱坐標(biāo)的關(guān)系擬合的變異期望值咽块，為用于排序的基因是黑色，其他基因是灰色

image.png

#降維
sc_cds <- reduceDimension(sc_cds)
#細(xì)胞排序
sc_cds <- orderCells(sc_cds)

4.繪圖展示

4.1 發(fā)育軌跡圖

library(ggsci)
p1 = plot_cell_trajectory(sc_cds)+ scale_color_nejm()
p2 = plot_cell_trajectory(sc_cds, color_by = 'Pseudotime') 
p3 = plot_cell_trajectory(sc_cds, color_by = 'celltype')  + scale_color_npg()
library(patchwork)
p2+p1/p3

image.png

4.2 經(jīng)典的擬時序熱圖

展示了一些基因是如何隨著時間軌跡的變化而漸變的，這個漸變不同于findmarkers哲银，是體現(xiàn)變化過程的扛吞，而不是直接給出差異表達(dá)的基因。
理論上應(yīng)該是一句代碼一行荆责，如果兩句代碼想放在同一行上就要用;隔開滥比。可以減少點(diǎn)Run或者按快捷鍵的次數(shù)做院。
選擇了q值最小的50個基因來畫熱圖盲泛。數(shù)量可以按需調(diào)整（head(50)）濒持，其他代碼不用改。

gene_to_cluster = diff_test_res %>% arrange(qval) %>% head(50) %>% pull(gene_short_name);head(gene_to_cluster)
plot_pseudotime_heatmap(sc_cds[gene_to_cluster,],
                        num_clusters = nlevels(Idents(scRNA)),   #是把基因聚成幾簇寺滚，推薦是由幾類細(xì)胞就設(shè)置幾柑营，用nlevels實現(xiàn)了自動化
                        show_rownames = TRUE,    #是把3種顏色切割成漸變的100種顏色
                        cores = 4,return_heatmap = TRUE,
                        hmcols = colorRampPalette(c("navy", "white", "firebrick3"))(100))

image.png

4.3 基因軌跡圖

用感興趣的基因給軌跡圖著色

gs = head(gene_to_cluster)     #感興趣的基因
plot_cell_trajectory(sc_cds,markers=gs,
                     use_color_gradient=T)

image.png

4.4 基因擬時序點(diǎn)圖

plot_genes_in_pseudotime(sc_cds[gs,],
                 color_by = "celltype",
                 nrow= 3, #6個基因所以排了3行，數(shù)量有變化時要改
                 ncol = NULL )

橫坐標(biāo)是按照pseudotime排好順序的村视。

image.png

改成2個B細(xì)胞群："naive B","plasma B"

改細(xì)胞群重新跑的時候官套，創(chuàng)建CellDataSet對象出bug

image.png

解決方法：更新igraph包
ps：另一個裝包error

install.packages("igraph")
Error in install.packages : Updating loaded packages

這種情況通常發(fā)生在某些需要更新的包當(dāng)前正在使用。需要重啟R蚁孔，在什么包都沒加載的時候再更新奶赔。
切換細(xì)胞子集重跑結(jié)果：

image.png

image.png

image.png

image.png

image.png

image.png