梗概
參照官網(wǎng)教程：http://cole-trapnell-lab.github.io/monocle-release/docs/#getting-started-with-monocle
或在R中運行以下命令：
browseVignettes("monocle")
關(guān)于Monocle
Monocle introduced the strategy of ordering single cells in pseudotime,
placing them along a trajectory corresponding to a biological process such as cell differentiation. Monocle learns this trajectory directly from the data, in either a fully unsupervised or a semi-supervised manner. It also performs differential gene expression and clustering to identify important genes and cell states.
Monocle是用于分析single-cell RNA-seq的R包的烁，能夠?qū)⒓?xì)胞按照模擬的時間順序進(jìn)行排列诈闺，顯示它們的發(fā)展軌跡如細(xì)胞分化等生物學(xué)過程买雾。Monocle從數(shù)據(jù)中用無監(jiān)督或半監(jiān)督（？）的方式直接獲得這個軌跡嗤军。除此之外晃危，它還能夠做差異表達(dá)分析和聚類來揭示重要的基因和細(xì)胞。
在發(fā)育過程中震叮，為了響應(yīng)刺激和生命苇瓣，細(xì)胞從一種功能性“狀態(tài)”轉(zhuǎn)變?yōu)榱硪环N功能性“狀態(tài)”偿乖。處于不同狀態(tài)的細(xì)胞表達(dá)不同的基因組哲嘲，產(chǎn)生蛋白質(zhì)和代謝物的動態(tài)重復(fù)序列眠副，從而完成它們的工作竣稽。當(dāng)細(xì)胞在狀態(tài)之間移動時，經(jīng)歷轉(zhuǎn)錄重組的過程娃弓，一些基因被沉默岛宦，另一些被新激活恋博。這些瞬態(tài)通常難以表征私恬，因為在更穩(wěn)定的端點狀態(tài)之間純化細(xì)胞可能是困難的或不可能的。單細(xì)胞RNA-Seq可以讓您在不需要純化的情況下查看這些狀態(tài)疫衩。但是闷煤，為此涮瞻，我們必須確定每個單元格的可能狀態(tài)范圍。
Monocle介紹了使用RNA-Seq進(jìn)行單細(xì)胞軌跡分析的策略近顷。Monocle不是通過實驗將細(xì)胞純化為離散狀態(tài)宁否，而是使用算法來學(xué)習(xí)每個細(xì)胞必須經(jīng)歷的基因表達(dá)變化的序列慕匠，作為動態(tài)生物過程的一部分。一旦它了解了基因表達(dá)變化的整體“軌跡”蓉媳，Monocle就可以將每個細(xì)胞放置在軌跡中的適當(dāng)位置。然后殴瘦，您可以使用Monocle的差異分析工具包來查找在軌跡過程中受到調(diào)控的基因号杠，如“ 查找作為假性函數(shù)的變化的基因 ”一節(jié)中所述姨蟋。如果該過程有多個結(jié)果，Monocle將重建“分支”軌跡悠砚。這些分支對應(yīng)于細(xì)胞“決策”堂飞，Monocle提供了強大的工具來識別受其影響的基因并參與制作它們绰筛。您可以在“分析單細(xì)胞軌跡中的分支”一節(jié)中查看如何分析分支。Monocle依賴于稱為反向圖嵌入的機器學(xué)習(xí)技術(shù)來構(gòu)建單細(xì)胞軌跡衡蚂。您可以 在Monocle理論背后的部分中閱讀更多關(guān)于Monocle方法的理論基礎(chǔ)骏庸，或參考小插圖末尾顯示的參考文獻(xiàn)
什么是偽時間具被？
Pseudotime是衡量單個細(xì)胞通過細(xì)胞分化等過程取得多少進(jìn)展的量度。在許多生物過程中补箍，細(xì)胞不能完美同步啸蜜。在諸如細(xì)胞分化的過程的單細(xì)胞表達(dá)研究中衬横，捕獲的細(xì)胞可以在進(jìn)展方面廣泛分布。也就是說遥诉，在完全同時捕獲的細(xì)胞群中，一些細(xì)胞可能很遠(yuǎn)霉翔，而其他細(xì)胞甚至可能尚未開始該過程苞笨。當(dāng)您想要了解細(xì)胞從一個狀態(tài)轉(zhuǎn)換到另一個狀態(tài)時發(fā)生的調(diào)節(jié)變化的順序時瀑凝，這種不同步會產(chǎn)生重大問題。跟蹤同時捕獲的細(xì)胞的表達(dá)產(chǎn)生非常壓縮的基因動力學(xué)感覺谚中，以及該基因的明顯變異性' 的表達(dá)會非常高宪塔。通過沿著學(xué)習(xí)軌跡根據(jù)其進(jìn)度對每個單元進(jìn)行排序囊拜，Monocle減輕了由于不同步引起的問題。Monocle不是根據(jù)時間跟蹤表達(dá)式的變化，而是根據(jù)軌跡上的進(jìn)度跟蹤變化蔽莱，我們稱之為pseudotime'戚长。Pseudotime是一個抽象的進(jìn)展單元：它只是一個單元格與軌跡起點之間的距離同廉，沿著最短路徑測量。軌跡的總長度是根據(jù)細(xì)胞從起始狀態(tài)移動到最終狀態(tài)時經(jīng)歷的轉(zhuǎn)錄變化的總量來定義的锅劝。
實驗記錄
偽時間分析實驗技術(shù)流程主要包括三個步驟故爵，而每一步都是一個機器學(xué)習(xí)任務(wù)：
step1：選擇輸入基因用于機器學(xué)習(xí)
這個過程稱為feature selection（特征選擇）隅津，這些基因?qū)壽E的形狀有著最重要的影響劲室。我們應(yīng)該要選擇的是最能反映細(xì)胞狀態(tài)的基因很洋。
如果直接通過Seurat輸出的一些重要的基因（比如每個cluster中的高FC值基因）作為輸入對象的話就能夠?qū)崿F(xiàn)一個“無監(jiān)督”分析隧枫∮贫猓或者也可以利用生物學(xué)知識手動選擇一些重要的基因進(jìn)行“半監(jiān)督”分析。
step2：數(shù)據(jù)降維
利用Reversed graph embedding算法將數(shù)據(jù)降維斤讥。沒有太懂湾趾。
相對于PCA來說這個算法更能夠反應(yīng)數(shù)據(jù)的內(nèi)部結(jié)構(gòu)（據(jù)說是這樣）
step3：將細(xì)胞按照偽時間排序
這個過程稱為manifold learning（流形學(xué)習(xí)）搀缠。Monocle利用軌跡來描述細(xì)胞如何從一個狀態(tài)轉(zhuǎn)換到另一個狀態(tài)。得到的是一個樹狀圖簸州，樹的根部或樹干表示的是細(xì)胞最初的狀態(tài)（如果有的話）歧譬，隨著細(xì)胞的變化就沿著分叉樹枝發(fā)展瑰步。一個細(xì)胞的偽時間值（pseudotime value）為它的位置沿著樹枝到根部的距離。
以下是正式實驗記錄：
step0：數(shù)據(jù)輸入
下載與安裝
實驗要求——
R version 3.4 or higher, Bioconductor version 3.5, and monocle 2.4.0 or higher to have access to the latest features.
查看R版本：
sessionInfo()
R version 3.4.4 (2018-03-15)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: OS X El Capitan 10.11.6
安裝bioconductor以及monocle
source("http://bioconductor.org/biocLite.R")
biocLite("monocle")
加載程輯包：
library(monocle)
將Seurat Object導(dǎo)入
Monocle包可以利用importCDS()命令將Seurat包中的數(shù)據(jù)導(dǎo)入，同時也能用exportCDS()將數(shù)據(jù)導(dǎo)出到其他包中袁滥。
由于之前的Seurat變量已經(jīng)被覆蓋题翻，故重新創(chuàng)建一個。
spleen_monocle <- CreateSeuratObject(raw.data = spleen.data, min.cells = 3, min.genes = 200, project = "10X_spleen")
spleen_monocle <- importCDS(spleen_monocle, import_all = TRUE)
查看數(shù)據(jù)：
spleen_monocle
CellDataSet (storageMode: environment)
assayData: 15655 features, 1959 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGAGGTGATAT AAACCTGAGTCATCCA ... TTTGTCAGTCGTCTTC (1959 total)
varLabels: nGene nUMI orig.ident Size_Factor
varMetadata: labelDescription
featureData
featureNames: RP11-34P13.7 FO538757.2 ... AC240274.1 (15655 total)
fvarLabels: gene_short_name
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
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
AAACCTGAGTCATCCA 1562 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對象相符偎窘。
（細(xì)胞名是一些barcode序列）
step1：選擇輸入基因用于機器學(xué)習(xí)
①脾臟數(shù)據(jù)subset的獲取
因為脾臟是在4個時間點經(jīng)過缺血處理的陌知，
spleen_monocle <- detectGenes(spleen_monocle)
查看數(shù)據(jù)，添加了一列num_cells_expressed
head(fData(spleen_monocle))
gene_short_name num_cells_expressed
RP11-34P13.7 RP11-34P13.7 3
FO538757.2 FO538757.2 214
AP006222.2 AP006222.2 157
RP4-669L17.10 RP4-669L17.10 10
RP11-206L10.9 RP11-206L10.9 25
LINC00115 LINC00115 35
print(pData(spleen_monocle))
nGene nUMI orig.ident percent.mito res.0.6 Size_Factor num_genes_expressed
AAACCTGAGGTGATAT 1072 3999 10X_spleen 0.0115086315 3 NA 1070
AAACCTGAGTCATCCA 1562 4640 10X_spleen 0.0370929480 5 NA 1560
AAACCTGCAGTGAGTG 995 3489 10X_spleen 0.0395528805 3 NA 995
AAACGGGAGACTGGGT 1704 7397 10X_spleen 0.0285250777 0 NA 1704
AAACGGGAGACTTGAA 976 3102 10X_spleen 0.0493230174 2 NA 976
AAACGGGAGATAGGAG 1236 4548 10X_spleen 0.0268249780 1 NA 1236
AAACGGGGTCCGAACC 2430 6863 10X_spleen 0.0379119277 6 NA 2425
AAACGGGTCAGGTAAA 1256 4750 10X_spleen 0.0381052632 1 NA 1256
AAACGGGTCGTCACGG 1887 8095 10X_spleen 0.0302730755 5 NA 1885
AAAGATGAGCGATGAC 1473 5729 10X_spleen 0.0185088179 4 NA 1471
…………
step2：數(shù)據(jù)降維
step3：將細(xì)胞按照偽時間排序
(此段內(nèi)容舍棄)將Seurat Object導(dǎo)入
Monocle包可以利用importCDS()命令將Seurat包中的數(shù)據(jù)導(dǎo)入，同時也能用exportCDS()將數(shù)據(jù)導(dǎo)出到其他包中把篓。
spleen_monocle <- importCDS(spleen, import_all = TRUE)
查看數(shù)據(jù)：
spleen_monocle
CellDataSet (storageMode: environment)
assayData: 15655 features, 1940 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGAGGTGATAT AAACCTGAGTCATCCA ... TTTGTCAGTCGTCTTC (1940 total)
varLabels: nGene nUMI ... Size_Factor (6 total)
varMetadata: labelDescription
featureData
featureNames: RP11-34P13.7 FO538757.2 ... AC240274.1 (15655 total)
fvarLabels: gene_short_name
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
15655個基因韧掩，1940個細(xì)胞窖铡，與之前的Seurat對象相符费彼。
選擇數(shù)據(jù)分布類型
需要根據(jù)實驗的數(shù)據(jù)類型選擇數(shù)據(jù)的分布，比如是相對表達(dá)量的數(shù)據(jù)/基于計數(shù)的數(shù)據(jù)（如UMI），官方認(rèn)為：“FPKM/TPM values are generally log-normally distributed, while UMIs or read counts are better modeled with the negative binomial.”即利用負(fù)二項式的模型會更加適合UMI計數(shù)的實驗數(shù)據(jù)虹钮。
spleen_M <- newCellDataSet(spleen_monocle, expressionFamily=negbinomial.size())
報錯：spleen_monocle不是矩陣.
Error in newCellDataSet(spleen_monocle, expressionFamily = negbinomial.size()) :
Error: argument cellData must be a matrix (either sparse from the Matrix package or dense)
In addition: Warning message:
In newCellDataSet(spleen_monocle, expressionFamily = negbinomial.size()) :
Warning: featureData must contain a column verbatim named 'gene_short_name' for certain functions
預(yù)估size factors和分散
spleen_monocle <- estimateSizeFactors(spleen_monocle)
spleen_monocle <- estimateDispersions(spleen_monocle)
removing 9 outliers.
比較奇怪芙粱，在輸入這兩個命令之后顯示去除了9個異常值春畔，但是查看spleen_monocle變量時返回：
CellDataSet (storageMode: environment)
assayData: 15655 features, 1940 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGAGGTGATAT AAACCTGAGTCATCCA ... TTTGTCAGTCGTCTTC (1940 total)
varLabels: nGene nUMI ... Size_Factor (6 total)
varMetadata: labelDescription
featureData
featureNames: RP11-34P13.7 FO538757.2 ... AC240274.1 (15655 total)
fvarLabels: gene_short_name
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
> dim(spleen_monocle)
Features Samples
15655 1940
仍然是15655個基因，1940個細(xì)胞振峻，與之前無異扣孟。因此采用了另一種輸入方式荣赶。
remove(spleen_monocle)
