|
文章來源于：sci666 概要  本文主要討論Seurat對象導(dǎo)入到Monocle中直接進(jìn)行分析的可行性，分兩種情況： 以下先進(jìn)行Monocle包的簡單介紹，再分這兩種情況進(jìn)行嘗試。 為什么要分這兩種情況進(jìn)行嘗試？ 1. Seurat包中也有將數(shù)據(jù)標(biāo)準(zhǔn)化的步驟，作者的建議是在Monocle中要再次進(jìn)行標(biāo)準(zhǔn)化，但是他自己也沒有嘗試過，所以不確定會怎么樣。 2.Seurat包中有個ScaleData的命令，目的是去除測序產(chǎn)生的批次效應(yīng)和技術(shù)噪音，但對于我們的數(shù)據(jù)（按不同時間缺血處理的脾臟，根據(jù)錐蟲感染小鼠的時間進(jìn)行測序），我們要觀察的就是這些不同時間批次之間的差別，有可能這個命令會將這個差別掩蓋了。因此如果直接輸入已經(jīng)聚類好的Seurat對象，也許會出現(xiàn)問題。
關(guān)于Monocle 
http://cole-trapnell-lab./monocle-release/ Introduction Monocle介紹了使用RNA-Seq進(jìn)行單細(xì)胞軌跡分析的策略，能夠?qū)⒓?xì)胞按照模擬的時間順序進(jìn)行排列，顯示它們的發(fā)展軌跡如細(xì)胞分化等生物學(xué)過程。Monocle從數(shù)據(jù)中用無監(jiān)督或半監(jiān)督學(xué)習(xí)獲得這個軌跡。 無監(jiān)督：利用Monocle的自己一套工具或Seurat生成一個基因列表 Monocle不是通過實驗將細(xì)胞純化為離散狀態(tài)，而是使用算法來學(xué)習(xí)每個細(xì)胞必須經(jīng)歷的基因表達(dá)變化的序列，作為動態(tài)生物過程的一部分。一旦它了解了基因表達(dá)變化的整體“軌跡”，Monocle就可以將每個細(xì)胞放置在軌跡中的適當(dāng)位置。然后，可以使用Monocle的差異分析工具包來查找在軌跡過程中受到調(diào)控的基因。如果該過程有多個結(jié)果，Monocle將重建“分支”軌跡。這些分支對應(yīng)于細(xì)胞“決策”，Monocle提供了強(qiáng)大的工具來識別受其影響的基因并參與制作它們。網(wǎng)站中也提供了分析分支的方法。Monocle依靠Reversed Graph Embedding的機(jī)器學(xué)習(xí)技術(shù)來構(gòu)建單細(xì)胞軌跡。 除了構(gòu)建單細(xì)胞軌跡之外，它還能夠做差異表達(dá)分析和聚類來揭示重要的基因和細(xì)胞。這與Seurat的功能相似。 Workflow以及與Seurat的異同 ①創(chuàng)建CellDataSet對象（下簡稱CDS對象） 若要創(chuàng)建新的CDS對象，則需要整理出3個輸入文件（基因-細(xì)胞表達(dá)矩陣、細(xì)胞-細(xì)胞特征注釋矩陣、基因-基因特征注釋矩陣），但方便的是，Monocle支持從Seurat中直接導(dǎo)入對象，通過 在創(chuàng)建之后，也會進(jìn)行數(shù)據(jù)過濾和標(biāo)準(zhǔn)化，不同的是Seurat是基于作圖的方法去篩選掉異常的細(xì)胞，而Monocle的過濾方法則是根據(jù)表達(dá)數(shù)據(jù)的數(shù)學(xué)關(guān)系來實現(xiàn)。 其中 ②通過已知的Marker基因分類細(xì)胞 方法一：通過一些現(xiàn)有的生物/醫(yī)學(xué)知識手動賦予它們細(xì)胞名，將這些細(xì)胞分為幾大類，然后再聚類細(xì)胞。 方法二：與Seurat包的分類細(xì)胞方法類似，篩選出差異表達(dá)基因用于聚類，然后再根據(jù)每一個cluster的marker賦予細(xì)胞名。 ③聚類細(xì)胞 采用的也是t-SNE算法。 ④將細(xì)胞按照偽時間的順序排列在軌跡上 Step1：選擇輸入基因用于機(jī)器學(xué)習(xí) ⑤差異表達(dá)分析 還沒細(xì)看
情況①：經(jīng)過清洗、標(biāo)準(zhǔn)化和聚類的Seurat對象導(dǎo)入 spleen clustered_spleen_monocle <- importCDS(spleen, import_all = TRUE) head(pData(clustered_spleen_monocle)) nGene nUMI orig.ident percent.mito res.0.6 Size_Factor AAACCTGAGGTGATAT 1072 3999 10X_spleen 0.01150863 3 NA AAACCTGAGTCATCCA 1562 4640 10X_spleen 0.03709295 5 NA AAACCTGCAGTGAGTG 995 3489 10X_spleen 0.03955288 3 NA AAACGGGAGACTGGGT1704 7397 10X_spleen 0.02852508 0 NA AAACGGGAGACTTGAA 976 3102 10X_spleen 0.04932302 2 NA AAACGGGAGATAGGAG1236 4548 10X_spleen 0.02682498 1 NA res.0.6為cluster的編號，將此列名更改為“Cluster”，方便后面使用。 names(pData(clustered_spleen_monocle))[names(pData(clustered_spleen_monocle))== “res.0.6”]=“Cluster” range(pData(clustered_spleen_monocle)$Cluster)
確認(rèn)此列的范圍為0到8，也就是共9個cluster。 計算SizeFactor和Dispersions
計算用于方便之后的分析使用。
clustered_spleen_monocle <- estimateSizeFactors(clustered_spleen_monocle)
跳過-數(shù)據(jù)清洗
由于此數(shù)據(jù)已經(jīng)經(jīng)過數(shù)據(jù)清洗，因此沒有必要在Monocle中再一次處理。
數(shù)據(jù)標(biāo)準(zhǔn)化 根據(jù)作者的建議，即使在Seurat包中已經(jīng)標(biāo)準(zhǔn)化處理過的數(shù)據(jù)，在轉(zhuǎn)化到Monocle中時仍然需要再一次進(jìn)行標(biāo)準(zhǔn)化。#將表達(dá)矩陣中所有值進(jìn)行l(wèi)og標(biāo)準(zhǔn)化 L <-log(exprs(clustered_spleen_monocle[expressed_genes,])) #將每個基因都標(biāo)準(zhǔn)化，melt方便作圖 library(reshape) #作圖，查看標(biāo)準(zhǔn)化的基因表達(dá)值的分布 qplot(value, geom = “density”, data = melted_dens_df) +
 利用細(xì)胞Marker分類細(xì)胞
首先，因為一種細(xì)胞下又可以細(xì)分成更加小的類別，因此在用marker給細(xì)胞類時，就要考慮到它們的對應(yīng)關(guān)系。如CD4基因所對應(yīng)的細(xì)胞為CD4+ T細(xì)胞，而CD4+T細(xì)胞屬于T細(xì)胞中的一種，我們就要告訴Monocle CD4+T細(xì)胞屬于T細(xì)胞的子集，讓它不要再分類的過程中將它們分成兩類。Monocle提供了一個將細(xì)胞分級的函數(shù)
cth <- newCellTypeHierarchy()
將marker和細(xì)胞對應(yīng)起來，以及排列好它們的從屬關(guān)系。
cth <- addCellType(cth, “T cell”, cth <- addCellType(cth, “CD4+ T cell”, cth <- addCellType(cth, “CD8+ T cell”, cth <- addCellType(cth, “B cell”, cth <- addCellType(cth, “Monocyte”, cth <- addCellType(cth, “Neutrophil”,
下面這一步將細(xì)胞進(jìn)行分類。
clustered_spleen_monocle <- classifyCells(clustered_spleen_monocle, cth, 0.1)
查看細(xì)胞分類的情況。
table(pData(clustered_spleen_monocle)$CellType)
 偽時間分析
Step1: 選擇定義細(xì)胞發(fā)展的基因
diff_test_res <- differentialGeneTest(clustered_spleen_monocle,fullModelFormulaStr = “~Cluster”) ordering_genes <- row.names (subset(diff_test_res, qval < 0.01)) Step2: 數(shù)據(jù)集降維clustered_spleen_monocle <- setOrderingFilter(clustered_spleen_monocle, ordering_genes)
plot_ordering_genes(clustered_spleen_monocle) 
clustered_spleen_monocle <– reduceDimension(clustered_spleen_monocle, max_components = 2, reduction_method = “DDRTree”)
Step3: 將細(xì)胞按照偽時間排序
clustered_spleen_monocle <- orderCells(clustered_spleen_monocle)
查看能用于上色區(qū)分的變量：
colnames(pData(clustered_spleen_monocle)) [6] “Size_Factor””CellType””Pseudotime” “State”
其實這一步按照正常情況下來說，是應(yīng)該最好有一個時間的變量（如Hour或者Time），用來區(qū)別不同時間批次處理產(chǎn)生的數(shù)據(jù)，子亮部分的數(shù)據(jù)就可以根據(jù)不同的寄生時間來給細(xì)胞上色，從而觀察隨著寄生時間的變化細(xì)胞狀態(tài)（發(fā)育/分化）的變化。而像脾臟數(shù)據(jù)這一種，雖然按照了4個時間點進(jìn)行了處理，但是數(shù)據(jù)并沒有按照不同時間點區(qū)分出來，因此只能夠根據(jù)細(xì)胞的分化的過程來確定哪些是原始狀態(tài)。
plot_cell_trajectory(clustered_spleen_monocle, color_by = “Cluster”) plot_cell_trajectory(clustered_spleen_monocle, color_by = “CellType”) plot_cell_trajectory(clustered_spleen_monocle, color_by = “State”)
 
這是一種樹狀圖，有三條細(xì)胞軌跡，表示細(xì)胞狀態(tài)主要分為3個階段，中間的數(shù)字1表示一個分叉。 上圖的細(xì)胞依據(jù)不同的Cluster進(jìn)行上色，根據(jù)之前的Seurat聚類分析，Cluster5（淺藍(lán)色）對應(yīng)的是中性粒細(xì)胞，在此圖位于上面那一分支的頂端；Cluster0（紅色）對應(yīng)的是B細(xì)胞，主要位于右邊一支的最頂端；左下角頂端的藍(lán)色有可能是NK細(xì)胞，但不確定。這樣看來比較適合做初始狀態(tài)的是右邊那一支。與下圖對比得到的結(jié)果也差不多。
plot_cell_trajectory(clustered_spleen_monocle, color_by = “CellType”) + facet_wrap(~CellType, nrow = 1)
上圖將每個細(xì)胞的分布軌跡分開顯示，比較明顯地看到B細(xì)胞集中的位置在右支頂端，然后集中為T細(xì)胞，中間摻雜一些中性粒細(xì)胞（也有可能沒分好）。但是大部分的細(xì)胞還是沒有分出來，這個結(jié)果還需要再處理一下。 由于Monocle不能分辨哪一條軌跡才是“樹根”，也就是不知道哪一個細(xì)胞狀態(tài)才是更初始的，可設(shè)定
head(pData(clustered_spleen_monocle)) nGene nUMI orig.ident percent.mito Cluster Size_Factor CellType AAACCTGAGGTGATAT 1072 3999 10X_spleen 0.01150863 3 0.8327676 T cell AAACCTGAGTCATCCA 1562 4640 10X_spleen 0.03709295 5 0.9661104 Ambiguous AAACCTGCAGTGAGTG 995 3489 10X_spleen 0.03955288 3 0.7269267 Monocyte AAACGGGAGACTGGGT1704 7397 10X_spleen 0.02852508 0 1.5411513 Ambiguous AAACGGGAGACTTGAA 976 3102 10X_spleen 0.04932302 2 0.6462960 B cell AAACGGGAGATAGGAG1236 4548 10X_spleen 0.02682498 1 0.9475674 Ambiguous 情況②：未經(jīng)過標(biāo)準(zhǔn)化處理的Seurat對象導(dǎo)入 創(chuàng)建CellDataSet對象（下簡稱CDS對象） 創(chuàng)建Seurat對象spleen_monocle，先去除一些測序質(zhì)量差的細(xì)胞：
spleen_monocle <– CreateSeuratObject(raw.data = spleen.data, min.cells = 3, min.genes = 200, project = “10X_spleen”) 從Monocle中導(dǎo)入Seurat對象
library(monocle) 查看數(shù)據(jù)：
spleen_monocle Annotation: head(fData(spleen_monocle)) gene_short_name RP11-34P13.7 RP11-34P13.7 FO538757.2 FO538757.2 AP006222.2 AP006222.2 RP4-669L17.10 RP4-669L17.10 RP11-206L10.9 RP11-206L10.9 LINC00115 LINC00115 head(pData(spleen_monocle)) nGene nUMI orig.ident Size_Factor AAACCTGAGGTGATAT 1072 3999 10X_spleen NA AAACCTGAGTCATCCA1562 4640 10X_spleen NA AAACCTGCAGTGAGTG 995 3489 10X_spleen NA AAACGGGAGACTGGGT 1704 7397 10X_spleen NA AAACGGGAGACTTGAA 976 3102 10X_spleen NA AAACGGGAGATAGGAG 1236 4548 10X_spleen NA 15655個基因，1959個細(xì)胞，與之前創(chuàng)建的的Seurat對象相符。 計算SizeFactor和Dispersions計算用于方便之后的分析使用。
spleen_monocle <- estimateSizeFactors(spleen_monocle) 數(shù)據(jù)清洗 根據(jù)前文所說的表達(dá)式，可以利用nUMI值進(jìn)行過濾
head(pData(spleen_monocle)) nGene nUMI orig.ident Size_Factor num_genes_expressed AAACCTGAGGTGATAT 1072 3999 10X_spleen 0.8207242 1070 AAACCTGAGTCATCCA1562 4640 10X_spleen 0.9521386 1560 AAACCTGCAGTGAGTG 995 3489 10X_spleen 0.7164140 995 AAACGGGAGACTGGGT1704 7397 10X_spleen 1.5188634 1704 AAACGGGAGACTTGAA 976 3102 10X_spleen 0.6369493 976 AAACGGGAGATAGGAG1236 4548 10X_spleen 0.9338638 1236
觀察到這里多了一個Size_Factor的列
#設(shè)置上下限： upper_bound_spleen <- 10^(mean(log10(pData(spleen_monocle)$nUMI)) #可視化 qplot(nUMI,data = pData(spleen_monocle), geom = “density”) + geom_vline(xintercept = lower_bound_spleen) + geom_vline(xintercept = upper_bound_spleen)
qplot_spleenmoncle_filter.jpeg
留下兩條豎線中間的部分：
spleen_monocle <- spleen_monocle[,pData(spleen_monocle)$nUMI > lower_bound_spleen & pData(spleen_monocle)$nUMI < upper_bound_spleen]
spleen_monocle CellDataSet (storageMode: environment) Annotation: |
|
|
