這次學(xué)習(xí)的是單細(xì)胞測序的第三個比較常用的包：Scater。本文主要參考教程：單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-15-利用scRNAseq包學(xué)習(xí)scater來走分析流程。
scater包的官方網(wǎng)站：here
這篇文章寫的些許混亂裳食，因?yàn)楣倬W(wǎng)里的很多函數(shù)與教程里的并不一樣，是因?yàn)橛行┖瘮?shù)是最新版scater新更新的疗琉，所以我在學(xué)習(xí)的時候有的部分按照官網(wǎng)的來雁竞，有的部分則按照教程里的來單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-15-利用scRNAseq包學(xué)習(xí)scater。但是實(shí)際上大同小異漾抬，主要的流程都是一樣的宿亡。
這個包像之前介紹的兩個包一樣，是為了分析單細(xì)胞測序分析而開發(fā)的纳令，也可以進(jìn)行質(zhì)控挽荠、可視化、下游分析前的預(yù)處理平绩。scater包的流程是這樣的：
簡單的來說圈匆，scater包可以做下面幾件事：
（1）QC
（2）從read數(shù)據(jù)里進(jìn)行轉(zhuǎn)錄本定量（根據(jù)偽時間排序）
（3）數(shù)據(jù)格式標(biāo)準(zhǔn)化
（4）可視化
（5）Seamless integration into the Bioconductor universe
（6）標(biāo)準(zhǔn)化數(shù)據(jù)
安裝scater包
>BiocManager::install("scater")
創(chuàng)建SingleCellExperiment對象
和之前介紹的Monocle2和Seurat包一樣，scater包也需要一個對象捏雌。
SingleCellExperiment這個對象跃赚，畫出來大概是這么個意思（圖片來自：單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-15-利用scRNAseq包學(xué)習(xí)scater)：
在導(dǎo)入數(shù)據(jù)之前，請將表達(dá)數(shù)據(jù)存為矩陣性湿。如果你的數(shù)據(jù)是10×Genomics的測序結(jié)果纬傲，可以用[DropletUtils]包的read10xCounts函數(shù)，它可以自動的生成SingleCellExperiment對象窘奏。這里我們用scater包自帶的示例數(shù)據(jù)來走一遍流程嘹锁。
首先，你要把數(shù)據(jù)導(dǎo)入進(jìn)來：
第一個要導(dǎo)入的是基因的表達(dá)量矩陣着裹，行是基因领猾，列是細(xì)胞：
> rm(list = ls())
> options(stringsAsFactors = F)
> data("sc_example_counts") #導(dǎo)入表達(dá)矩陣
> dim(sc_example_counts) #看一下表達(dá)矩陣?yán)锒嗌傩卸嗌倭?[1] 2000 40
> sc_example_counts[1:3,1:3] # 簡單的看一下表達(dá)矩陣長什么樣
Cell_001 Cell_002 Cell_003
Gene_0001 0 123 2
Gene_0002 575 65 3
Gene_0003 0 0 0
第二，導(dǎo)入細(xì)胞信息骇扇，注釋信息可以是細(xì)胞名稱摔竿，組織來源，收集樣品的時間點(diǎn)少孝，處理?xiàng)l件等等, 行是細(xì)胞继低，列是屬性：
> data("sc_example_cell_info")
> sc_example_cell_info[1:3,1:3] # 樣本注釋信息
Cell Mutation_Status Cell_Cycle
Cell_001 Cell_001 positive S
Cell_002 Cell_002 positive G0
Cell_003 Cell_003 negative G1
有了上述的文件，就可以構(gòu)建對象了：
> example_sce <- SingleCellExperiment(
assays = list(counts = sc_example_counts),
colData = sc_example_cell_info)
#瞅一眼對象里有啥
> example_sce
class: SingleCellExperiment
dim: 2000 40
metadata(0):
assays(1): counts
rownames(2000): Gene_0001 Gene_0002 ... Gene_1999 Gene_2000
rowData names(0):
colnames(40): Cell_001 Cell_002 ... Cell_039 Cell_040
colData names(4): Cell Mutation_Status Cell_Cycle Treatment
reducedDimNames(0):
spikeNames(0):
上面在構(gòu)建對象的時候稍走，我們將表達(dá)矩陣起了一個名字是“counts”袁翁，是為了方便有的時候你如果需要把矩陣?yán)锏腸ounts提取出來柴底，用counts函數(shù)比較方便：
> str(counts(example_sce))
int [1:2000, 1:40] 0 575 0 0 0 0 0 0 416 12 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:2000] "Gene_0001" "Gene_0002" "Gene_0003" "Gene_0004" ...
..$ : chr [1:40] "Cell_001" "Cell_002" "Cell_003" "Cell_004" ...
你可以通過下面的方法來新增一列meatdata：
> example_sce$whee <- sample(LETTERS, ncol(example_sce), replace=TRUE)
> colData(example_sce)
DataFrame with 40 rows and 5 columns
Cell Mutation_Status Cell_Cycle Treatment whee
<character> <character> <character> <character> <character>
Cell_001 Cell_001 positive S treat1 Z
Cell_002 Cell_002 positive G0 treat1 I
Cell_003 Cell_003 negative G1 treat1 W
Cell_004 Cell_004 negative S treat1 O
Cell_005 Cell_005 negative G1 treat2 I
... ... ... ... ... ...
Cell_036 Cell_036 negative G0 treat1 C
Cell_037 Cell_037 negative G0 treat1 E
Cell_038 Cell_038 negative G0 treat2 D
Cell_039 Cell_039 negative G1 treat1 F
Cell_040 Cell_040 negative G0 treat2 F
新增一行的方法：
#比如我想給每一行加上一個隨機(jī)數(shù)，隨機(jī)數(shù)生成的函數(shù)是runif()
> rowData(example_sce)$stuff <- runif(nrow(example_sce))
> rowData(example_sce)
DataFrame with 2000 rows and 1 column
stuff
<numeric>
Gene_0001 0.0491068968549371
Gene_0002 0.611523448256776
Gene_0003 0.282845608890057
Gene_0004 0.0250186347402632
Gene_0005 0.182038959115744
... ...
Gene_1996 0.633494263514876
Gene_1997 0.813403354492038
Gene_1998 0.987525492208079
Gene_1999 0.204807848902419
Gene_2000 0.369947757339105
還有一些其他的函數(shù)粱胜，以后如果有需要可以用的上柄驻，比如（摘自：單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-15-利用scRNAseq包學(xué)習(xí)scater)：
（1）如果你的數(shù)據(jù)里有spike-in control，可以用isSpike 對spike-in操作
（2）sizeFactors 是進(jìn)行標(biāo)準(zhǔn)化時對細(xì)胞文庫大小計算的結(jié)果
（3）reducedDim 對降維結(jié)果（reduced dimensionality results）操作
其他導(dǎo)入數(shù)據(jù)的方法
如果你的count矩陣是存在一個CSV文件里的焙压，可以用read.table() 來讀取鸿脓。
如果是大的數(shù)據(jù)集，可以用"Matrix"包里的readSparseCounts()函數(shù)讀取涯曲。
或者野哭，如果你的數(shù)據(jù)來自10X Genomics測序，你可以用 read10xCounts函數(shù)讀取幻件，它會自動產(chǎn)生一個sparse矩陣拨黔。
另外，kallisto和Salmon的轉(zhuǎn)錄本豐度數(shù)據(jù)可以使用tximeta包來進(jìn)行讀取傲武。
質(zhì)量控制
scater可以從三個水平上進(jìn)行質(zhì)量控制（QC）：
（1）過濾細(xì)胞
（2）過濾基因
（3）實(shí)驗(yàn)變量的質(zhì)控
細(xì)胞水平QC
細(xì)胞水平指標(biāo)的計算是利用perCellQCMetrics()函數(shù)進(jìn)行的蓉驹，它包括：
（1）sum: 細(xì)胞里的總count數(shù)（文庫大小）
（2）detected: 細(xì)胞里含有count的基因數(shù)
（3）subsets_X_percent: count的百分比（某個基因集）
NOTE這里需要安裝最新版的scater(1.14.4)揪利，如果是老版的scater，這一步是不能完成的
由于我的破電腦用了N種方法也沒安裝上最新版的scater狠持，所以質(zhì)控這些步驟我只把代碼放了上來疟位，自己并沒有親自運(yùn)行一遍，如果有同學(xué)可以安裝喘垂，最好要試一遍質(zhì)控的過程甜刻，因?yàn)橄裣乱徊絧erCellQCMetrics函數(shù)是最新版里新增的功能:
library(scater)
per.cell <- perCellQCMetrics(example_sce,
subsets=list(Mito=grep("mt-", rownames(example_sce))))
summary(per.cell$sum)
summary(per.cell$detected)
summary(per.cell$subsets_Mito_percent)
colData(example_sce) <- cbind(colData(example_sce), per.cell)
教程里并沒有用這個函數(shù)，而是用了另外一種：calculateQCMetrics
這個函數(shù)對細(xì)胞和基因都計算了很多種的指標(biāo)正勒，分別存在colData和rowData里得院，默認(rèn)情況下對count值進(jìn)行計算，你也可以通過參數(shù)exprs_values進(jìn)行修改：
> example_sce <- calculateQCMetrics(example_sce)
> colnames(colData(example_sce)) #列是樣品
[1] "Cell" "Mutation_Status" "Cell_Cycle"
[4] "Treatment" "whee" "is_cell_control"
[7] "total_features_by_counts" "log10_total_features_by_counts" "total_counts"
[10] "log10_total_counts" "pct_counts_in_top_50_features" "pct_counts_in_top_100_features"
[13] "pct_counts_in_top_200_features" "pct_counts_in_top_500_features"
> colnames(rowData(example_sce)) #行是基因
[1] "stuff" "is_feature_control" "mean_counts" #平均表達(dá)量 "log10_mean_counts"
[5] "n_cells_by_counts"#多少個細(xì)胞表達(dá)了該基因 "pct_dropout_by_counts"#基因在細(xì)胞中表達(dá)量為0的細(xì)胞數(shù)占比 "total_counts" "log10_total_counts"
如果你的測序里添加了ERCC章贞，那么這里還可以設(shè)置對照：
> example_sce <- calculateQCMetrics(example_sce,
feature_controls = list(ERCC = 1:20, mito = 500:1000),
cell_controls = list(empty = 1:5, damaged = 31:40))
>
> all_col_qc <- colnames(colData(example_sce))
> all_col_qc <- all_col_qc[grep("ERCC", all_col_qc)]
> all_col_qc #這里就是ERCC相關(guān)的信息了
[1] "total_features_by_counts_ERCC" "log10_total_features_by_counts_ERCC"
[3] "total_counts_ERCC" "log10_total_counts_ERCC"
[5] "pct_counts_ERCC" "pct_counts_in_top_50_features_ERCC"
[7] "pct_counts_in_top_100_features_ERCC" "pct_counts_in_top_200_features_ERCC"
[9] "pct_counts_in_top_500_features_ERCC"
可視化
（一）細(xì)胞水平上的質(zhì)控可以畫一個總基因表達(dá)量 比上feature control的百分比祥绞，如果一個基因在總基因表達(dá)量上的比例多，在ERCC里少鸭限，就是正常的細(xì)胞蜕径，反之則不正常。
> plotColData(example_sce, x = "total_features_by_counts",
y = "pct_counts_feature_control", colour = "Mutation_Status") +
theme(legend.position = "top") +
stat_smooth(method = "lm", se = FALSE, size = 2, fullrange = TRUE)
（二）plotScater（單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-15-利用scRNAseq包學(xué)習(xí)scater）
先在表達(dá)量最高的基因中選一部分（默認(rèn)是500個）败京，然后從高到低累加兜喻，看看它們對每個細(xì)胞文庫的貢獻(xiàn)值如何。它將不同細(xì)胞的不同基因表達(dá)分布繪制出來赡麦，就像芯片數(shù)據(jù)或bulk轉(zhuǎn)錄組數(shù)據(jù)每個樣本做一個箱線圖朴皆，來看基因表達(dá)量分布帕识。但是由于單細(xì)胞數(shù)據(jù)樣本太多，沒辦法全畫出箱線圖遂铡。這個圖在處理不同批次細(xì)胞時就可以很清楚地看到它們之間的差異肮疗。
為了看不同細(xì)胞的表達(dá)量差異，可以利用colData中的數(shù)據(jù)將細(xì)胞分類忧便；默認(rèn)使用counts值族吻，但是只要在assays存在的表達(dá)矩陣，就可以在exprs_values參數(shù)中自定義:
> plotScater(example_sce, block1 = "Mutation_Status", block2 = "Treatment",
colour_by = "Cell_Cycle", nfeatures = 300, exprs_values = "counts")
（三）QC結(jié)果兩兩比較（單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-15-利用scRNAseq包學(xué)習(xí)scater）
# 比較feature的屬性
> plotRowData(example_sce, x = "n_cells_by_counts", y = "mean_counts")
# 比較樣本的屬性
> p1 <- plotColData(example_sce, x = "total_counts",
y = "total_features_by_counts")
> p2 <- plotColData(example_sce, x = "pct_counts_feature_control",
y = "total_features_by_counts")
> p3 <- plotColData(example_sce, x = "pct_counts_feature_control",
y = "pct_counts_in_top_50_features")
# 最后組合
> multiplot(p1, p2, p3, cols = 3)
過濾低質(zhì)量細(xì)胞
舉個栗子珠增，我們只保留有10萬個counts和500個基因表達(dá)的細(xì)胞：
> keep.total <- example_sce$total_counts > 1e5
> keep.n <- example_sce$total_features_by_counts > 500
> filtered <- example_sce[,keep.total & keep.n]
> dim(filtered)
[1] 1973 37
isOutlier 函數(shù)提供了一個更為合適的方法來選擇閾值超歌。它定義距離中位數(shù)一定數(shù)量的中位數(shù)絕對偏差(MADs)的閾值。超過這個值的細(xì)胞認(rèn)為是離群的蒂教，從而被過濾掉：
> keep.total <- isOutlier(example_sce$total_counts, nmads=3,
type="lower", log=TRUE)
> filtered <- example_sce[,keep.total]
使用quickPerCellQC()函數(shù)還可以更方便地檢測幾個常見指標(biāo)的異常值巍举。該方法使用總count數(shù)、檢測到的基因數(shù)和某個基因集的count百分比(如線粒體基因凝垛、spik -in轉(zhuǎn)錄本)來確定要舍棄哪些細(xì)胞：
qc.stats <- quickPerCellQC(per.cell, percent_subsets="subsets_Mito_percent")
colSums(as.matrix(qc.stats))
filtered <- example_sce[,!qc.stats$discard]
isOutlier的方法需要根據(jù)你的實(shí)驗(yàn)需求進(jìn)行調(diào)整懊悯，如測序深度、spike-in的加入梦皮、細(xì)胞類型炭分。
你也可以對基因進(jìn)行過濾：
> keep_feature <- nexprs(example_sce, byrow=TRUE) >= 4
> example_sce <- example_sce[keep_feature,]
> dim(example_sce)
[1] 1753 40
基因水平的質(zhì)控
指標(biāo)的定義
使用的perFeatureQCMetrics()函數(shù)包括下面幾個指標(biāo)：
（1） mean: 所有的細(xì)胞之間基因的平均count
（2）detected: 對于每一個基因都是非零的細(xì)胞百分比
（3）subsets_Y_ratio: 某一個子集的細(xì)胞和總細(xì)胞的平均count的比值
# 這里舉例：前10個孔是空孔
per.feat <- perFeatureQCMetrics(example_sce, subsets=list(Empty=1:10))
summary(per.feat$mean)
summary(per.feat$detected)
summary(per.feat$subsets_Empty_ratio)
calculateAverage()函數(shù)提供了更精確的平均值計算，它根據(jù)取平均值之前的相對文庫大小來調(diào)整計數(shù)：
#這里我用的是過濾前的對象來進(jìn)行演示的
> ave <- calculateAverage(example_sce)
> summary(ave)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.009 5.529 67.730 197.135 190.286 12428.798
我們還可以直接計算出表達(dá)基因的細(xì)胞數(shù)量：
> summary(nexprs(example_sce, byrow=TRUE))
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.00 6.00 12.00 14.66 21.00 40.00
可視化：
在默認(rèn)值下剑肯，只顯示前50個基因捧毛。行是基因，橙色的線(bar)代表該基因在每一個細(xì)胞里的表達(dá)量让网。圓圈代表這個基因表達(dá)量的中位數(shù)呀忧。
> plotHighestExprs(example_sce, exprs_values = "counts")
下面這個圖是對于每一個基因都是非零的細(xì)胞百分比：
> plotExprsFreqVsMean(example_sce)
根據(jù)行來取子集
從SingleCellExperiment對象里移除一些基因，比如我們把在所有細(xì)胞中都不表達(dá)的基因刪除：
> keep_feature <- rowSums(counts(example_sce) > 0) > 0
> example_sce <- example_sce[keep_feature,]
> dim(example_sce)
[1] 1973 40
也可以用現(xiàn)有的基因注釋進(jìn)行過濾溃睹，比如線粒體基因而账，通常我們對這些基因不感興趣，對我們了解細(xì)胞異質(zhì)性也沒有多大幫助因篇。
變量水平的質(zhì)控
用getVarianceExplained()函數(shù)(經(jīng)過標(biāo)準(zhǔn)化之后泞辐，請參閱下面的內(nèi)容)。
example_sce <- logNormCounts(example_sce) # see below.
vars <- getVarianceExplained(example_sce,
variables=c("tissue", "total mRNA mol", "sex", "age"))
head(vars)
plotExplanatoryVariables(vars)
計算表達(dá)值
標(biāo)準(zhǔn)化文庫大小差異
第一種方法：
最常用的函數(shù)是logNormCounts()惜犀，結(jié)果儲存在"logcounts"里铛碑。它是用每一個count除以自身的size factor，加一個 pseudo-count 和log-transforming虽界。
example_sce <- logNormCounts(example_sce)
assayNames(example_sce)
summary(librarySizeFactors(example_sce))
第二種方法：
用calculateCPM()函數(shù)汽烦，結(jié)果儲存在"cpm"中。相關(guān)函數(shù)還包括：calculateTPM()和calculateFPKM()莉御。
cpm(example_sce) <- calculateCPM(example_sce)
assay(example_sce, "normed") <- normalizeCounts(example_sce,
size_factors=runif(ncol(example_sce)), pseudo_count=1.5)
上面兩種方法是官網(wǎng)的代碼撇吞，我用的是下面這個俗冻，是教程里用的代碼：
> cpm(example_sce) <- calculateCPM(example_sce)
> example_sce <- normalize(example_sce)
> assayNames(example_sce)
[1] "counts" "cpm" "logcounts"
> plotExplanatoryVariables(example_sce)
也可以選擇一部分進(jìn)行畫圖：
> plotExplanatoryVariables(example_sce,
variables = c("total_features_by_counts", "total_counts",
"Mutation_Status", "Treatment", "Cell_Cycle"))
可視化表達(dá)值
plotExpression()函數(shù)可以畫某一個基因子集的表達(dá)值牍颈，默認(rèn)使用"logcounts"里的數(shù)據(jù)迄薄，也可以通過添加exprs_values參數(shù)來修改。
> plotExpression(example_sce, rownames(example_sce)[1:6])
你也可以按照分組的情況畫基因表達(dá)情況煮岁，比如在這個包自帶的數(shù)據(jù)info里讥蔽，有Mutation_Status，Cell_cycle的分組情況画机，你也可以根據(jù)這些來畫圖：
> plotExpression(example_sce, rownames(example_sce)[1:6],
x = "Mutation_Status", exprs_values = "logcounts")
> plotExpression(example_sce, rownames(example_sce)[1:6],
x = "Cell_Cycle", exprs_values = "logcounts")
你還可以同時將兩個變量一起畫出來：
> plotExpression(example_sce, rownames(example_sce)[1:6],
colour_by = "Cell_Cycle", shape_by = "Mutation_Status",
size_by = "Gene_0002")
降維
（一）PCA
最常用的是PCA降維的方法冶伞，PCA結(jié)果將存放在SingleCellExperiment對象的reducedDims里：
> example_sce <- runPCA(example_sce)
> str(reducedDim(example_sce, "PCA"))
num [1:40, 1:2] -11.35 -1.41 2.01 -9.8 6.35 ...
- attr(*, "dimnames")=List of 2
..$ : chr [1:40] "Cell_001" "Cell_002" "Cell_003" "Cell_004" ...
..$ : chr [1:2] "PC1" "PC2"
- attr(*, "percentVar")= num [1:2] 0.0679 0.0422
默認(rèn)情況下，runPCA()函數(shù)用第一個PCs（聚類）里變化量最大的500個基因來進(jìn)行計算步氏。你也可以用subset_row參數(shù)來調(diào)整你想要進(jìn)行計算的基因數(shù)响禽，或者基因子集。用name參數(shù)可以修改在reducedDims slot里的名稱：
> example_sce <- runPCA(example_sce, name="PCA2",
subset_row=rownames(example_sce)[1:1000],
ncomponents=25)
（二）其他降維方法(tSNE, UMAP)
> set.seed(1000)
> example_sce <- runTSNE(example_sce, perplexity=10)
> head(reducedDim(example_sce, "TSNE"))
[,1] [,2]
Cell_001 3.7408876 -0.8718345
Cell_002 0.9961570 -1.0279169
Cell_003 0.3328361 1.1045595
Cell_004 4.9745275 0.2320463
Cell_005 -4.1067009 0.1909082
Cell_006 0.5936069 -2.1781319
這里官方網(wǎng)站推薦多用幾個不同的隨機(jī)種子和perplexity值（表示近鄰細(xì)胞的數(shù)量）來跑tSNE荚醒，以確保你的結(jié)果是比較穩(wěn)定的芋类。
一個更常見的模式是用你之前跑的PCA結(jié)果，作為tSNE的輸入界阁，進(jìn)行分析侯繁。這有助于提高運(yùn)行速度，并且降低噪音泡躯，只集中在幾個主要的變量上巫击。下面這個代碼是用的PCA前10個dimension的結(jié)果來進(jìn)行tSNE計算：
> set.seed(1000)
> example_sce <- runTSNE(example_sce, perplexity=50,
dimred="PCA", n_dimred=10)
> head(reducedDim(example_sce, "TSNE"))
你也可以用UMAP來計算：
> example_sce <- runUMAP(example_sce)
> head(reducedDim(example_sce, "UMAP"))
[,1] [,2]
Cell_001 2.1841452 -0.4125833
Cell_002 0.2940562 1.0461869
Cell_003 -0.5211948 1.2969395
Cell_004 2.1369329 0.1383032
Cell_005 -2.0110685 -0.4748272
Cell_006 1.0558627 0.7567186
你還可以用diffusion maps來降維：
> example_sce <- runDiffusionMap(example_sce)
Registered S3 method overwritten by 'xts':
method from
as.zoo.xts zoo
> plotDiffusionMap(example_sce, colour_by = "Gene_0001", size_by = "Gene_1000")
降維的可視化
#根據(jù)Mutation_Status來畫
> plotReducedDim(example_sce, use_dimred = "PCA", colour_by = "Mutation_Status")
還可以利用表達(dá)量添加顏色、大小分組（參考單細(xì)胞轉(zhuǎn)錄組學(xué)習(xí)筆記-15-利用scRNAseq包學(xué)習(xí)scater）：
> # 看特定基因在PCA過程中起到的作用
> plotReducedDim(example_sce, use_dimred = "PCA",
colour_by = "Gene_1000", size_by = "Gene_0500")
如果你想同時畫2個以上成分的可視化圖的時候精续，可以用下面的代碼：
> example_sce <- runPCA(example_sce, ncomponents=20)
> plotPCA(example_sce, ncomponents = 4, colour_by = "Mutation_Status")
