Aaron T. L. Lun1, Davis J. McCarthy2,3 and John C. Marioni1,2,4

1Cancer Research UK Cambridge Institute, Li Ka Shing Centre, Robinson Way, Cambridge CB2 0RE, United Kingdom
2EMBL European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SD, United Kingdom
3St Vincent's Institute of Medical Research, 41 Victoria Parade, Fitzroy, Victoria 3065, Australia
4Wellcome Trust Sanger Institute, Wellcome Genome Campus, Hinxton, Cambridge CB10 1SA, United Kingdom

2019-05-03

簡介

在這個工作流程中，我們檢查了來自小鼠大腦細(xì)胞類型研究的異構(gòu)數(shù)據(jù)集(Zeisel et al. 2015)。這包括大約3000個不同類型的細(xì)胞幻赚，如少突膠質(zhì)細(xì)胞、小膠質(zhì)細(xì)胞和神經(jīng)元尿孔。使用Fluidigm C1微流體系統(tǒng)(2014年)分離單個細(xì)胞，使用基于umi的方案對每個細(xì)胞進(jìn)行文庫制備沾凄。測序后焚廊，通過計算每個基因的UMIs數(shù)目來量化表達(dá)。所有內(nèi)源性基因淋袖、線粒體基因和spikein轉(zhuǎn)錄本的計數(shù)數(shù)據(jù)（ count data ）可從http://linnarssonlab.org/cortex獲得鸿市。

library(BiocFileCache)
bfc <- BiocFileCache("raw_data", ask = FALSE)
base.url <- file.path("https://storage.googleapis.com",
    "linnarsson-lab-www-blobs/blobs/cortex")
mRNA.path <- bfcrpath(bfc, file.path(base.url, 
    "expression_mRNA_17-Aug-2014.txt"))
mito.path <- bfcrpath(bfc, file.path(base.url, 
    "expression_mito_17-Aug-2014.txt"))
spike.path <- bfcrpath(bfc, file.path(base.url, 
    "expression_spikes_17-Aug-2014.txt"))

格式化數(shù)據(jù)

計數(shù)數(shù)據(jù)分布在多個文件中，因此需要進(jìn)行一些工作來將它們合并到單個矩陣中即碗。我們定義了一個簡單的實用函數(shù)，用于從每個文件加載數(shù)據(jù)陌凳。(我們強(qiáng)調(diào)此函數(shù)僅與當(dāng)前數(shù)據(jù)集相關(guān)剥懒，不應(yīng)用于其他數(shù)據(jù)集。如果所有的計數(shù)都在一個單獨的文件中合敦，并且與元數(shù)據(jù)分離初橘，則通常不需要這種工作。

readFormat <- function(infile) { 
    # First column is empty.
    metadata <- read.delim(infile, stringsAsFactors=FALSE, header=FALSE, nrow=10)[,-1] 
    rownames(metadata) <- metadata[,1]
    metadata <- metadata[,-1]
    metadata <- as.data.frame(t(metadata))

    # First column after row names is some useless filler.
    counts <- read.delim(infile, stringsAsFactors=FALSE, 
        header=FALSE, row.names=1, skip=11)[,-1] 
    counts <- as.matrix(counts)
    return(list(metadata=metadata, counts=counts))
}

利用這個函數(shù)充岛，我們讀取內(nèi)源性基因保檐、ERCC spik -in轉(zhuǎn)錄本和線粒體基因的計數(shù)。

endo.data <- readFormat(mRNA.path)
spike.data <- readFormat(spike.path)
mito.data <- readFormat(mito.path)

我們還需要重新排列線粒體數(shù)據(jù)的列崔梗，因為順序與其他文件不一致夜只。

m <- match(endo.data$metadata$cell_id, mito.data$metadata$cell_id)
mito.data$metadata <- mito.data$metadata[m,]
mito.data$counts <- mito.data$counts[,m]

在這個特定的數(shù)據(jù)集中，一些基因由多個行表示蒜魄，這些行對應(yīng)于可選的基因組位置扔亥。為了便于解釋，我們將與單個基因?qū)?yīng)的所有行計數(shù)相加谈为。

raw.names <- sub("_loc[0-9]+$", "", rownames(endo.data$counts))
new.counts <- rowsum(endo.data$counts, group=raw.names, reorder=FALSE)
endo.data$counts <- new.counts

然后這些計數(shù)被合并成一個單獨的矩陣來構(gòu)建一個單獨的SingleCellExperiment對象旅挤。為了方便起見，所有cell的元數(shù)據(jù)都存儲在同一個對象中伞鲫，供以后訪問粘茄。

library(SingleCellExperiment)
all.counts <- rbind(endo.data$counts, mito.data$counts, spike.data$counts)
sce <- SingleCellExperiment(list(counts=all.counts), colData=endo.data$metadata)
dim(sce)

我們添加了基于基因的注釋來標(biāo)識對應(yīng)于每一類特性的行。我們還將為每一行確定Ensembl標(biāo)識符秕脓。

# Specifying the nature of each row.
nrows <- c(nrow(endo.data$counts), nrow(mito.data$counts), nrow(spike.data$counts))
is.spike <- rep(c(FALSE, FALSE, TRUE), nrows)
is.mito <- rep(c(FALSE, TRUE, FALSE), nrows)
isSpike(sce, "Spike") <- is.spike

# Adding Ensembl IDs.
library(org.Mm.eg.db)
ensembl <- mapIds(org.Mm.eg.db, keys=rownames(sce), keytype="SYMBOL", column="ENSEMBL")
rowData(sce)$ENSEMBL <- ensembl

sce

cell QC

該研究的最初作者已經(jīng)在數(shù)據(jù)發(fā)表之前移除了低質(zhì)量的細(xì)胞柒瓣。盡管如此儒搭，我們使用scater (McCarthy et al. 2017)計算一些質(zhì)量控制指標(biāo)來檢查剩余的細(xì)胞是否令人滿意。

library(scater)
sce <- calculateQCMetrics(sce, feature_controls=list(Mt=is.mito)) 

我們檢查了QC指標(biāo)在所有單元中的分布(圖1)嘹朗。這里的庫大小至少比416B數(shù)據(jù)集中觀察到的要小一個數(shù)量級师妙。這與使用UMI count而不是read count是一致的，因為每個轉(zhuǎn)錄分子只能產(chǎn)生一個UMI count屹培，但在片段化后可以產(chǎn)生多個reads默穴。此外，與416B數(shù)據(jù)集相比褪秀，spik - In比例變化更大蓄诽。這可能反映了當(dāng)存在多種細(xì)胞類型時，每個細(xì)胞內(nèi)源性RNA總量的更大變異性媒吗。

par(mfrow=c(2,2), mar=c(5.1, 4.1, 0.1, 0.1))
hist(sce$total_counts/1e3, xlab="Library sizes (thousands)", main="", 
    breaks=20, col="grey80", ylab="Number of cells")
hist(sce$total_features_by_counts, xlab="Number of expressed genes", main="", 
    breaks=20, col="grey80", ylab="Number of cells")
hist(sce$pct_counts_Spike, xlab="ERCC proportion (%)",
    ylab="Number of cells", breaks=20, main="", col="grey80")
hist(sce$pct_counts_Mt, xlab="Mitochondrial proportion (%)", 
    ylab="Number of cells", breaks=20, main="", col="grey80")

圖1:質(zhì)控指標(biāo)的直方圖仑氛，包括文庫的大小、表達(dá)基因的數(shù)量和UMIs在大腦數(shù)據(jù)集中所有細(xì)胞的spike-in 轉(zhuǎn)錄或線粒體基因中的比例闸英。

我們刪除了庫大小和表示的特性數(shù)量的小異常值锯岖，并刪除了大異常值。再一次甫何，有了線粒體轉(zhuǎn)錄的存在出吹，就意味著我們不必使用線粒體的比例。

libsize.drop <- isOutlier(sce$total_counts, nmads=3, type="lower", log=TRUE)
feature.drop <- isOutlier(sce$total_features_by_counts, nmads=3, type="lower", log=TRUE)
spike.drop <- isOutlier(sce$pct_counts_Spike, nmads=3, type="higher")

然后通過組合所有指標(biāo)的過濾器來移除低質(zhì)量的細(xì)胞辙喂。大部分細(xì)胞被保留捶牢，這表明原來的質(zhì)量控制程序一般是適當(dāng)?shù)摹?/p>

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

我們可以進(jìn)一步改進(jìn)我們的細(xì)胞過濾程序，在離群點設(shè)置一個或多個已知因素的批次巍耗，例如秋麸，mouse/plate。如前所述炬太，這將避免MAD的膨脹灸蟆，提高功率來移除低質(zhì)量的電池。但是娄琉，為了簡單起見次乓，我們不會這樣做，因為已經(jīng)執(zhí)行了足夠的質(zhì)量控制孽水。

細(xì)胞周期分類

cyclone(Scialdone et al. 2015)對大腦數(shù)據(jù)集的應(yīng)用表明票腰，大多數(shù)細(xì)胞處于G1期(圖2)。這需要使用Ensembl標(biāo)識符來匹配預(yù)定義的分類器女气。

library(scran)
mm.pairs <- readRDS(system.file("exdata", "mouse_cycle_markers.rds", package="scran"))
assignments <- cyclone(sce, mm.pairs, gene.names=rowData(sce)$ENSEMBL)
table(assignments$phase)

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

圖2:在大腦數(shù)據(jù)集上應(yīng)用基于成對分類器的細(xì)胞周期階段得分杏慰，其中每個點代表一個細(xì)胞。

然而，由于訓(xùn)練數(shù)據(jù)集和測試數(shù)據(jù)集之間的差異缘滥，這個結(jié)果的完整性需要謹(jǐn)慎轰胁。該分類器在C1數(shù)據(jù)上進(jìn)行訓(xùn)練，并考慮到協(xié)議中的偏差朝扼。大腦數(shù)據(jù)集使用UMI計數(shù)赃阀，它有一組不同的偏差，例如擎颖，只有3 '端覆蓋榛斯，沒有長度偏差，沒有放大噪聲搂捧。此外驮俗，許多神經(jīng)元細(xì)胞類型預(yù)計位于G0靜止期，這與細(xì)胞周期的其他階段不同(Coller, Sang, and Roberts 2006)允跑。cyclone一般會將這些細(xì)胞分配到訓(xùn)練集中已知的最近的階段王凑，即G1。

基因水平的度量

圖3顯示了大腦數(shù)據(jù)集中細(xì)胞群中表達(dá)最高的基因聋丝。這主要是由spik -in轉(zhuǎn)錄本占據(jù)索烹，反映了使用的spik -in濃度跨越整個表達(dá)范圍。與預(yù)期的一樣弱睦，也有一些組成性表達(dá)基因术荤。

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

圖3:分配給大腦數(shù)據(jù)集中最豐富的50個特征的總計數(shù)的百分比。對于每個特性每篷，每個條表示分配給單個cell的特性的百分比，而圓圈表示所有單元格的平均值端圈。條形圖由每個cell中表示的特征的總數(shù)來著色焦读，而圓形圖則根據(jù)該特征是否被標(biāo)記為控制特征來著色。

通過計算所有細(xì)胞的平均計數(shù)來量化基因豐度(圖4)舱权。如前所述矗晃，UMI計數(shù)通常低于read計數(shù)。

ave.counts <- calcAverage(sce, use_size_factors=FALSE)
hist(log10(ave.counts), breaks=100, main="", col="grey",
    xlab=expression(Log[10]~"average count"))

我們將平均計數(shù)保存到singlecellexper對象中供以后使用宴倍。我們還刪除了平均計數(shù)為零的基因张症，因為這意味著它們在任何細(xì)胞中都不表達(dá)。

rowData(sce)$ave.count <- ave.counts
to.keep <- ave.counts > 0
sce <- sce[to.keep,]
summary(to.keep)
##    Mode   FALSE    TRUE 
## logical       2   19894

Normalization of cell-specific biases

我們使用computeSumFactors()函數(shù)和額外的預(yù)聚類步驟(Lun鸵贬、Bach和Marioni 2016)對內(nèi)源性基因進(jìn)行標(biāo)準(zhǔn)化俗他。每個集群中的單元單獨歸一化，并重新調(diào)整大小因子以使其在集群之間具有可比性阔逼。這就避免了假定大多數(shù)基因在整個種群中都是non-DE——在成對的簇之間只需要non-DE的多數(shù)兆衅。然后進(jìn)行標(biāo)準(zhǔn)化，以確保不同集群中的細(xì)胞大小因子具有可比性。

• 我們使用平均計數(shù)閾值0.1來定義在標(biāo)準(zhǔn)化期間使用的高豐度基因羡亩。這低于min.mean=1的默認(rèn)閾值摩疑，反映了UMI計數(shù)通常小于讀計數(shù)的事實。
• 我們通過快速降維來加速聚類畏铆，然后在pc上對單元進(jìn)行聚類雷袋。這就是BSPARAM=參數(shù)的目的，它指示quickCluster()使用PCA的近似算法辞居。這個近似依賴于隨機(jī)初始化楷怒，所以我們需要設(shè)置隨機(jī)種子的可重復(fù)性。

library(BiocSingular)
set.seed(1000)
clusters <- quickCluster(sce, use.ranks=FALSE, BSPARAM=IrlbaParam())
sce <- computeSumFactors(sce, cluster=clusters, min.mean=0.1)
summary(sizeFactors(sce))

與416B分析相比速侈，每個細(xì)胞的總計數(shù)和大小因子之間的趨勢出現(xiàn)了更多的散點(圖5)率寡。這與不同類型細(xì)胞之間DE的增加相一致，這影響了庫大小標(biāo)準(zhǔn)化的準(zhǔn)確性(Robinson和Oshlack 2010)倚搬。相比之下冶共，大小因子是根據(jù)中位數(shù)比值來估計的，并且對細(xì)胞間DE的存在更可靠每界。

plot(sizeFactors(sce), sce$total_counts/1e3, log="xy",
    ylab="Library size (thousands)", xlab="Size factor")

圖5:來自池的大小因子捅僵，與大腦數(shù)據(jù)集中所有細(xì)胞的庫大小作圖。坐標(biāo)軸以對數(shù)刻度表示眨层。

還計算特定于spik -in集的大小因子庙楚，如前所述。

sce <- computeSpikeFactors(sce, type="Spike", general.use=FALSE)

最后趴樱，使用適當(dāng)?shù)拇笮∫蜃佑嬎忝總€內(nèi)源性基因或spik -in轉(zhuǎn)錄本的歸一化對數(shù)表達(dá)值馒闷。

sce <- normalize(sce)

建模和消除技術(shù)噪音

我們通過使用trendVar()函數(shù)擬合spik -in轉(zhuǎn)錄本的均值-方差趨勢，對技術(shù)噪聲進(jìn)行建模叁征。

var.fit <- trendVar(sce, parametric=TRUE, loess.args=list(span=0.4))
var.out <- decomposeVar(sce, var.fit)

圖6顯示了趨勢與技術(shù)方差的精確擬合纳账。技術(shù)和總方差也比416B數(shù)據(jù)集中的小得多。這是由于UMIs的使用捺疼，減少了可變PCR擴(kuò)增引起的噪聲(Islam et al. 2014)疏虫。此外，內(nèi)源基因的變異量始終低于內(nèi)源基因的變異量啤呼。這反映了不同類型細(xì)胞間基因表達(dá)的異質(zhì)性卧秘。

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

圖6:大腦數(shù)據(jù)集中所有細(xì)胞的歸一化對數(shù)表達(dá)值相對于每個基因平均值的方差。藍(lán)線代表的是在技術(shù)差異中官扣，spik -in轉(zhuǎn)錄本的平均值相關(guān)趨勢(也用紅點表示)翅敌。
我們檢查具有最大生物成分的基因的表達(dá)值分布，以確保它們不會受到離群值的驅(qū)動(圖7)醇锚。

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

圖7:大腦數(shù)據(jù)集中前10位HVGs的標(biāo)準(zhǔn)化日志表達(dá)式值的小提琴圖哼御。對于每個基因坯临，每個點表示單個細(xì)胞的對數(shù)表達(dá)值。

最后恋昼，我們使用PCA對表達(dá)式值進(jìn)行降噪處理看靠，為每個去除技術(shù)噪聲的單元格生成一組坐標(biāo)。在()中設(shè)置BSPARAM=IrlbaParam()將使用來自irlba包的方法執(zhí)行一個近似奇異值分解(SVD)液肌。這比在大型數(shù)據(jù)集上的精確算法要快得多挟炬，而且不會損失很多精度。近似算法包含隨機(jī)初始化嗦哆，因此我們設(shè)置種子以保證可重復(fù)性谤祖。

set.seed(1000)
sce <- denoisePCA(sce, technical=var.fit$trend, BSPARAM=IrlbaParam())
ncol(reducedDim(sce, "PCA"))

• 理論上，我們應(yīng)該在每個細(xì)胞的原板上進(jìn)行阻擋老速。然而粥喜，每個培養(yǎng)皿中只有20-40個細(xì)胞，而且種群的異質(zhì)性也很強(qiáng)橘券。這意味著我們不能假設(shè)每個平板上采樣的細(xì)胞類型的分布是相同的额湘。因此，為了避免倒退出潛在的生物學(xué)旁舰，我們不會在這個分析中阻止任何因素锋华。

• denoisePCA()中PCs保留的上限由max指定。max.rank=參數(shù),這是默認(rèn)設(shè)置為100箭窜，以確保近似的SVD快速運行毯焕。雖然默認(rèn)設(shè)置對于降維通常是令人滿意的，但rank可能更適合于極端異質(zhì)的人群磺樱。

降維數(shù)據(jù)挖掘

我們在去噪后的PCs上進(jìn)行降維纳猫，檢查是否有子結(jié)構(gòu)。在圖8的t-SNE圖(Van der Maaten和Hinton 2008)中竹捉，細(xì)胞分裂成清晰的集群续担，對應(yīng)不同的亞群。這與大腦中多種細(xì)胞類型的存在是一致的活孩。我們增加了復(fù)雜性，以犧牲局部規(guī)模來支持整體結(jié)構(gòu)的可視化乖仇。

set.seed(1000)
sce <- runTSNE(sce, use_dimred="PCA", perplexity=50)
tsne1 <- plotTSNE(sce, colour_by="Neurod6") + fontsize
tsne2 <- plotTSNE(sce, colour_by="Mog") + fontsize
multiplot(tsne1, tsne2, cols=2)

圖8:由去噪的腦數(shù)據(jù)集PCs構(gòu)建的t-SNE圖憾儒。每個點代表一個細(xì)胞，并根據(jù)其Neurod6(左)或Mog(右)的表達(dá)進(jìn)行著色乃沙。

PCA圖在將細(xì)胞分成許多不同的簇方面效果較差(圖9)起趾，這是因為前兩個PCs是由特定亞群之間的強(qiáng)烈差異驅(qū)動的，這降低了其他一些亞群之間更細(xì)微差異的分辨率警儒。盡管如此训裆，仍然可以看到一些子結(jié)構(gòu)眶根。

pca1 <- plotReducedDim(sce, use_dimred="PCA", colour_by="Neurod6") + fontsize
pca2 <- plotReducedDim(sce, use_dimred="PCA", colour_by="Mog") + fontsize
multiplot(pca1, pca2, cols=2)

圖9:用去噪的大腦數(shù)據(jù)集PCs構(gòu)建的PCA圖。每個點代表一個細(xì)胞边琉，并根據(jù)其神經(jīng)d6(左)或Mog(右)的表達(dá)進(jìn)行著色属百。

對于這兩種方法，我們根據(jù)特定基因的表達(dá)為每個細(xì)胞上色变姨。這是在低維空間中可視化表達(dá)式變化的有用策略族扰。如果選擇的基因是特定細(xì)胞類型的已知標(biāo)記，它也可以用來描述每個簇定欧。例如渔呵，Mog可用于識別與少突膠質(zhì)細(xì)胞相對應(yīng)的簇。

聚類

降維坐標(biāo)用于將細(xì)胞聚集成假定的亞群體砍鸠。我們通過構(gòu)建共享最近鄰圖(shared-nearest-neighbour扩氢，Xu和Su 2015)來做到這一點，其中單元是共享最近鄰的單元之間形成的節(jié)點和邊爷辱。然后使用igraph包中的方法將集群定義為圖中高度連接的細(xì)胞群落录豺。這比形成成對距離矩陣對大量細(xì)胞的層次聚類更有效。

snn.gr <- buildSNNGraph(sce, use.dimred="PCA")
cluster.out <- igraph::cluster_walktrap(snn.gr)
my.clusters <- cluster.out$membership
table(my.clusters)
## my.clusters
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16 
## 103 200 234 191 194 355 161 658 384  85  90  16 209  38  45  26

我們在圖10中可視化了t-SNE圖中所有cell的集群分配托嚣。相鄰的cell通常被分配到相同的集群中巩检，這表明該集群過程得到了正確的應(yīng)用。

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

圖10:腦數(shù)據(jù)集去噪PCs的t-SNE圖示启。每個點表示一個cell兢哭，并根據(jù)其指定的集群標(biāo)識進(jìn)行著色。

另一種方法是使用基于圖形(graph-based)的可視化夫嗓，如強(qiáng)制導(dǎo)向的布局(force-directed layouts,圖11)迟螺。它們很有吸引力，因為它們直接表示集群期間使用的關(guān)系舍咖。然而矩父，對于大型圖，收斂速度往往較慢排霉，因此可能需要對niter=進(jìn)行一些修改窍株，以確保結(jié)果是穩(wěn)定的。

set.seed(2000)
reducedDim(sce, "force") <- igraph::layout_with_fr(snn.gr, niter=5000)
plotReducedDim(sce, colour_by="cluster", use_dimred="force")

非常異構(gòu)的數(shù)據(jù)集可能會在第一輪集群中產(chǎn)生一些大型集群攻柠。它可以是有用的重復(fù)方差建模球订，去噪和聚類只使用每個初始集群中的細(xì)胞。這可以通過根據(jù)my.clusters的特定級別設(shè)置sce來實現(xiàn), 并在子集上重新應(yīng)用相關(guān)函數(shù)瑰钮。這樣做可能集中在一組不同的基因上冒滩，這些基因定義了一個初始簇內(nèi)的異質(zhì)性，而不是那些定義初始簇之間差異的基因浪谴。這將允許以更高的分辨率探索每個集群中的精細(xì)結(jié)構(gòu)开睡。不過因苹，為了簡單起見，我們將只使用與此工作流中清除子種群相對應(yīng)的廣泛集群篇恒。

• igraph包中提供了許多不同的聚類方法扶檐。我們發(fā)現(xiàn)Walktrap算法通常是一個不錯的默認(rèn)選擇(Yang, heimer, Tessone 2016)，盡管我們鼓勵用戶嘗試不同的算法婚度。

• 在buildSNNGraph中減少鄰居k的數(shù)量會降低圖的連接性蘸秘。這通常會導(dǎo)致形成更小的集群(Xu和Su 2015)，如果需要更大的分辨率蝗茁，這可能是可取的醋虏。

• 注意，我們不運行l(wèi)ibrary(igraph)哮翘，而是使用igraph::從包中提取方法颈嚼。這是因為igraph包含一個規(guī)范化方法，它將覆蓋來自scater的對應(yīng)方法饭寺，從而導(dǎo)致一些不尋常的bug阻课。

模塊化評分為社區(qū)檢測方法提供了一個全局的集群性能度量。簡單地說艰匙，它將簇內(nèi)邊緣的數(shù)量與隨機(jī)邊緣的空模型下的期望數(shù)量進(jìn)行了比較限煞。較高的模塊性分?jǐn)?shù)(接近最大值1)表明被檢測到的簇的內(nèi)部邊緣比較豐富，簇之間的邊緣相對較少员凝。

igraph::modularity(cluster.out)

我們通過檢驗每對簇的邊的總權(quán)值來進(jìn)一步研究簇署驻。對于每一對，觀察到的總權(quán)重將與空模型下的預(yù)期值進(jìn)行比較健霹，類似于模塊化計算旺上。大多數(shù)集群包含的內(nèi)部鏈接比預(yù)期的多(圖12)，而集群之間的鏈接比預(yù)期的少糖埋。這表明我們成功地將細(xì)胞聚集成高度連接的社區(qū)宣吱。

mod.out <- clusterModularity(snn.gr, my.clusters, get.values=TRUE)
ratio <- mod.out$observed/mod.out$expected
lratio <- log10(ratio + 1)

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

圖12:同一集群或不同集群中節(jié)點之間的總重量與隨機(jī)鏈接的空模型下的總重量之比log10的熱圖。

為了總結(jié)集群之間的關(guān)系瞳别，我們使用觀察到的總權(quán)重與期望總權(quán)重的比值來構(gòu)建跨集群的圖征候。基于集群的圖可以使用力指向的布局來可視化祟敛，以識別高度互連的“集群的集群”倍奢。這類似于Wolf等人(2017)提出的“圖的抽象”策略。

cluster.gr <- igraph::graph_from_adjacency_matrix(ratio, 
    mode="undirected", weighted=TRUE, diag=FALSE)
plot(cluster.gr, edge.width=igraph::E(cluster.gr)$weight*10)  

圖13:基于不同集群中節(jié)點間觀察到的總權(quán)值與期望總權(quán)值之比垒棋，力導(dǎo)向的布局顯示了集群之間的關(guān)系。一對簇之間的邊緣厚度與對應(yīng)的比例成正比痪宰。

• 我們不使用剪影系數(shù)來評估聚類的大型數(shù)據(jù)集叼架。這是因為cluster::silhouette需要構(gòu)建一個距離矩陣畔裕，當(dāng)涉及到許多細(xì)胞時，這可能是不可行的乖订。

+從技術(shù)上講扮饶，模塊化分?jǐn)?shù)是通過從期望總權(quán)重中減去觀測值得到的。我們使用這個比率乍构，因為它保證是正的甜无，并且不會因為每個簇中細(xì)胞數(shù)量的不同而出現(xiàn)大小上的差異。

marker genes

我們使用帶有direction="up"的findmarker函數(shù)來識別每個集群中上調(diào)的標(biāo)記基因哥遮。如前所述岂丘，我們重點關(guān)注上調(diào)基因，因為這些基因可以快速在異質(zhì)人群中提供細(xì)胞類型的陽性識別眠饮。我們研究了集群4的表奥帘，其中報告了集群4和其他每個集群之間的日志倍變化。為每個集群提供相同的輸出仪召，以便識別區(qū)分集群的基因寨蹋。

markers <- findMarkers(sce, my.clusters, direction="up")
marker.set <- markers[["4"]]
head(marker.set[,1:8], 10) # only first 8 columns, for brevity

## DataFrame with 10 rows and 8 columns
##               Top               p.value                   FDR
##         <integer>             <numeric>             <numeric>
## Snap25          1  2.6436828195046e-260 5.24427360905115e-256
## Mllt11          1 2.64835618578792e-189 1.05070883314949e-185
## Gad1            1 6.21481316278381e-174 1.54104060887682e-170
## Atp1a3          1 2.42819725183299e-171 5.35201654273444e-168
## Celf4           1  2.1901645927764e-161 2.89641966846033e-158
## Rcan2           1  4.64612038923461e-99  1.35536897295951e-96
## Synpr           1  4.34116077754636e-75  5.66550041738068e-73
## Slc32a1         1  1.84158661555892e-65  1.63818626425301e-63
## Ndrg4           2 2.35144138055514e-242 2.33227713330369e-238
## Stmn3           2 6.36954591686139e-195  3.1588170588194e-191
##                     logFC.1           logFC.2          logFC.3
##                   <numeric>         <numeric>        <numeric>
## Snap25      1.2501919230692  0.67652736125569 3.68263627661648
## Mllt11     1.54667484840171  1.09590198758324 2.99539448469482
## Gad1       3.91523916973606  3.65238977769234 3.98295747045627
## Atp1a3  0.00896823664721724 0.962530545122529 3.26123555545746
## Celf4      0.44643142748766 0.740188939705178 2.71368354840561
## Rcan2      2.20233550942963  1.69523449324012 1.85629710950148
## Synpr      3.14964318830122  2.80066779182528 3.23587929298924
## Slc32a1    1.64312498454162  1.64133992421794 1.74089829213631
## Ndrg4      1.00144003064511  1.20871495924897 3.67086401750027
## Stmn3      1.77997011413237  1.01476517703874 4.20227407178467
##                    logFC.5            logFC.6
##                  <numeric>          <numeric>
## Snap25  -0.298729905898302  0.717172380132666
## Mllt11   0.446025722578893  0.495996491073405
## Gad1       4.1518488345932   4.13283704263783
## Atp1a3   0.549478199572573 -0.211309991972772
## Celf4   -0.100730536941031  0.120453982728159
## Rcan2     1.15269900728503   2.30454808040103
## Synpr     3.04938017855904   3.25850885734223
## Slc32a1   1.74737321385146    1.7380206484717
## Ndrg4      0.3445857243819  0.754609477952267
## Stmn3    0.546398537343761  0.524932745209836

圖14顯示，與部分或全部其他簇相比扔茅，簇4細(xì)胞中的大多數(shù)頂標(biāo)記都具有很強(qiáng)的DE已旧。我們可以使用這些標(biāo)記在獨立細(xì)胞群的驗證研究中識別來自cluster 4的細(xì)胞≌倌龋快速查看標(biāo)記表明运褪，cluster 4代表基于Gad1和Slc6a1表達(dá)的間神經(jīng)元(Zeng et al. 2012)，與cluster 11中緊密相關(guān)的細(xì)胞不同萤晴，后者具有較高的Synpr表達(dá)吐句。

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

圖14:以平均值為中心的熱圖，以及大腦數(shù)據(jù)集中cluster 4頂部標(biāo)記集的規(guī)范化日志表達(dá)式值店读。列顏色表示每個cell格被分配到的集群嗦枢，如圖例所示。

另一種可視化方法是直接繪制所有其他集群的log變化(圖15)屯断。這種方法更簡潔文虏，在涉及許多包含不同數(shù)量cell的集群的情況下非常有用。

logFCs <- as.matrix(marker.set[1:50,-(1:3)])
colnames(logFCs) <- sub("logFC.", "", colnames(logFCs))

library(pheatmap)
max.lfc <- max(abs(range(logFCs)))
pheatmap(logFCs, breaks=seq(-5, 5, length.out=101))

圖15:與大腦數(shù)據(jù)集中的其他集群相比殖演，第4集群頂部標(biāo)記集的表達(dá)變化的熱度圖氧秘。
我們保存候選標(biāo)記基因的列表以供進(jìn)一步檢查，使用壓縮來減少文件大小趴久。

gzout <- gzfile("brain_marker_1.tsv.gz", open="wb")
write.table(marker.set, file=gzout, sep="\t", quote=FALSE, col.names=NA)
close(gzout)

結(jié)論

完成基本分析后丸相，我們將帶有相關(guān)數(shù)據(jù)的singlecellexper搽對象保存到文件中。這一點在這里尤其重要彼棍，因為大腦數(shù)據(jù)集非常大灭忠。如果要執(zhí)行進(jìn)一步的分析膳算，就不方便重復(fù)上面描述的所有預(yù)處理步驟。

saveRDS(file="brain_data.rds", sce)

在這個工作流程中使用的所有軟件包都可以從綜合的R檔案網(wǎng)絡(luò)(https://cran.r-project.org)或Bioconductor項目(http://bioconductor.org)獲得弛作。下面顯示了使用的包的具體版本號涕蜂，以及R安裝的版本。

sessionInfo()

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.2 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] pheatmap_1.0.12                       
##  [2] scRNAseq_1.9.0                        
##  [3] org.Mm.eg.db_3.8.2                    
##  [4] readxl_1.3.1                          
##  [5] TxDb.Mmusculus.UCSC.mm10.ensGene_3.4.0
##  [6] GenomicFeatures_1.36.0                
##  [7] Matrix_1.2-17                         
##  [8] edgeR_3.26.0                          
##  [9] limma_3.40.0                          
## [10] BiocNeighbors_1.2.0                   
## [11] TENxBrainData_1.3.0                   
## [12] HDF5Array_1.12.0                      
## [13] rhdf5_2.28.0                          
## [14] Rtsne_0.15                            
## [15] BiocSingular_1.0.0                    
## [16] scran_1.12.0                          
## [17] scater_1.12.0                         
## [18] ggplot2_3.1.1                         
## [19] SingleCellExperiment_1.6.0            
## [20] SummarizedExperiment_1.14.0           
## [21] DelayedArray_0.10.0                   
## [22] BiocParallel_1.18.0                   
## [23] matrixStats_0.54.0                    
## [24] GenomicRanges_1.36.0                  
## [25] GenomeInfoDb_1.20.0                   
## [26] org.Hs.eg.db_3.8.2                    
## [27] AnnotationDbi_1.46.0                  
## [28] IRanges_2.18.0                        
## [29] S4Vectors_0.22.0                      
## [30] Biobase_2.44.0                        
## [31] BiocGenerics_0.30.0                   
## [32] BiocFileCache_1.8.0                   
## [33] dbplyr_1.4.0                          
## [34] knitr_1.22                            
## [35] BiocStyle_2.12.0                      
## 
## loaded via a namespace (and not attached):
##   [1] tidyselect_0.2.5              RSQLite_2.1.1                
##   [3] trimcluster_0.1-2.1           grid_3.6.0                   
##   [5] ranger_0.11.2                 munsell_0.5.0                
##   [7] destiny_2.14.0                codetools_0.2-16             
##   [9] statmod_1.4.30                sROC_0.1-2                   
##  [11] withr_2.1.2                   batchelor_1.0.0              
##  [13] colorspace_1.4-1              highr_0.8                    
##  [15] robustbase_0.93-4             vcd_1.4-4                    
##  [17] VIM_4.8.0                     TTR_0.23-4                   
##  [19] labeling_0.3                  cvTools_0.3.2                
##  [21] GenomeInfoDbData_1.2.1        bit64_0.9-7                  
##  [23] xfun_0.6                      ggthemes_4.1.1               
##  [25] diptest_0.75-7                R6_2.4.0                     
##  [27] ggbeeswarm_0.6.0              robCompositions_2.1.0        
##  [29] rsvd_1.0.0                    RcppEigen_0.3.3.5.0          
##  [31] locfit_1.5-9.1                flexmix_2.3-15               
##  [33] mvoutlier_2.0.9               reshape_0.8.8                
##  [35] bitops_1.0-6                  assertthat_0.2.1             
##  [37] promises_1.0.1                scales_1.0.0                 
##  [39] nnet_7.3-12                   beeswarm_0.2.3               
##  [41] gtable_0.3.0                  beachmat_2.0.0               
##  [43] processx_3.3.0                rlang_0.3.4                  
##  [45] scatterplot3d_0.3-41          splines_3.6.0                
##  [47] rtracklayer_1.44.0            lazyeval_0.2.2               
##  [49] BiocManager_1.30.4            yaml_2.2.0                   
##  [51] abind_1.4-5                   httpuv_1.5.1                 
##  [53] tools_3.6.0                   bookdown_0.9                 
##  [55] zCompositions_1.2.0           RColorBrewer_1.1-2           
##  [57] proxy_0.4-23                  dynamicTreeCut_1.63-1        
##  [59] Rcpp_1.0.1                    plyr_1.8.4                   
##  [61] progress_1.2.0                zlibbioc_1.30.0              
##  [63] purrr_0.3.2                   RCurl_1.95-4.12              
##  [65] ps_1.3.0                      prettyunits_1.0.2            
##  [67] viridis_0.5.1                 cowplot_0.9.4                
##  [69] zoo_1.8-5                     haven_2.1.0                  
##  [71] cluster_2.0.9                 magrittr_1.5                 
##  [73] data.table_1.12.2             openxlsx_4.1.0               
##  [75] lmtest_0.9-37                 truncnorm_1.0-8              
##  [77] mvtnorm_1.0-10                hms_0.4.2                    
##  [79] mime_0.6                      evaluate_0.13                
##  [81] xtable_1.8-4                  smoother_1.1                 
##  [83] XML_3.98-1.19                 rio_0.5.16                   
##  [85] mclust_5.4.3                  gridExtra_2.3                
##  [87] compiler_3.6.0                biomaRt_2.40.0               
##  [89] tibble_2.1.1                  crayon_1.3.4                 
##  [91] htmltools_0.3.6               pcaPP_1.9-73                 
##  [93] later_0.8.0                   tidyr_0.8.3                  
##  [95] rrcov_1.4-7                   DBI_1.0.0                    
##  [97] ExperimentHub_1.10.0          fpc_2.1-11.2                 
##  [99] MASS_7.3-51.4                 rappdirs_0.3.1               
## [101] boot_1.3-22                   car_3.0-2                    
## [103] sgeostat_1.0-27               igraph_1.2.4.1               
## [105] forcats_0.4.0                 pkgconfig_2.0.2              
## [107] GenomicAlignments_1.20.0      foreign_0.8-71               
## [109] laeken_0.5.0                  sp_1.3-1                     
## [111] vipor_0.4.5                   dqrng_0.2.0                  
## [113] XVector_0.24.0                NADA_1.6-1                   
## [115] stringr_1.4.0                 callr_3.2.0                  
## [117] digest_0.6.18                 pls_2.7-1                    
## [119] Biostrings_2.52.0             rmarkdown_1.12               
## [121] cellranger_1.1.0              DelayedMatrixStats_1.6.0     
## [123] kernlab_0.9-27                curl_3.3                     
## [125] modeltools_0.2-22             shiny_1.3.2                  
## [127] Rsamtools_2.0.0               Rhdf5lib_1.6.0               
## [129] carData_3.0-2                 viridisLite_0.3.0            
## [131] pillar_1.3.1                  GGally_1.4.0                 
## [133] lattice_0.20-38               survival_2.44-1.1            
## [135] httr_1.4.0                    DEoptimR_1.0-8               
## [137] interactiveDisplayBase_1.22.0 glue_1.3.1                   
## [139] xts_0.11-2                    zip_2.0.1                    
## [141] prabclus_2.2-7                simpleSingleCell_1.8.0       
## [143] bit_1.1-14                    class_7.3-15                 
## [145] stringi_1.4.3                 blob_1.1.1                   
## [147] AnnotationHub_2.16.0          memoise_1.1.0                
## [149] dplyr_0.8.0.1                 irlba_2.3.3                  
## [151] e1071_1.7-1

References

Coller, H. A., L. Sang, and J. M. Roberts. 2006. “A new description of cellular quiescence.” PLoS Biol. 4 (3):e83.

Islam, S., A. Zeisel, S. Joost, G. La Manno, P. Zajac, M. Kasper, P. Lonnerberg, and S. Linnarsson. 2014. “Quantitative single-cell RNA-seq with unique molecular identifiers.” Nat. Methods 11 (2):163–66.

Lun, A. T., K. Bach, and J. C. Marioni. 2016. “Pooling across cells to normalize single-cell RNA sequencing data with many zero counts.” Genome Biol. 17 (April):75.

McCarthy, D. J., K. R. Campbell, A. T. Lun, and Q. F. Wills. 2017. “Scater: pre-processing, quality control, normalization and visualization of single-cell RNA-seq data in R.” Bioinformatics 33 (8):1179–86.

Pollen, A. A., T. J. Nowakowski, J. Shuga, X. Wang, A. A. Leyrat, J. H. Lui, N. Li, et al. 2014. “Low-coverage single-cell mRNA sequencing reveals cellular heterogeneity and activated signaling pathways in developing cerebral cortex.” Nat. Biotechnol. 32 (10):1053–8.

Robinson, M. D., and A. Oshlack. 2010. “A scaling normalization method for differential expression analysis of RNA-seq data.” Genome Biol. 11 (3):R25.

Scialdone, A., K. N. Natarajan, L. R. Saraiva, V. Proserpio, S. A. Teichmann, O. Stegle, J. C. Marioni, and F. Buettner. 2015. “Computational assignment of cell-cycle stage from single-cell transcriptome data.” Methods 85 (September):54–61.

Van der Maaten, L., and G. Hinton. 2008. “Visualizing Data Using T-SNE.” J. Mach. Learn. Res. 9 (2579-2605):85.

Wolf, F. Alexander, Fiona Hamey, Mireya Plass, Jordi Solana, Joakim S. Dahlin, Berthold Gottgens, Nikolaus Rajewsky, Lukas Simon, and Fabian J. Theis. 2017. “Graph Abstraction Reconciles Clustering with Trajectory Inference Through a Topology Preserving Map of Single Cells.” bioRxiv. Cold Spring Harbor Laboratory. https://doi.org/10.1101/208819.

Xu, C., and Z. Su. 2015. “Identification of cell types from single-cell transcriptomes using a novel clustering method.” Bioinformatics 31 (12):1974–80.

Yang, Z., R. Algesheimer, and C. J. Tessone. 2016. “A comparative analysis of community detection algorithms on artificial networks.” Sci. Rep. 6 (August):30750.

Zeisel, A., A. B. Munoz-Manchado, S. Codeluppi, P. Lonnerberg, G. La Manno, A. Jureus, S. Marques, et al. 2015. “Brain structure. Cell types in the mouse cortex and hippocampus revealed by single-cell RNA-seq.” Science 347 (6226):1138–42.

Zeng, H., E. H. Shen, J. G. Hohmann, S. W. Oh, A. Bernard, J. J. Royall, K. J. Glattfelder, et al. 2012. “Large-scale cellular-resolution gene profiling in human neocortex reveals species-specific molecular signatures.” Cell 149 (2):483–96.

url