hello，大家好，要過禮拜了，但是今天還是要好好工作，今天我們來分享一個(gè)針對空間轉(zhuǎn)錄組的聚類算法，目前大多數(shù)多空間轉(zhuǎn)錄組的聚類還是基于找高變基因降維聚類，但是文獻(xiàn)中采用這個(gè)方法的很少，我們來看看新方法有什么改進(jìn)。文獻(xiàn)在Spatial transcriptomics at subspot resolution with BayesSpace,發(fā)表于Nature Biotechnology（IF54，相當(dāng)高了）。

單細(xì)胞轉(zhuǎn)錄組測序（scRNA-seq）實(shí)現(xiàn)了高通量和高分辨率的基因表達(dá)譜分析，但由于樣品制備過程中組織被分離，空間信息并沒有被保留。近年來，新興起的空間轉(zhuǎn)錄組學(xué)技術(shù)能夠在保留空間背景的同時(shí)全面檢測轉(zhuǎn)錄組譜，使人們對轉(zhuǎn)錄組表達(dá)的空間位置有了更深的認(rèn)知，并且為生物學(xué)功能和病理提供重要的見解。

目前，利用空間轉(zhuǎn)錄組學(xué)，人們已經(jīng)對腫瘤異質(zhì)性、腦功能等不同領(lǐng)域產(chǎn)生了更加深刻的認(rèn)識(shí)。但這些空間基因表達(dá)平臺(tái)仍然存在分辨率的技術(shù)限制。雖然后續(xù)開發(fā)的結(jié)合熒光原位雜交（FISH）技術(shù)的空間轉(zhuǎn)錄組測序方法以及其他方法，如Slide-seq和ZipSeq等提供了更高的分辨率，但大多數(shù)方法仍然存在通量較低、敏感性較低、依賴某種特定條件或不能夠被廣泛使用等限制問題。

此外，在測序技術(shù)發(fā)展的同時(shí)，還需要新的統(tǒng)計(jì)方法來分析空間基因表達(dá)數(shù)據(jù)，從而有效地利用現(xiàn)有的空間信息進(jìn)行科學(xué)研究。其中，聚類分群是此類數(shù)據(jù)分析的重要一步，也是下游深入分析的基礎(chǔ)，后續(xù)研究諸如細(xì)胞類型或組織注釋、差異表達(dá)等均依賴有效的分群。有學(xué)者提出了一種隱馬爾科夫隨機(jī)場模型（HMRF），通過對基因表達(dá)和空間鄰域結(jié)構(gòu)聯(lián)合建模，將低分辨率原位雜交數(shù)據(jù)聚類到不同的空間域，選擇空間差異表達(dá)基因進(jìn)行聚類。但是，現(xiàn)有的大多數(shù)空間基因表達(dá)數(shù)據(jù)分析方法往往依賴于非空間scRNA-seq數(shù)據(jù)的聚類方法，難以精確反映空間轉(zhuǎn)錄組所提供的信息。

近日，美國西雅圖弗雷德·哈欽森癌癥研究中心、華盛頓大學(xué)等單位的聯(lián)合研究團(tuán)隊(duì)在Nature Biotechnology發(fā)表了題為“Spatial transcriptomics at subspot resolution with BayesSpace”的文章。該研究提出了BayesSpace算法，這是一種利用空間轉(zhuǎn)錄組數(shù)據(jù)中的鄰域結(jié)構(gòu)來增加子點(diǎn)級別分辨率的計(jì)算方法，通過使用貝葉斯統(tǒng)計(jì)來實(shí)現(xiàn)超分辨率圖像分析。研究人員利用BayesSpace進(jìn)行空間聚類改進(jìn)了對空間分布的組織域的識(shí)別，提高了基因表達(dá)圖譜的分辨率，并可以重現(xiàn)接近單細(xì)胞分辨率的真實(shí)空間結(jié)構(gòu)。

圖片.png

主要研究內(nèi)容

BayesSpace基本原理

文章介紹，BayesSpace通過對低維的基因表達(dá)矩陣進(jìn)行建模，并通過空間先驗(yàn)知識(shí)誘導(dǎo)真實(shí)的鄰近點(diǎn)聚集，以此進(jìn)行推廣，從而實(shí)現(xiàn)空間聚類。這種方法來自于此前開發(fā)的用于圖像分析和微陣列數(shù)據(jù)的空間統(tǒng)計(jì)方法。與已有的方法相比，BayesSpace允許對聚類結(jié)構(gòu)和錯(cuò)誤項(xiàng)進(jìn)行更靈活的調(diào)整和規(guī)范。

圖片.png

圖1. BayesSpace基本原理，

BayesSpace空間聚類提高了對腦組織中已知皮層的認(rèn)識(shí)

為了檢測BayesSpace的性能，研究人員使用Maynard等人公開發(fā)表的12個(gè)背外側(cè)前額葉皮層 （DLPFC） 樣本的Visium空間表達(dá)譜數(shù)據(jù)，以及每個(gè)樣本的6個(gè)皮質(zhì)層和白質(zhì)的手工注釋，這些是作為R包spatialLIBD的一部分。借助該數(shù)據(jù)集，研究人員評估了BayesSpace識(shí)別不同皮層特定表達(dá)輪廓的能力，并將其性能與其他空間和非空間聚類方法進(jìn)行比較。

結(jié)果顯示，BayesSpace大大優(yōu)于原始的spatialLIBD聚類分區(qū)，以及為空間轉(zhuǎn)錄組數(shù)據(jù)開發(fā)的其他所有非空間聚類算法和空間聚類方法。除BayesSpace外，大多數(shù)聚類分區(qū)表現(xiàn)出大量的噪聲，且聚類之間缺乏明確的空間分隔。相比之下，BayesSpace利用了空間信息平滑數(shù)據(jù)，并提供不同的集群分層。此外，BayesSpace的運(yùn)行時(shí)和內(nèi)存占用與其他空間聚類方法相當(dāng)。

圖片.png

圖2. BayesSpace性能評估及比較

BayesSpace能夠鑒定出易被其他方法遺漏的組織結(jié)構(gòu)

隨后，研究人員還使用BayesSpace分析了由Thrane等人首次注釋和描述的黑色素瘤空間轉(zhuǎn)錄組樣本，這些數(shù)據(jù)包括了人工標(biāo)注識(shí)別的黑色素瘤、間質(zhì)和淋巴組織的區(qū)域。結(jié)果顯示，利用BayesSpace得到的4個(gè)空間聚類與手工標(biāo)注的組織類型相吻合。此外，BayesSpace增強(qiáng)的空間聚類提供了更高分辨率的組織類型圖，例如增強(qiáng)識(shí)別的腫瘤邊緣的淋巴樣區(qū)域和可能的免疫浸潤腫瘤區(qū)域。這些區(qū)域在最初的分辨率下是無法識(shí)別的，并且這些區(qū)域在很大程度上也沒有被其他聚類方法識(shí)別。

差異表達(dá)分析結(jié)果提示，淋巴區(qū)域有一個(gè)獨(dú)特的表達(dá)譜，其淋巴細(xì)胞標(biāo)志物CD52和MS4A1的表達(dá)升高，黑色素瘤標(biāo)志物MCAM和SPP1的表達(dá)相對于周圍腫瘤邊界的表達(dá)降低，四個(gè)聚類間的增強(qiáng)分辨率差異表達(dá)分析強(qiáng)調(diào)了基因表達(dá)的額外空間變異。

圖片.png

圖3. BayesSpace鑒定出黑色素瘤樣本中的腫瘤近端淋巴組織結(jié)構(gòu)

BayesSpace能夠區(qū)分浸潤性導(dǎo)管癌的瘤內(nèi)異質(zhì)性

接下來，研究人員進(jìn)一步分析了浸潤性導(dǎo)管癌（乳腺癌的一種病理分型）的組織切片以確定聚類分群的生物學(xué)相關(guān)性。每一張組織切片，病理學(xué)家注釋了主要的浸潤性癌、原位癌和良性增生的區(qū)域，從中可以得到每個(gè)點(diǎn)的真實(shí)標(biāo)簽。結(jié)果顯示，BayesSpace聚類與組織病理學(xué)注釋基本一致。

此外，如果沒有蘇木精和伊紅（H&E）染色或免疫熒光染色作為腫瘤標(biāo)志物，腫瘤-基質(zhì)界面在組織學(xué)上不能被完全描述。BayesSpace增強(qiáng)的聚類可以識(shí)別出組織內(nèi)的異質(zhì)性，并得到關(guān)鍵腫瘤標(biāo)記基因的明確數(shù)據(jù)支持。也就是說，已知腫瘤標(biāo)記基因的空間表達(dá)模式和這些聚類之間的差異表達(dá)分析與臨床和組織病理學(xué)注釋基本一致。例如，在整個(gè)腫瘤群體中觀察到ERBB2和ESR1基因的高表達(dá)水平；非腫瘤細(xì)胞群1、7和10的特征是免疫基因的表達(dá)，如PTPRC（白細(xì)胞共同抗原CD45） 的高表達(dá)。這些空間表達(dá)模式表明，侵襲性腫瘤存在明顯瘤內(nèi)轉(zhuǎn)錄異質(zhì)性，目前的組織病理學(xué)分析方法無法識(shí)別這些差異，表明空間轉(zhuǎn)錄組數(shù)據(jù)相對于單純免疫熒光的優(yōu)勢。

圖片.png

圖4. BayesSpace區(qū)分浸潤性導(dǎo)管癌的瘤內(nèi)異質(zhì)性

研究總結(jié)

綜上所述，該研究首次報(bào)道了BayesSpace這一基于空間轉(zhuǎn)錄組模型的聚類方法。BayesSpace使用t分布錯(cuò)誤模型來識(shí)別空間聚類，這些聚類對于技術(shù)噪聲引起的離群值的存在更加魯棒。隨后的應(yīng)用及分析結(jié)果證明，BayesSpace在識(shí)別具有相似表達(dá)譜的空間群體和提高空間轉(zhuǎn)錄組的分辨率方面有較高的實(shí)用性。

我們來看看代碼

準(zhǔn)備數(shù)據(jù)

BayesSpace supports three ways of loading a SingleCellExperiment for analysis.（這個(gè)是R語言常見的數(shù)據(jù)結(jié)構(gòu)）

Visium datasets processed with Space Ranger can be loaded directly via the readVisium() function. This function takes only the path to the Space Ranger output directory (containing the spatial/ and filtered_feature_bc_matrix/ subdirectories) and returns a SingleCellExperiment.（直接讀取10XSpaceranger的文件）.

sce <- readVisium("path/to/spaceranger/outs/")

其次，為 BayesSpace manuscript分析的所有數(shù)據(jù)集都可以通過 getRDS() 函數(shù)輕松訪問。 這個(gè)函數(shù)有兩個(gè)參數(shù)——數(shù)據(jù)集的名稱和數(shù)據(jù)集中樣本的名稱。

melanoma <- getRDS(dataset="2018_thrane_melanoma", sample="ST_mel1_rep2")

也可以從計(jì)數(shù)矩陣和行列數(shù)據(jù)表手動(dòng)構(gòu)建 SingleCellExperiment 對象。 BayesSpace 僅要求將點(diǎn)數(shù)組坐標(biāo)作為 colData 中名為 row 和 col 的列提供。 （注意，Visium 數(shù)據(jù)集的增強(qiáng)還需要組織圖像中每個(gè)點(diǎn)的像素坐標(biāo)，但在這種情況下，數(shù)據(jù)集應(yīng)使用 readVisium() 加載，它會(huì)自動(dòng)加載這些數(shù)據(jù)。）

library(Matrix)

rowData <- read.csv("path/to/rowData.csv", stringsAsFactors=FALSE)
colData <- read.csv("path/to/colData.csv", stringsAsFactors=FALSE, row.names=1)
counts <- read.csv("path/to/counts.csv.gz",
                   row.names=1, check.names=F, stringsAsFactors=FALSE))

sce <- SingleCellExperiment(assays=list(counts=as(counts, "dgCMatrix")),
                            rowData=rowData,
                            colData=colData)

前處理

BayesSpace 需要最少的數(shù)據(jù)預(yù)處理，提供了一個(gè)輔助函數(shù)來自動(dòng)化它。

spatialPreprocess() 對計(jì)數(shù)矩陣進(jìn)行對數(shù)歸一化，并對頂部 n.HVGs 高度可變的基因執(zhí)行 PCA，保留頂部 n.PCs 主成分。 此外，空間排序平臺(tái)作為meta數(shù)據(jù)添加到SingleCellExperiment 中，用于下游分析。 如果不想重新運(yùn)行 PCA，運(yùn)行帶有skip.PCA=TRUE 標(biāo)志的 spatialPreprocess() 只會(huì)添加 BayesSpace 需要的meta數(shù)據(jù)。

在這里，省略了對數(shù)歸一化，因?yàn)樗型ㄟ^ getRDS() 可用的數(shù)據(jù)集都已經(jīng)包含對數(shù)歸一化計(jì)數(shù)。

set.seed(102)
melanoma <- spatialPreprocess(melanoma, platform="ST", 
                              n.PCs=7, n.HVGs=2000, log.normalize=FALSE)

聚類

Selecting the number of clusters

我們可以使用 qTune() 和 qPlot() 函數(shù)來幫助選擇 q，即我們分析中要使用的聚類數(shù)。

• qTune() 為 q 的多個(gè)指定值（默認(rèn)為 3 到 7）運(yùn)行 BayesSpace 聚類算法，并計(jì)算它們的平均偽對數(shù)似然。 它接受spatialCluster() 的任何參數(shù)。
• qPlot() 將偽對數(shù)似然繪制為 q 的函數(shù)； 我們建議在該圖的肘部周圍選擇一個(gè) q。

melanoma <- qTune(melanoma, qs=seq(2, 10), platform="ST", d=7)
qPlot(melanoma)

圖片.png

Clustering with BayesSpace

spatialCluster() 函數(shù)對spot進(jìn)行聚類，并將預(yù)測的聚類標(biāo)簽添加到SingleCellExperiment。 通常，建議至少運(yùn)行 10,000 次迭代 (nrep=10000)，但為了運(yùn)行示例中使用了 1,000 次迭代。 （注意必須設(shè)置隨機(jī)種子才能使結(jié)果可重現(xiàn)。）

set.seed(149)
melanoma <- spatialCluster(melanoma, q=4, platform="ST", d=7,
                           init.method="mclust", model="t", gamma=2,
                           nrep=1000, burn.in=100,
                           save.chain=TRUE)

mclust 初始化 (cluster.init) 和 BayesSpace 集群分配 (spatial.cluster) 現(xiàn)在都可以在 SingleCellExperiment 的 colData 中使用。

head(colData(melanoma))
#> DataFrame with 6 rows and 5 columns
#>            row       col sizeFactor cluster.init spatial.cluster
#>      <integer> <integer>  <numeric>    <numeric>       <numeric>
#> 7x15         7        15   0.795588            1               1
#> 7x16         7        16   0.307304            1               1
#> 7x17         7        17   0.331247            2               2
#> 7x18         7        18   0.420747            3               2
#> 8x13         8        13   0.255453            1               1
#> 8x14         8        14   1.473439            1               1

可視化

clusterPlot(melanoma)

圖片.png

由于 clusterPlot() 返回一個(gè) ggplot 對象，因此可以通過組合熟悉的 ggplot2 函數(shù)對其進(jìn)行自定義。 此外，參數(shù)調(diào)色板設(shè)置用于每個(gè)簇的顏色，clusterPlot() 將附加參數(shù)傳遞給 geom_polygon()，例如大小或顏色，以控制斑點(diǎn)邊界。

clusterPlot(melanoma, palette=c("purple", "red", "blue", "yellow"), color="black") +
  theme_bw() +
  xlab("Column") +
  ylab("Row") +
  labs(fill="BayesSpace\ncluster", title="Spatial clustering of ST_mel1_rep2")

圖片.png

Enhanced resolution（加強(qiáng)精度）

Clustering at enhanced resolution

spatialEnhance() 函數(shù)將提高主成分的分辨率，并將這些 PC 以及子點(diǎn)分辨率的預(yù)測聚類標(biāo)簽添加到新的 SingleCellExperiment。 與我們上面對 spatialCluster() 的演示一樣，我們在本示例中使用的迭代次數(shù) (nrep=1000) 比我們在實(shí)踐中推薦的次數(shù) (nrep=100000 或更大) 少。

melanoma.enhanced <- spatialEnhance(melanoma, q=4, platform="ST", d=7,
                                    model="t", gamma=2,
                                    jitter_prior=0.3, jitter_scale=3.5,
                                    nrep=1000, burn.in=100,
                                    save.chain=TRUE)

The enhanced SingleCellExperiment includes an index to the parent spot in the original sce (spot.idx), along with an index to the subspot. It adds the offsets to the original spot coordinates, and provides the enhanced cluster label (spatial.cluster).

clusterPlot(melanoma.enhanced)

圖片.png

Enhancing the resolution of gene expression

ayesSpace 對基因表達(dá)矩陣的主要成分進(jìn)行運(yùn)算，因此 spatialEnhance() 計(jì)算增強(qiáng)分辨率的 PC 向量。增強(qiáng)的基因表達(dá)不是直接計(jì)算的，而是使用回歸算法估算的。對于每個(gè)基因，訓(xùn)練使用每個(gè)點(diǎn)的 PC 向量的模型來預(yù)測點(diǎn)級別的基因表達(dá)，并使用擬合模型來預(yù)測來自subspot PC 的subspot表達(dá)。

基因表達(dá)增強(qiáng)在enhanceFeatures() 函數(shù)中實(shí)現(xiàn)。 BayesSpace 默認(rèn)使用 xgboost 預(yù)測表達(dá)式，但也可以通過模型參數(shù)使用線性和狄利克雷回歸。使用 xgboost 時(shí)，建議通過將 nrounds 參數(shù)設(shè)置為 0 來自動(dòng)調(diào)整它，盡管這是以增加運(yùn)行時(shí)間為代價(jià)的（實(shí)際上比預(yù)先指定的 nrounds 慢約 4 倍）。

EnhanceFeatures() 可用于估算所有基因或感興趣基因子集的subspot水平表達(dá)。在這里，我們將通過增強(qiáng)四種標(biāo)記基因的表達(dá)來證明：PMEL（黑色素瘤）、CD2（T 細(xì)胞）、CD19（B 細(xì)胞）和 COL1A1（成纖維細(xì)胞）。

markers <- c("PMEL", "CD2", "CD19", "COL1A1")
melanoma.enhanced <- enhanceFeatures(melanoma.enhanced, melanoma,
                                     feature_names=markers,
                                     nrounds=0)

By default, log-normalized expression (logcounts(sce)) is imputed, although other assays or arbitrary feature matrices can be specified.

logcounts(melanoma.enhanced)[markers, 1:5]
#>        subspot_1.1 subspot_2.1 subspot_3.1 subspot_4.1 subspot_5.1
#> PMEL     2.6428437   1.8550344   2.4704804   3.1827958   2.3572879
#> CD2      0.3489273   0.6066852   0.2315192   0.2210583   0.3489273
#> CD19     0.6170074   0.6528370   0.4164957   0.2558892   0.6259792
#> COL1A1   0.0000000   2.9053805   1.1940085   0.1023711   0.8157212

Diagnostic measures from each predictive model, such as rmse when using xgboost, are added to the rowData of the enhanced dataset.

rowData(melanoma.enhanced)[markers, ]
#> DataFrame with 4 rows and 4 columns
#>                gene_id   gene_name    is.HVG enhanceFeatures.rmse
#>            <character> <character> <logical>            <numeric>
#> PMEL   ENSG00000185664        PMEL      TRUE             0.804628
#> CD2    ENSG00000116824         CD2      TRUE             0.614575
#> CD19   ENSG00000177455        CD19      TRUE             0.697328
#> COL1A1 ENSG00000108821      COL1A1      TRUE             0.704845

可視化enhanced gene expression

featurePlot(melanoma.enhanced, "PMEL")

圖片.png

enhanced.plots <- purrr::map(markers, function(x) featurePlot(melanoma.enhanced, x))
patchwork::wrap_plots(enhanced.plots, ncol=2)

圖片.png

spot.plots <- purrr::map(markers, function(x) featurePlot(melanoma, x))
patchwork::wrap_plots(c(enhanced.plots, spot.plots), ncol=4)

圖片.png

方法真的不錯(cuò)，最后，大家周末愉快

生活很好，有你更好