前一陣子練習(xí)了幾個(gè)單細(xì)胞測序分析筝闹，然而對于所用的R包并沒有深入的學(xué)習(xí)。作為強(qiáng)迫癥晚期患者划提，想知其然，還想知其所以然邢享，所以打算仔細(xì)的學(xué)習(xí)這些R包的使用鹏往。這篇文章將深入的學(xué)習(xí)Monocle2這個(gè)包，官方網(wǎng)站http://cole-trapnell-lab.github.io/monocle-release/docs/#introduction骇塘，我主要是翻譯一下這個(gè)包的使用說明和注意事項(xiàng)伊履，一般來說如果你想深入了解一個(gè)包，看官網(wǎng)是最靠譜的方法了款违。里面一些“廢話”我就不翻譯了唐瀑，只翻譯必要的。以便以后自己分析數(shù)據(jù)時(shí)使用插爹，同時(shí)也分享給懶得看英文的童鞋們～

Abstract

單細(xì)胞基因表達(dá)研究使得我們可以研究在復(fù)雜的生物過程和高度異質(zhì)性細(xì)胞群里的轉(zhuǎn)錄調(diào)節(jié)哄辣，這些研究有助于基因的發(fā)現(xiàn)，鑒定某些細(xì)胞的類型赠尾，或者標(biāo)記生物過程的中間態(tài)力穗，和兩種細(xì)胞“命運(yùn)”的分支。在許多單細(xì)胞研究中气嫁，細(xì)胞個(gè)體的基因調(diào)節(jié)的執(zhí)行過程并不同步当窗。事實(shí)上，每個(gè)細(xì)胞都是轉(zhuǎn)錄過程中的一個(gè)“瞬間”寸宵。Monocle包提供了一個(gè)工具用來分析單細(xì)胞測序結(jié)果超全。Monocle引入了一個(gè)“偽時(shí)間”排序單個(gè)細(xì)胞的策略，利用細(xì)胞個(gè)體的異步性把細(xì)胞按照一個(gè)軌道排列邓馒，比如說細(xì)胞分化。Monocle利用機(jī)器學(xué)習(xí)的技術(shù)(Reversed Graph Embedding)來對細(xì)胞進(jìn)行排序蛾坯，這個(gè)方法能夠可靠而準(zhǔn)確地解決復(fù)雜的生物學(xué)過程光酣。Monocle也可以對細(xì)胞進(jìn)行聚類分析。Monocle還可以在聚類后進(jìn)行差異分析脉课，可以鑒定在不同階段和不同細(xì)胞類型的差異基因救军。Monocle是專門為單細(xì)胞測序分析設(shè)計(jì)的包，但是也可以用來做其他的分析倘零。Monocle2可以分析非常大的單細(xì)胞測序數(shù)據(jù)（數(shù)萬個(gè)細(xì)胞甚至更多）唱遭。

Introduction

這篇文章主要提供一個(gè)基于Monocle的單細(xì)胞測序分析的流程。
Monocle可以幫助你進(jìn)行三個(gè)主要類型的分析：

• 聚類呈驶，分類拷泽，細(xì)胞計(jì)數(shù)。Monocle可以幫助你鑒定新的細(xì)胞類型。
• 構(gòu)建單細(xì)胞軌跡司致。在發(fā)育拆吆、疾病以及整個(gè)生命過程里，細(xì)胞從一個(gè)狀態(tài)轉(zhuǎn)成另一個(gè)狀態(tài)脂矫。Monocle可以幫助你發(fā)現(xiàn)這些轉(zhuǎn)變枣耀。
• 差異表達(dá)分析。

在使用Monocle前庭再，先來看一下如何安裝Monocle捞奕。

安裝 Monocle

Monocle運(yùn)行在R語言環(huán)境下。你需要R3.4版本或者更高的版本拄轻、Bioconductor V3.5版本颅围， 在這個(gè)網(wǎng)址安裝Bioconductor:

> source("http://bioconductor.org/biocLite.R")
> biocLite()

安裝了Bioconductor，就可以準(zhǔn)備安裝Monocle和所有需要的依賴包了哺眯。

> biocLite("monocle")

檢查你的Monocle是否安裝正確了谷浅，打開你的R，并且打如下代碼：

> library(monocle) #實(shí)際上就是調(diào)用一下

安裝最近的版本：
我們推薦你使用Bioconductor安裝我們最近的版本奶卓。你也可以通過GitHub安裝一疯，在R里輸入如下代碼：

> install.packages("devtools")
> devtools::install_github("cole-trapnell-lab/monocle-release@develop")

有時(shí)我們添加了一些features需要你安裝額外的包。如果你用上述代碼時(shí)有報(bào)錯(cuò)夺姑，你可以安裝報(bào)錯(cuò)中提示你缺失的包來解決問題墩邀。比如說：

> biocLite(c("DDRTree", "pheatmap"))

分析protocol

你不需要知道Monocle包所有的分析功能，有些步驟你可以有很多種方法去做盏浙。流程分為多個(gè)步驟眉睹，我們會標(biāo)記出來哪些步驟你可以用其他的方法來做。像這樣的標(biāo)記：

流程步驟

Workflow一覽（這里先簡要的介紹有哪些步驟废膘，下面會具體的把每一步的都詳細(xì)的講解）

（一）儲存數(shù)據(jù)到一個(gè)CellDataSet對象

第一步就是load你的數(shù)據(jù)到Monocle的 CellDataSet里:

> pd <- new("AnnotatedDataFrame", data = sample_sheet)
> fd <- new("AnnotatedDataFrame", data = gene_annotation)
> cds <- newCellDataSet(expr_matrix, phenoData = pd, featureData = fd)

（二）用已知marker將細(xì)胞分類

利用你對關(guān)鍵marker基因的知識將細(xì)胞分類（NOTE:當(dāng)然有時(shí)候你的marker不是很特異或者很難分的很清晰竹海，后面會詳細(xì)講解分4種情況來進(jìn)行細(xì)胞分類）：

> cth <- newCellTypeHierarchy()
> MYF5_id <- row.names(subset(fData(cds), gene_short_name == "MYF5"))
> ANPEP_id <- row.names(subset(fData(cds), gene_short_name == "ANPEP"))
> cth <- addCellType(cth, "Myoblast", classify_func =
    function(x) { x[MYF5_id,] >= 1 })
> cth <- addCellType(cth, "Fibroblast", classify_func =
    function(x) { x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 } )
> cds <- classifyCells(cds, cth, 0.1)

（三）細(xì)胞聚類

> cds <- clusterCells(cds)

（四）沿軌跡在偽時(shí)間內(nèi)對細(xì)胞進(jìn)行排序

現(xiàn)在，把你的細(xì)胞按照你所研究的生物過程來排序丐黄，例如分化斋配，重編程，或者棉衣應(yīng)答灌闺。

> disp_table <- dispersionTable(cds)
> ordering_genes <- subset(disp_table, mean_expression >= 0.1)
> cds <- setOrderingFilter(cds, ordering_genes)
> cds <- reduceDimension(cds)
> cds <- orderCells(cds)

（五）差異分析

這一步就有很多種方法尋找差異基因和處理批次之間的效應(yīng)艰争。比如說：

> diff_test_res <- differentialGeneTest(cds,
    fullModelFormulaStr = "~Media")
> sig_genes <- subset(diff_test_res, qval < 0.1)

開始Monocle（詳細(xì)講解）：

Monocle采用Cufflinks或者包來計(jì)算基因表達(dá)矩陣。Monocle可以處理相對表達(dá)值（比如FPKM或者TPM）桂对，也可以處理絕對轉(zhuǎn)錄數(shù)（UMI實(shí)驗(yàn)）甩卓。Monocle還可以直接使用CellRanger（專門用于10×Genomics單細(xì)胞測序）處理得到的轉(zhuǎn)錄矩陣。Monocle還可以很好地處理來自其他RNA-Seq工作流程(如sci-RNA-Seq)和Biorad ddSEQ等儀器的數(shù)據(jù)蕉斜。盡管Monocle可以用原始counts分析逾柿，但除非你按照基因長度對counts進(jìn)行標(biāo)準(zhǔn)化缀棍，否則它們不會與表達(dá)值成正比，因此有些Monocle函數(shù)可能會產(chǎn)生無意義的結(jié)果鹿寻。如果你沒有UMI counts睦柴，我們建議你加載FPKM或TPM值進(jìn)行后續(xù)分析，而不是原始counts數(shù)毡熏。

（一）CellDataSet

Monocle把單細(xì)胞表達(dá)數(shù)據(jù)存放在CellDataSet對象里坦敌。這個(gè)對象源自Bioconductor的ExpressionSet對象，該對象要求三個(gè)輸入文件：

• exprs：是一個(gè)數(shù)值型的表達(dá)矩陣痢法，行是基因狱窘，列是細(xì)胞。
• phenoData：是一個(gè)AnnotatedDataFrame 對象,行是細(xì)胞财搁，列是細(xì)胞屬性（比如細(xì)胞類型蘸炸，培養(yǎng)條件，培養(yǎng)天數(shù)尖奔，等等） 搭儒。
• featureData：是一個(gè)AnnotatedDataFrame對象, 行是features（基因），列是基因?qū)傩裕ū热鏱iotype,GC含量提茁，等等）淹禾。

表達(dá)矩陣必須滿足如下要求：
（1）列數(shù)與phenoData行數(shù)一致。
（2）行數(shù)與featureData的data frame行數(shù)一致茴扁。
另外：
（3）phenoData對象的行名應(yīng)該與表達(dá)矩陣的列名一致铃岔。
（4）featureData對象的行名應(yīng)該與表達(dá)矩陣的行名一致。
（5）featureData的其中一列應(yīng)該命名為"gene_short_name"

你可以按照下面的步驟創(chuàng)建一個(gè)CellDataSet對象：

> HSMM_expr_matrix <- read.table("fpkm_matrix.txt")
> HSMM_sample_sheet <- read.delim("cell_sample_sheet.txt")
> HSMM_gene_annotation <- read.delim("gene_annotations.txt")

一旦上面這些table導(dǎo)入之后峭火，你可以創(chuàng)建cellDataSet對象：

> pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
> fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)
> HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
    phenoData = pd, featureData = fd)

這時(shí)將創(chuàng)建一個(gè)CellDataSet對象毁习，其中表達(dá)值以FPKM表示（Cufflinks方法計(jì)算）。默認(rèn)情況下卖丸，Monocle假定你的表達(dá)數(shù)據(jù)以轉(zhuǎn)錄本count數(shù)為單位纺且，并使用負(fù)二項(xiàng)式模型進(jìn)行下游的差異表達(dá)分析步驟。然而稍浆，如果你使用的是相對表達(dá)值例如TPM或者FPKM數(shù)據(jù)隆檀，請參照下面的內(nèi)容，如何“告訴”Monocle怎樣進(jìn)行下游的步驟粹湃。

如果你的數(shù)據(jù)里有UMI數(shù)據(jù)，在創(chuàng)建cellDataSet對象時(shí)不應(yīng)該標(biāo)準(zhǔn)化你的數(shù)據(jù)泉坐。你也不應(yīng)該用UMI的count去轉(zhuǎn)換相對豐度（比如FPKM/TPM）为鳄，你也不應(yīng)該用下面會提到的relative2abs()。monocle內(nèi)部會做所有需要的標(biāo)準(zhǔn)化步驟腕让。如果你自己標(biāo)準(zhǔn)化可能會對monocle某些關(guān)鍵步驟中斷孤钦。

（二）從其他包輸出文件導(dǎo)入/導(dǎo)出數(shù)據(jù)

Monocle可以將"Seurat"包的Seurat對象歧斟，和"Scater"包的SCESets轉(zhuǎn)換成CellDataSet對象，以便于Monocle的使用偏形。也可以處理Scran的SCESets静袖。要做轉(zhuǎn)換，用importCDS()功能：

# 'data_to_be_imported' 可以是一個(gè)Seurat 對象俊扭，也可以是一個(gè)SCESet
> importCDS(data_to_be_imported)
# 如果我們設(shè)置參數(shù)'import_all' 為 TRUE队橙，我們就可以把Seurat對象或者SCESet所有的slots導(dǎo)入進(jìn)來，默認(rèn)值是FALSE萨惑，只導(dǎo)入最小的dataset捐康。 
> importCDS(data_to_be_imported, import_all = TRUE)

Monocle還可以導(dǎo)出CellDataSets到Seurat和scater包，用exportCDS()功能：

> lung <- load_lung()
# To convert to Seurat object
> lung_seurat <- exportCDS(lung, 'Seurat')
# To convert to SCESet
> lung_SCESet <- exportCDS(lung, 'Scater')

（三）給你的數(shù)據(jù)選擇一個(gè)合適的distribution（分布）（RequiredＳ拱＝庾堋！）

Monocle可以處理相對表達(dá)數(shù)據(jù)和基于counts的測定(如UMIs)姐仅。一般來說花枫，它最適合轉(zhuǎn)錄本count數(shù)據(jù)，特別是UMI數(shù)據(jù)掏膏。無論你的數(shù)據(jù)類型是什么劳翰，為其指定適當(dāng)?shù)膁istribution是非常重要的。FPKM/TPM值通常是對數(shù)正態(tài)分布壤追，而UMIs或counts計(jì)數(shù)最好使用負(fù)二項(xiàng)式來建模磕道。處理counts數(shù)據(jù)，將負(fù)二項(xiàng)分布指定為newCellDataSet的expressionFamily參數(shù):

> HSMM <- newCellDataSet(count_matrix,
                phenoData = pd,
                featureData = fd,
                expressionFamily=negbinomial.size())

expressionFamily參數(shù)有幾個(gè)選項(xiàng)可以選擇：

官網(wǎng)說明：請根據(jù)你的情況選擇正確的expressionFamily參數(shù)選項(xiàng)行冰，如果選擇錯(cuò)誤溺蕉，可能會導(dǎo)致錯(cuò)誤的結(jié)果。如果你先用了relative2abs()處理了你的原始數(shù)據(jù)為FPKM/TPM悼做，你仍然可以使用negative binomial這個(gè)選項(xiàng)（但首先你用relative2abs()將相對表達(dá)量轉(zhuǎn)換為轉(zhuǎn)錄本counts）疯特。會給你更準(zhǔn)確的結(jié)果。

（四）處理大數(shù)據(jù)（Recommended）

一些單細(xì)胞RNA-Seq實(shí)驗(yàn)來自數(shù)萬個(gè)或更多細(xì)胞的測量結(jié)果肛走。隨著儀器的改進(jìn)和成本的降低漓雅，實(shí)驗(yàn)將變得越來越龐大，也越來越復(fù)雜朽色，包含有許多條件谭期、對照和重復(fù)老玛。一個(gè)包含50,000個(gè)細(xì)胞的表達(dá)矩陣，并且每個(gè)細(xì)胞包含有對人類基因組中25,000多個(gè)基因的測量就會占用大量的內(nèi)存。然而圆雁，由于目前的protocol通常不會捕獲每個(gè)細(xì)胞中的全部甚至大部分mRNA分子腊满，因此表達(dá)矩陣的許多基因?yàn)榱闾羯纭Ｊ褂孟∈杈仃嚕╯parse matrices）可以幫助你在普通的計(jì)算機(jī)上處理大型數(shù)據(jù)。我們通常建議大多數(shù)用戶使用sparseMatrices赵讯，因?yàn)樗梢约铀僭S多計(jì)算。

在sparse格式下處理你的數(shù)據(jù)耿眉，可以給Monocle提供一個(gè)sparse matrix:

> HSMM <- newCellDataSet(as(umi_matrix, "sparseMatrix"),
                phenoData = pd,
                featureData = fd,
                lowerDetectionLimit = 0.5,
                expressionFamily = negbinomial.size())

NOTE：CellRanger輸出的一系列RNA-Seq pipelines边翼，已經(jīng)是sparseMatrix的格式，所以你在做這一步的時(shí)候就不應(yīng)該再用as.matrix()來轉(zhuǎn)換格式了鸣剪，否則會占用你大量的內(nèi)存组底。
如果你的是10X Genomics單細(xì)胞測序數(shù)據(jù)（利用cellrangerRkit）處理的，你可以用下面的代碼導(dǎo)入你的數(shù)據(jù)到Monocle：

> cellranger_pipestance_path <- "/path/to/your/pipeline/output/directory"
> gbm <- load_cellranger_matrix(cellranger_pipestance_path)
> fd <- fData(gbm)

# 這里的2是隨便寫的西傀，你需要根據(jù)你的數(shù)據(jù)來改相應(yīng)的數(shù)字
# 這個(gè)“2”你應(yīng)該寫列數(shù)斤寇，這個(gè)列數(shù)與你的featureData的gene short names相一致。

> colnames(fd)[2] <- "gene_short_name"
> gbm_cds <- newCellDataSet(exprs(gbm),
                  phenoData = new("AnnotatedDataFrame", data = pData(gbm)),
                  featureData = new("AnnotatedDataFrame", data = fd),
                  lowerDetectionLimit = 0.5,
                  expressionFamily = negbinomial.size())

（五）將TPM/FPKM值轉(zhuǎn)換為mRNA counts（Alternative）

如果你的單細(xì)胞測序?qū)嶒?yàn)里有spike-in control拥褂，你可以將這些測量值轉(zhuǎn)換為每個(gè)細(xì)胞的mRNA的值 (RPC)娘锁。RPC值通常比FPKM或TPM值分析起來更容易。事實(shí)上饺鹃，即使沒有在實(shí)驗(yàn)中加入spike-in control莫秆，也有方法可以將FPKM或TPM值轉(zhuǎn)換為RPC值。Monocle 2包含一個(gè)名為Census的算法可以執(zhí)行這種轉(zhuǎn)換悔详。您可以在創(chuàng)建CellDataSet對象之前使用relative2abs()函數(shù)將其轉(zhuǎn)換為RPC值镊屎，如下所示:

> pd <- new("AnnotatedDataFrame", data = HSMM_sample_sheet)
> fd <- new("AnnotatedDataFrame", data = HSMM_gene_annotation)

# 首先用相對表達(dá)水平創(chuàng)建一個(gè)CellDataSet 
> HSMM <- newCellDataSet(as.matrix(HSMM_expr_matrix),
                phenoData = pd,
                featureData = fd,
                lowerDetectionLimit = 0.1,
                expressionFamily = tobit(Lower = 0.1))

# 接下來，轉(zhuǎn)換成RNA的counts
> rpc_matrix <- relative2abs(HSMM, method = "num_genes")

# 現(xiàn)在茄螃，可以用RNA counts創(chuàng)建一個(gè)新的CellDataSet
> HSMM <- newCellDataSet(as(as.matrix(rpc_matrix), "sparseMatrix"),
                phenoData = pd,
                featureData = fd,
                lowerDetectionLimit = 0.5,
                expressionFamily = negbinomial.size())

TPM/FPKM數(shù)據(jù)和mRNA count數(shù)據(jù)使用不同的newCellDataSet參數(shù)缝驳。由于我們使用了Census的mRNA count值，我們已經(jīng)修改了lowerDetectionLimit這個(gè)參數(shù)（從0.1改成0.5）來反應(yīng)新度量的表達(dá)量归苍。重要的一點(diǎn)是用狱，我們還修改了expressionFamily參數(shù)的選項(xiàng)為negbinomial()，告訴Monocle用負(fù)二項(xiàng)分布來進(jìn)行下游統(tǒng)計(jì)分析拼弃。如果你不修改這兩個(gè)參數(shù)值夏伊，后續(xù)可能會產(chǎn)生問題，所以如果你使用Census counts吻氧，不要忘記修改這兩個(gè)參數(shù)Ｄ缬恰！盯孙！

（六）估計(jì)size factor和分散度（RequiredＢ成！Ｕ穸琛）

我們還要用兩個(gè)函數(shù)來pre-calculate一些信息歌溉。size facotr幫助我們標(biāo)準(zhǔn)化細(xì)胞之間的mRNA的差異。分散度值可以幫助我們進(jìn)行后續(xù)的差異分析报账。
如果你在上面的步驟里研底，CellDataSet用的是negbinomial() or negbinomial.size()，那么你只需使用estimateSizeFactors() and estimateDispersions()透罢。

> HSMM <- estimateSizeFactors(HSMM)
> HSMM <- estimateDispersions(HSMM)

（七）過濾低質(zhì)量的細(xì)胞（Recommended）

任何單細(xì)胞測序分析的第一步都是要鑒定低質(zhì)量文庫榜晦。大部分單細(xì)胞處理過程中都會包含一些死細(xì)胞，或者空孔羽圃。另外去掉doublets也是很重要的：就是有兩個(gè)甚至更多的細(xì)胞的文庫乾胶。這些細(xì)胞會破壞下游步驟，例如偽時(shí)間排序和聚類朽寞。這部分舉例一個(gè)比較經(jīng)典的質(zhì)控步驟识窿。要看多少細(xì)胞里表達(dá)某個(gè)基因，或者看一個(gè)細(xì)胞里表達(dá)多少基因是非常容易的脑融。Monocle提供了一個(gè)非常簡單的函數(shù)來計(jì)算：

> HSMM <- detectGenes(HSMM, min_expr = 0.1)
> print(head(fData(HSMM)))

> HSMM <- detectGenes(HSMM, min_expr = 0.1)
> print(head(fData(HSMM)))
> expressed_genes <- row.names(subset(fData(HSMM),
    num_cells_expressed >= 10))

現(xiàn)在這個(gè)expressed_genes 里的基因都是至少在10個(gè)細(xì)胞里有表達(dá)的（這里官網(wǎng)寫錯(cuò)啦喻频，他寫成50了）。我們將在后續(xù)生物過程分析里用這個(gè)表肘迎。你也可以根據(jù)自己的需要進(jìn)行篩選甥温。

你的單細(xì)胞測序?qū)嶒?yàn)過程中也可以用成像來進(jìn)行質(zhì)控（捕獲細(xì)胞之后，裂解細(xì)胞之前）妓布。成像法可以評估你的細(xì)胞姻蚓，確定你的文庫里不會包含空孔，或者細(xì)胞碎片匣沼。一些儀器在捕獲細(xì)胞時(shí)可能會使得你的一個(gè)孔里有多于一個(gè)的細(xì)胞狰挡。你應(yīng)該把這些文庫去除掉。特別是空孔和細(xì)胞碎片的孔释涛。檢查每個(gè)細(xì)胞的RNA-seq文庫的sequence是否達(dá)到可接受的程度也是一個(gè)過濾的方法加叁。

然而，對于測序的“深度”并沒有一個(gè)最低標(biāo)準(zhǔn)枢贿。但是運(yùn)用你的判斷：如果一個(gè)細(xì)胞測序只有幾千個(gè)reads的話殉农，是不太可能產(chǎn)生有意義的結(jié)果的。

在CellDataSet對象里局荚，存放每個(gè)細(xì)胞評分的data在phenoData里超凳。評分屬性作為一列。 你可以過濾那些不符合你要求的細(xì)胞耀态。你也可以用FastQC來過濾細(xì)胞轮傍。這類軟件可以鑒別RNA-seq文庫里嚴(yán)重降解的RNA，也可以鑒別包含有核糖體首装，線粒體或者其他類型RNA污染的文庫创夜。
在我們的例子中，HSMM數(shù)據(jù)庫包含的評分列：

> print(head(pData(HSMM)))

如果你用的是RPC值來測定表達(dá)量仙逻，也可以根據(jù)樣本表達(dá)量進(jìn)行過濾：

> pData(HSMM)$Total_mRNAs <- Matrix::colSums(exprs(HSMM))

> HSMM <- HSMM[,pData(HSMM)$Total_mRNAs < 1e6]

> upper_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) +
            2*sd(log10(pData(HSMM)$Total_mRNAs)))
> lower_bound <- 10^(mean(log10(pData(HSMM)$Total_mRNAs)) -
            2*sd(log10(pData(HSMM)$Total_mRNAs)))

> qplot(Total_mRNAs, data = pData(HSMM), color = Hours, geom =
"density") +
> geom_vline(xintercept = lower_bound) +
> geom_vline(xintercept = upper_bound)

> HSMM <- HSMM[,pData(HSMM)$Total_mRNAs > lower_bound &
      pData(HSMM)$Total_mRNAs < upper_bound]
> HSMM <- detectGenes(HSMM, min_expr = 0.1)

我們還進(jìn)行了一步過濾驰吓，去掉那些mRNA水平非常低或者非常高的細(xì)胞涧尿。通常來講，doublets或者triplets是正常細(xì)胞mRNA的兩倍檬贰，所以要把這些細(xì)胞也要排除掉姑廉。如果你沒有條件通過成像法檢查你的細(xì)胞，那么你可以在這一步來過濾翁涤。一般設(shè)置的“門檻”是10000和40000 mRNA桥言。你可以根據(jù)你的實(shí)驗(yàn)需要調(diào)整這個(gè)門檻。一旦你排除掉不符合你要求的細(xì)胞葵礼，你應(yīng)該驗(yàn)證一下儲存在CelDataSet里的表達(dá)量大致是一個(gè)lognormal的分布：

# 先把表達(dá)矩陣用log轉(zhuǎn)換一下
> L <- log(exprs(HSMM[expressed_genes,]))

# 標(biāo)準(zhǔn)化每個(gè)基因号阿，然后合并你的data
> melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L))))

# 畫標(biāo)準(zhǔn)化后的基因表達(dá)值的分布
> qplot(value, geom = "density", data = melted_dens_df) +
stat_function(fun = dnorm, size = 0.5, color = 'red') +
xlab("Standardized log(FPKM)") +
ylab("Density")

（二）細(xì)胞分類和細(xì)胞計(jì)數(shù)

單細(xì)胞實(shí)驗(yàn)經(jīng)常用復(fù)雜混合的多細(xì)胞類型的樣品。從組織里分離出來的細(xì)胞往往含有兩個(gè)鸳粉、三個(gè)扔涧、甚至更多不同類型的細(xì)胞。在這種情況下赁严，最好是根據(jù)細(xì)胞類型的特異性marker將細(xì)胞分類扰柠。在我們的這個(gè)例子里，培養(yǎng)的細(xì)胞是從肌肉樣品來源的（包含成纖維細(xì)胞）疼约，用于原代細(xì)胞培養(yǎng)卤档。成肌細(xì)胞表達(dá)一些成纖維細(xì)胞沒有的基因。選擇特異性表達(dá)的基因程剥，比如說劝枣，利用高水平的MYF5可以排除成纖維細(xì)胞。另外织鲸，成纖維細(xì)胞表達(dá)高水平的ANPEP(CD13)舔腾，而成肌細(xì)胞不表達(dá)。

有幾種細(xì)胞分類的方法：

（1）按照細(xì)胞類型將細(xì)胞分類 （Recommended）

Monocle提供了一個(gè)簡單的系統(tǒng)來按照你選擇的基因來標(biāo)記細(xì)胞搂擦。比如說：你可以提供每一個(gè)類型細(xì)胞的function稳诚。這些functions可以是每一個(gè)細(xì)胞的表達(dá)data的input，返回TRUE來告訴Monocle這個(gè)細(xì)胞滿足你定義的這個(gè)function瀑踢。你可以用一個(gè)function來返回成肌細(xì)胞特異表達(dá)基因的細(xì)胞TRUR扳还，另一個(gè)function為成纖維細(xì)胞特異基因，等等橱夭。下面是個(gè)例子：

# 根據(jù)基因名字找到其在表達(dá)矩陣的ID
> MYF5_id <- row.names(subset(fData(HSMM), gene_short_name == "MYF5"))
> ANPEP_id <- row.names(subset(fData(HSMM),
   gene_short_name == "ANPEP"))
# 這里選取的基因取決于自己的單細(xì)胞實(shí)驗(yàn)設(shè)計(jì)
> cth <- newCellTypeHierarchy()
> cth <- addCellType(cth, "Myoblast", classify_func =
   function(x) { x[MYF5_id,] >= 1 })
> cth <- addCellType(cth, "Fibroblast", classify_func = function(x)
{ x[MYF5_id,] < 1 & x[ANPEP_id,] > 1 })

> HSMM <- classifyCells(HSMM, cth, 0.1)

classifyCells 功能將所有的gating function去判斷每一個(gè)細(xì)胞是否復(fù)合標(biāo)準(zhǔn)氨距，從而進(jìn)行分類，然后會返回 CellDataSet對象里新的一列,稱為 CellType棘劣。我們現(xiàn)在可以計(jì)數(shù)每一個(gè)細(xì)胞類型里都有多少個(gè)細(xì)胞俏让。

> table(pData(HSMM)$CellType)

| Fibroblast | Myoblast | Unknown |
| 56 | 85 | 121 |

還可以畫圖來展示細(xì)胞數(shù)：

> pie <- ggplot(pData(HSMM),
> aes(x = factor(1), fill = factor(CellType))) + geom_bar(width = 1)
> pie + coord_polar(theta = "y") +
> theme(axis.title.x = element_blank(), axis.title.y = element_blank())

你會發(fā)現(xiàn)很多細(xì)胞被列為"Unknown"。主要是因?yàn)樵诖蟛糠值膯渭?xì)胞測序?qū)嶒?yàn)中，mRNA捕獲效率低首昔，比如一個(gè)細(xì)胞表達(dá)很低水平的MYF5寡喝，但是我們很不幸沒有捕獲到它。 當(dāng)一個(gè)細(xì)胞不滿足你定下的任何一個(gè)function標(biāo)準(zhǔn)勒奇，就會標(biāo)記為 "Unknown"拘荡。如果它滿足多個(gè)標(biāo)準(zhǔn)，會被標(biāo)記為"Ambiguous"撬陵。你可以排除這些細(xì)胞，但是你會丟掉你的很多data网缝，在這個(gè)例子里巨税，你會損失大概一半的細(xì)胞。

（2）不根據(jù)marker基因?qū)⒓?xì)胞分類 （Alternative）

Monocle提供一種算法粉臊，你可以用來劃分你的"Unknown"的細(xì)胞草添。這種算法可以通過clusterCells來實(shí)現(xiàn), 根據(jù)整體的表達(dá)量來把細(xì)胞分組。比方說扼仲，如果你的細(xì)胞表達(dá)大量的成肌細(xì)胞基因远寸，但是缺少M(fèi)YF5，我們?nèi)匀豢梢园阉R別為成肌細(xì)胞屠凶。 clustercell可以以無監(jiān)督的方式使用驰后，也可以以半監(jiān)督的方式使用。我們先來看看無監(jiān)督模式：

第一步是確定哪些基因被用來劃分細(xì)胞群矗愧。我們可以用所有的基因灶芝，但是很多基因沒有足夠高的表達(dá)量。加上這些基因只會增加背景噪音唉韭。我們可以根據(jù)平均表達(dá)水平過濾基因夜涕，我們還可以另外選取那些在細(xì)胞之間不經(jīng)常變動(dòng)的基因。這些基因能夠?qū)?xì)胞狀態(tài)提供很有效的信息属愤。

> disp_table <- dispersionTable(HSMM)
> unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
> HSMM <- setOrderingFilter(HSMM, unsup_clustering_genes$gene_id)
> plot_ordering_genes(HSMM)

setOrderingFilter函數(shù)會標(biāo)記基因（在后續(xù)clusterCells細(xì)胞聚類所用的基因）女器。 盡管我們也可以提供一個(gè)我們自己寫的基因列表。plot_ordering_genes功能根據(jù)細(xì)胞之間的平局表達(dá)量來顯示一個(gè)基因的變異性 住诸。紅線顯示的是Monocle基于這種關(guān)系的預(yù)期值驾胆。我們把用于細(xì)胞分群的基因標(biāo)記為黑點(diǎn)，其他基因標(biāo)記為灰點(diǎn)只壳。

下面我們嘗試給細(xì)胞聚類：

# HSMM@auxClusteringData[["tSNE"]]$variance_explained <- NULL
# 這里看看基因的表達(dá)量和基因的變異度之間的關(guān)系
> plot_pc_variance_explained(HSMM, return_all = F) # norm_method='log'

> HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 6,
                reduction_method = 'tSNE', verbose = T)
> HSMM <- clusterCells(HSMM, num_clusters = 2)
# 這里先用tSNE的聚類方法處理HSMM數(shù)據(jù)集俏拱，并可視化展示
> plot_cell_clusters(HSMM, 1, 2, color = "CellType",
    markers = c("MYF5", "ANPEP"))

Monocle使用t-SNE來給細(xì)胞聚類，使用的方法與Seurat包非常相似吼句。

上圖里的根據(jù)我們提供的gating functions被標(biāo)記為成肌細(xì)胞的用綠色表示锅必，成纖維細(xì)胞標(biāo)記為紅色。不表達(dá)任何marker的細(xì)胞標(biāo)記為藍(lán)色。在許多實(shí)驗(yàn)中搞隐，細(xì)胞可以被很清楚的劃分驹愚，但是在這個(gè)實(shí)驗(yàn)里，細(xì)胞聚類并不容易劣纲。綠色和紅色的細(xì)胞并沒有清晰的分界線逢捺，是由于成肌細(xì)胞和混合進(jìn)來的間質(zhì)的成纖維細(xì)胞表達(dá)許多同樣的基因（在這個(gè)培養(yǎng)條件下）。還有其他的一些因素可以影響分群癞季。實(shí)驗(yàn)中的一個(gè)影響因素源于實(shí)驗(yàn)設(shè)計(jì)劫瞳。為了讓成肌細(xì)胞分化，我們將培養(yǎng)基從高絲裂原生長培養(yǎng)基(GM)轉(zhuǎn)換為低絲裂原分化培養(yǎng)基(DM)绷柒。也許這些細(xì)胞是基于它們所培養(yǎng)的培養(yǎng)基而聚集在一起的?可以試一下用培養(yǎng)基來分群：

> plot_cell_clusters(HSMM, 1, 2, color = "Media")

Monocle允許我們?nèi)サ簟安桓信d趣”來源變異的影響志于，從而降低這些因素對聚類的影響。你可以在clusterCells里添加一個(gè)參數(shù)residualModelFormulaStr废睦，這個(gè)參數(shù)接受R字符串模型伺绽，指定要在聚類之前減去的效果。

# 我們假設(shè)只有2種細(xì)胞類型嗜湃，所以在做聚類的時(shí)候可以把這個(gè)參數(shù)添加進(jìn)去奈应，這樣可以去除無關(guān)變量的干擾。
> HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 2,
           reduction_method = 'tSNE',
           residualModelFormulaStr = "~Media + num_genes_expressed",
           verbose = T)
> HSMM <- clusterCells(HSMM, num_clusters = 2)
> plot_cell_clusters(HSMM, 1, 2, color = "CellType")

現(xiàn)在紅色和綠色的細(xì)胞就能被分離開了

既然我們已經(jīng)解釋了一些不想要的變異來源购披，我們準(zhǔn)備再嘗試一次用無監(jiān)督聚類來分類細(xì)胞：

>HSMM <- clusterCells(HSMM, num_clusters = 2)
>plot_cell_clusters(HSMM, 1, 2, color = "Cluster") +
    facet_wrap(~CellType)

現(xiàn)在杖挣，大部分的成肌細(xì)胞在一個(gè)群里，大部分成纖維細(xì)胞在另一個(gè)群里刚陡，然而我們?nèi)匀豢吹接行┘?xì)胞在兩個(gè)群里都有程梦，這可能是因?yàn)槲覀兊膍arker基因和CellTypeHierarchy的functions缺乏特異性，但也可能是因?yàn)榫垲惒粔蚶硐腴佘榱伺懦笳哂旄剑屛覀儑L試以其半監(jiān)督模式運(yùn)行clusterCells。

（3）根據(jù)marker基因聚類細(xì)胞（Recommended）

首先哥童，我們選擇一個(gè)不同的基因列表來對細(xì)胞聚類挺份。之前我們是選取了那些高表達(dá)并且高變異性的基因。現(xiàn)在我們選取那些與我們的marker基因共變的基因贮懈。在某種意義上匀泊，我們建立一個(gè)大的基因列表作為marker，以便于即使一個(gè)細(xì)胞不表達(dá)MYF5朵你，但是它可能表達(dá)其他的基因各聘，而被識別為成肌細(xì)胞。這種聚類也叫“半監(jiān)督聚類”抡医。

> marker_diff <- markerDiffTable(HSMM[expressed_genes,],
            cth,
            residualModelFormulaStr = "~Media + num_genes_expressed",
            cores = 1)

markerDiffTable功能可以把所有細(xì)胞按照你提供的functions進(jìn)行分類躲因。 然后它會在識別出不同類型間差異表達(dá)的基因之前去掉所有的"Unknown" 和"Ambiguous" functions早敬。該函數(shù)返回一個(gè)data frame，你可以用這個(gè)去挑選用來聚類的基因大脉。通常情況下搞监，最好選擇每種細(xì)胞類型最具特異性的前10或20個(gè)基因。這確保了聚類基因不會被某一種細(xì)胞類型的marker所主導(dǎo)镰矿。Monocle提供了一個(gè)方便的功能琐驴，可以根據(jù)每種基因的表達(dá)受限制程度對基因進(jìn)行排序。

# 對每個(gè)基因增加了pval和qval兩列信息秤标，挑選出那些在不同media培養(yǎng)條件下顯著差異表達(dá)的基因
> candidate_clustering_genes <-
    row.names(subset(marker_diff, qval < 0.01))
# 計(jì)算這些基因在不同的celltype的specificity值
> marker_spec <-
  calculateMarkerSpecificity(HSMM[candidate_clustering_genes,], cth)
> head(selectTopMarkers(marker_spec, 3))

上面最后一行顯示了成肌細(xì)胞和成纖維細(xì)胞的top3的marker gene绝淡。 "specificity"評分范圍可以從0到1。它越接近1苍姜，對細(xì)胞類型的限制就越多够委。你可以用它為已知的細(xì)胞類型定義新的marker，也可以挑選出那些你用來純化新發(fā)現(xiàn)的細(xì)胞類型的基因怖现。這對于后續(xù)的實(shí)驗(yàn)非常有價(jià)值。
要將這些細(xì)胞聚在一起玉罐，我們將為每種細(xì)胞類型選擇前500個(gè)標(biāo)記:

> semisup_clustering_genes <- unique(selectTopMarkers(marker_spec, 500)$gene_id)
> HSMM <- setOrderingFilter(HSMM, semisup_clustering_genes)
> plot_ordering_genes(HSMM)

需要注意的是屈嗤，我們有一小組基因，它們其中一些在實(shí)驗(yàn)中并沒有特別高的表達(dá)吊输。然而饶号，它們對于區(qū)分表達(dá)MYF5和表達(dá)ANPEP的細(xì)胞非常有用。我們已經(jīng)將它們在聚類中使用季蚂，但是即使我們不這樣做茫船，我們?nèi)匀豢梢詫⑺鼈冎苯犹峁┙oclusterCells來使用它們。

# 重新挑選基因扭屁，用黑色點(diǎn)基因來進(jìn)行聚類算谈。
> plot_pc_variance_explained(HSMM, return_all = F)

> HSMM <- reduceDimension(HSMM, max_components = 2, num_dim = 3,
  norm_method = 'log',
  reduction_method = 'tSNE',
  residualModelFormulaStr = "~Media + num_genes_expressed",
  verbose = T)
> HSMM <- clusterCells(HSMM, num_clusters = 2)
> plot_cell_clusters(HSMM, 1, 2, color = "CellType")

（4）Imputing cell type （Alternative）（這個(gè)title不知道怎么翻譯才到位）

需要注意，我們已經(jīng)減少了成肌細(xì)胞群中“污染”的成纖維細(xì)胞的數(shù)量料滥。但是那些“unknown”的細(xì)胞怎么處理?如果你用CellTypeHierarcy來提供clusterCells然眼，Monocle將使用它作為一個(gè)整體進(jìn)行分類，而不僅僅是對單個(gè)細(xì)胞處理葵腹。clustercell計(jì)算每個(gè)細(xì)胞類型在每個(gè)聚類里的頻率高每。當(dāng)一個(gè)聚類由某個(gè)細(xì)胞類型的10%(在本例中為10%)以上組成時(shí)，這個(gè)聚類中的所有細(xì)胞都設(shè)置為該類型践宴。如果一個(gè)聚類里由多個(gè)細(xì)胞類型組成鲸匿，那么整個(gè)集群將被標(biāo)記為“Ambiguous”。如果沒有細(xì)胞類型符合上述的threshold阻肩，則這個(gè)聚類被標(biāo)記為“Unknown”带欢。因此，Monocle可以幫助你歸類每個(gè)細(xì)胞的類型，即使沒有現(xiàn)存的marker洪囤。

> HSMM <- clusterCells(HSMM,
              num_clusters = 2,
              frequency_thresh = 0.1,
              cell_type_hierarchy = cth)
> plot_cell_clusters(HSMM, 1, 2, color = "CellType",
    markers = c("MYF5", "ANPEP"))

可以看出徒坡，根據(jù)MYF5的表達(dá)可以清晰的把細(xì)胞聚類。在兩個(gè)聚類里都有一些細(xì)胞表達(dá)ANPEP瘤缩，但是在成肌細(xì)胞群里還表達(dá)MYF5喇完。這并不令人驚訝，因?yàn)锳NPEP并不是成纖維細(xì)胞的一個(gè)特別特異性的marker基因剥啤〗跸總之，我們可以成功的把這些細(xì)胞分群：

> pie <- ggplot(pData(HSMM), aes(x = factor(1), fill =
    factor(CellType))) +
   geom_bar(width = 1)
> pie + coord_polar(theta = "y") +
theme(axis.title.x = element_blank(), axis.title.y = element_blank())

構(gòu)建單細(xì)胞軌跡

在發(fā)育過程中府怯，細(xì)胞應(yīng)對外界的刺激刻诊，從一種“狀態(tài)”轉(zhuǎn)變?yōu)榱硪环N“狀態(tài)”。細(xì)胞在不同的“狀態(tài)”下表達(dá)著不同的基因牺丙，產(chǎn)生不同的蛋白行使各自的功能则涯。隨著細(xì)胞在不同“狀態(tài)”間轉(zhuǎn)換，其自身也經(jīng)歷著轉(zhuǎn)錄水平的重構(gòu)冲簿，一些基因沉默粟判，而一些基因被激活。這些瞬時(shí)的“狀態(tài)”很難被鑒定峦剔，因?yàn)樵诟€(wěn)定的兩個(gè)端點(diǎn)狀態(tài)之間純化細(xì)胞幾乎是不可能的档礁。單細(xì)胞RNA-Seq可以讓你在不需要純化細(xì)胞的情況下看到這些狀態(tài)。然而吝沫，為了做到這一點(diǎn)呻澜，我們必須確定每個(gè)細(xì)胞可能的“狀態(tài)”范圍。

Monocle不是通過實(shí)驗(yàn)將細(xì)胞純化成離散的狀態(tài)惨险，而是使用一種算法來“學(xué)習(xí)”每個(gè)細(xì)胞必須經(jīng)歷的基因表達(dá)變化的順序羹幸，這是動(dòng)態(tài)生物學(xué)過程的一部分。一旦它了解了基因表達(dá)變化的整體“軌跡”辫愉，Monocle就可以將每個(gè)細(xì)胞置于其軌跡的適當(dāng)位置睹欲。你就可以使用Monocle的differential analysis toolkit來查找在軌跡過程中基因的調(diào)節(jié)。

如果該分析過程有多個(gè)結(jié)果一屋，Monocle將重建一個(gè)“分支”軌跡窘疮。這些分支對應(yīng)著細(xì)胞的“命運(yùn)決定”，Monocle提供了強(qiáng)大的工具來識別受影響的基因冀墨。后面我們會講到如何分析分支闸衫。Monocle依靠一種稱為“反向圖嵌入”的機(jī)器學(xué)習(xí)技術(shù)來構(gòu)建單細(xì)胞軌跡。

什么是偽時(shí)間诽嘉？
偽時(shí)間是衡量一個(gè)細(xì)胞在生物過程中（比如分化）進(jìn)程的一個(gè)方法蔚出。在許多生物學(xué)過程中弟翘，細(xì)胞并不是完全同步進(jìn)行的。在細(xì)胞分化等過程的單細(xì)胞表達(dá)研究中骄酗，捕獲的細(xì)胞分布在不同的過程里稀余。也就是說，即使在同一時(shí)間捕獲的細(xì)胞群中趋翻，有些細(xì)胞可能已經(jīng)進(jìn)入某個(gè)生物過程很久了睛琳，而有些細(xì)胞甚至還沒有開始這個(gè)過程。當(dāng)您想要了解細(xì)胞從一種狀態(tài)過渡到另一種狀態(tài)時(shí)所發(fā)生的調(diào)節(jié)變化的順序時(shí)踏烙，這種異步性會產(chǎn)生很大的問題师骗。跟蹤同時(shí)捕獲的細(xì)胞的基因表達(dá)會產(chǎn)生一個(gè)非常壓縮的基因動(dòng)力學(xué)，該基因表達(dá)的變異性將非常高讨惩。Monocle根據(jù)每個(gè)細(xì)胞在軌跡上的進(jìn)展對其進(jìn)行排序辟癌，從而降低了異步帶來的問題。Monocle不是隨時(shí)間來追蹤變化荐捻，而是沿著軌跡來追蹤的黍少，我們稱之為偽時(shí)間。偽時(shí)間是一個(gè)抽象的進(jìn)程單位：它只是一個(gè)細(xì)胞到軌跡起點(diǎn)的距離处面，沿著最短路徑測量厂置。軌跡的總長度是根據(jù)一個(gè)細(xì)胞從起始狀態(tài)移動(dòng)到結(jié)束狀態(tài)所經(jīng)歷的總轉(zhuǎn)錄水平的變化量來定義的。

分析流程：

在我們深入了解沿軌跡排列細(xì)胞的細(xì)節(jié)之前鸳君，先來了解一下Monocle在做什么。有三個(gè)主要步驟患蹂，每個(gè)步驟都涉及一個(gè)重要的機(jī)器學(xué)習(xí)或颊。

Step 1: 選擇定義過程的基因

推斷一個(gè)單細(xì)胞軌跡是一個(gè)機(jī)器學(xué)習(xí)的過程。第一步就是選擇基因传于，Monocle將把這些基因作為其機(jī)器學(xué)習(xí)方法的input囱挑。這個(gè)過程叫做“feature selection”，它對軌跡的形狀有很大的影響沼溜。在單細(xì)胞RNA-Seq中平挑，低水平表達(dá)的基因通常非常嘈雜，但某些基因包含著有關(guān)細(xì)胞狀態(tài)的重要信息系草。Monocle通過檢查這些基因在細(xì)胞群之間的表達(dá)模式來對細(xì)胞進(jìn)行排序通熄。Monocle尋找那些以“有趣的”(即不只是嘈雜)方式變化的基因，并利用它們來構(gòu)建數(shù)據(jù)找都。Monocle提供了多種工具來選擇基因唇辨，并產(chǎn)生一個(gè)準(zhǔn)確、具有生物學(xué)意義的軌跡能耻。你可以用這些工具進(jìn)行完全的“無監(jiān)督”分析赏枚，或者你也可以用你的專業(yè)知識選擇基因亡驰，從而形成Monocle的軌跡。我們稱為這種方式為“半監(jiān)督”饿幅。

Step 2: 數(shù)據(jù)降維

一旦我們選擇了用來給細(xì)胞排序的基因凡辱，Monocle就會對數(shù)據(jù)進(jìn)行降維處理。Monocle使用最近開發(fā)的一種稱為反向圖嵌入的算法來降低數(shù)據(jù)的維數(shù)。

Step 3:按偽時(shí)間排序細(xì)胞

將表達(dá)數(shù)據(jù)投射到低維空間后鞍陨，Monocle就準(zhǔn)備學(xué)習(xí)描述細(xì)胞如何從一種狀態(tài)過渡到另一種狀態(tài)的軌跡侠仇。Monocle假設(shè)軌跡是“樹狀”結(jié)構(gòu)的，一端是“根”续徽，另一端是“葉”。Monocle的工作就是盡可能地把最佳的“樹狀結(jié)構(gòu)”與數(shù)據(jù)相匹配亲澡。這項(xiàng)任務(wù)被稱為“manifold learning”钦扭。一個(gè)細(xì)胞從生物過程的開始階段（從根部開始），沿著“主干”前進(jìn)床绪，直到它到達(dá)第一個(gè)分支(如果有的話)客情。然后，這個(gè)細(xì)胞必須選擇一條分支癞己，并且沿著它“走”得越來越遠(yuǎn)膀斋，直到到達(dá)一片葉子。一個(gè)細(xì)胞的偽時(shí)間值是它從這片“葉子”返回到“根”的距離痹雅。
（下面是對以上三個(gè)步驟的詳細(xì)介紹）

Trajectory step 1: 選擇定義細(xì)胞過程的基因

比如：我們必須首先要決定用哪些基因來定義細(xì)胞的肌生成過程（這里官網(wǎng)開始用他的例子來具體的說明）仰担。我們想要的是一組基因，它們的表達(dá)量隨著我們研究的生物過程的進(jìn)展而增加(或減少)绩社。
理想情況下摔蓝，我們想盡可能少的使用系統(tǒng)生物學(xué)的先驗(yàn)知識。我們還想從數(shù)據(jù)中發(fā)現(xiàn)一些重要的排序基因愉耙，而不是僅僅依賴于文獻(xiàn)和教科書贮尉，因?yàn)槟菢拥脑捒赡軙谂判蛑幸肫睢Ｎ覀兺ǔＭ扑]一種更成熟的方法朴沿，稱為“dpFeature”猜谚。
分離排序基因的一種有效方法是比較生物過程開始時(shí)和結(jié)束時(shí)收集的細(xì)胞，并尋找差異表達(dá)基因赌渣。下面的代碼就是為了找到從生長培養(yǎng)基到分化培養(yǎng)基的轉(zhuǎn)換中產(chǎn)生差異表達(dá)的基因：

> diff_test_res <- differentialGeneTest(HSMM_myo[expressed_genes,],
              fullModelFormulaStr = "~Media")
> ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))

根據(jù)時(shí)間點(diǎn)的差異分析來選擇基因通常是非常有效的魏铅，但是如果我們沒有時(shí)間排列的數(shù)據(jù)該怎么辦呢?如果這些細(xì)胞在我們所研究的生物過程中是不同步的，Monocle通臣嵛撸可以從同時(shí)捕獲的單個(gè)細(xì)胞群中重建它們的軌跡沦零。)
一旦我們有了用于排序的基因id列表，我們需要在HSMM對象中設(shè)置它們货岭，因?yàn)榻酉聛淼膸讉€(gè)函數(shù)都將用到它們：

> HSMM_myo <- setOrderingFilter(HSMM_myo, ordering_genes)
> plot_ordering_genes(HSMM_myo)

Trajectory step 2: 降維

接下來我們把空間降到2維路操，便于可視化：

> HSMM_myo <- reduceDimension(HSMM_myo, max_components = 2,
    method = 'DDRTree')

Trajectory step 3:按照軌跡排序細(xì)胞

> HSMM_myo <- orderCells(HSMM_myo)

可視化：

> plot_cell_trajectory(HSMM_myo, color_by = "Hours")

軌跡是樹狀結(jié)構(gòu)的疾渴。我們可以看到，在0時(shí)間點(diǎn)收集的細(xì)胞位于樹的一個(gè)尖端附近屯仗，而其他細(xì)胞分布在兩個(gè)“分支”上搞坝。Monocle并不知道哪條分支作為“開始”，所以我們經(jīng)常需要使用root_state參數(shù)調(diào)用ordercell來指定“起點(diǎn)”魁袜。首先桩撮，我們繪制軌跡，這次按“狀態(tài)”給細(xì)胞上色:

> plot_cell_trajectory(HSMM_myo, color_by = "State")

“State”是Monocle對象的一部分峰弹。下面的函數(shù)“手動(dòng)”的定義細(xì)胞的“State”店量，然后將其傳遞給ordercell：

#實(shí)際上就是自己寫個(gè)函數(shù)
> GM_state <- function(cds){
  if (length(unique(pData(cds)$State)) > 1){
    T0_counts <- table(pData(cds)$State, pData(cds)$Hours)[,"0"]
    return(as.numeric(names(T0_counts)[which
          (T0_counts == max(T0_counts))]))
  } else {
    return (1)
  }
}
> HSMM_myo <- orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
> plot_cell_trajectory(HSMM_myo, color_by = "Pseudotime")

如果你畫出來的“樹”里有大量的"state“，那么很難確定每個(gè)”state“落在”樹“的位置鞠呈。你可以將每一個(gè)"state"單獨(dú)畫出來融师，比如說：

> plot_cell_trajectory(HSMM_myo, color_by = "State") +
    facet_wrap(~State, nrow = 1)

如果你沒有一系列的時(shí)間點(diǎn)，你可以根據(jù)marker基因的表達(dá)來設(shè)置”樹根“的位置蚁吝，例如旱爆，在這個(gè)實(shí)驗(yàn)中，一個(gè)高度增殖的祖細(xì)胞群產(chǎn)生兩種有絲分裂后的細(xì)胞窘茁。所以“樹根”所在位置的細(xì)胞應(yīng)該有高水平的增殖marker表達(dá)怀伦。我們可以使用jitter plot來找出快速增殖所對應(yīng)的"state":

> blast_genes <- row.names(subset(fData(HSMM_myo),
#比如這里作者就選了3個(gè)marker基因
> gene_short_name %in% c("CCNB2", "MYOD1", "MYOG")))
> plot_genes_jitter(HSMM_myo[blast_genes,],
    grouping = "State",
    min_expr = 0.1)

為了確認(rèn)排序是正確的，我們可以選擇幾個(gè)肌源性進(jìn)展的marker來畫一下圖：

> HSMM_expressed_genes <-  row.names(subset(fData(HSMM_myo),
num_cells_expressed >= 10))
> HSMM_filtered <- HSMM_myo[HSMM_expressed_genes,]
> my_genes <- row.names(subset(fData(HSMM_filtered),
          gene_short_name %in% c("CDK1", "MEF2C", "MYH3")))
> cds_subset <- HSMM_filtered[my_genes,]
> plot_genes_in_pseudotime(cds_subset, color_by = "Hours")

其他的方法選擇排序基因

（1）根據(jù)聚類的差異基因排序（Recommended）

我們建議使用“dpFeature”的非監(jiān)督過程選擇基因山林。
為了使用dpFeature房待，我們首先要挑選出至少在5%的細(xì)胞里都表達(dá)的基因：

> HSMM_myo <- detectGenes(HSMM_myo, min_expr = 0.1)
> fData(HSMM_myo)$use_for_ordering <-
    fData(HSMM_myo)$num_cells_expressed > 0.05 * ncol(HSMM_myo)

然后，我們將執(zhí)行PCA分析驼抹，以確定每個(gè)主成分的variance桑孩。我們可以根據(jù)散點(diǎn)圖來確定需要多少個(gè)pca維度。

> plot_pc_variance_explained(HSMM_myo, return_all = F)

然后砂蔽，我們將使用t-SNE降維（運(yùn)行上圖位于頂部的variance）洼怔，并將它們進(jìn)一步投射到二維上署惯。

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

然后在二維t-SNE空間中進(jìn)行密度peak聚類來鑒定聚類左驾。densityPeak算法根據(jù)每個(gè)細(xì)胞的局部密度(Ρ)和兩個(gè)細(xì)胞之間最近的距離(Δ)進(jìn)行細(xì)胞聚類。默認(rèn)情況下,clusterCells選擇Ρ和Δ定義閾值的95%极谊。我們還可以設(shè)置一些想要進(jìn)行聚類的數(shù)量(n)诡右。默認(rèn)設(shè)置通常可以比較好的聚類轻猖。

> HSMM_myo <- clusterCells(HSMM_myo, verbose = F)

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

我們也可以根據(jù)自定義的閾值重新運(yùn)行聚類：

> HSMM_myo <- clusterCells(HSMM_myo,
                 rho_threshold = 2,
                 delta_threshold = 4,
                 skip_rho_sigma = T,
                 verbose = F)

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

在確定聚類是make sense的情況下帆吻，就可以提取差異基因了：

> clustering_DEG_genes <-
    differentialGeneTest(HSMM_myo[HSMM_expressed_genes,],
          fullModelFormulaStr = '~Cluster',
          cores = 1)

我們選取top1000作為排序基因：

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

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

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

> HSMM_myo <-
    orderCells(HSMM_myo)

> HSMM_myo <-
    orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))

> plot_cell_trajectory(HSMM_myo, color_by = "Hours")

（2）選擇在細(xì)胞間高度變化的基因用來排序（Alternative）

變化很大的基因通常為鑒定細(xì)胞亞群或沿軌跡對細(xì)胞進(jìn)行排序提供了大量信息。在RNA-Seq中咙边，一個(gè)基因的variance通常取決于它的均值猜煮，所以要小心地根據(jù)它們的variance選擇基因次员。

> disp_table <- dispersionTable(HSMM_myo)
> ordering_genes <- subset(disp_table,
                  mean_expression >= 0.5 &
                  dispersion_empirical >= 1 * dispersion_fit)$gene_id

（3）用已知marker基因進(jìn)行排序（Alternative）

無監(jiān)督排序是為了避免在分析中引入偏差。然而王带，非監(jiān)督的機(jī)器學(xué)習(xí)有時(shí)會focus在一個(gè)feature上淑蔚，而這并不是你實(shí)驗(yàn)的重點(diǎn)。例如愕撰，當(dāng)你使用無監(jiān)督學(xué)習(xí)時(shí)刹衫，每個(gè)細(xì)胞在細(xì)胞周期中的位置對軌跡的形狀有很大的影響。但是搞挣，如果你想研究的生物過程與周期無關(guān)該怎么辦呢?這時(shí)就需要用Monocle的“半監(jiān)督”排序模式带迟。
以半監(jiān)督的方式對細(xì)胞進(jìn)行排序非常簡單。首先使用CellTypeHierchy定義marker基因囱桨，這與我們之前使用它進(jìn)行細(xì)胞類型分類的方式非常相似仓犬。然后，用它來選擇與這些marker共變的排序基因蝇摸。最后婶肩，根據(jù)這些基因?qū)?xì)胞進(jìn)行排序。因此貌夕，無監(jiān)督和半監(jiān)督排序之間的唯一區(qū)別就是我們用來排序的基因律歼。
正如我們之前看到的，成肌細(xì)胞退出細(xì)胞周期開始分化啡专，然后通過一系列調(diào)控事件使得肌肉收縮所需的一些關(guān)鍵肌肉特異性蛋白的表達(dá)险毁。我們可以用cyclin B2 (CCNB2)標(biāo)記在細(xì)胞周期內(nèi)的細(xì)胞，并識別肌管们童，因?yàn)檫@些細(xì)胞表達(dá)了高水平的肌球蛋白重鏈3 (MYH3)畔况。

> CCNB2_id <-
    row.names(subset(fData(HSMM_myo), gene_short_name == "CCNB2"))
> MYH3_id <-
    row.names(subset(fData(HSMM_myo), gene_short_name == "MYH3"))

> cth <- newCellTypeHierarchy()
> cth <- addCellType(cth,
                   "Cycling myoblast",
                   classify_func = function(x) { x[CCNB2_id,] >= 1 })
> cth <- addCellType(cth,
                   "Myotube",
                   classify_func = function(x) { x[MYH3_id,] >= 1 })
> cth <- addCellType(cth,
                   "Reserve cell",
                   classify_func =
                function(x) { x[MYH3_id,] == 0 & x[CCNB2_id,] == 0 })

> HSMM_myo <- classifyCells(HSMM_myo, cth)

現(xiàn)在我們選擇與上面的兩個(gè)基因共變的一些基因：

> marker_diff <- markerDiffTable(HSMM_myo[HSMM_expressed_genes,],
                       cth,
                       cores = 1)
#下面這行代碼也可以寫成這樣：semisup_clustering_genes <- row.names(subset(marker_diff, qval < 0.05))
> semisup_clustering_genes <-
    row.names(marker_diff)[order(marker_diff$qval)][1:1000]

使用top1000個(gè)基因排序產(chǎn)生的軌跡與我們用非監(jiān)督方法得到的軌跡非常相似，但它更“干凈”一些:

> HSMM_myo <- setOrderingFilter(HSMM_myo, semisup_clustering_genes)
> HSMM_myo <- reduceDimension(HSMM_myo, max_components = 2,
    method = 'DDRTree', norm_method = 'log')
> HSMM_myo <- orderCells(HSMM_myo)
> HSMM_myo <- orderCells(HSMM_myo, root_state = GM_state(HSMM_myo))
> plot_cell_trajectory(HSMM_myo, color_by = "CellType") +
    theme(legend.position = "right")

為了確認(rèn)排序是正確的慧库，我們可以選擇幾個(gè)肌生成過程的marker跷跪。在這個(gè)實(shí)驗(yàn)中，一個(gè)分支對應(yīng)于成功融合形成肌管的細(xì)胞齐板，另一個(gè)分支對應(yīng)于未能完全分化的細(xì)胞〕痴埃現(xiàn)在我們排除后者，只把前者畫出來：

> HSMM_filtered <- HSMM_myo[HSMM_expressed_genes,]

> my_genes <- row.names(subset(fData(HSMM_filtered),
                gene_short_name %in% c("CDK1", "MEF2C", "MYH3")))

> cds_subset <- HSMM_filtered[my_genes,]
> plot_genes_branched_pseudotime(cds_subset,
                       branch_point = 1,
                       color_by = "Hours",
                       ncol = 1)

差異分析

差異分析是RNA-Seq分析中常見的任務(wù)甘磨。Monocle可以幫助你找到不同細(xì)胞組間差異表達(dá)的基因橡羞，并評估這些變化的顯著性。首先需要把你的細(xì)胞分到到兩個(gè)或更多的組里济舆。這些組由每個(gè)CellDataSet的phenoData的列定義卿泽。

（一）最基本的差異分析

對人類基因組中的所有基因進(jìn)行差異表達(dá)分析需要相當(dāng)長的時(shí)間。對于我們舉的例子這樣大的數(shù)據(jù)集(其中包含數(shù)百個(gè)細(xì)胞)滋觉，在單個(gè)CPU上進(jìn)行分析可能需要幾個(gè)小時(shí)签夭。所以先選擇一小組我們已知在肌生成中很重要的基因來試一下:

> marker_genes <- row.names(subset(fData(HSMM_myo),
                   gene_short_name %in% c("MEF2C", "MEF2D", "MYF5",
                                          "ANPEP", "PDGFRA","MYOG",
                                          "TPM1",  "TPM2",  "MYH2",
                                          "MYH3",  "NCAM1", "TNNT1",
                                          "TNNT2", "TNNC1", "CDK1",
                                          "CDK2",  "CCNB1", "CCNB2",
                                          "CCND1", "CCNA1", "ID1")))

在成肌細(xì)胞數(shù)據(jù)中齐邦，實(shí)驗(yàn)開始時(shí)收集的細(xì)胞在“生長培養(yǎng)基”(GM)中培養(yǎng)，以防止它們分化第租。在收集細(xì)胞后侄旬，其余細(xì)胞轉(zhuǎn)入“分化培養(yǎng)基”(DM)以促進(jìn)分化。首先來用Monocle尋找一下哪些基因可能受上述培養(yǎng)基轉(zhuǎn)換的影響：

> diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,],
                                      fullModelFormulaStr = "~Media")

# 根據(jù) FDR < 10%的基因
> sig_genes <- subset(diff_test_res, qval < 0.1)
> sig_genes[,c("gene_short_name", "pval", "qval")]

Monocle還提供了一些簡單的方法來繪制一小組基因的表達(dá)煌妈，這些基因是根據(jù)你在差異分析中使用的因素分組的儡羔。這有助于你可視化上述的差異。

> MYOG_ID1 <- HSMM_myo[row.names(subset(fData(HSMM_myo),
              gene_short_name %in% c("MYOG", "CCNB2"))),]
> plot_genes_jitter(MYOG_ID1, grouping = "Media", ncol= 2)

這里要注意的是璧诵，我們可以控制在上圖里怎么樣排布你的基因汰蜘，包括行數(shù)和列數(shù)。請參閱
plot_genes_jitter手冊之宿。

（二）根據(jù)細(xì)胞類型和狀態(tài)尋找差異基因

在動(dòng)態(tài)的生物過程中族操，例如分化，細(xì)胞可能呈現(xiàn)出不同的中間態(tài)或最終態(tài)比被。我們之前根據(jù)幾個(gè)關(guān)鍵的marker來區(qū)分成肌細(xì)胞和污染的成纖維細(xì)胞∩眩現(xiàn)在我們嘗試其他幾個(gè)也可能區(qū)分成纖維細(xì)胞和成肌細(xì)胞的基因：

> to_be_tested <- row.names(subset(fData(HSMM),
              gene_short_name %in% c("UBC", "NCAM1", "ANPEP")))
> cds_subset <- HSMM[to_be_tested,]

#下面這個(gè)differentialGeneTest函數(shù)是為了驗(yàn)證我們根據(jù)細(xì)胞類型所尋找的差異基因是不是具有顯著差異性
> diff_test_res <- differentialGeneTest(cds_subset,
                    fullModelFormulaStr = "~CellType")
> diff_test_res[,c("gene_short_name", "pval", "qval")]

> plot_genes_jitter(cds_subset,
                  grouping = "CellType",
                  color_by = "CellType",
                  nrow= 1,
                  ncol = NULL,
                  plot_trend = TRUE)

實(shí)際上Monocle的差異分析步驟非常的靈活多變，你可以用pData表里任何存在的列來進(jìn)行分析等缀。比如說枷莉，你之前用的是clusterCells來分類你的細(xì)胞，那么在差異分析時(shí)尺迂，你就可以用pData里Cluster這列來尋找差異基因笤妙。

（三）根據(jù)偽時(shí)間功能尋找差異基因

Monocle的主要工作是將細(xì)胞按照生物過程(如細(xì)胞分化)的順序進(jìn)行排列。這樣你就可以分析細(xì)胞里找到隨著生物過程進(jìn)展而改變的基因噪裕。比如說蹲盘，你可以發(fā)現(xiàn)當(dāng)細(xì)胞“成熟”時(shí)有哪些基因顯著上調(diào)。我們看一下一組對肌生成很重要的基因：

> to_be_tested <- row.names(subset(fData(HSMM),
> gene_short_name %in% c("MYH3", "MEF2C", "CCNB2", "TNNT1")))
> cds_subset <- HSMM_myo[to_be_tested,]

> diff_test_res <- differentialGeneTest(cds_subset,
fullModelFormulaStr = "~sm.ns(Pseudotime)")

sm.ns函數(shù)是通過表達(dá)值來擬合一個(gè)曲線膳音，以幫助它將基因表達(dá)中的變化描述為一個(gè)隨生物過程變化的函數(shù)召衔。

> diff_test_res[,c("gene_short_name", "pval", "qval")]
> plot_genes_in_pseudotime(cds_subset, color_by = "Hours")

（四）根據(jù)偽時(shí)間表達(dá)pattern聚類基因（熱圖）

Monocle提供了一個(gè)簡單的方法可視化所有偽時(shí)間依賴的基因。
plot_pseudotime_heatmap接受 CellDataSet 對象（通常包含一小組差異基因）祭陷。然后把這些基因聚類苍凛，并用pheatmap包畫熱圖。

> diff_test_res <- differentialGeneTest(HSMM_myo[marker_genes,],
              fullModelFormulaStr = "~sm.ns(Pseudotime)")
> sig_gene_names <- row.names(subset(diff_test_res, qval < 0.1))
> plot_pseudotime_heatmap(HSMM_myo[sig_gene_names,],
                num_clusters = 3,
                cores = 1,
                show_rownames = T)

多因素差異分析

Monocle可以在多個(gè)因素存在的情況下進(jìn)行差異分析颗胡，可以幫助你減去一些不必要的影響因素毫深。在下面的簡單例子中吩坝，Monocle測試在成肌細(xì)胞和成纖維細(xì)胞之間的差異表達(dá)的3個(gè)基因毒姨，同時(shí)減去Hours的影響（Hours指每個(gè)細(xì)胞收集的日期）。

> to_be_tested <-
    row.names(subset(fData(HSMM),
        gene_short_name %in% c("TPM1", "MYH3", "CCNB2", "GAPDH")))

> cds_subset <- HSMM[to_be_tested,]

> diff_test_res <- differentialGeneTest(cds_subset,
                        fullModelFormulaStr = "~CellType + Hours",
                        reducedModelFormulaStr = "~Hours")
> diff_test_res[,c("gene_short_name", "pval", "qval")]
> plot_genes_jitter(cds_subset,
    grouping = "Hours", color_by = "CellType", plot_trend = TRUE) +
        facet_wrap( ~ feature_label, scales= "free_y")

單細(xì)胞軌跡的“分支”分析

通常钉寝，單細(xì)胞的軌跡里有分支弧呐。分支的產(chǎn)生是因?yàn)榧?xì)胞有可選擇的基因表達(dá)模式闸迷。比如在發(fā)育過程中，當(dāng)細(xì)胞做出“命運(yùn)”的選擇時(shí)俘枫，分支就會出現(xiàn)在軌跡中:一個(gè)發(fā)育譜系沿著一條路徑前進(jìn)腥沽，而另一個(gè)譜系產(chǎn)生第二條路徑。Monocle可以分析這些分支事件鸠蚪。用BEAM函數(shù)今阳。

> lung <- load_lung()
> plot_cell_trajectory(lung, color_by = "Time")

BEAM的輸入對象是排序好的CellDataSet對象。你需要用orderCells 和軌跡分支點(diǎn)的名字來排序你的CellDataSet對象茅信。它會返回一個(gè)每個(gè)基因的顯著值的table盾舌。

> BEAM_res <- BEAM(lung, branch_point = 1, cores = 1)
> BEAM_res <- BEAM_res[order(BEAM_res$qval),]
> BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]

你可以用一種特殊的熱圖，將所有“分支依賴”的基因可視化蘸鲸。該熱圖顯示的是同一時(shí)間點(diǎn)兩個(gè)譜系的變化妖谴。熱圖的列是偽時(shí)間的點(diǎn)，行是基因酌摇。從熱圖中間往右讀膝舅，是偽時(shí)間的一個(gè)譜系；往左是另一個(gè)譜系窑多∪韵。基因是被按照等級聚類的，所以你看到的基因表達(dá)模式和譜系的表達(dá)模式是非常相似的埂息。

> plot_genes_branched_heatmap(lung[row.names(subset(BEAM_res,
                                          qval < 1e-4)),],
                                          branch_point = 1,
                                          num_clusters = 4,
                                          cores = 1,
                                          use_gene_short_name = T,
                                          show_rownames = T)

我們還可以用plot_genes_branched_pseudotime函數(shù)單獨(dú)畫幾個(gè)基因（細(xì)胞命運(yùn)決定的marker）：

> lung_genes <- row.names(subset(fData(lung),
          gene_short_name %in% c("Ccnd2", "Sftpb", "Pdpn")))
> plot_genes_branched_pseudotime(lung[lung_genes,],
                       branch_point = 1,
                       color_by = "Time",
                       ncol = 1)

附錄

計(jì)算單細(xì)胞表達(dá)值

用Monocle之前琳轿，你必須計(jì)算每一個(gè)基因在每個(gè)細(xì)胞里的表達(dá)。有很多種方法耿芹，我們推薦Cufflinks崭篡，你也可以用RSEM, eXpress, Sailfish等其他的包來計(jì)算。這里我們展示一個(gè)簡單的流程利用TopHat和Cufflinks組合來計(jì)算
首先你必須有reads矩陣吧秕，如果你是雙端測序琉闪，每個(gè)細(xì)胞應(yīng)該是兩個(gè)文件，像這樣：CELL_TXX_YYY.RZ.fastq.gz（XX是你收集細(xì)胞的時(shí)間點(diǎn)砸彬，YY是你在準(zhǔn)備文庫的時(shí)候96孔板颠毙，Z一般是1或者2，表示雙端測序時(shí)left mate或者right mate砂碉。所以例如CELL_T24_A01.R1.fastq.gz 意思是細(xì)胞在24小時(shí)收集的在A01孔里的文庫left mate文件）

用TopHat比對reads（這個(gè)有點(diǎn)過時(shí)了）

tophat -o CELL_T24_A01_thout -G GENCODE.gtf bowtie-hg19-idx CELL_T24_
    A01.R1.fastq.gz CELL_T24_A01.R2.fastq.gz
tophat -o CELL_T24_A02_thout -G GENCODE.gtf bowtie-hg19-idx CELL_T24_
    A02.R1.fastq.gz CELL_T24_A02.R2.fastq.gz
tophat -o CELL_T24_A03_thout -G GENCODE.gtf bowtie-hg19-idx CELL_T24_
    A03.R1.fastq.gz CELL_T24_A03.R2.fastq.gz

用Cufflinks計(jì)算基因表達(dá)

cuffquant -o CELL_T24_A01_cuffquant_out GENCODE.gtf
    CELL_T24_A01_thout/accepted_hits.bam
cuffquant -o CELL_T24_A02_cuffquant_out GENCODE.gtf
    CELL_T24_A02_thout/accepted_hits.bam
cuffquant -o CELL_T24_A03_cuffquant_out GENCODE.gtf
    CELL_T24_A03_thout/accepted_hits.bam

接下來把所有的表達(dá)矩陣merge到一個(gè)文件里：

> cuffnorm --use-sample-sheet -o sc_expr_out GENCODE.gtf
    sample_sheet.txt

這里的--use-sample-sheet前提是你要有一個(gè)表蛀蜜，里面是所有表達(dá)矩陣文件的名字，比如這樣：

sample_name                                   group
CELL_T24_A01_cuffquant_out/abundances.cxb     T24_A01
CELL_T24_A02_cuffquant_out/abundances.cxb     T24_A02
CELL_T24_A03_cuffquant_out/abundances.cxb     T24_A03

然后就可以快樂的用Monocle分析了～

總結(jié)：
后面還有一些數(shù)學(xué)公式解釋Monocle里面一些函數(shù)的原理增蹭，我就沒看了滴某，數(shù)學(xué)也不好，就不誤導(dǎo)別人了■荩總體來說户誓，過了一遍原始官方說明書，感覺就是覺得“非常細(xì)”幕侠，我覺得如果有時(shí)間的同學(xué)在用到每一個(gè)非常重要的包帝美，最好還是閱讀一下官方說明書。像我這樣成天白天做實(shí)驗(yàn)晤硕，只有晚上才有時(shí)間學(xué)習(xí)生信的人來說悼潭，有時(shí)為了偷懶，就直接copy別人的代碼了舞箍。但是還是想多了解一下里面的用法女责，以及“這一步到底是在干嘛？”總之這次還是有些收獲的创译，之后還會深入學(xué)習(xí)其他包（有時(shí)間的情況下）抵知。生信的東西就是積少成多，我這個(gè)小白就感觸頗深软族，今年夏天的時(shí)候我連“R包”是什么都不知道刷喜，但現(xiàn)在至少知道些東西了～（雖然還是非常的菜。立砸。掖疮。）好好學(xué)習(xí)，天天向上～共勉之颗祝！