1 簡(jiǎn)介

1.1 scRNA-seq

• 2009年由湯富酬等人首次公開
• 測(cè)量每個(gè)基因在細(xì)胞群中表達(dá)水平的分布
• 允許研究在轉(zhuǎn)錄組中特定于細(xì)胞的變化的非常重要的生物學(xué)問題膏萧，例如細(xì)胞類型識(shí)別绣张，細(xì)胞反應(yīng)的異質(zhì)性焊刹，基因表達(dá)的隨機(jī)性，細(xì)胞內(nèi)基因調(diào)控網(wǎng)絡(luò)的推斷
• bulk和單細(xì)胞RNA-seq之間的主要區(qū)別在于劫映，每個(gè)測(cè)序文庫代表單個(gè)細(xì)胞蕴侧，而不是細(xì)胞群體。因此判帮，必須特別注意比較不同細(xì)胞（測(cè)序文庫）的結(jié)果。
• 目前單細(xì)胞測(cè)序主要由于轉(zhuǎn)錄物的起始量較低溉箕，在擴(kuò)增（1 million fold）的時(shí)候會(huì)產(chǎn)生許多問題晦墙。目前，提高轉(zhuǎn)錄物捕獲效率和降低擴(kuò)增偏差是研究的活躍領(lǐng)域肴茄∩纬可以通過適當(dāng)?shù)臍w一化緩解其中的一些問題。

1.2 scRNA-seq的技術(shù)

• CEL-seq (Hashimshony et al. 2012)
• CEL-seq2 (Hashimshony et al. 2016)
• Drop-seq (Macosko et al. 2015)
• InDrop-seq (Klein et al. 2015)
• MARS-seq (Jaitin et al. 2014)
• SCRB-seq (Soumillon et al. 2014)
• Seq-well (Gierahn et al. 2017)
• Smart-seq (Picelli et al. 2014)
• Smart-seq2 (Picelli et al. 2014)
• SMARTer
• STRT-seq (Islam et al. 2013)

2 分析流程

2.1 讀入數(shù)據(jù)

本次使用的數(shù)據(jù)集（Lun et al.2017 ）該數(shù)據(jù)集包含兩塊416B細(xì)胞板（永生的小鼠骨髓祖細(xì)胞系）寡痰，使用Smart-seq2進(jìn)行處理（Picelli et al抗楔，2014）。在制備文庫之前拦坠，還將spike-in RNA添加到每個(gè)細(xì)胞的裂解物中连躏。進(jìn)行高通量測(cè)序，并通過計(jì)數(shù)映射到其外顯子區(qū)域的讀數(shù)總數(shù)來量化每個(gè)基因的表達(dá)贞滨。類似地入热，通過計(jì)數(shù)映射到刺入?yún)⒖夹蛄械淖x數(shù)的數(shù)目來測(cè)量每個(gè)刺入轉(zhuǎn)錄本的量。

首先使用BiocFileCache包下載兩個(gè)416B細(xì)胞板的數(shù)據(jù)晓铆，并將其保存在指定路徑中勺良，隨后使用read.delim讀取并使用SingleCellExperiment轉(zhuǎn)化為單細(xì)胞分析中的對(duì)象，它可以儲(chǔ)存各種信息骄噪，包括表達(dá)矩陣尚困，基因注釋信息，樣本分組信息链蕊，文庫信息和降維信息等絕大部分信息尾组，極大地簡(jiǎn)化了我們分析時(shí)的操作忙芒。

library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
lun.zip <- bfcrpath(bfc, 
                    file.path("https://www.ebi.ac.uk/arrayexpress/files",
                              "E-MTAB-5522/E-MTAB-5522.processed.1.zip"))
lun.sdrf <- bfcrpath(bfc, 
                     file.path("https://www.ebi.ac.uk/arrayexpress/files",
                               "E-MTAB-5522/E-MTAB-5522.sdrf.txt"))
unzip(lun.zip, exdir='E:/simpleSingcell/simpleSingcell')

plate1 <- read.delim(file.path('E:/simpleSingcell/simpleSingcell', "counts_Calero_20160113.tsv"), 
                     header=TRUE, row.names=1, check.names=FALSE)
plate2 <- read.delim(file.path('E:/simpleSingcell/simpleSingcell', "counts_Calero_20160325.tsv"), 
                     header=TRUE, row.names=1, check.names=FALSE)

gene.lengths <- plate1$Length # First column is the gene length.
plate1 <- as.matrix(plate1[,-1]) # Discarding gene length (as it is not a cell).
plate2 <- as.matrix(plate2[,-1])
rbind(Plate1=dim(plate1), Plate2=dim(plate2))

stopifnot(identical(rownames(plate1), rownames(plate2)))
all.counts <- cbind(plate1, plate2)

library(SingleCellExperiment)
sce <- SingleCellExperiment(list(counts=all.counts))
rowData(sce)$GeneLength <- gene.lengths
sce
##查看sce這個(gè)SingleCellExperiment，它是一個(gè)S4對(duì)象讳侨，我們可以使用rowData，colData奏属，assay等命令獲取對(duì)象內(nèi)的主要信息跨跨，seurat中的CreateSeuratObject也是一個(gè)S4對(duì)象，不過命令稍有不同

#class: SingleCellExperiment 
#dim: 46703 192 
#metadata(0):
#assays(1): counts
#rownames(46703): ENSMUSG00000102693 ENSMUSG00000064842 ... SIRV7 CBFB-MYH11-mcherry
#rowData names(1): GeneLength
#colnames(192): SLX-9555.N701_S502.C89V9ANXX.s_1.r_1 SLX-9555.N701_S503.C89V9ANXX.s_1.r_1 ...
  SLX-11312.N712_S508.H5H5YBBXX.s_8.r_1 SLX-11312.N712_S517.H5H5YBBXX.s_8.r_1
#colData names(0):
#reducedDimNames(0):
#spikeNames(0):

2.2基因注釋

由于本例中存在ERCC基因囱皿，我們需要將這些基因單獨(dú)保存在對(duì)象中勇婴，以便于后續(xù)步驟提取嘱腥；不僅如此耕渴，我們還需要獲取基因的ensembl id，symbol id齿兔，染色體定位等信息橱脸，儲(chǔ)存在rowData中，將樣本的信息保存在colData中分苇。

#首先獲取表達(dá)矩陣中的ERCC信息
isSpike(sce, "ERCC") <- grepl("^ERCC", rownames(sce))
summary(isSpike(sce, "ERCC"))

#Mode   FALSE    TRUE 
#logical   46611      92

#本例中不僅包含ERCC添诉，作者還添加了SIRV信息，我們將其刪掉
is.sirv <- grepl("^SIRV", rownames(sce))
sce <- sce[!is.sirv,] 
summary(is.sirv)

#Mode   FALSE    TRUE 
#logical   46696       7 

#讀取樣本信息医寿，并將其添加到colData中栏赴，主要包含批次以及實(shí)驗(yàn)處理等信息
metadata <- read.delim(lun.sdrf, check.names=FALSE, header=TRUE)
m <- match(colnames(sce), metadata[["Source Name"]]) # Enforcing identical order.
stopifnot(all(!is.na(m))) # Checking that nothing's missing.
metadata <- metadata[m,]
head(colnames(metadata))

#[1] "Source Name"                "Comment[ENA_SAMPLE]"        "Comment[BioSD_SAMPLE]"     
#[4] "Characteristics[organism]"  "Characteristics[cell line]" "Characteristics[cell type]"

colData(sce)$Plate <- factor(metadata[["Factor Value[block]"]])
pheno <- metadata[["Factor Value[phenotype]"]]
levels(pheno) <- c("induced", "control")
colData(sce)$Oncogene <- pheno
table(colData(sce)$Oncogene, colData(sce)$Plate)

#         20160113 20160325
# induced       48       48
#   control       48       48

#我們可以使用getBMFeatureAnnos使用調(diào)取biomaRt包來對(duì)基因進(jìn)行注釋
library(scater)
sce <- getBMFeatureAnnos(
  sce,
  filters = "ensembl_gene_id", 
  attributes = c(
    "ensembl_gene_id",
    "external_gene_name",
    "chromosome_name",
    "start_position",
    "end_position"
  ), 
  biomart = "ENSEMBL_MART_ENSEMBL", 
  dataset = "mmusculus_gene_ensembl",
  host = "www.ensembl.org"
)

rownames(sce) <- uniquifyFeatureNames(rowData(sce)$ensembl_gene_id, rowData(sce)$external_gene_name)
rowData(sce)$ensembl_gene_id <- rownames(sce)
head(rownames(sce))

#[1] "4933401J01Rik" "Gm26206"       "Xkr4"          "Gm18956"       "Gm37180"       "Gm37363" 

mito <- which(rowData(sce)$chromosome_name=="MT")

##  [1] "Plate"                          "Oncogene"                      
##  [3] "is_cell_control"                "total_features_by_counts"      
##  [5] "log10_total_features_by_counts" "total_counts"                  
##  [7] "log10_total_counts"             "pct_counts_in_top_50_features" 
##  [9] "pct_counts_in_top_100_features" "pct_counts_in_top_200_features"

這些指標(biāo)的分布如圖所示，按致癌基因誘導(dǎo)狀態(tài)和板進(jìn)行了分層靖秩。目的是去除具有低文庫大小须眷，少量表達(dá)特征和高ERCC（或線粒體）比例的低質(zhì)量細(xì)胞。這樣的細(xì)胞可以干擾下游分析沟突，例如通過形成使結(jié)果解釋復(fù)雜化的不同簇花颗。

2.3 檢測(cè)離群值

為這些指標(biāo)選擇一個(gè)閾值并不容易，因?yàn)樗鼈兊闹等Q于實(shí)驗(yàn)方案事扭。不管細(xì)胞的質(zhì)量如何捎稚，測(cè)序到更大的深度將導(dǎo)致更多的讀取和更多表達(dá)的特征。同樣求橄，在方案中使用更多的刺入RNA將導(dǎo)致更高的ERCC比例今野。為了獲得自適應(yīng)閾值，我們假設(shè)大多數(shù)數(shù)據(jù)集均由高質(zhì)量單元格組成罐农，并確定對(duì)于各種QC指標(biāo)而言異常的單元格条霜。

離群值是根據(jù)與所有單元格中每個(gè)指標(biāo)的中值的中值絕對(duì)偏差（MAD）定義的。刪除的庫大小比中位數(shù)庫大小大3個(gè)MAD的單元格涵亏。對(duì)數(shù)轉(zhuǎn)換可提高小值時(shí)的分辨率宰睡，尤其是當(dāng)原始值的MAD等于或大于中值時(shí)蒲凶。還刪除了表達(dá)基因的對(duì)數(shù)轉(zhuǎn)化數(shù)量低于中位數(shù)3 MAD的細(xì)胞。

#首先使用calculateQCMetrics函數(shù)計(jì)算這些質(zhì)量控制指標(biāo)拆内，隨后做圖展示
sce <- calculateQCMetrics(sce, feature_controls=list(Mt=mito))
head(colnames(colData(sce)), 10)

sce$PlateOnco <- paste0(sce$Oncogene, ".", sce$Plate)
multiplot(
  plotColData(sce, y="total_counts", x="PlateOnco"),
  plotColData(sce, y="total_features_by_counts", x="PlateOnco"),
  plotColData(sce, y="pct_counts_ERCC", x="PlateOnco"),
  plotColData(sce, y="pct_counts_Mt", x="PlateOnco"),
  cols=2)

文庫大小旋圆，基因表達(dá)以及ERCC轉(zhuǎn)錄本或線粒體基因的比例

par(mfrow=c(1,3))
plot(sce$total_features_by_counts, sce$total_counts/1e6, xlab="Number of expressed genes",
     ylab="Library size (millions)")
plot(sce$total_features_by_counts, sce$pct_counts_ERCC, xlab="Number of expressed genes",
     ylab="ERCC proportion (%)")
plot(sce$total_features_by_counts, sce$pct_counts_Mt, xlab="Number of expressed genes",
     ylab="Mitochondrial proportion (%)")

每個(gè)QC指標(biāo)與所表達(dá)特征的總數(shù)的散點(diǎn)圖

#使用isOutlier計(jì)算大于3mad的離群值
libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", 
    log=TRUE, batch=sce$PlateOnco)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", 
    log=TRUE, batch=sce$PlateOnco)
spike.drop <- isOutlier(sce$pct_counts_ERCC, nmads=3, type="higher",
    batch=sce$PlateOnco)

該batch=參數(shù)可確保在指定批次的每個(gè)水平內(nèi)識(shí)別異常值，防止不感興趣的因素影響結(jié)果麸恍。

按上述條件將僅保留高質(zhì)量單元格灵巧。檢查每個(gè)過濾器除去的細(xì)胞數(shù)以及保留的細(xì)胞總數(shù)。刪除大量的單元格（> 10％）可能表明數(shù)據(jù)質(zhì)量出現(xiàn)整體問題抹沪。

keep <- !(libsize.drop | feature.drop | spike.drop)
data.frame(ByLibSize=sum(libsize.drop), ByFeature=sum(feature.drop),
    BySpike=sum(spike.drop), Remaining=sum(keep))

##   ByLibSize ByFeature BySpike Remaining
## 1         5         4       6       183

然后刻肄，將SingleCellExperiment保留高質(zhì)量細(xì)胞，還將原始對(duì)象保存到文件中以供以后使用融欧。

sce$PassQC <- keep
saveRDS(sce, file="416B_preQC.rds")
sce <- sce[,keep]
dim(sce)

## [1] 46696   183

還可以使用PCA進(jìn)行離群值的檢測(cè)敏弃，在運(yùn)行PCA時(shí)，使用detect_outliers函數(shù)檢測(cè)離群值

sce_tmp <- runPCA(sce, use_coldata=TRUE, detect_outliers=TRUE)
table(sce_tmp  $outlier)

FALSE  TRUE 
  182     1 

還存在一些R包可以供我們質(zhì)控時(shí)使用噪馏，cellity使用支持向量機(jī)從scRNA-seq中識(shí)別低質(zhì)量細(xì)胞麦到。

2.4 細(xì)胞周期識(shí)別

基于基因表達(dá)數(shù)據(jù)將細(xì)胞分類為細(xì)胞周期階段。使用訓(xùn)練數(shù)據(jù)集逝薪，為每對(duì)基因計(jì)算兩個(gè)基因之間表達(dá)差異的跡象隅要。選擇跨細(xì)胞周期階段具有符號(hào)變化的對(duì)作為標(biāo)記。然后董济，可以根據(jù)每個(gè)標(biāo)記對(duì)的觀察符號(hào)是否與一個(gè)或另一個(gè)相一致步清，將測(cè)試數(shù)據(jù)集中的單元格劃分為適當(dāng)?shù)南啵@種方法是在scran包的cyclone函數(shù)中實(shí)現(xiàn)的虏肾。

#使用scran中的cyclone確定細(xì)胞周期廓啊，目前scran自帶人和小鼠的數(shù)據(jù)可以直接調(diào)用
set.seed(100)
library(scran)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", 
                                package="scran"))
assignments <- cyclone(sce, mm.pairs, gene.names=rowData(sce)$ensembl_gene_id)

plot(assignments$score$G1, assignments$score$G2M, 
    xlab="G1 score", ylab="G2/M score", pch=16)

如果G1分?jǐn)?shù)高于0.5且大于G2 / M分?jǐn)?shù)，則將細(xì)胞分類為處于G1階段封豪。如果G2 / M得分高于0.5且高于G1得分谴轮，則處于G2 / M階段；如果兩個(gè)分?jǐn)?shù)均未超過0.5吹埠，則為S階段第步。在此，絕大多數(shù)細(xì)胞被分類為處于G1期缘琅。我們將這些分配保存到SingleCellExperiment對(duì)象中以供以后使用粘都。

sce$phases <- assignments$phases
table(sce$phases)

## 
##  G1 G2M   S 
##  98  62  23

2.5 基因表達(dá)質(zhì)控

檢查了表達(dá)量最高的基因的身份。通常應(yīng)以組成型表達(dá)的轉(zhuǎn)錄本為主，例如核糖體或線粒體蛋白的轉(zhuǎn)錄本。如果其他類別的特征與預(yù)期的生物學(xué)不一致堡掏，則可能會(huì)引起關(guān)注宋梧。例如堆生，包含許多ERCC轉(zhuǎn)錄本的集合表明在文庫制備過程中添加了太多spike-in RNA专缠，而缺少核糖體蛋白和/或它們的假基因則表明存在次佳的比對(duì)。

fontsize <- theme(axis.text=element_text(size=12), axis.title=element_text(size=16))
plotHighestExprs(sce, n=50) + fontsize

低豐度基因存在問題淑仆，因?yàn)榱慊蚪咏愕挠?jì)數(shù)沒有太多信息可用于可靠的統(tǒng)計(jì)推斷涝婉。在涉及假設(shè)檢驗(yàn)的應(yīng)用中，這些基因通常不能提供足夠的證據(jù)來拒絕無效假設(shè)糯景，但它們?nèi)詴?huì)增加多重檢驗(yàn)校正的嚴(yán)重性嘁圈。計(jì)數(shù)的離散性也可能會(huì)干擾統(tǒng)計(jì)程序，例如通過破壞連續(xù)逼近的準(zhǔn)確性蟀淮。因此，在應(yīng)用下游方法之前钞澳，通常會(huì)在許多RNA-seq分析流程中刪除低豐度基因怠惶。

#首先做一個(gè)直方圖觀察基因表達(dá)的分布情況
ave.counts <- calcAverage(sce, use_size_factors=FALSE)
hist(log10(ave.counts), breaks=100, main="", col="grey80", 
    xlab=expression(Log[10]~"average count"))

基因平均分布的直方圖

可以將最小閾值應(yīng)用于此值，以濾除表達(dá)水平較低的基因轧粟。

#表達(dá)每個(gè)基因的細(xì)胞數(shù)量策治，這與大多數(shù)基因的平均計(jì)數(shù)密切相關(guān)
num.cells <- nexprs(sce, byrow=TRUE)
smoothScatter(log10(ave.counts), num.cells, ylab="Number of cells", 
    xlab=expression(Log[10]~"average count"))

#在這里我們過濾掉了平均表達(dá)量為0的基因，代表在所有細(xì)胞中都沒有表達(dá)的基因兰吟。
to.keep <- num.cells > 0
sce <- sce[to.keep,]
summary(to.keep)

##    Mode   FALSE    TRUE 
## logical   22833   23863

2.6 標(biāo)準(zhǔn)化表達(dá)量

讀取count數(shù)據(jù)受細(xì)胞之間捕獲效率和測(cè)序深度的差異的影響通惫。在進(jìn)行下游定量分析之前，需要進(jìn)行標(biāo)準(zhǔn)化以消除這些細(xì)胞特異性偏差混蔼。這通常是通過假設(shè)大多數(shù)基因在細(xì)胞之間不差異表達(dá)（DE）來完成的履腋。假設(shè)兩個(gè)細(xì)胞之間跨非DE多數(shù)基因的計(jì)數(shù)大小的任何系統(tǒng)性差異都代表偏差，并通過縮放去除惭嚣。更具體地說遵湖，將計(jì)算“size factor”，該大小因子表示每個(gè)庫中應(yīng)縮放計(jì)數(shù)的程度延旧。

#首先快速聚類，如果細(xì)胞數(shù)量較少可以不使用聚類
set.seed(1000)
clusters <- quickCluster(sce, use.ranks=FALSE, BSPARAM=IrlbaParam())
table(clusters)
#我們使用computeSumFactors計(jì)算每個(gè)細(xì)胞文庫大小
sce <- computeSumFactors(sce,min.mean=0.1槽地，clusters=clusters)
summary(sizeFactors(sce))

#以散點(diǎn)圖的形式展示size factor與文庫之間的相關(guān)性
plot(sce$total_counts/1e6, sizeFactors(sce), log="xy",
     xlab="Library size (millions)", ylab="Size factor",
     col=c("red", "black")[sce$Oncogene], pch=16)
legend("bottomright", col=c("red", "black"), pch=16, cex=1.2,
       legend=levels(sce$Oncogene))

size factor與文庫之間的相關(guān)性

由于ERCC基因并不是細(xì)胞本身存在的迁沫，所以計(jì)算size factor時(shí)并不適用于ERCC基因，而scater專門為ERCC開發(fā)了computeSpikeFactors函數(shù)來對(duì)ERCC進(jìn)行標(biāo)準(zhǔn)化

#general.use表示是否應(yīng)用到全局
sce <- computeSpikeFactors(sce, type="ERCC", general.use=FALSE)

隨后使用normalize對(duì)表達(dá)矩陣進(jìn)行反卷積標(biāo)準(zhǔn)化捌蚊，并在同時(shí)取對(duì)數(shù)降維集畅，對(duì)數(shù)轉(zhuǎn)換是有用的，因?yàn)檫@意味著值的任何差異都直接表示單元格之間表達(dá)的對(duì)數(shù)2倍變化逢勾。這通常比覆蓋范圍的絕對(duì)差異更相關(guān)牡整，后者需要在總體豐度的背景下進(jìn)行解釋。對(duì)數(shù)轉(zhuǎn)換還提供了方差穩(wěn)定化的一些措施溺拱，因此具有大方差的高豐度基因不會(huì)主導(dǎo)下游分析逃贝。數(shù)據(jù)保存在logcounts中

sce <- normalize(sce)

assayNames(sce)
##[1] "counts"    "logcounts"

seurat方法

在seurat中首先使用NormalizeData進(jìn)行全局縮放歸一化方法“ LogNormalize”谣辞，該方法將每個(gè)單元格的特征表達(dá)式測(cè)量結(jié)果與總表達(dá)式進(jìn)行歸一化，再乘以比例因子（默認(rèn)為10,000）沐扳，然后對(duì)結(jié)果進(jìn)行對(duì)數(shù)轉(zhuǎn)換泥从。

suerat_object<- NormalizeData(suerat_object, normalization.method = "LogNormalize", scale.factor = 10000)

接下來，我們應(yīng)用線性變換沪摄，這是像PCA這樣的降維技術(shù)之前的標(biāo)準(zhǔn)預(yù)處理步驟躯嫉。對(duì)應(yīng)ScaleData函數(shù)

• 移動(dòng)每個(gè)基因的表達(dá)，從而使整個(gè)細(xì)胞的平均表達(dá)為0
• 縮放每個(gè)基因的表達(dá)杨拐，從而使細(xì)胞之間的差異為1
• 此步驟在下游分析中具有相同的權(quán)重祈餐，因此高表達(dá)的基因不會(huì)占主導(dǎo)地位
• 結(jié)果存儲(chǔ)在 pbmc[["RNA"]]@scale.data

all.genes <- rownames(suerat_object)
suerat_object<- ScaleData(suerat_object, features = all.genes)

2.7 方差建模篩選HVG

真正的生物異質(zhì)性或無趣的技術(shù)噪音可能會(huì)驅(qū)動(dòng)跨基因觀察到的表達(dá)值的差異。為了區(qū)分這兩種可能性哄陶，我們需要對(duì)每個(gè)基因的表達(dá)值差異的技術(shù)成分進(jìn)行建模帆阳。我們使用一組ERCC轉(zhuǎn)錄本來完成，這些轉(zhuǎn)錄本以相同的數(shù)量添加到每個(gè)單元格中屋吨。因此蜒谤，ERCC轉(zhuǎn)錄本不應(yīng)表現(xiàn)出生物學(xué)變異性，即其表達(dá)的任何變異都應(yīng)是技術(shù)來源至扰。

我們使用該trendVar函數(shù)將均值相關(guān)趨勢(shì)擬合為ERCC轉(zhuǎn)錄本的對(duì)數(shù)表達(dá)值的方差鳍徽。block參數(shù)將對(duì)設(shè)定批次，以確保各板之間的技術(shù)差異不會(huì)擴(kuò)大差異敢课。給定一個(gè)基因的平均豐度阶祭，然后將趨勢(shì)的擬合值用作該基因的技術(shù)成分的估計(jì)值。最后翎猛，通過decomposeVar函數(shù)從每個(gè)基因的總方差中減去技術(shù)成分胖翰，來計(jì)算方差的生物學(xué)成分。

var.fit <- trendVar(sce, parametric=TRUE, block=sce$Plate,
    loess.args=list(span=0.3))
var.out <- decomposeVar(sce, var.fit)
head(var.out)

## DataFrame with 6 rows and 6 columns
##                                   mean               total
##                              <numeric>           <numeric>
## ENSMUSG00000103377 0.00807160215928894   0.011921865486065
## ENSMUSG00000103147  0.0346526072192529  0.0722196162535234
## ENSMUSG00000103161 0.00519472222570747 0.00493857699521053
## ENSMUSG00000102331   0.018666093059853   0.032923591860573
## ENSMUSG00000102948   0.059057000132083  0.0881371257735823
## Rp1                 0.0970243712569606    0.45233813529556
##                                    bio               tech            p.value
##                              <numeric>          <numeric>          <numeric>
## ENSMUSG00000103377 -0.0238255786088717 0.0357474440949367                  1
## ENSMUSG00000103147 -0.0812680860584481  0.153487702311972  0.999999999992144
## ENSMUSG00000103161 -0.0180705438722202 0.0230091208674307                  1
## ENSMUSG00000102331 -0.0497487337065681  0.082672325567141  0.999999999998056
## ENSMUSG00000102948  -0.173441452696662  0.261578578470245                  1
## Rp1                 0.0226096722909625  0.429728463004597 0.0354980966384924
##                                  FDR
##                            <numeric>
## ENSMUSG00000103377                 1
## ENSMUSG00000103147                 1
## ENSMUSG00000103161                 1
## ENSMUSG00000102331                 1
## ENSMUSG00000102948                 1
## Rp1                0.153727758280855

隨后將擬合出來的結(jié)果進(jìn)行展示

plot(var.out$mean, var.out$total, pch=16, cex=0.6, xlab="Mean log-expression", 
    ylab="Variance of log-expression")
curve(var.fit$trend(x), col="dodgerblue", lwd=2, add=TRUE)
cur.spike <- isSpike(sce)
points(var.out$mean[cur.spike], var.out$total[cur.spike], col="red", pch=16)

藍(lán)色線代表擬合線切厘，紅色表示ERCC轉(zhuǎn)錄本

隨后檢查具有最大生物學(xué)成分的基因的表達(dá)值分布萨咳。這樣可確保方差估計(jì)不受一個(gè)或兩個(gè)異常值的影響。

chosen.genes <- order(var.out$bio, decreasing=TRUE)[1:10]
plotExpression(sce, features=rownames(var.out)[chosen.genes]) + fontsize

生物學(xué)差異最大的10個(gè)基因的小提琴圖

• 當(dāng)我們的測(cè)序樣本中沒有添加spike-in基因的時(shí)候疫稿，我們就需要使用全部基因來建立擬合趨勢(shì)培他。

new.trend <- makeTechTrend(x=sce)

fit <- trendVar(sce, use.spikes=FALSE, loess.args=list(span=0.05))
plot(fit$mean, fit$var, pch=16)
curve(fit$trend(x), col="dodgerblue", add=TRUE)
curve(new.trend(x), col="red", add=TRUE)

#使用全部基因擬合出來的趨勢(shì)尋找HVG
fit$trend <- new.trend # overwrite trend.
dec <- decomposeVar(fit=fit) # use per-gene variance estimates in 'fit'.
top.dec <- dec[order(dec$bio, decreasing=TRUE),] 
head(top.dec)

2.8 去除批次效應(yīng)

limma

如前所述，數(shù)據(jù)是在兩個(gè)板上收集的遗座。板之間的加工中不可控的微小差異會(huì)導(dǎo)致批次效應(yīng)舀凛，即不同板上細(xì)胞之間表達(dá)的系統(tǒng)差異。這種差異我們并不關(guān)心途蒋，可以通過removeBatchEffect實(shí)現(xiàn)猛遍。這消除了批次的作用，同時(shí)考慮了致癌基因誘導(dǎo)的作用。

library(limma)
assay(sce, "corrected") <- removeBatchEffect(logcounts(sce), 
    design=model.matrix(~sce$Oncogene), batch=sce$Plate)
assayNames(sce)

## [1] "counts"    "logcounts" "corrected"

當(dāng)我們的數(shù)據(jù)類型相對(duì)簡(jiǎn)單的時(shí)候可以使用removeBatchEffect懊烤，它是假定細(xì)胞群的組成在各批次之間是已知的或相同的梯醒。如果我們想要合并更加復(fù)雜的數(shù)據(jù)集的時(shí)候，由于scRNA-seq的特性腌紧，有一些方法專門用于處理單細(xì)胞測(cè)序的批次效應(yīng)茸习，例如MNN方法，seurat的CCA方法等壁肋。

MNN去除批次效應(yīng)

使用fastMNN函數(shù)應(yīng)用于消除批次效應(yīng)号胚。為了減少計(jì)算工作量和技術(shù)噪聲，將所有數(shù)據(jù)投影到低維空間中浸遗。然后在該低維空間中執(zhí)行MNN的識(shí)別和校正向量的計(jì)算猫胁。該函數(shù)返回一個(gè)SingleCellExperiment包含低維校正值的對(duì)象，用于下游分析（如聚類或可視化）跛锌。

#首先對(duì)兩個(gè)樣本的基因名取交集
universe <- intersect(rownames(sce1), rownames(sce2))
#計(jì)算兩個(gè)樣本平均生物學(xué)差異
mean.bio <- (sce1[universe,"bio"] + sce2[universe,"bio"])/2
chosen <- universe[mean.bio > 0]
length(chosen)

重新縮放每個(gè)批次杜漠，以調(diào)整批次之間測(cè)序深度的差異。multiBatchNorm在調(diào)整大小因數(shù)以實(shí)現(xiàn)SingleCellExperiment對(duì)象之間覆蓋率的系統(tǒng)差異之后察净，該函數(shù)將重新計(jì)算對(duì)數(shù)歸一化的表達(dá)式值。（請(qǐng)注意盼樟，先前計(jì)算的大小的因素僅除去細(xì)胞之間的偏差內(nèi)的單個(gè)批次）氢卡。這通過去除的批次之間的技術(shù)差異一方面提高了校正的質(zhì)量。

rescaled <- batchelor::multiBatchNorm(
    sce.gse85241[universe,], 
    sce.gse81076[universe,]
)
rescaled.gse85241 <- rescaled[[1]]
rescaled.gse81076 <- rescaled[[2]]

set.seed(100) 
unc.gse81076 <- logcounts(rescaled.gse81076)[chosen,]
unc.gse85241 <- logcounts(rescaled.gse85241)[chosen,]

mnn.out <- batchelor::fastMNN(
    GSE81076=unc.gse81076, GSE85241=unc.gse85241,
    k=20, d=50, BSPARAM=IrlbaParam(deferred=TRUE)
)
mnn.out

2.8 降維

我們使用denoisePCA進(jìn)行PCA線性降維晨缴，使用technical參數(shù)設(shè)定技術(shù)噪音译秦，denoisePCA可以同時(shí)進(jìn)行PC的篩選，刪除不重要的PC以便于后續(xù)的分析击碗。

sce <- denoisePCA(sce, technical=var.out, assay.type="corrected")
dim(reducedDim(sce, "PCA")) 

## [1] 183  24

在降維之后筑悴，我們可以在低維空間對(duì)結(jié)果進(jìn)行可視化。

#查看前三個(gè)PC的情況
plotReducedDim(sce, use_dimred="PCA", ncomponents=3, 
    colour_by="Oncogene") + fontsize

按照處理?xiàng)l件劃分細(xì)胞

相比之下稍途，我們使用批次信息對(duì)細(xì)胞進(jìn)行分組阁吝，可以發(fā)現(xiàn)并沒有明顯的細(xì)胞分離，這證明我們的批次更正步驟已成功械拍。

plotReducedDim(sce, use_dimred="PCA", ncomponents=3, 
    colour_by="Plate") + fontsize

使用批次信息進(jìn)行分組

請(qǐng)注意突勇，plotReducedDim將使用已經(jīng)存儲(chǔ)在主成分分析sce的denoisePCA。這使我們能夠快速生成具有不同美感的新圖坷虑，而無需重復(fù)整個(gè)PCA計(jì)算甲馋。同樣，plotPCA將使用現(xiàn)有結(jié)果（如果可用）迄损，否則將重新計(jì)算定躏。用戶應(yīng)設(shè)置rerun=TRUE為在存在現(xiàn)有結(jié)果的情況下強(qiáng)制重新計(jì)算PC。

• 在PCA的基礎(chǔ)上，我們繼續(xù)使用t-SNE進(jìn)行非線性降維痊远， t- SNE可以比PCA更好地分離更多種群體中的細(xì)胞垮抗。這是因?yàn)榍罢呖梢灾苯硬东@高維空間中的非線性關(guān)系，而后者必須在線性軸上表示它們拗引。但是借宵，這種改進(jìn)是以更多的計(jì)算工作為代價(jià)的。

set.seed(100)
out5 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=5),
    colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 5")

set.seed(100)
out10 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=10),
    colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 10")

set.seed(100)
out20 <- plotTSNE(sce, run_args=list(use_dimred="PCA", perplexity=20),
    colour_by="Oncogene") + fontsize + ggtitle("Perplexity = 20")

multiplot(out5, out10, out20, cols=3)

使用不同的perplexity參數(shù)的t-SNE結(jié)果

t-SNE是一種隨機(jī)方法矾削，因此用戶應(yīng)多次運(yùn)行該算法以確保結(jié)果具有代表性壤玫。腳本應(yīng)通過set.seed設(shè)置一個(gè)種子，以確保所選結(jié)果可再現(xiàn)哼凯。還建議測(cè)試perplexity參數(shù)的不同設(shè)置欲间，因?yàn)檫@會(huì)影響低維空間中點(diǎn)的分布。

在這里断部，我們runTSNE以20的perplexity進(jìn)行調(diào)用猎贴，以將t-SNE結(jié)果存儲(chǔ)在SingleCellExperiment對(duì)象中。這避免了每當(dāng)我們想使用創(chuàng)建新圖時(shí)都重復(fù)計(jì)算plotTSNE蝴光，因?yàn)閷⑹褂么鎯?chǔ)的結(jié)果她渴。同樣，用戶可以設(shè)置rerun=TRUE為在存在存儲(chǔ)結(jié)果的情況下強(qiáng)制重新計(jì)算蔑祟。

set.seed(100)
sce <- runTSNE(sce, use_dimred="PCA", perplexity=20)
reducedDimNames(sce)

## [1] "PCA"  "TSNE"

2.9 細(xì)胞聚類

使用層次聚類劃分亞群

去噪的對(duì)數(shù)表達(dá)值用于將細(xì)胞聚類為推定的亞群趁耗。具體來說，我們使用Ward.D2對(duì)單元之間的歐式距離進(jìn)行層次聚類疆虚，以使每個(gè)聚類中的總方差最小苛败。這將產(chǎn)生一個(gè)樹狀圖，該樹狀圖將跨所選基因具有相似表達(dá)模式的細(xì)胞分組在一起径簿。

pcs <- reducedDim(sce, "PCA")
my.dist <- dist(pcs)
my.tree <- hclust(my.dist, method="ward.D2")

通過對(duì)樹狀圖應(yīng)用動(dòng)態(tài)樹切割來明確定義聚類罢屈。這利用了樹狀圖中的分支形狀來完善聚類定義，并且比cutree復(fù)雜樹狀圖更合適篇亭。通過在中手動(dòng)指定cutHeight缠捌，cutreeDynamic函數(shù)可以更好地控制亞群的劃分。我們還設(shè)置minClusterSize了一個(gè)低于默認(rèn)值20的值暗赶，以避免劃分虛假的小亞群鄙币。

library(dynamicTreeCut)
my.clusters <- unname(cutreeDynamic(my.tree, distM=as.matrix(my.dist), 
    minClusterSize=10, verbose=0))

隨后使用不同條件對(duì)聚類結(jié)果進(jìn)行檢查

#使用批次信息
table(my.clusters, sce$Plate)

#使用實(shí)驗(yàn)處理信息
table(my.clusters, sce$Oncogene)

在t-SNE圖中展示亞群的劃分，通常比較接近的細(xì)胞會(huì)被劃分在同一亞群中蹂随。

sce$cluster <- factor(my.clusters)
plotTSNE(sce, colour_by="cluster") + fontsize

亞群劃分之后的t-SNE圖

在聚類之后十嘿，我們使用輪廓寬度檢查聚類的分離性。理想情況下岳锁，每個(gè)亞群將包含許多具有較大正寬度的單元绩衷，這表明該群集與其他群集完全分開。

library(cluster)
clust.col <- scater:::.get_palette("tableau10medium") # hidden scater colours
sil <- silhouette(my.clusters, dist = my.dist)
sil.cols <- clust.col[ifelse(sil[,3] > 0, sil[,1], sil[,2])]
sil.cols <- sil.cols[order(-sil[,1], sil[,3])]
plot(sil, main = paste(length(unique(my.clusters)), "clusters"), 
    border=sil.cols, col=sil.cols, do.col.sort=FALSE) 

每個(gè)亞群中單元的輪廓寬度的條形圖

使用圖聚類劃分亞群

我們建立一個(gè)共享的最近鄰圖，并使用Walktrap算法來識(shí)別聚類咳燕。

snn.gr <- buildSNNGraph(sce, use.dimred="PCA")
clusters <- igraph::cluster_walktrap(snn.gr)
sce$Cluster <- factor(clusters$membership)
table(sce$Cluster)

## 1  2  3  4 
![Rplot15.png](https://upload-images.jianshu.io/upload_images/18112569-4555c77eff365d98.png?imageMogr2/auto-orient/strip%7CimageView2/2/w/1240)
##53 13 37 80 

我們查看觀察到的邊緣權(quán)重與預(yù)期邊緣權(quán)重之比勿决，以確認(rèn)群集是模塊化的。（我們沒有看模塊性評(píng)分本身招盲，因?yàn)槟K性評(píng)分在群集之間的變化幅度不同低缩，并且難以解釋。）下圖指出曹货，大多數(shù)群集都很好地分離開了咆繁，幾乎沒有強(qiáng)烈的相關(guān)亞群。

cluster.mod <- clusterModularity(snn.gr, sce$Cluster, get.values=TRUE)
log.ratio <- log2(cluster.mod$observed/cluster.mod$expected + 1)

library(pheatmap)
pheatmap(log.ratio, cluster_rows=FALSE, cluster_cols=FALSE, 
    color=colorRampPalette(c("white", "blue"))(100))

模塊之間的權(quán)重占比

隨后在t-SNE中檢查聚類情況

plotTSNE(sce, colour_by="Cluster")

基于圖聚類最后亞群劃分結(jié)果

使用一致性聚類對(duì)亞群進(jìn)行劃分

讓我們SC3對(duì)數(shù)據(jù)進(jìn)行聚類顶籽。SC3的優(yōu)點(diǎn)是它可以直接攝取SingleCellExperiment對(duì)象玩般，當(dāng)我們不清楚聚類條目K時(shí)，SC3可以估計(jì)多個(gè)群集

library(SC3) # BiocManager::install('SC3')
sce <- sc3_estimate_k(sce) # 先預(yù)估一下聚類亞群
sce@metadata$sc3$k_estimation # 預(yù)估13個(gè)亞群
rowData(sce)$feature_symbol <- rownames(rowData(sce))

# 接下來正式運(yùn)行礼饱，kn參數(shù)表示的預(yù)估聚類數(shù)
# 我們這里自定義為5組
kn <- 5 
start=Sys.time()
sce <- sc3(sce, ks = kn, biology = TRUE) 
end=Sys.time()
(dur=end-start)


# 會(huì)將聚類結(jié)果放入表型信息(sce@colData)中去坏为，默認(rèn)叫sc3_cluster，這里人為改個(gè)名稱 
sc3_cluster="sc3_5_clusters"
# 最后進(jìn)行可視化比較之前tSNE的Kmeans聚類和SC3的聚類的一致性
sc3_plot_consensus(sce, k = kn, show_pdata = c("cluster",'sc3_5_clusters'))

SC3聚類情況

2.10 Marker基因的篩選

我們使用findMarkers來尋找每個(gè)亞群中的Marker基因

top.markers <- rownames(marker.set)[marker.set$Top <= 10]
plotHeatmap(sce, features=top.markers, columns=order(sce$cluster), 
    colour_columns_by=c("cluster", "Plate", "Oncogene"),
    cluster_cols=FALSE, center=TRUE, symmetric=TRUE, zlim=c(-5, 5)) 

marker基因展示情況