代碼分析 | 單細胞轉錄組clustering詳解

我們在單細胞轉錄組分析中最為常用的聚類可視化即為tSNE和UMAP（Hemberg-lab單細胞轉錄組數據分析（十二）- Scater單細胞表達譜tSNE可視化），不過非線性可視化方法（例如t-SNE）通常會擾亂數據中的全局結構。diffusion maps能夠很好的表達局部和全局結構，但無法針對二維（2D）可視化進行優化，因為它們不是專門為可視化設計的。（其實即便如此，我現在也習慣性應用diffusion maps，可以說明更多的信息，雖然有時候確實是很難看。。。）

芬蘭CSC-IT科學中心主講的生物信息課程(https://www.csc.fi/web/training/-/scrnaseq)視頻，官網上還提供了練習素材以及詳細代碼，今天就來練習一下單細胞數據clustering的過程。

在本教程中，我們將研究不同的聚類scRNA-seq數據集方法以表徵不同的細胞亞群。

所需加載包

suppressMessages(require(tidyverse))
suppressMessages(require(Seurat))
suppressMessages(require(cowplot))
suppressMessages(require(scater))
suppressMessages(require(scran))
suppressMessages(require(igraph))

數據集

本教程使用的是來自發育中的小鼠胚胎的小細胞數據集[1]。該數據集已進行了預處理、創建了SingleCellExperiment對象並對細胞進行了注釋（cellassign：用於腫瘤微環境分析的單細胞注釋工具（9月Nature））。

數據下載連結：https://github.com/NBISweden/excelerate-scRNAseq/tree/master/session-clustering/session-clustering_files。

#加載表達矩陣
deng <- readRDS("session-clustering_files/deng-reads.rds")
deng

#查看細胞類型注釋
table(colData(deng)$cell_type2)

特徵選擇

第一步是決定在細胞聚類中使用哪些基因（Hemberg-lab單細胞轉錄組數據分析（十）- Scater基因評估和過濾）。單細胞RNA-seq可以分析細胞中的大量基因，大多數基因的表達不足以提供有意義的信號並且通常受到技術噪音的幹擾，其中可能會添加一些不必要的信號，從而使生物變異模糊，而基因過濾還可以加快下遊分析的計算時間。

首先看一下所有基因的表達平均值和方差。哪些基因似乎不太重要，哪些可能是技術噪音？

#計算整個細胞的基因平均值
gene_mean <- rowMeans(counts(deng))

#計算整個細胞的基因方差
gene_var <- rowVars(counts(deng))

#ggplot plot
library(tidyverse)
gene_stat_df <- tibble(gene_mean,gene_var)
ggplot(data=gene_stat_df ,aes(x=log(gene_mean), y=log(gene_var))) + geom_point(size=0.5) + theme_classic()

濾除低豐度基因

低豐度基因大多無信息，不能代表數據的生物學差異。它們通常是由技術噪音（例如dropout）導致，在下遊分析中的存在通常會導致準確性降低，因為它們會干擾將要使用的某些統計模型，並且會毫無意義地增加計算時間，這在處理非常大的數據時可能至關重要。

abundant_genes <- gene_mean >= 0.5 #去除低豐度基因
# 繪製基因表達分布
hist(log10(gene_mean), breaks=100, main="", col="grey80",
xlab=expression(Log[10]~"average count"))
abline(v=log10(0.5), col="red", lwd=2, lty=2)

#在SingleCellExperiment Object中去除低豐度基因
deng <- deng[abundant_genes,]
dim(deng)

過濾在很少細胞中表達的基因

我們還可以過濾少量細胞中的某些基因。此過程將刪除一些在一兩個細胞中高度表達的異常基因。這些基因不需要進一步分析，因為它們主要來自人為的不規則擴增。值得注意的是，當分析的目的是檢測數據中非常稀有的亞群時，我們可能不希望進行此過程。

#計算每個非零表達基因的數量
numcells <- nexprs(deng, byrow=TRUE)

#過濾不到5個細胞中檢測到的基因
numcells2 <- numcells >= 5
deng <- deng[numcells2,]
dim(deng)

檢測高變的基因

highly variable gene（HVG）假設基因在細胞之間的表達差異很大，則其中的一些差異是由於細胞之間的生物學差異而不是技術噪音引起的。但是，基因的平均表達與不同細胞reads計數的差異之間存在正相關關係。如果僅保留高變異基因，將保留許多在每個細胞中表達但不能代表生物學變異的高表達管家基因。因此為了正確識別HVG，必須糾正此關係。

我們可以使用下列任一方法來確定高變基因。

選項1 ：將變異係數建模為平均值的函數。

# out <- technicalCV2(deng, spike.type=NA, assay.type= "counts")
# out$genes <- rownames(deng)
# out$HVG <- (out$FDR<0.05)
# as_tibble(out)
#
# # 繪製高變基因
# ggplot(data = out) + geom_point(aes(x=log2(mean), y=log2(cv2), color=HVG), size=0.5) + geom_point(aes(x=log2(mean), y=log2(trend)), color="red", size=0.1)
#
# ## 將HVG存入metadata
# metadata(deng)$hvg_genes <- rownames(deng)[out$HVG]

選項2 ：根據平均數對生物成分的方差建模。

首先，我們估算每個基因表達的差異，然後將差異分解為生物學和技術成分。然後將HVGs鑑定為具有最大生物學成分的基因。這避免了優先排序對由於技術因素（例如在RNA捕獲和文庫製備過程中的採樣噪聲）而高度可變的基因。詳細內容查看scran vignette：https://bioconductor.org/packages/devel/bioc/vignettes/scran/inst/doc/scran.html#5_variance_modelling

fit <- trendVar(deng, parametric=TRUE, use.spikes = FALSE)
dec <- decomposeVar(deng, fit)
dec$HVG <- (dec$FDR<0.00001)
hvg_genes <- rownames(dec[dec$FDR < 0.00001, ])

# plot highly variable genes
plot(dec$mean, dec$total, pch=16, cex=0.6, xlab="Mean log-expression",
ylab="Variance of log-expression")
o <- order(dec$mean)
lines(dec$mean[o], dec$tech[o], col="dodgerblue", lwd=2)
points(dec$mean[dec$HVG], dec$total[dec$HVG], col="red", pch=16)

## save the decomposed variance table and hvg_genes into metadata for safekeeping
metadata(deng)$hvg_genes <- hvg_genes
metadata(deng)$dec_var <- dec

降維

由於高頻率的噪聲（技術和生物學的）和大量的維度（即基因）導致聚類在計算上具有一定的困難，因此需要通過應用降維方法（例如PCA、tSNE和UMAP）解決這些問題。

#PCA (選擇組分數量)
deng <- runPCA(deng, method = "irlba",
ncomponents = 30,
feature_set = metadata(deng)$hvg_genes)

#繪製不同PC對變異解釋的差異
X <- attributes(deng@reducedDims$PCA)
plot(X$percentVar~c(1:30), type="b", lwd=1, ylab="Percentage variance" , xlab="PCs" , bty="l" , pch=1)

畫PCA plot(PC1 vs. PC2) ：

plotReducedDim(deng, "PCA", colour_by = "cell_type2")

繪製tSNE圖時需要注意tSNE是隨機方法。每次運行它時都會得到略有不同的結果。為了方便起見可以設置隨機種子，這樣就可以獲得相同的結果。

#tSNE
deng <- runTSNE(deng,
perplexity = 30,
feature_set = metadata(deng)$hvg_genes,
set.seed = 1)

plotReducedDim(deng, "TSNE", colour_by = "cell_type2")

聚類層次聚類

#計算距離（默認值：歐幾裡得距離）
distance_eucledian <- dist(t(logcounts(deng)))

#使用ward linkage執行分層聚類
ward_hclust_eucledian <- hclust(distance_eucledian,method = "ward.D2")
plot(ward_hclust_eucledian, main = "dist = eucledian, Ward linkage")

我反正是啥也看不見。。。但是看起來層次分類蠻明顯的。。。

現在，設置樹狀圖參數生成10個亞群，並在PCA圖上繪製聚類標籤。

#設置樹狀圖參數生成10個亞群
cluster_hclust <- cutree(ward_hclust_eucledian,k = 10)
colData(deng)$cluster_hclust <- factor(cluster_hclust)

plot_grid(plotReducedDim(deng, "PCA", colour_by = "cluster_hclust"),
plotReducedDim(deng, "PCA", colour_by = "cell_type2"))

在tSNE圖上繪製聚類標籤。

plot_grid(plotReducedDim(deng, "TSNE", colour_by = "cluster_hclust"),
plotReducedDim(deng, "TSNE", colour_by = "cell_type2"))

我們嘗試另一種距離參數。常用的距離參數是1-correlation。

# 計算距離（1 - correlation）
C <- cor(logcounts(deng))

# Run clustering based on the correlations, where the distance will
# be 1-correlation, e.g. higher distance with lower correlation.
distance_corr <- as.dist(1-C)

#使用ward linkage進行分層聚類
ward_hclust_corr <- hclust(distance_corr,method="ward.D2")
plot(ward_hclust_corr, main = "dist = 1-corr, Ward linkage")

再次設置樹狀圖參數生成10個亞群，並在PCA圖上繪製聚類標籤。

#設置樹狀圖參數生成10個亞群
cluster_hclust <- cutree(ward_hclust_corr,k = 10)
colData(deng)$cluster_hclust <- factor(cluster_hclust)

plot_grid(plotReducedDim(deng, "PCA", colour_by = "cluster_hclust"),
plotReducedDim(deng, "PCA", colour_by = "cell_type2"))

除了更改距離度量，我們還可以更改連結方法 —— 使用完全連結，而不是使用Ward的方法。

# 計算距離 (default: Eucledian distance)
distance_eucledian <- dist(t(logcounts(deng)))

#使用ward linkage進行分層聚類
comp_hclust_eucledian <- hclust(distance_eucledian,method = "complete")
plot(comp_hclust_eucledian, main = "dist = eucledian, complete linkage")

再一次，切割樹狀圖以生成10個細胞亞群，並在PCA圖上繪製亞群標籤。

#設置樹狀圖參數以生成10個細胞亞群
cluster_hclust <- cutree(comp_hclust_eucledian,k = 10)
colData(deng)$cluster_hclust <- factor(cluster_hclust)

plot_grid(plotReducedDim(deng, "PCA", colour_by = "cluster_hclust"),
plotReducedDim(deng, "PCA", colour_by = "cell_type2"))

tSNE + Kmeans

# 在tSNE坐標上執行kmeans算法
deng_kmeans <- kmeans(x = deng@reducedDims$TSNE,centers = 10)
TSNE_kmeans <- factor(deng_kmeans$cluster)
colData(deng)$TSNE_kmeans <- TSNE_kmeans
#Compare with ground truth
plot_grid(plotTSNE(deng, colour_by = "TSNE_kmeans"),
plotTSNE(deng, colour_by = "cell_type2"))

Graph Based Clustering

#k=5
#The k parameter defines the number of closest cells to look for each cells
SNNgraph_k5 <- buildSNNGraph(x = deng, k=5)
SNNcluster_k5 <- cluster_louvain(SNNgraph_k5)
colData(deng)$SNNcluster_k5 <- factor(SNNcluster_k5$membership)
p5<- plotPCA(deng, colour_by="SNNcluster_k5")+ guides(fill=guide_legend(ncol=2))

# k30
SNNgraph_k30 <- buildSNNGraph(x = deng, k=30)
SNNcluster_k30 <- cluster_louvain(SNNgraph_k30)
colData(deng)$SNNcluster_k30 <- factor(SNNcluster_k30$membership)
p30 <- plotPCA(deng, colour_by="SNNcluster_k30")

#plot the different clustering.
plot_grid(p5+ guides(fill=guide_legend(ncol=1)),p30)

Session info

sessionInfo()

## R version 3.5.3 (2019-03-11)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17763)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=English_United States.1252
## [2] LC_CTYPE=English_United States.1252
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C
## [5] LC_TIME=English_United States.1252
##
## attached base packages:
## [1] parallel stats4 stats graphics grDevices utils datasets
## [8] methods base
##

[1]：https://science.sciencemag.org/content/343/6167/193

撰文：Tiger
校對：生信寶典

你可能還想看

往期精品(點擊圖片直達文字對應教程)

後臺回復「生信寶典福利第一波」或點擊閱讀原文獲取教程合集