參考：
測試所用的數(shù)據(jù)鏈接: https://pan.baidu.com/s/1TUysqHHbXrWGi7QGt-6L1g 提取碼: uqgh
Monocle[7]是一個基于pseudotime分析的軟件包掸宛。它通過從一系列差異表達(dá)的基因來描述細(xì)胞間的距離遠(yuǎn)近耳幢，由此得到類似樹狀的分叉圖象來表示細(xì)胞間的差異性半夷。 提供了聚類塞椎，pseudotime, 差異分析等多種功能悦陋，該項(xiàng)目的網(wǎng)址如下
主要內(nèi)容
1. 讀取數(shù)據(jù)
2.預(yù)處理過濾細(xì)胞
3. 細(xì)胞分群/聚類
? ? step1：dispersionTable()
? ? step2：plot_pc_variance_explained() 耗時
? ? step3: 進(jìn)行降維和聚類
4.篩選基因(用于擬分析時序)
? ? 4.1 一般過濾標(biāo)準(zhǔn)
? ? 4.2 使用差異表達(dá)基因BEAM分析
5.推斷發(fā)育軌跡
? ? 5.1. 選合適基因
? ? 5.2 降維
? ? 5.3 細(xì)胞排序
小結(jié)：
主要介紹使用該R包進(jìn)行pseudotime分析的步驟.
Tips : pseudotime分析本質(zhì)上就是輸入基因及其表達(dá)量矩陣巫俺，計算擬時序的值。并且選取一些基因剖淀，用于降維可視化纯蛾。聚類只是降維后顏色填充而已.
1. 讀取數(shù)據(jù)
monocle包有很多種讀取數(shù)據(jù)的方式 ，這里只展示其中三種：（I）直接讀取;（II）來源于seurat2 包; （III）來源于seurat3包祷蝌。
Tips : 由于版本更替較快茅撞，seurat2 可以直接導(dǎo)入monocle2.但是seurat3卻要手動寫代碼導(dǎo)入.
library(monocle)
library(DDRTree)
library(pheatmap)
library(data.table)
library(plyr)
下面以10x 數(shù)據(jù)為例，最好通過seurat對象直接讀入monocle2 方便一些.
（1）直接讀染揠： 當(dāng)有三個文件米丘，expr_matrix ; sample_info ;gene_annotation.
expr_matrix <- read.delim(nUMI.summary.file)
sample_info <- data.frame(sampleID=colnames(expr_matrix))
gene_annotation <- data.frame(symbol=rownames(expr_matrix))
pd <- new("AnnotatedDataFrame", data = sample_info)
fd <- new("AnnotatedDataFrame", data = gene_annotation)
d <- newCellDataSet(as.matrix(expr_matrix), phenoData = pd, featureData = fd)
（2）來源于seurat2 包，數(shù)據(jù)如下：
## 切換到工作目錄糊啡，讀入三個文本文件.
pbmc.data <- Read10X(data.dir = ".")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "10XPBMC")
## seurat2 可以直接讀入
sce <- importCDS(pbmc)
（3） 來源于seurat3包：
差別之處拄查，不能直接用importCDS 導(dǎo)入.
## 由seurat3 導(dǎo)入函數(shù)newimport 替代上面的importCDS功能.
newimport <- function(otherCDS, import_all = FALSE) {
if(class(otherCDS)[1] == 'Seurat') {
requireNamespace("Seurat")
data <- otherCDS@assays$RNA@counts
if(class(data) == "data.frame") {
data <- as(as.matrix(data), "sparseMatrix")
}
pd <- tryCatch( {
pd <- new("AnnotatedDataFrame", data = otherCDS@meta.data)
pd
},
#warning = function(w) { },
error = function(e) {
pData <- data.frame(cell_id = colnames(data), row.names = colnames(data))
pd <- new("AnnotatedDataFrame", data = pData)
message("This Seurat object doesn't provide any meta data");
pd
})
# remove filtered cells from Seurat
if(length(setdiff(colnames(data), rownames(pd))) > 0) {
data <- data[, rownames(pd)]
}
fData <- data.frame(gene_short_name = row.names(data), row.names = row.names(data))
fd <- new("AnnotatedDataFrame", data = fData)
lowerDetectionLimit <- 0
if(all(data == floor(data))) {
expressionFamily <- negbinomial.size()
} else if(any(data < 0)){
expressionFamily <- uninormal()
} else {
expressionFamily <- tobit()
}
valid_data <- data[, row.names(pd)]
monocle_cds <- newCellDataSet(data,
phenoData = pd,
featureData = fd,
lowerDetectionLimit=lowerDetectionLimit,
expressionFamily=expressionFamily)
if(import_all) {
if("Monocle" %in% names(otherCDS@misc)) {
otherCDS@misc$Monocle@auxClusteringData$seurat <- NULL
otherCDS@misc$Monocle@auxClusteringData$scran <- NULL
monocle_cds <- otherCDS@misc$Monocle
mist_list <- otherCDS
} else {
# mist_list <- list(ident = ident)
mist_list <- otherCDS
}
} else {
mist_list <- list()
}
if(1==1) {
var.genes <- setOrderingFilter(monocle_cds, otherCDS@assays$RNA@var.features)
}
monocle_cds@auxClusteringData$seurat <- mist_list
} else if (class(otherCDS)[1] == 'SCESet') {
requireNamespace("scater")
message('Converting the exprs data in log scale back to original scale ...')
data <- 2^otherCDS@assayData$exprs - otherCDS@logExprsOffset
fd <- otherCDS@featureData
pd <- otherCDS@phenoData
experimentData = otherCDS@experimentData
if("is.expr" %in% slotNames(otherCDS))
lowerDetectionLimit <- otherCDS@is.expr
else
lowerDetectionLimit <- 1
if(all(data == floor(data))) {
expressionFamily <- negbinomial.size()
} else if(any(data < 0)){
expressionFamily <- uninormal()
} else {
expressionFamily <- tobit()
}
if(import_all) {
# mist_list <- list(iotherCDS@sc3,
# otherCDS@reducedDimension)
mist_list <- otherCDS
} else {
mist_list <- list()
}
monocle_cds <- newCellDataSet(data,
phenoData = pd,
featureData = fd,
lowerDetectionLimit=lowerDetectionLimit,
expressionFamily=expressionFamily)
# monocle_cds@auxClusteringData$sc3 <- otherCDS@sc3
# monocle_cds@auxOrderingData$scran <- mist_list
monocle_cds@auxOrderingData$scran <- mist_list
} else {
stop('the object type you want to export to is not supported yet')
}
return(monocle_cds)
}
導(dǎo)入monocle2
## 讀入seurat3
pbmc.data <- Read10X(data.dir = ".")
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "10XPBMC")
## 讀入monocle2
sce <- newimport(pbmc)
本文采用第1種方式，讀入文件
library(Seurat)
library(monocle)
expression_matrix <- readRDS("packer_embryo_expression.rds")
cell_metadata <- readRDS("packer_embryo_colData.rds")
gene_annotation <- readRDS("packer_embryo_rowData.rds")
## 數(shù)據(jù)框需要先保存為AnnotatedDataFrame 格式棚蓄，才可以創(chuàng)建newCellDataSet 對象
# Error: CellDataSet 'phenoData' is class 'data.frame' but should be or extend 'AnnotatedDataFrame'
pd <- new("AnnotatedDataFrame", data = cell_metadata)
fd <- new("AnnotatedDataFrame", data = gene_annotation)
sce <- newCellDataSet(expression_matrix,
phenoData = pd,
featureData = fd,
lowerDetectionLimit = 0.1,
expressionFamily = VGAM::negbinomial.size()) # 默認(rèn)參數(shù)
查看sce 類型：
> sce
CellDataSet (storageMode: environment)
assayData: 20222 features, 6188 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGCAAGACGTG-300.1.1 AAACCTGGTGTGAATA-300.1.1 ...
TTTGTCAAGTACACCT-b02 (6188 total)
varLabels: cell n.umi ... State (21 total)
varMetadata: labelDescription
featureData
featureNames: WBGene00010957 WBGene00010958 ... WBGene00007064 (20222 total)
fvarLabels: id gene_short_name num_cells_expressed use_for_ordering
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
查看一下數(shù)據(jù)讀入到monocle2里面堕扶，它的表達(dá)分布類型：
sce@expressionFamily
## Family: negbinomial.size
##
## Negative-binomial distribution with size known
##
## Links: loge(mu)
## Mean: mu
## Variance: mu * (1 + mu / size) for NB-2
? 注意到構(gòu)建CDS對象過程中有一個參數(shù)是：expressionFamily ，它是選擇了一個數(shù)據(jù)分布梭依，例如FPKM/TPM 值是log-正態(tài)分布的稍算；UMIs和原始count值用負(fù)二項(xiàng)分布模擬的效果更好。負(fù)二項(xiàng)分布有兩種方法役拴，這里選用了negbinomial.size糊探，另外一種negbinomial稍微更準(zhǔn)確一點(diǎn)，但速度大打折扣，它主要針對非常小的數(shù)據(jù)集科平，如下表：
由于輸入10x genomic 數(shù)count 類似褥紫，所以表達(dá)分布為 negbinomial.size()
| Family function | 數(shù)據(jù)類型 | 備注 |
|---|---|---|
| negbinomial.size() | 原始計數(shù)值, 比如UMIs counts | Negative binomial distribution with fixed variance |
| negbinomial() | 原始計數(shù)值, 比如UMIs counts | 比negbinomial.size()精確一點(diǎn), 但是非常慢. |
| tobit() | FPKM, TPM | Tobits are truncated normal distributions. 注意這點(diǎn)不是log-transformed FPKM/TPM |
| gaussianff() | log-transformed FPKM/TPMs, Ct values from single-cell qPCR | 符合高期分布的數(shù)值 |
2.預(yù)處理
和處理RNA-seq數(shù)據(jù)一樣，monocle會有一些預(yù)處理的步驟瞪慧，這包括估計size factors和dispersions髓考，以及去除一些質(zhì)量比較差的細(xì)胞。
第一行代碼用于評估每個細(xì)胞的sizefactor, 用于后續(xù)歸一化弃酌；
第二行代碼計算基因的方差,這一步可能會提醒我們有多少的outlier的細(xì)胞氨菇。所以下一步就是要去除掉這些細(xì)胞。
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
# Removing 143 outliers
通過pData和fData兩個函數(shù)矢腻，可以查看基因和細(xì)胞的相關(guān)信息门驾，示意如下
- 查看細(xì)胞注釋信息：pData.
- 查看基因注釋信息：fData
過濾細(xì)胞
在準(zhǔn)備樣品細(xì)胞的過程中，不可避免地會引入一下死的細(xì)胞多柑，或者說是doublets(兩個細(xì)胞被標(biāo)記成一個細(xì)胞)。 這些細(xì)胞對下游分析的影響會比較大楣责，它們很有可能會占去很大分析權(quán)重竣灌。所以我們需要把它們?nèi)コ簟?去除它們的辦法就是設(shè)置一個比較合理的表達(dá)總值的上下限，然后把處于上下限外圍的細(xì)胞都去除掉秆麸。
- 首先我們來查看一下低表達(dá)的基因初嘹。
### detectGenes計算每一個細(xì)胞表達(dá)的基因個數(shù).
sce <- detectGenes(sce, min_expr = 0.1)
### 保留在10個細(xì)胞中表達(dá)的基因，當(dāng)做表達(dá)基因.
expressed_genes <- row.names(subset(fData(sce),
num_cells_expressed >= 10))
我們將表達(dá)值在10以上的基因保存為expressed_genes沮趣。下面的步驟會用到保存好的表達(dá)的基因屯烦。
-
接下來我們要去除掉死細(xì)胞或者doublets了。
head(pData(sce))
一般的房铭，我們是沒有數(shù)據(jù)告訴我們哪些細(xì)胞是死的或者空的驻龟，哪些細(xì)胞是doublets的。但是通常而方言缸匪，我們都會拿到RPC(reads per cell)的計數(shù)翁狐，我們可以人為的設(shè)定一下表達(dá)總值的上下限。
pData(sce)$Total_mRNAs <- Matrix::colSums(exprs(sce))
sce <- sce[, pData(sce)$Total_mRNAs < 1e6 ]
##設(shè)置上下限：
upper_bound <- 10^(mean(log10(pData(sce)$Total_mRNAs)) +
2*sd(log10(pData(sce)$Total_mRNAs)))
lower_bound <- 10^(mean(log10(pData(sce)$Total_mRNAs)) -
2*sd(log10(pData(sce)$Total_mRNAs)))
## 可視化
qplot(Total_mRNAs, data = pData(sce), geom = "density") +
geom_vline(xintercept = lower_bound) +
geom_vline(xintercept = upper_bound) + xlim(0, 6000)+
theme_classic()
可視化RPC(reads per cell) 分布
過濾不合格細(xì)胞
## 保留下通過標(biāo)準(zhǔn)的細(xì)胞
sce <- sce[,pData(sce)$Total_mRNAs > lower_bound &
pData(sce)$Total_mRNAs < upper_bound]
## 新的sce 數(shù)據(jù),評估細(xì)胞表達(dá)基因數(shù).
sce <- detectGenes(sce, min_expr = 0.1)
> sce
CellDataSet (storageMode: environment)
assayData: 20222 features, 5943 samples
element names: exprs
protocolData: none
phenoData
sampleNames: AAACCTGCAAGACGTG-300.1.1 AAACCTGGTGTGAATA-300.1.1
... TTTGTCAAGTACACCT-b02 (5943 total)
varLabels: cell n.umi ... Total_mRNAs (20 total)
varMetadata: labelDescription
featureData
featureNames: WBGene00010957 WBGene00010958 ... WBGene00007064
(20222 total)
fvarLabels: id gene_short_name num_cells_expressed
fvarMetadata: labelDescription
experimentData: use 'experimentData(object)'
Annotation:
細(xì)胞數(shù)目有6188 變成現(xiàn)在的5943個.
檢查過濾細(xì)胞的效果凌蔬，過濾后的值基本上是符合lognormal分布的露懒。
# Log-transform each value in the expression matrix.
L <- log(exprs(sce[expressed_genes,]))
# Standardize each gene, so that they are all on the same scale,
# Then melt the data with plyr so we can plot it easily
melted_dens_df <- melt(Matrix::t(scale(Matrix::t(L))))
# Plot the distribution of the standardized gene expression values.
qplot(value, geom = "density", data = melted_dens_df) +
stat_function(fun = dnorm, size = 0.5, color = 'red') +
xlab("Standardized log(counts)") +
ylab("Density")
理想的圖：
得到過濾后的矩陣，可以分別進(jìn)行聚類和擬時序分析.
3. 細(xì)胞分群/聚類
不使用marker 基因的細(xì)胞分類方法
使用函數(shù)clusterCells()砂心，根據(jù)整體的表達(dá)量對細(xì)胞進(jìn)行分組懈词。例如，細(xì)胞表達(dá)了大量的與成肌細(xì)胞相關(guān)的基因辩诞，但就是沒有成肌細(xì)胞的marker--MYF5 坎弯，我們依然可以判斷這個細(xì)胞屬于成肌細(xì)胞。
step1：dispersionTable()
首先就是判斷使用哪些基因進(jìn)行細(xì)胞分群。當(dāng)然荞怒，可以使用全部基因洒琢，但這會摻雜很多表達(dá)量不高而檢測不出來的基因，反而會增加噪音褐桌。挑有差異的衰抑，挑表達(dá)量不太低的
disp_table <- dispersionTable(sce)
unsup_clustering_genes <- subset(disp_table,
mean_expression >= 0.1) ## 數(shù)值因?qū)嶒?yàn)而不同
sce <- setOrderingFilter(sce, unsup_clustering_genes$gene_id)
plot_ordering_genes(sce)
其中, setOrderingFilter函數(shù)設(shè)置了將來會用于分群的細(xì)胞。這一步對于分群非常關(guān)鍵荧嵌，所以可以試著使用不同的方法來過濾基因呛踊，以期達(dá)到滿意的效果。
plot_ordering_genes根據(jù)這些基因的平均表達(dá)水平展示了基因的表達(dá)差異程度啦撮，紅線表示Monocle基于這種關(guān)系對分散的期望谭网。 上圖中實(shí)黑的為會使用的基因，灰色的是過濾掉的基因赃春。
step2：plot_pc_variance_explained() 耗時
然后選一下主成分
plot_pc_variance_explained(cds, return_all = F,max_components = 100) # norm_method='log'
上圖與Seurat中使用到的PCElbowPlot類似愉择，我們可以看到大約多少的PC將會被我們使用。 可以看到大約在15的時候就已經(jīng)到的平原區(qū)了织中。下面試著分群锥涕。
step3: 進(jìn)行降維和聚類
根據(jù)上面??的圖，選擇合適的主成分?jǐn)?shù)量（這個很主觀狭吼，可以多試幾次）层坠，這里選前15個成分（需要注意的 使用主成分會影響后面結(jié)果）
關(guān)于處理批次效應(yīng)：例如在芯片數(shù)據(jù)中經(jīng)常會利用SVA的combat函數(shù)。
磨平批次效應(yīng)實(shí)際上就是去掉各個組的前幾個主成分
# 進(jìn)行降維 刁笙，設(shè)置降維維度15.
cds <- reduceDimension(sce, max_components = 2, norm_method = 'log',
num_dim = 15, reduction_method = "tSNE", verbose = TRUE)
# Remove noise by PCA ...
# Reduce dimension by tSNE ...
### 如果需要去除批次效應(yīng)(batch 和 幾群細(xì)胞基因表達(dá)量數(shù)目)
# cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
# reduction_method = 'tSNE',
# residualModelFormulaStr = "~batch + #num_genes_expressed",
# verbose = T)
# 進(jìn)行聚類 （指定cluster為4的時候破花，只能畫出來3個,為什么？）
cds <- clusterCells(cds, num_clusters = 4)
# Distance cutoff calculated to 4.826512
## 查看聚類效果
head(cds@phenoData@data)
# 查看tSNE圖
plot_cell_clusters(cds,x=1,y=2, color_by = "Cluster")
# 畫具體基因的表達(dá)分布圖
#plot_cell_clusters(HSMM, 1, 2, color = "Cluster",
# + markers = c("CDC20", "GABPB2"))
4.篩選基因(用于擬分析時序)
4.1 一般過濾標(biāo)準(zhǔn)
篩選基因沒有一個固定標(biāo)準(zhǔn)疲吸，可以根據(jù)自己的需要靈活選擇座每，這里只展示其中1種 .
disp_table <- dispersionTable(sce)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
## 得到用于擬時序分析基因
ordering_genes <- row.names(unsup_clustering_genes)
4.2 使用差異表達(dá)基因
對于差異表達(dá)分析，可以像RNA-seq一樣磅氨，直接到不同屬性的樣品進(jìn)行分析
fullModelFormulaStr ： 分組進(jìn)行差異分析的變量.
我們可以看到尺栖，只要有不同的pData屬性，我們都可以做差異表達(dá)分析烦租。比如State：
## 差異分析非常耗時延赌，運(yùn)行時間在幾分鐘至十幾分鐘不等
diff_test_res <- differentialGeneTest(sce[expressed_genes, ],
fullModelFormulaStr = "Cluster")
sig_genes <- subset(diff_test_res, qval < 0.1)
sig_genes[1:6,c("gene_short_name", "pval", "qval")]
## 得到用于擬時序分析基因
ordering_genes <- row.names(sig_genes)
下面是可選部分??
再比如 Pseudotime
diff_test_res <- differentialGeneTest(sce[expressed_genes, ],
fullModelFormulaStr = "~sm.ns(Pseudotime)")
sig_genes <- subset(diff_test_res, qval < 0.1)
sig_genes[1:6,c("gene_short_name", "pval", "qval")]
作圖（注意要將基因名變成character）
cg=as.character(head(sig_genes$gene_short_name))
# 普通圖
plot_genes_jitter(sce[cg,],
grouping = "State", ncol= 2)
# 還能上色
plot_genes_jitter(sce[cg,],
grouping = "State",
color_by = "State",
nrow= 3,
ncol = NULL )
還可以繪制熱圖
sig_genes <- sig_genes[order(sig_genes$qval), ]
plot_pseudotime_heatmap(sce[row.names(sig_genes)[1:20],],
num_clusters = 3,
cores = 1,
show_rownames = T)
得到類似圖
BEAM分析
BEAM分析的目的是比較分枝點(diǎn)與分枝末端的差異。
BEAM_res <- BEAM(sce[expressed_genes, ], branch_point = 1, cores = 4)
BEAM_res <- BEAM_res[order(BEAM_res$qval),]
BEAM_res <- BEAM_res[,c("gene_short_name", "pval", "qval")]
plot_genes_branched_heatmap(sce[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(sce[rownames(subset(BEAM_res, qval < 1e-8)), ],
branch_point = 1, color_by = "orig.ident",
ncol = 3)
上面是可選部分 ??.
5.推斷發(fā)育軌跡
用上面過濾的基因進(jìn)行降維叉橱， 以方便下一步將細(xì)胞分布到不同的狀態(tài)之間挫以，可視化擬時序得分.
5.1: choosing genes that define progress
5.2: reducing the dimensionality of the data
5.3: ordering the cells in pseudotime
5.1. 選合適基因
得到上面兩種過濾的基因集：unsup_clustering_genes 或者 diff_test_res，選取合適基因進(jìn)行排序.
## 以差異基因?yàn)槔?## ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)
5.2 降維
默認(rèn)使用DDRTree的方法進(jìn)行降維分析窃祝，代碼如下
cds <- reduceDimension(
cds,
max_components = 2,
method = 'DDRTree')
5.3 細(xì)胞排序
然后就是細(xì)胞排序 掐松，代碼如下：
cds <- orderCells(cds)
運(yùn)行之后，對于每個細(xì)胞，會給出一個pseudotime值大磺，示意如下
可視化代碼：
- 默認(rèn)用"state" 填充抡句。應(yīng)該是cluster 含義.
plot_cell_trajectory(sce)
- 降維之后，在二維空間展示細(xì)胞pseudotime的分布杠愧，可以看到是一個樹狀結(jié)構(gòu)待榔，除了上述方法外，還可以根據(jù)pseudotime的值給細(xì)胞賦顏色流济，代碼如下
plot_cell_trajectory(sce, color_by = "Pseudotime")
其實(shí)就是將fData中對應(yīng)的列設(shè)置為顏色锐锣，如果想要觀察不同細(xì)胞亞型的分布，可以在fData中新增一列細(xì)胞對應(yīng)的cluster ID, 然后用這一類來設(shè)置顏色绳瘟。
對于pseudotime分析雕憔，我們需要明白它的基本輸入就是一張基因在細(xì)胞中表達(dá)量的表格，與細(xì)胞的聚類結(jié)果無關(guān)糖声，只不過在可視化的時候根據(jù)聚類的結(jié)果填充了顏色而已斤彼。
另外，plot_genes_in_pseudotime 可以對基因在不同細(xì)胞中的表達(dá)量變化進(jìn)行繪圖
plot_genes_in_pseudotime(cds[cg,],
color_by = "State")
小結(jié)：
【W(wǎng)orkflow以及與Seurat的異同】
-
①創(chuàng)建CellDataSet對象（下簡稱CDS對象）
若要創(chuàng)建新的CDS對象蘸泻，則需要整理出3個輸入文件（基因-細(xì)胞表達(dá)矩陣畅卓、細(xì)胞-細(xì)胞特征注釋矩陣、基因-基因特征注釋矩陣）蟋恬，但方便的是，Monocle支持從Seurat中直接導(dǎo)入對象趁冈，通過importCDS命令實(shí)現(xiàn)歼争。
在創(chuàng)建之后，也會進(jìn)行數(shù)據(jù)過濾和標(biāo)準(zhǔn)化渗勘，不同的是Seurat是基于作圖的方法去篩選掉異常的細(xì)胞沐绒，而Monocle的過濾方法則是根據(jù)表達(dá)數(shù)據(jù)的數(shù)學(xué)關(guān)系來實(shí)現(xiàn)。
上限：image.png
其中X為基因表達(dá)總數(shù), n 為細(xì)胞數(shù)旺坠，sd為表達(dá)水平方差
大概就是以一個大致的細(xì)胞表達(dá)水平為基準(zhǔn)乔遮，表達(dá)量太高的跟太低的去除掉。
②通過已知的Marker基因分類細(xì)胞
方法一：通過一些現(xiàn)有的生物/醫(yī)學(xué)知識手動賦予它們細(xì)胞名取刃，將這些細(xì)胞分為幾大類蹋肮，然后再聚類細(xì)胞。
方法二：與Seurat包的分類細(xì)胞方法類似璧疗，篩選出差異表達(dá)基因用于聚類坯辩，然后再根據(jù)每一個cluster的marker賦予細(xì)胞名。③聚類細(xì)胞
采用的也是t-SNE算法崩侠。-
④將細(xì)胞按照偽時間的順序排列在軌跡上
Step1：選擇輸入基因用于機(jī)器學(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ù)monocle網(wǎng)站說是這樣）Step3：將細(xì)胞按照偽時間排序
這個過程稱為manifold learning（流形學(xué)習(xí)）句灌。Monocle利用軌跡來描述細(xì)胞如何從一個狀態(tài)轉(zhuǎn)換到另一個狀態(tài)。得到的是一個樹狀圖阵赠，樹的根部或樹干表示的是細(xì)胞最初的狀態(tài)（如果有的話）涯塔，隨著細(xì)胞的變化就沿著分叉樹枝發(fā)展。一個細(xì)胞的偽時間值（pseudotime value）為它的位置沿著樹枝到根部的距離清蚀。
分析scRNA-seq的數(shù)據(jù)的關(guān)鍵在于如何對細(xì)胞進(jìn)行cluster匕荸。這其中有很多的算法，而之后的降維分析比如tSNE其實(shí)主要還是為了數(shù)據(jù)圖形化顯示方便枷邪。在細(xì)胞分群之后榛搔，差異表達(dá)分析其實(shí)與第三章的RNA-seq并無二致，我們只需要對需要比較的因素做到心中有數(shù)即可东揣。
