聚類可視化

NGS系列文章包括NGS基礎(chǔ)、轉(zhuǎn)錄組分析?（Nature重磅綜述|關(guān)于RNA-seq你想知道的全在這）谨设、ChIP-seq分析?（ChIP-seq基本分析流程）蜓谋、單細(xì)胞測序分析?(重磅綜述：三萬字長文讀懂單細(xì)胞RNA測序分析的最佳實踐教程 （原理担忧、代碼和評述）)馁启、DNA甲基化分析岗喉、重測序分析云稚、GEO數(shù)據(jù)挖掘（典型醫(yī)學(xué)設(shè)計實驗GEO數(shù)據(jù)分析 (step-by-step) - Limma差異分析、火山圖沈堡、功能富集）等內(nèi)容。

我們在單細(xì)胞轉(zhuǎn)錄組分析中最為常用的聚類可視化即為tSNE和UMAP（Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（十二）- Scater單細(xì)胞表達(dá)譜tSNE可視化）燕雁，不過非線性可視化方法（例如t-SNE）通常會擾亂數(shù)據(jù)中的全局結(jié)構(gòu)诞丽。diffusion maps能夠很好的表達(dá)局部和全局結(jié)構(gòu)，但無法針對二維（2D）可視化進(jìn)行優(yōu)化拐格，因為它們不是專門為可視化設(shè)計的僧免。（其實即便如此，我現(xiàn)在也習(xí)慣性應(yīng)用diffusion maps捏浊，可以說明更多的信息懂衩，雖然有時候確實是很難看。金踪。浊洞。）

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

?

在本教程中，我們將研究不同的聚類scRNA-seq數(shù)據(jù)集方法以表征不同的細(xì)胞亞群靶瘸。

所需加載包

suppressMessages(require(tidyverse))

suppressMessages(require(Seurat))

suppressMessages(require(cowplot))

suppressMessages(require(scater))

suppressMessages(require(scran))

suppressMessages(require(igraph))

數(shù)據(jù)集

本教程使用的是來自發(fā)育中的小鼠胚胎的小細(xì)胞數(shù)據(jù)集[1]苫亦。該數(shù)據(jù)集已進(jìn)行了預(yù)處理、創(chuàng)建了SingleCellExperiment對象并對細(xì)胞進(jìn)行了注釋（cellassign：用于腫瘤微環(huán)境分析的單細(xì)胞注釋工具（9月Nature））怨咪。

數(shù)據(jù)下載鏈接：https://github.com/NBISweden/excelerate-scRNAseq/tree/master/session-clustering/session-clustering_files屋剑。

#加載表達(dá)矩陣

deng <- readRDS("session-clustering_files/deng-reads.rds")

deng

?

#查看細(xì)胞類型注釋

table(colData(deng)$cell_type2)

?

特征選擇

第一步是決定在細(xì)胞聚類中使用哪些基因（Hemberg-lab單細(xì)胞轉(zhuǎn)錄組數(shù)據(jù)分析（十）- Scater基因評估和過濾）。單細(xì)胞RNA-seq可以分析細(xì)胞中的大量基因诗眨，大多數(shù)基因的表達(dá)不足以提供有意義的信號并且通常受到技術(shù)噪音的干擾唉匾，其中可能會添加一些不必要的信號，從而使生物變異模糊匠楚，而基因過濾還可以加快下游分析的計算時間肄鸽。

首先看一下所有基因的表達(dá)平均值和方差。哪些基因似乎不太重要油啤，哪些可能是技術(shù)噪音典徘？

#計算整個細(xì)胞的基因平均值

gene_mean <- rowMeans(counts(deng))

#計算整個細(xì)胞的基因方差

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()

?

濾除低豐度基因

低豐度基因大多無信息，不能代表數(shù)據(jù)的生物學(xué)差異益咬。它們通常是由技術(shù)噪音（例如dropout）導(dǎo)致逮诲，在下游分析中的存在通常會導(dǎo)致準(zhǔn)確性降低帜平，因為它們會干擾將要使用的某些統(tǒng)計模型，并且會毫無意義地增加計算時間梅鹦，這在處理非常大的數(shù)據(jù)時可能至關(guān)重要裆甩。

abundant_genes <- gene_mean >= 0.5 #去除低豐度基因

# 繪制基因表達(dá)分布

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)

?

過濾在很少細(xì)胞中表達(dá)的基因

我們還可以過濾少量細(xì)胞中的某些基因。此過程將刪除一些在一兩個細(xì)胞中高度表達(dá)的異称胨簦基因嗤栓。這些基因不需要進(jìn)一步分析，因為它們主要來自人為的不規(guī)則擴(kuò)增箍邮。值得注意的是茉帅，當(dāng)分析的目的是檢測數(shù)據(jù)中非常稀有的亞群時，我們可能不希望進(jìn)行此過程锭弊。

#計算每個非零表達(dá)基因的數(shù)量

numcells <- nexprs(deng, byrow=TRUE)

#過濾不到5個細(xì)胞中檢測到的基因

numcells2 <- numcells >= 5

deng <- deng[numcells2,]

dim(deng)

檢測高變的基因

highly variable gene（HVG）假設(shè)基因在細(xì)胞之間的表達(dá)差異很大堪澎，則其中的一些差異是由于細(xì)胞之間的生物學(xué)差異而不是技術(shù)噪音引起的。但是味滞，基因的平均表達(dá)與不同細(xì)胞reads計數(shù)的差異之間存在正相關(guān)關(guān)系樱蛤。如果僅保留高變異基因，將保留許多在每個細(xì)胞中表達(dá)但不能代表生物學(xué)變異的高表達(dá)管家基因剑鞍。因此為了正確識別HVG昨凡，必須糾正此關(guān)系。

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

選項1：將變異系數(shù)建模為平均值的函數(shù)土匀。

# 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：根據(jù)平均數(shù)對生物成分的方差建模。

首先形用，我們估算每個基因表達(dá)的差異就轧，然后將差異分解為生物學(xué)和技術(shù)成分。然后將HVGs鑒定為具有最大生物學(xué)成分的基因田度。這避免了優(yōu)先排序?qū)τ捎诩夹g(shù)因素（例如在RNA捕獲和文庫制備過程中的采樣噪聲）而高度可變的基因妒御。詳細(xì)內(nèi)容查看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

降維

由于高頻率的噪聲（技術(shù)和生物學(xué)的）和大量的維度（即基因）導(dǎo)致聚類在計算上具有一定的困難，因此需要通過應(yīng)用降維方法（例如PCA镇饺、tSNE和UMAP）解決這些問題乎莉。

#PCA (選擇組分?jǐn)?shù)量)

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是隨機(jī)方法。每次運(yùn)行它時都會得到略有不同的結(jié)果奸笤。為了方便起見可以設(shè)置隨機(jī)種子惋啃，這樣就可以獲得相同的結(jié)果。

#tSNE

deng <- runTSNE(deng,

? ? ? ? ? ? ? perplexity = 30,

? ? ? ? ? ? ? feature_set = metadata(deng)$hvg_genes,

? ? ? ? ? ? ? set.seed = 1)

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

?

聚類

層次聚類

#計算距離（默認(rèn)值：歐幾里得距離）

distance_eucledian <- dist(t(logcounts(deng)))

#使用ward linkage執(zhí)行分層聚類

ward_hclust_eucledian <- hclust(distance_eucledian,method = "ward.D2")

plot(ward_hclust_eucledian, main = "dist = eucledian, Ward linkage")

?

我反正是啥也看不見监右。边灭。。但是看起來層次分類蠻明顯的健盒。绒瘦。称簿。

現(xiàn)在，設(shè)置樹狀圖參數(shù)生成10個亞群惰帽，并在PCA圖上繪制聚類標(biāo)簽憨降。

#設(shè)置樹狀圖參數(shù)生成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圖上繪制聚類標(biāo)簽。

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

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

?

我們嘗試另一種距離參數(shù)该酗。常用的距離參數(shù)是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進(jìn)行分層聚類

ward_hclust_corr <- hclust(distance_corr,method="ward.D2")

plot(ward_hclust_corr, main = "dist = 1-corr, Ward linkage")

?

再次設(shè)置樹狀圖參數(shù)生成10個亞群，并在PCA圖上繪制聚類標(biāo)簽呜魄。

#設(shè)置樹狀圖參數(shù)生成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進(jìn)行分層聚類

comp_hclust_eucledian <- hclust(distance_eucledian,method = "complete")

plot(comp_hclust_eucledian, main = "dist = eucledian, complete linkage")

?

再一次，切割樹狀圖以生成10個細(xì)胞亞群膳殷，并在PCA圖上繪制亞群標(biāo)簽操骡。

#設(shè)置樹狀圖參數(shù)以生成10個細(xì)胞亞群

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坐標(biāo)上執(zhí)行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