聚類可視化

?

NGS系列文章包括NGS基礎、轉錄組分析?（Nature重磅綜述|關于RNA-seq你想知道的全在這）、ChIP-seq分析?（ChIP-seq基本分析流程）、單細胞測序分析?(重磅綜述：三萬字長文讀懂單細胞RNA測序分析的最佳實踐教程 （原理、代碼和評述）)、DNA甲基化分析、重測序分析、GEO數據挖掘（典型醫學設計實驗GEO數據分析 (step-by-step) - Limma差異分析、火山圖、功能富集）等內容。

?

?

?

我們在單細胞轉錄組分析中最為常用的聚類可視化即為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

創作挑戰賽新人創作獎勵來咯，堅持創作打卡瓜分現金大獎

總結

以上是生活随笔為你收集整理的代码分析 | 单细胞转录组clustering详解的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。