這里是佳奧！我們繼續(xù)學(xué)習(xí)最重要的三個(gè)R包。

1 scater Package(calculateQCMetrics已停用，用addPerCellQC替代)

加載R包

rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(scater))
suppressMessages(library(scRNAseq))

創(chuàng)建測(cè)試數(shù)據(jù)集

library(scRNAseq)
## ----- Load Example Data -----
fluidigm<-ReprocessedFluidigmData()
data(fluidigm)

# Set assay to RSEM estimated counts
assay(fluidigm) <- assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4]
table(rowSums(ct)==0)
# 這里使用原始表達(dá)矩陣，所以有很多基因在所有細(xì)胞均無(wú)表達(dá)量，即表現(xiàn)為沒(méi)有被檢測(cè)到，這樣的基因是需要過(guò)濾掉的。

pheno_data <- as.data.frame(colData(fluidigm))
## 這里需要把Pollen的表達(dá)矩陣做成我們的 scater 要求的對(duì)象
#data("sc_example_counts")
#data("sc_example_cell_info") 
# 也可以嘗試該R包自帶的數(shù)據(jù)集。
# https://bioconductor.org/packages/release/bioc/vignettes/scater/inst/doc/vignette-intro.R
sce <- SingleCellExperiment(
    assays = list(counts = ct), 
    colData = pheno_data
    )
sce
# 后面所有的分析都是基于 sce 這個(gè)變量
# 是一個(gè) SingleCellExperiment 對(duì)象，被很多單細(xì)胞R包采用。
> sce
class: SingleCellExperiment 
dim: 26255 130 
metadata(0):
assays(1): counts
rownames(26255): A1BG A1BG-AS1 ... ZZEF1 ZZZ3
rowData names(0):
colnames(130): SRR1275356 SRR1274090 ... SRR1275366 SRR1275261
colData names(28): NREADS NALIGNED ... Cluster1 Cluster2
reducedDimNames(0):
mainExpName: NULL
altExpNames(0):

探究SingleCellExperiment對(duì)象——sce

##sce是S4對(duì)象，與getGEO返回的對(duì)象相似

> str(sce,max.level=2)

##羅列了所有屬性以及下一級(jí)菜單的屬性
> str(sce)
Formal class 'SingleCellExperiment' [package "SingleCellExperiment"] with 9 slots
  ..@ int_elementMetadata:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. ..@ rownames       : NULL
  .. .. ..@ nrows          : int 26255
  .. .. ..@ elementType    : chr "ANY"
  .. .. ..@ elementMetadata: NULL
  .. .. ..@ metadata       : list()
  .. .. ..@ listData       :List of 1
  .. .. .. ..$ rowPairs:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
  .. .. .. .. .. ..@ rownames       : NULL
  .. .. .. .. .. ..@ nrows          : int 26255
  .. .. .. .. .. ..@ elementType    : chr "ANY"
  .. .. .. .. .. ..@ elementMetadata: NULL
  .. .. .. .. .. ..@ metadata       : list()
  .. .. .. .. .. ..@ listData       : Named list()
  ..@ int_colData        :Formal class 'DFrame' [package "S4Vectors"] with 6 slots
 略

一些質(zhì)量控制

這里稍作講解，包括基因?qū)用婕凹?xì)胞層面，但是并不設(shè)計(jì)實(shí)驗(yàn)本身引入的干擾因素。

首先看基因?qū)用娴倪^(guò)濾

使用 calculateQCMetrics函數(shù)作用于 sce 那個(gè)單細(xì)胞數(shù)據(jù)對(duì)象后，就可以 用 rowData(object) 可以查看各個(gè)基因統(tǒng)計(jì)指標(biāo)：

• mean_counts: 平均表達(dá)量
• log10_mean_counts: 歸一化 log10-scale.
• pct_dropout_by_counts: 該基因丟失率。
• n_cells_by_counts: 多少個(gè)細(xì)胞表達(dá)了該基因

上面那些指標(biāo)可以用來(lái)過(guò)濾，也可以自己重新再次計(jì)算一下那些統(tǒng)計(jì)學(xué)指標(biāo)。

主要是過(guò)濾低表達(dá)量基因，還有 線粒體基因 和 ERCC spike-ins 的控制。

exprs(sce) <- log2( calculateCPM( sce ) + 1)
## 只有運(yùn)行了下面的函數(shù)后才有各式各樣的過(guò)濾指標(biāo)
genes <- rownames(rowData(sce))
genes[grepl('^MT-',genes)]
genes[grepl('^ERCC-',genes)]
# 比較不幸的是，這個(gè)測(cè)試數(shù)據(jù)里面既沒(méi)有線粒體基因，有沒(méi)有ERCC序列。
# > genes[grepl('^MT-',genes)]
# character(0)
# > genes[grepl('^ERCC-',genes)]
# character(0)

sce <- calculateQCMetrics(sce, 
        feature_controls = list(ERCC = grep('^ERCC',genes)))

sce <- addPerCellQC(sce, 
        subsets = list(ERCC = grep('^ERCC',genes)))

# 也沒(méi)有定義啥細(xì)胞是屬于control組，所以這里并不需要完全follow教程
# example_sce <- calculateQCMetrics(example_sce, 
# feature_controls = list(ERCC = 1:20, mito = 500:1000),
# cell_controls = list(empty = 1:5, damaged = 31:40))

查看信息

tmp <- as.data.frame(rowData(sce))
colnames(tmp)
head(tmp)

目前只過(guò)濾那些在所有細(xì)胞都沒(méi)有表達(dá)的基因, 但是這個(gè)過(guò)濾條件可以自行調(diào)整?？梢钥吹交驍?shù)量大幅度減少。

keep_feature <- rowSums(exprs(sce) > 0) > 0
table(keep_feature)
sce <- sce[keep_feature,]
sce

然后看細(xì)胞層面的過(guò)濾

用 colData(object) 可以查看各個(gè)樣本統(tǒng)計(jì)情況

• total_counts: total number of counts for the cell (aka ‘library size’)

• log10_total_counts: total_counts on the log10-scale

• total_features: the number of features for the cell that have expression above the detection limit (default detection limit is zero)

每個(gè)數(shù)據(jù)集都有適合自己的過(guò)濾閾值，不要局限于教程

tmp <- as.data.frame(colData(sce))
colnames(tmp) 
# 可以發(fā)現(xiàn)細(xì)胞質(zhì)控屬性非常多，可以說(shuō)是實(shí)用至極。
## 比如看每個(gè)樣本測(cè)到的基因數(shù)量
tf <- sce$total_features_by_counts 
boxplot(tf)
fivenum(tf)

還有很多其它指標(biāo)，比如哪些細(xì)胞的ERCC含量過(guò)高，或者某些其高表達(dá)量基因占比太高，等等。

一些可視化

可視化函數(shù)非常多，所以這個(gè)R包才會(huì)被引用那么廣泛，里面包括如下：

首先是一些細(xì)胞距離情況

• plotPCA: produce a principal components plot for the cells.
• plotTSNE: produce a t-distributed stochastic neighbour embedding (reduced dimension) plot for the cells.
• plotDiffusionMap: produce a diffusion map (reduced dimension) plot for the cells.
• plotMDS: produce a multi-dimensional scaling plot for the cells.
• plotReducedDim: plot a reduced-dimension representation of the cells.

然后是一些表達(dá)量相關(guān)的:

• plotExpression: plot expression levels for a defined set of features.
• plotPlatePosition: plot cells in their position on a plate, coloured by cell metadata and QC metrics or feature expression level.
• plotColData: plot cell metadata and QC metrics.
• plotRowData: plot feature metadata and QC metrics.

首先可視化表達(dá)量

可以添加一些細(xì)胞屬性，展示如下：

# 挑選一些細(xì)胞屬性（臨床信息） 來(lái)進(jìn)行可視化展示。
colnames(tmp)[25:28]
# 主要展現(xiàn)下面這些基因
rownames(sce)[1:6]

展示一些基因在不同細(xì)胞分類的表達(dá)區(qū)別

plotExpression(sce, rownames(sce)[1:6],
               x = "Biological_Condition", 
               exprs_values = "logcounts") 

散點(diǎn)圖，展示兩個(gè)基因的相關(guān)性，這里批量展示6對(duì)相關(guān)性

plotExpression(sce, rownames(sce)[1:6],
               x = rownames(sce)[1])

還可以更復(fù)雜

plotExpression(sce, rownames(sce)[1:6],
               colour_by = "Biological_Condition", 
               shape_by = "Cluster1", 
               size_by = rownames(sce)[1])

還可以可視化細(xì)胞距離分布

這里可以學(xué)習(xí)PCA或者tSNE的使用

sce <- runPCA(sce)
# 這里并沒(méi)有進(jìn)行任何基因的挑選，就直接進(jìn)行了PCA，與 Seurat包不一樣。
reducedDimNames(sce)

PCA分布圖上面添加臨床信息，同樣的可以看到GW16混在了GW21里面

plotReducedDim(sce, use_dimred = "PCA", 
               colour_by = "Biological_Condition" )

PCA分布圖上面添加表達(dá)量信息, 但這個(gè)基本上不怎么用

plotReducedDim(sce, use_dimred = "PCA", 
               colour_by = rownames(sce)[1], 
               size_by = rownames(sce)[2])

最原始的

plotPCA(sce)

其它降維算法結(jié)果的可視化

## ----plot-pca-feature-controls-----
# 這里在真實(shí)場(chǎng)景中應(yīng)用比較多。
sce2 <- runPCA(sce, 
               feature_set = rowData(sce)$is_feature_control)
plotPCA(sce2)

僅僅是選取前20個(gè)PC, 可以看到繪圖時(shí)候細(xì)胞分布并沒(méi)有太大區(qū)別

## ----plot-pca-4comp-colby-shapeby-----
sce <- runPCA(sce, ncomponents=20)
plotPCA(sce,  colour_by = "Biological_Condition" )

但是可以挑選指定的PC來(lái)可視化, 這里選擇的是第4個(gè)PC

plotPCA(sce, ncomponents = 4,  
        colour_by = "Biological_Condition" )

意義不大的展示

## ----plot-pca-4comp-colby-sizeby-exprs-----
plotPCA(sce,
        colour_by = rownames(sce)[1], 
        size_by = rownames(sce)[2])

tSNE可視化

## ----plot-tsne-1comp-colby-sizeby-exprs------
# Perplexity of 10 just chosen here arbitrarily. 
set.seed(1000)
# 這里的這個(gè) perplexity 參數(shù)很重要
sce <- runTSNE(sce, perplexity=30)
plotTSNE(sce, 
         colour_by = "Biological_Condition" )

## ----plot-tsne-from-pca------
set.seed(1000)
sce <- runTSNE(sce, perplexity=10, use_dimred="PCA", n_dimred=10)
plotTSNE(sce,  
         colour_by = "Biological_Condition" )

這一步會(huì)調(diào)用destiny

## ----plot-difmap-1comp-colby-sizeby-exprs-------
sce <- runDiffusionMap(sce)
plotDiffusionMap(sce,  
                 colour_by = "Biological_Condition" )

關(guān)于SC3包

雖然 SC3包 跟 scater 聯(lián)系很緊密，但畢竟不屬于其知識(shí)點(diǎn)，這里就簡(jiǎn)單運(yùn)行即可。

預(yù)估亞群, 這里顯示有24個(gè), 比較多

library(SC3) # BiocManager::install('SC3')
sce <- sc3_estimate_k(sce)
metadata(sce)$sc3$k_estimation
rowData(sce)$feature_symbol <- rownames(rowData(sce))

一步運(yùn)行sc3的所有分析, 相當(dāng)耗費(fèi)時(shí)間

這里kn表示的預(yù)估聚類數(shù), 考慮到數(shù)據(jù)集是已知的，我們強(qiáng)行設(shè)置為4組, 具體數(shù)據(jù)要具體考慮。

# 耗費(fèi)時(shí)間
kn <- 4 ## 這里可以選擇 3:5 看多種分類結(jié)果。
sc3_cluster <- "sc3_4_clusters"
Sys.time()
sce <- sc3(sce, ks = kn, biology = TRUE)
Sys.time()

可視化展示部分, kn就是聚類數(shù)

熱圖: 比較先驗(yàn)分類和SC3的聚類的一致性

sc3_plot_consensus(sce, k = kn, show_pdata = c("Biological_Condition",sc3_cluster))

展示表達(dá)量信息

sc3_plot_expression(sce, k = kn, show_pdata =  c("Biological_Condition",sc3_cluster))

展示可能的標(biāo)記基因

sc3_plot_markers(sce, k = kn, show_pdata =  c("Biological_Condition",sc3_cluster))

在PCA上展示SC3的聚類結(jié)果

plotPCA(sce, colour_by =  sc3_cluster )
# sc3_interactive(sce)

顯示運(yùn)行環(huán)境

sessionInfo()

2 Seurat Package（最新4.X版本改了很多函數(shù)，需要重新查看參考文檔）

參考文檔

https://satijalab.org/seurat/articles/pbmc3k_tutorial.html

加載R包

rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(Seurat))
# 加載R包，請(qǐng)注意R包版本，可能會(huì)有莫名其妙的版本錯(cuò)誤
# 單細(xì)胞轉(zhuǎn)錄組領(lǐng)域發(fā)展太快，不同版本的 同一個(gè)R包差異很大。

單細(xì)胞轉(zhuǎn)錄組領(lǐng)域發(fā)展太快，不同版本的 同一個(gè)R包差異很大，這個(gè)Seurat包也不例外，大部分的函數(shù)都改了，還專門(mén)有一個(gè)對(duì)照表格供大家學(xué)習(xí)：

 https://satijalab.org/seurat/essential_commands.html 

也有新的基于3X的教程：（但是現(xiàn)在的版本已經(jīng)是4.1.1了......）

https://satijalab.org/seurat/pbmc3k_tutorial.html 

創(chuàng)建測(cè)試數(shù)據(jù)集

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)
fluidigm<-ReprocessedFluidigmData()
# Set assay to RSEM estimated counts
assay(fluidigm) <- assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4]

簡(jiǎn)單看看表達(dá)矩陣的性質(zhì)，主要是基因數(shù)量，細(xì)胞數(shù)量;以及每個(gè)細(xì)胞表達(dá)基因的數(shù)量，和每個(gè)基因在多少個(gè)細(xì)胞里面表達(dá)。

fivenum(apply(ct,1,function(x) sum(x>0) ))
boxplot(apply(ct,1,function(x) sum(x>0) ))
fivenum(apply(ct,2,function(x) sum(x>0) ))
hist(apply(ct,2,function(x) sum(x>0) ))

names(metadata(fluidigm))
meta <- as.data.frame(colData(fluidigm))
counts <- ct

檢測(cè)了 counts 和 meta 兩個(gè)變量，后面需要使用

identical(rownames(meta),colnames(counts))

這里需要把Pollen的表達(dá)矩陣做成我們的Seurat要求的對(duì)象

Pollen <- CreateSeuratObject(counts = counts, 
                             meta.data =meta,
                             min.cells = 3, 
                             min.genes = 200, 
                             project = "Pollen")
Pollen
## 后續(xù)所有的分析都基于這個(gè) Pollen 變量，是一個(gè)對(duì)象
# An object of class Seurat in project Pollen 

> Pollen
An object of class Seurat 
14656 features across 130 samples within 1 assay 
Active assay: RNA (14656 features, 0 variable features)

檢查表達(dá)矩陣

將Pollen賦值給sce，目的是代碼復(fù)用

sce <- Pollen

為原信息增加線粒體基因的比例，如果線粒體基因所占比例過(guò)高，意味著這可能是死細(xì)胞

由于包的更新，這里運(yùn)行會(huì)提示sce沒(méi)有名為raw.data的對(duì)象，需要進(jìn)行修改。

mito.genes <- grep(pattern = "^MT-", x = rownames(x = sce@meta.data), value = TRUE)
# 恰好這個(gè)例子的表達(dá)矩陣?yán)锩鏇](méi)有線粒體基因
percent.mito <- Matrix::colSums(sce@raw.data[mito.genes, ]) / Matrix::colSums(sce@raw.data)
## 也可以加入很多其它屬性，比如 ERCC 等。

# AddMetaData adds columns to object@meta.data, and is a great place to stash QC stats
sce <- AddMetaData(object = sce, metadata = percent.mito,
                   col.name = "percent.mito")

這里繪圖，可以指定分組，前提是這個(gè)分組變量存在于meta信息里面

這里的例子是：'Biological_Condition'

VlnPlot(object = sce, features.plot = c("nGene", "nUMI", "percent.mito"), group.by = 'Biological_Condition', nCol = 3)
GenePlot(object = sce, gene1 = "nUMI", gene2 = "nGene")

可以看看高表達(dá)量基因是哪些

tail(sort(Matrix::rowSums(sce@raw.data)))
## 散點(diǎn)圖可視化任意兩個(gè)基因的一些屬性（通常是細(xì)胞的度量）
GenePlot(object = sce, gene1 = "SOX11", gene2 = "EEF1A1")
# 散點(diǎn)圖可視化任意兩個(gè)細(xì)胞的一些屬性（通常是基因的度量）

CellPlot(sce,sce@cell.names[3],sce@cell.names[4],do.ident = FALSE)

表達(dá)矩陣的歸一化

起初sce對(duì)象里面的data就是原始表達(dá)矩陣

identical(sce@raw.data,sce@data)
sce <- NormalizeData(object = sce, 
                     normalization.method = "LogNormalize", 
                     scale.factor = 10000,
                     display.progress = F)

經(jīng)過(guò)了歸一化，sce對(duì)象里面的data被改變。

identical(sce@raw.data,sce@data)

尋找波動(dòng)比較明顯的基因，后續(xù)用這些基因而非全部基因進(jìn)行分析，主要為了降低計(jì)算量。

FindVariableGenes函數(shù)當(dāng)前版本無(wú)法調(diào)用

sce <- FindVariableGenes(object = sce, mean.function = ExpMean, dispersion.function = LogVMR, 
                         x.low.cutoff = 0.0125, 
                         x.high.cutoff = 3, 
                         y.cutoff = 0.5)
# 通過(guò)調(diào)整參數(shù)可以得到不同數(shù)量的 var.genes
length(sce@var.genes)

對(duì)歸一化后的矩陣進(jìn)行去除混雜因素和降維

對(duì)矩陣進(jìn)行回歸建模，以及scale,主要是為了去除一些文庫(kù)大小，線粒體基因含量，ERCC含量等因素。

sce <- ScaleData(object = sce, 
                 vars.to.regress = c("nUMI"),##去除文庫(kù)大小影響
                 display.progress = F)

現(xiàn)在sce對(duì)象的 sce@scale.data 也有了數(shù)值

運(yùn)行PCA進(jìn)行線性降維，這里僅僅是挑選高變化的基因組成的表達(dá)矩陣進(jìn)行PCA分析。

sce <- RunPCA(object = sce, 
              pc.genes = sce@var.genes, 
              do.print = TRUE, 
              pcs.print = 1:5, 
              genes.print = 5)
sce@dr

這樣就能拿到PC的基因的重要性占比情況。

tmp <- sce@dr$pca@gene.loadings
VizPCA( sce, pcs.use = 1:2)

PCAPlot(sce, dim.1 = 1, dim.2 = 2,
        group.by = 'Biological_Condition')
sce <- ProjectPCA(sce, do.print = FALSE)

因?yàn)榧?xì)胞數(shù)量不多，所以可以全部畫(huà)出來(lái)

PCHeatmap(object = sce, 
          pc.use = 1, 
          cells.use = ncol(sce@data), 
          do.balanced = TRUE, 
          label.columns = FALSE)

PCHeatmap(object = sce, 
          pc.use = 1:10, 
          cells.use = ncol(sce@data), 
          do.balanced = TRUE, 
          label.columns = FALSE)

基于PCA結(jié)果看看細(xì)胞如何分群

重點(diǎn): 需要搞懂這里的 resolution 參數(shù)，而且降維算法可以選PCA或者ICA ， 分群算法也可以選擇。

sce <- FindClusters(object = sce, 
                    reduction.type = "pca", 
                    dims.use = 1:10, force.recalc = T,
                    resolution = 0.9, print.output = 0,
                    save.SNN = TRUE)
PrintFindClustersParams(sce)
table(sce@meta.data$res.0.9)

看單細(xì)胞分群后的tSNE圖

跟前面的 RunPCA 函數(shù)功能差不多，都是為了降維。這里為了節(jié)省計(jì)算量，首先使用PCA的線性降維結(jié)果，再進(jìn)行tSNE

sce <- RunTSNE(object = sce, 
               dims.use = 1:10, 
               do.fast = TRUE, 
               perplexity=10)
# note that you can set do.label=T to help label individual clusters
TSNEPlot(object = sce)

可以看到，雖然說(shuō)有4類細(xì)胞，但是 GW16和GW21沒(méi)有區(qū)分開(kāi)來(lái)，需要探索參數(shù)。

table(meta$Biological_Condition)
table(meta$Biological_Condition,sce@meta.data$res.0.9)
TSNEPlot(object = sce,group.by = 'Biological_Condition')

對(duì)每個(gè)分類都尋找其marker基因

# 下面的代碼是需要適應(yīng)性修改，因?yàn)椴煌臄?shù)據(jù)集分組不一樣，本次是3組，所以演示3組后的代碼。
markers_df <- FindMarkers(object = sce, ident.1 = 0, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features.plot =markers_genes, use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = sce, 
            features.plot =markers_genes, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne")

markers_df <- FindMarkers(object = sce, ident.1 = 1, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features.plot =markers_genes, use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = sce, 
            features.plot =markers_genes, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne")

markers_df <- FindMarkers(object = sce, ident.1 = 2, min.pct = 0.25)
print(x = head(markers_df))
markers_genes =  rownames(head(x = markers_df, n = 5))
VlnPlot(object = sce, features.plot =markers_genes, use.raw = TRUE, y.log = TRUE)
FeaturePlot(object = sce, 
            features.plot =markers_genes, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne") 

展現(xiàn)各個(gè)分類的marker基因的表達(dá)情況

sce.markers <- FindAllMarkers(object = sce, only.pos = TRUE, min.pct = 0.25, 
                              thresh.use = 0.25)
DT::datatable(sce.markers)
library(dplyr)
sce.markers %>% group_by(cluster) %>% top_n(2, avg_logFC)
top10 <- sce.markers %>% group_by(cluster) %>% top_n(10, avg_logFC)

# setting slim.col.label to TRUE will print just the cluster IDS instead of# every cell name
DoHeatmap(object = sce, genes.use = top10$gene, slim.col.label = TRUE, remove.key = TRUE)

FeaturePlot(object = sce, 
            features.plot =top10$gene, 
            cols.use = c("grey", "blue"), 
            reduction.use = "tsne")

Seurat總結(jié)

counts矩陣進(jìn)來(lái)后被包裝為對(duì)象，方便操作。

然后一定要經(jīng)過(guò) NormalizeData 和 ScaleData 的操作

函數(shù) FindVariableGenes 可以挑選適合進(jìn)行下游分析的基因集。

函數(shù) RunPCA 和 RunTSNE 進(jìn)行降維

函數(shù) FindClusters 直接就分群了，非常方便
函數(shù) FindAllMarkers 可以對(duì)分群后各個(gè)亞群找標(biāo)志基因。

函數(shù) FeaturePlot 可以展示不同基因在所有細(xì)胞的表達(dá)量
函數(shù) VlnPlot 可以展示不同基因在不同分群的表達(dá)量差異情況
函數(shù) DoHeatmap 可以選定基因集后繪制熱圖

顯示運(yùn)行環(huán)境

sessionInfo()

3 monocle Package

引言

教程，當(dāng)然是以官網(wǎng)為主：

http://cole-trapnell-lab.github.io/monocle-release/docs/
http://cole-trapnell-lab.github.io/monocle-release/tutorials/

if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
if (!requireNamespace("monocle"))
    BiocManager::install("monocle")

加載R包

rm(list = ls()) # clear the environment
#load all the necessary libraries
options(warn=-1) # turn off warning message globally
suppressMessages(library(monocle))

創(chuàng)建數(shù)據(jù)集

后續(xù)分析的前提就是將數(shù)據(jù)構(gòu)建成monocle需要的對(duì)象。

因此這里先介紹一下monocle需要的用來(lái)構(gòu)建 CellDataSet 對(duì)象的三個(gè)數(shù)據(jù)集

• 表達(dá)量矩陣exprs:數(shù)值矩陣 行名是基因, 列名是細(xì)胞編號(hào).
• 細(xì)胞的表型信息phenoData: 第一列是細(xì)胞編號(hào)，其他列是細(xì)胞的相關(guān)信息
• 基因注釋featureData: 第一列是基因編號(hào), 其他列是基因?qū)?yīng)的信息

并且這三個(gè)數(shù)據(jù)集要滿足如下要求:

表達(dá)量矩陣必須：

• 保證它的列數(shù)等于phenoData的行數(shù)
• 保證它的行數(shù)等于featureData的行數(shù)

而且

• phenoData的行名需要和表達(dá)矩陣的列名匹配
• featureData和表達(dá)矩陣的行名要匹配
• featureData至少要有一列"gene_short_name", 就是基因的symbol

如下是幾個(gè)例子

測(cè)試數(shù)據(jù)集

該R包團(tuán)隊(duì)給出的例子有點(diǎn)可怕。

需要使用網(wǎng)絡(luò)公共數(shù)據(jù)，取決于網(wǎng)速，是浙江大學(xué)郭老師的40萬(wàn)小鼠單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)。

第一個(gè)rds文件大小接近3G，太難下載了，即使僥幸下載，一般人的電腦也沒(méi)辦法處理它。

下載地址:

http://trapnell-lab.gs.washington.edu/public_share/MCA_merged_mat.rds
http://trapnell-lab.gs.washington.edu/public_share/MCA_All-batch-removed-assignments.csv

第二個(gè)csv文件也有 26.5M

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE72857

MCA <- readRDS("./MCA_merged_mat.rds")
cell.meta.data <- read.csv("./MCA_All-batch-removed-assignments.csv", row.names = 1)

39855 x 405191 sparse Matrix of class "dgCMatrix"

class(MCA)
## 這是一個(gè)對(duì)象，需要仔細(xì)理解，這里略

overlapping_cells <- intersect(row.names(cell.meta.data), colnames(MCA))
gene_ann <- data.frame(gene_short_name = row.names(MCA), row.names = row.names(MCA))

pd <- new("AnnotatedDataFrame",
          data=cell.meta.data[overlapping_cells, ])
fd <- new("AnnotatedDataFrame",
          data=gene_ann)

下面構(gòu)造了一個(gè)MCA數(shù)據(jù)集的子集，基于：overlapping_cells

MCA_cds <- newCellDataSet(
  MCA[, overlapping_cells], 
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=1)

save(MCA_cds,file = './MCA_cds_monocle_example.Rdata')   

但是數(shù)據(jù)集實(shí)在是太大，我這里只能載入前面步驟保存的小數(shù)據(jù)集。

## 這個(gè)文件很大，我放在了Rmd路徑下面。
load(file = './MCA_cds_monocle_example.Rdata')

MCA_cds
# 是一個(gè) CellDataSet 對(duì)象， 有著24萬(wàn)的單細(xì)胞，演示軟件用法實(shí)在是太難

這里需要學(xué)習(xí) CellDataSet 對(duì)象，就跟之前GEO視頻教程教大家的 ExpressionSet 對(duì)象類似，可以看到這里的數(shù)據(jù)集仍然是24萬(wàn)細(xì)胞。

scRNAseq R包中的數(shù)據(jù)集

我們選擇 scRNAseq 這個(gè)R包。

這個(gè)包內(nèi)置的是 Pollen et al. 2014 數(shù)據(jù)集，人類單細(xì)胞細(xì)胞，分成4類，分別是 pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，還有 “GW16” and “GW21” ，“GW21+3” 這種孕期細(xì)胞。

library(scRNAseq)
## ----- Load Example Data -----
data(fluidigm)
# Set assay to RSEM estimated counts
assay(fluidigm)  <-  assays(fluidigm)$rsem_counts
ct <- floor(assays(fluidigm)$rsem_counts)
ct[1:4,1:4] 
sample_ann <- as.data.frame(colData(fluidigm))

準(zhǔn)備Monocle對(duì)象需要的phenotype data 和 feature data 以及表達(dá)矩陣，從 scRNAseq 這個(gè)R包里面提取這三種數(shù)據(jù)

gene_ann <- data.frame(
  gene_short_name = row.names(ct), 
  row.names = row.names(ct)
)

pd <- new("AnnotatedDataFrame",
          data=sample_ann)
fd <- new("AnnotatedDataFrame",
          data=gene_ann)

構(gòu)建Monocle后續(xù)分析的所用對(duì)象，主要是根據(jù)包的說(shuō)明書(shū)，仔細(xì)探索其需要的構(gòu)建對(duì)象的必備元素

注意點(diǎn): 因?yàn)楸磉_(dá)矩陣是counts值，所以注意 expressionFamily 參數(shù)

sc_cds <- newCellDataSet(
  ct, 
  phenoData = pd,
  featureData =fd,
  expressionFamily = negbinomial.size(),
  lowerDetectionLimit=1)
sc_cds

下面是monocle對(duì)新構(gòu)建的CellDataSet 對(duì)象的標(biāo)準(zhǔn)操作

library(dplyr)
colnames(phenoData(sc_cds)@data)
## 為了電腦的健康，我這里選擇小數(shù)據(jù)集。 
sc_cds <- estimateSizeFactors(sc_cds)
sc_cds <- estimateDispersions(sc_cds)

本地加載RPKM數(shù)據(jù)

你也可以從本地RPKM值文件讀入R語(yǔ)言后構(gòu)造 CellDataSet 對(duì)象，下面是簡(jiǎn)單的例子：

#do not run
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")
# Once these tables are loaded, you can create the CellDataSet object like this:

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)

值得注意的是，因?yàn)閙onocle和前面我們講解的scater，還有Seurat，它們基于的對(duì)象都不一樣，所以monocle干脆提供了轉(zhuǎn)換函數(shù)：

# 加入你把上面的 HSMM 賦值給 lung ，然后使用函數(shù)進(jìn)行轉(zhuǎn)換：
lung  <-  HSMM
# To convert to Seurat object
lung_seurat <- exportCDS(lung, 'Seurat')

# To convert to SCESet
lung_SCESet <- exportCDS(lung, 'Scater')

直接讀取10X結(jié)果

因?yàn)?0X實(shí)在是太流行了，所以monocle跟Seurat一樣，也提供了直接讀取10X上游分析結(jié)果的接口函數(shù)（其實(shí)是使用另外一個(gè)R包），因?yàn)楸疚臄?shù)據(jù)來(lái)源于smart-seq2，所以并不演示下面的代碼：

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

fd <- fData(gbm)

# The number 2 is picked arbitrarily in the line below.
# Where "2" is placed you should place the column number that corresponds to your
# featureData's 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())

假設(shè)使用monocle3版本

具體教程以官網(wǎng)為主：

http://cole-trapnell-lab.github.io/monocle-release/monocle3/  

這里仍然是演示monocle2，所以下面代碼不要運(yùn)行，僅僅作為演示而已。

devtools::install_github("cole-trapnell-lab/monocle-release", ref="monocle3_alpha")
library(monocle)
cds <- sc_cds
cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
cds <- preprocessCDS(cds, num_dim = 20)
cds <- reduceDimension(cds, reduction_method = 'UMAP')
cds <- partitionCells(cds)
cds <- learnGraph(cds,  RGE_method = 'SimplePPT')
plot_cell_trajectory(cds,
                     color_by = "cell_type2") +
                     scale_color_manual(values = cell_type_color)

首先是質(zhì)控

這里通常也是對(duì)基因 和 細(xì)胞進(jìn)行質(zhì)控，質(zhì)控指標(biāo)需要根據(jù)項(xiàng)目來(lái)進(jìn)行具體探索，這里只是演示一下用法。

cds=sc_cds
cds
## 起初是： 26255 features, 130 samples 
cds <- detectGenes(cds, min_expr = 0.1)
print(head(fData(cds)))
expressed_genes <- row.names(subset(fData(cds),
                                    num_cells_expressed >= 5))
length(expressed_genes)
cds <- cds[expressed_genes,]
cds
# 過(guò)濾基因后是：assayData: 13385 features, 130 samples 
print(head(pData(cds)))
tmp=pData(cds)
fivenum(tmp[,1])
fivenum(tmp[,30])
# 這里并不需要過(guò)濾細(xì)胞，如果想過(guò)濾，就自己摸索閾值，然后挑選細(xì)胞即可。
valid_cells <- row.names(pData(cds) )
cds <- cds[,valid_cells]
cds 

聚類

單細(xì)胞轉(zhuǎn)錄組最重要的就是把細(xì)胞分群啦，這里可供選擇的算法非常多，我們首先演示PCA結(jié)果。

# 并不是所有的基因都有作用，所以先進(jìn)行挑選，合適的基因用來(lái)進(jìn)行聚類。
disp_table <- dispersionTable(cds)
unsup_clustering_genes <- subset(disp_table, mean_expression >= 0.1)
cds <- setOrderingFilter(cds, unsup_clustering_genes$gene_id)
plot_ordering_genes(cds) 
plot_pc_variance_explained(cds, return_all = F) # norm_method='log'

# 其中 num_dim 參數(shù)選擇基于上面的PCA圖
cds <- reduceDimension(cds, max_components = 2, num_dim = 6,
                        reduction_method = 'tSNE', verbose = T)
cds <- clusterCells(cds, num_clusters = 4) 
plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")
table(pData(cds)$Biological_Condition)

可以看到，GW21 也是被打散在其它分組里面。

值得注意的是：這里選擇不同的PC個(gè)數(shù)進(jìn)行后續(xù)分析，是會(huì)影響聚類結(jié)果的。具體可以參考：

https://davetang.org/muse/2017/10/01/getting-started-monocle/

排除干擾因素后聚類

跟前面的質(zhì)控步驟一樣，所謂的干擾因素，也是看自己對(duì)數(shù)據(jù)集的認(rèn)識(shí)情況來(lái)自己摸索的，比如我們這里

tmp=pData(cds)
fivenum(tmp[,1])
fivenum(tmp[,30])
colnames(tmp)
# 放在 residualModelFormulaStr 里面的是需要去除的
cds <- reduceDimension(cds, max_components = 2, num_dim = 2,
                        reduction_method = 'tSNE',
                        residualModelFormulaStr = "~Biological_Condition + num_genes_expressed",
                        verbose = T)
cds <- clusterCells(cds, num_clusters = 2)
plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")

## 上面去除了Biological_Condition，所以接下來(lái)聚類它們就被打散了。

cds <- reduceDimension(cds, max_components = 2, num_dim = 2,
                        reduction_method = 'tSNE',
                        residualModelFormulaStr = "~NREADS + num_genes_expressed",
                        verbose = T)
cds <- clusterCells(cds, num_clusters = 2)
plot_cell_clusters(cds, 1, 2, color = "Biological_Condition")

差異分析

這個(gè)是轉(zhuǎn)錄組數(shù)據(jù)的常規(guī)分析了，在單細(xì)胞轉(zhuǎn)錄組領(lǐng)域也是如此，monocle這個(gè)包提供 differentialGeneTest 函數(shù)來(lái)做差異分析，作用就是挑選那些在某些類別細(xì)胞里面高表達(dá)的基因，假設(shè)其為那一組細(xì)胞的marker基因。

在我們的例子里面可以是已知的細(xì)胞分類，或者是自己推斷的聚類結(jié)果。

這一步時(shí)間會(huì)比較久，也可以不給定分組信息，使用reduced model來(lái)計(jì)算各個(gè)基因的差異與否，是另外一種算法了。

Sys.time()
diff_test_res <- differentialGeneTest(cds,
                                      fullModelFormulaStr = "~Biological_Condition")
Sys.time()
# 可以看到運(yùn)行耗時(shí)

# Select genes that are significant at an FDR < 10%
sig_genes <- subset(diff_test_res, qval < 0.1)
head(sig_genes[,c("gene_short_name", "pval", "qval")] )

cg=as.character(head(sig_genes$gene_short_name))
plot_genes_jitter(cds[cg,], 
                  grouping = "Biological_Condition", ncol= 2)

下面是不同的可視化參數(shù)的效果

plot_genes_jitter(cds[cg,],
                  grouping = "Biological_Condition",
                  color_by = "Biological_Condition",
                  nrow= 3,
                  ncol = NULL )

這個(gè)時(shí)候可以考慮去選用多個(gè)差異分析R包來(lái)進(jìn)行不同的比較。

尋找marker基因

使用起來(lái)有點(diǎn)復(fù)雜，需要預(yù)先給定好分組情況，具體自行看說(shuō)明書(shū)。3個(gè)R包其實(shí)是共通的。

推斷發(fā)育軌跡

前面介紹的monocle的功能都只能說(shuō)是中規(guī)中矩，而這個(gè)推斷發(fā)育軌跡才是monocle的拿手好戲，也是它榮升為3大R包的核心競(jìng)爭(zhēng)力。

第一步: 挑選合適的基因。有多個(gè)方法，例如提供已知的基因集，這里選取統(tǒng)計(jì)學(xué)顯著的差異基因列表

ordering_genes <- row.names (subset(diff_test_res, qval < 0.01))
cds <- setOrderingFilter(cds, ordering_genes)
plot_ordering_genes(cds)
cds <- reduceDimension(cds, max_components = 2,
                            method = 'DDRTree')

第二步: 降維。降維的目的是為了更好的展示數(shù)據(jù)。函數(shù)里提供了很多種方法, 不同方法的最后展示的圖都不太一樣, 其中“DDRTree”是Monocle2使用的默認(rèn)方法

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

第三步: 對(duì)細(xì)胞進(jìn)行排序

cds <- orderCells(cds)

最后兩個(gè)可視化函數(shù)，對(duì)結(jié)果進(jìn)行可視化

plot_cell_trajectory(cds, color_by = "Biological_Condition")  

可以很明顯看到細(xì)胞的發(fā)育軌跡，正好對(duì)應(yīng) pluripotent stem cells 分化而成的 neural progenitor cells (“NPC”) ，還有 “GW16” and “GW21” ，“GW21+3” 這種孕期細(xì)胞。

plot_genes_in_pseudotime可以展現(xiàn)marker基因，本例子隨便選取了6個(gè)差異表達(dá)基因。

plot_genes_in_pseudotime(cds[cg,],
                         color_by = "Biological_Condition")

最開(kāi)始挑選合適基因，除了我們演示的找統(tǒng)計(jì)學(xué)顯著的差異表達(dá)基因這個(gè)方法外，還可以于已知的標(biāo)記基因，主要是基于生物學(xué)背景知識(shí)。

如果是已知基因列表，就需要自己讀取外包文件，導(dǎo)入R里面來(lái)分析。

顯示運(yùn)行環(huán)境

sessionInfo()

三個(gè)包介紹結(jié)束。

下一篇我們回到文章數(shù)據(jù)，使用上述介紹的R包進(jìn)行分析。

我們下一篇再見(jiàn)！