|
文章來(lái)源于：sci666 概要  本文主要討論Seurat對(duì)象導(dǎo)入到Monocle中直接進(jìn)行分析的可行性,，分兩種情況： 以下先進(jìn)行Monocle包的簡(jiǎn)單介紹，再分這兩種情況進(jìn)行嘗試,。 為什么要分這兩種情況進(jìn)行嘗試,？ 1. Seurat包中也有將數(shù)據(jù)標(biāo)準(zhǔn)化的步驟，作者的建議是在Monocle中要再次進(jìn)行標(biāo)準(zhǔn)化,，但是他自己也沒(méi)有嘗試過(guò),，所以不確定會(huì)怎么樣。 2.Seurat包中有個(gè)ScaleData的命令,，目的是去除測(cè)序產(chǎn)生的批次效應(yīng)和技術(shù)噪音,，但對(duì)于我們的數(shù)據(jù)（按不同時(shí)間缺血處理的脾臟，根據(jù)錐蟲(chóng)感染小鼠的時(shí)間進(jìn)行測(cè)序）,，我們要觀察的就是這些不同時(shí)間批次之間的差別,，有可能這個(gè)命令會(huì)將這個(gè)差別掩蓋了。因此如果直接輸入已經(jīng)聚類好的Seurat對(duì)象,，也許會(huì)出現(xiàn)問(wèn)題,。
關(guān)于Monocle 
http://cole-trapnell-lab./monocle-release/ Introduction Monocle介紹了使用RNA-Seq進(jìn)行單細(xì)胞軌跡分析的策略，能夠?qū)⒓?xì)胞按照模擬的時(shí)間順序進(jìn)行排列,，顯示它們的發(fā)展軌跡如細(xì)胞分化等生物學(xué)過(guò)程,。Monocle從數(shù)據(jù)中用無(wú)監(jiān)督或半監(jiān)督學(xué)習(xí)獲得這個(gè)軌跡。 無(wú)監(jiān)督：利用Monocle的自己一套工具或Seurat生成一個(gè)基因列表 Monocle不是通過(guò)實(shí)驗(yàn)將細(xì)胞純化為離散狀態(tài),，而是使用算法來(lái)學(xué)習(xí)每個(gè)細(xì)胞必須經(jīng)歷的基因表達(dá)變化的序列，作為動(dòng)態(tài)生物過(guò)程的一部分,。一旦它了解了基因表達(dá)變化的整體“軌跡”,，Monocle就可以將每個(gè)細(xì)胞放置在軌跡中的適當(dāng)位置。然后,，可以使用Monocle的差異分析工具包來(lái)查找在軌跡過(guò)程中受到調(diào)控的基因,。如果該過(guò)程有多個(gè)結(jié)果，Monocle將重建“分支”軌跡,。這些分支對(duì)應(yīng)于細(xì)胞“決策”,，Monocle提供了強(qiáng)大的工具來(lái)識(shí)別受其影響的基因并參與制作它們。網(wǎng)站中也提供了分析分支的方法,。Monocle依靠Reversed Graph Embedding的機(jī)器學(xué)習(xí)技術(shù)來(lái)構(gòu)建單細(xì)胞軌跡,。 除了構(gòu)建單細(xì)胞軌跡之外，它還能夠做差異表達(dá)分析和聚類來(lái)揭示重要的基因和細(xì)胞,。這與Seurat的功能相似,。 Workflow以及與Seurat的異同 ①創(chuàng)建CellDataSet對(duì)象（下簡(jiǎn)稱CDS對(duì)象） 若要?jiǎng)?chuàng)建新的CDS對(duì)象,，則需要整理出3個(gè)輸入文件（基因-細(xì)胞表達(dá)矩陣,、細(xì)胞-細(xì)胞特征注釋矩陣、基因-基因特征注釋矩陣）,，但方便的是,，Monocle支持從Seurat中直接導(dǎo)入對(duì)象,，通過(guò) 在創(chuàng)建之后,，也會(huì)進(jìn)行數(shù)據(jù)過(guò)濾和標(biāo)準(zhǔn)化,，不同的是Seurat是基于作圖的方法去篩選掉異常的細(xì)胞，而Monocle的過(guò)濾方法則是根據(jù)表達(dá)數(shù)據(jù)的數(shù)學(xué)關(guān)系來(lái)實(shí)現(xiàn),。 其中 ②通過(guò)已知的Marker基因分類細(xì)胞 方法一：通過(guò)一些現(xiàn)有的生物/醫(yī)學(xué)知識(shí)手動(dòng)賦予它們細(xì)胞名，將這些細(xì)胞分為幾大類,，然后再聚類細(xì)胞,。 方法二：與Seurat包的分類細(xì)胞方法類似，篩選出差異表達(dá)基因用于聚類,，然后再根據(jù)每一個(gè)cluster的marker賦予細(xì)胞名,。 ③聚類細(xì)胞 采用的也是t-SNE算法。 ④將細(xì)胞按照偽時(shí)間的順序排列在軌跡上 Step1：選擇輸入基因用于機(jī)器學(xué)習(xí) ⑤差異表達(dá)分析 還沒(méi)細(xì)看
情況①：經(jīng)過(guò)清洗,、標(biāo)準(zhǔn)化和聚類的Seurat對(duì)象導(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的編號(hào),，將此列名更改為“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個(gè)cluster。 計(jì)算SizeFactor和Dispersions
計(jì)算用于方便之后的分析使用,。
clustered_spleen_monocle <- estimateSizeFactors(clustered_spleen_monocle)
跳過(guò)-數(shù)據(jù)清洗
由于此數(shù)據(jù)已經(jīng)經(jīng)過(guò)數(shù)據(jù)清洗,，因此沒(méi)有必要在Monocle中再一次處理。
數(shù)據(jù)標(biāo)準(zhǔn)化 根據(jù)作者的建議,，即使在Seurat包中已經(jīng)標(biāo)準(zhǔn)化處理過(guò)的數(shù)據(jù),，在轉(zhuǎn)化到Monocle中時(shí)仍然需要再一次進(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,])) #將每個(gè)基因都標(biāo)準(zhǔn)化,，melt方便作圖 library(reshape) #作圖,，查看標(biāo)準(zhǔn)化的基因表達(dá)值的分布 qplot(value, geom = “density”, data = melted_dens_df) +
 利用細(xì)胞Marker分類細(xì)胞
首先，因?yàn)橐环N細(xì)胞下又可以細(xì)分成更加小的類別,，因此在用marker給細(xì)胞類時(shí),，就要考慮到它們的對(duì)應(yīng)關(guān)系。如CD4基因所對(duì)應(yīng)的細(xì)胞為CD4+ T細(xì)胞,，而CD4+T細(xì)胞屬于T細(xì)胞中的一種,，我們就要告訴Monocle CD4+T細(xì)胞屬于T細(xì)胞的子集，讓它不要再分類的過(guò)程中將它們分成兩類,。Monocle提供了一個(gè)將細(xì)胞分級(jí)的函數(shù)
cth <- newCellTypeHierarchy()
將marker和細(xì)胞對(duì)應(yīng)起來(lái),，以及排列好它們的從屬關(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)
 偽時(shí)間分析
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ì)胞按照偽時(shí)間排序
clustered_spleen_monocle <- orderCells(clustered_spleen_monocle)
查看能用于上色區(qū)分的變量：
colnames(pData(clustered_spleen_monocle)) [6] “Size_Factor””CellType””Pseudotime” “State”
其實(shí)這一步按照正常情況下來(lái)說(shuō)，是應(yīng)該最好有一個(gè)時(shí)間的變量（如Hour或者Time）,，用來(lái)區(qū)別不同時(shí)間批次處理產(chǎn)生的數(shù)據(jù),，子亮部分的數(shù)據(jù)就可以根據(jù)不同的寄生時(shí)間來(lái)給細(xì)胞上色，從而觀察隨著寄生時(shí)間的變化細(xì)胞狀態(tài)（發(fā)育/分化）的變化,。而像脾臟數(shù)據(jù)這一種,，雖然按照了4個(gè)時(shí)間點(diǎn)進(jìn)行了處理，但是數(shù)據(jù)并沒(méi)有按照不同時(shí)間點(diǎn)區(qū)分出來(lái),，因此只能夠根據(jù)細(xì)胞的分化的過(guò)程來(lái)確定哪些是原始狀態(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”)
 
這是一種樹(shù)狀圖，有三條細(xì)胞軌跡，表示細(xì)胞狀態(tài)主要分為3個(gè)階段,，中間的數(shù)字1表示一個(gè)分叉,。 上圖的細(xì)胞依據(jù)不同的Cluster進(jìn)行上色，根據(jù)之前的Seurat聚類分析,，Cluster5（淺藍(lán)色）對(duì)應(yīng)的是中性粒細(xì)胞,，在此圖位于上面那一分支的頂端；Cluster0（紅色）對(duì)應(yīng)的是B細(xì)胞,，主要位于右邊一支的最頂端,；左下角頂端的藍(lán)色有可能是NK細(xì)胞，但不確定,。這樣看來(lái)比較適合做初始狀態(tài)的是右邊那一支。與下圖對(duì)比得到的結(jié)果也差不多,。
plot_cell_trajectory(clustered_spleen_monocle, color_by = “CellType”) + facet_wrap(~CellType, nrow = 1)
上圖將每個(gè)細(xì)胞的分布軌跡分開(kāi)顯示,，比較明顯地看到B細(xì)胞集中的位置在右支頂端，然后集中為T細(xì)胞,，中間摻雜一些中性粒細(xì)胞（也有可能沒(méi)分好）,。但是大部分的細(xì)胞還是沒(méi)有分出來(lái)，這個(gè)結(jié)果還需要再處理一下,。 由于Monocle不能分辨哪一條軌跡才是“樹(shù)根”,，也就是不知道哪一個(gè)細(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)過(guò)標(biāo)準(zhǔn)化處理的Seurat對(duì)象導(dǎo)入 創(chuàng)建CellDataSet對(duì)象（下簡(jiǎn)稱CDS對(duì)象） 創(chuàng)建Seurat對(duì)象spleen_monocle，先去除一些測(cè)序質(zhì)量差的細(xì)胞：
spleen_monocle <– CreateSeuratObject(raw.data = spleen.data, min.cells = 3, min.genes = 200, project = “10X_spleen”) 從Monocle中導(dǎo)入Seurat對(duì)象
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個(gè)基因，1959個(gè)細(xì)胞,，與之前創(chuàng)建的的Seurat對(duì)象相符,。 計(jì)算SizeFactor和Dispersions計(jì)算用于方便之后的分析使用。
spleen_monocle <- estimateSizeFactors(spleen_monocle) 數(shù)據(jù)清洗 根據(jù)前文所說(shuō)的表達(dá)式,，可以利用nUMI值進(jìn)行過(guò)濾
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
觀察到這里多了一個(gè)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: |
|
|
來(lái)自： 創(chuàng)客小組 > 《sci666》
