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

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

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

此外，在測序技術(shù)發(fā)展的同時尤蒿，還需要新的統(tǒng)計方法來分析空間基因表達(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)隊在Nature Biotechnology發(fā)表了題為“Spatial transcriptomics at subspot resolution with BayesSpace”的文章。該研究提出了BayesSpace算法齿兔，這是一種利用空間轉(zhuǎn)錄組數(shù)據(jù)中的鄰域結(jié)構(gòu)來增加子點級別分辨率的計算方法橱脸，通過使用貝葉斯統(tǒng)計來實現(xiàn)超分辨率圖像分析础米。研究人員利用BayesSpace進(jìn)行空間聚類改進(jìn)了對空間分布的組織域的識別，提高了基因表達(dá)圖譜的分辨率添诉，并可以重現(xiàn)接近單細(xì)胞分辨率的真實空間結(jié)構(gòu)屁桑。

圖片.png

主要研究內(nèi)容

BayesSpace基本原理

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

圖片.png

圖1. BayesSpace基本原理准浴，

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

為了檢測BayesSpace的性能，研究人員使用Maynard等人公開發(fā)表的12個背外側(cè)前額葉皮層 （DLPFC） 樣本的Visium空間表達(dá)譜數(shù)據(jù)捎稚，以及每個樣本的6個皮質(zhì)層和白質(zhì)的手工注釋乐横，這些是作為R包spatialLIBD的一部分。借助該數(shù)據(jù)集今野，研究人員評估了BayesSpace識別不同皮層特定表達(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)行時和內(nèi)存占用與其他空間聚類方法相當(dāng)灵巧。

圖片.png

圖2. BayesSpace性能評估及比較

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

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

差異表達(dá)分析結(jié)果提示隅要，淋巴區(qū)域有一個獨特的表達(dá)譜，其淋巴細(xì)胞標(biāo)志物CD52和MS4A1的表達(dá)升高董济，黑色素瘤標(biāo)志物MCAM和SPP1的表達(dá)相對于周圍腫瘤邊界的表達(dá)降低步清，四個聚類間的增強(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ū)域谴轮，從中可以得到每個點的真實標(biāo)簽。結(jié)果顯示吹埠，BayesSpace聚類與組織病理學(xué)注釋基本一致第步。

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

圖片.png

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

研究總結(jié)

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

我們來看看代碼

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

BayesSpace supports three ways of loading a SingleCellExperiment for analysis.（這個是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ù)輕松訪問策治。 這個函數(shù)有兩個參數(shù)——數(shù)據(jù)集的名稱和數(shù)據(jù)集中樣本的名稱脓魏。

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

也可以從計數(shù)矩陣和行列數(shù)據(jù)表手動構(gòu)建 SingleCellExperiment 對象。 BayesSpace 僅要求將點數(shù)組坐標(biāo)作為 colData 中名為 row 和 col 的列提供通惫。 （注意茂翔，Visium 數(shù)據(jù)集的增強(qiáng)還需要組織圖像中每個點的像素坐標(biāo)，但在這種情況下履腋，數(shù)據(jù)集應(yīng)使用 readVisium() 加載珊燎，它會自動加載這些數(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ù)處理遵湖，提供了一個輔助函數(shù)來自動化它悔政。

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

在這里，省略了對數(shù)歸一化牡整，因為所有通過 getRDS() 可用的數(shù)據(jù)集都已經(jīng)包含對數(shù)歸一化計數(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 的多個指定值（默認(rèn)為 3 到 7）運(yùn)行 BayesSpace 聚類算法谣辞，并計算它們的平均偽對數(shù)似然。 它接受spatialCluster() 的任何參數(shù)沐扳。
• qPlot() 將偽對數(shù)似然繪制為 q 的函數(shù)泥从； 我們建議在該圖的肘部周圍選擇一個 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() 返回一個 ggplot 對象帆阳，因此可以通過組合熟悉的 ggplot2 函數(shù)對其進(jìn)行自定義。 此外屋吨，參數(shù)調(diào)色板設(shè)置用于每個簇的顏色蜒谤，clusterPlot() 將附加參數(shù)傳遞給 geom_polygon()山宾，例如大小或顏色，以控制斑點邊界芭逝。

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 以及子點分辨率的預(yù)測聚類標(biāo)簽添加到新的 SingleCellExperiment渊胸。 與我們上面對 spatialCluster() 的演示一樣旬盯，我們在本示例中使用的迭代次數(shù) (nrep=1000) 比我們在實踐中推薦的次數(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() 計算增強(qiáng)分辨率的 PC 向量胖翰。增強(qiáng)的基因表達(dá)不是直接計算的，而是使用回歸算法估算的切厘。對于每個基因萨咳，訓(xùn)練使用每個點的 PC 向量的模型來預(yù)測點級別的基因表達(dá)，并使用擬合模型來預(yù)測來自subspot PC 的subspot表達(dá)疫稿。

基因表達(dá)增強(qiáng)在enhanceFeatures() 函數(shù)中實現(xiàn)培他。 BayesSpace 默認(rèn)使用 xgboost 預(yù)測表達(dá)式，但也可以通過模型參數(shù)使用線性和狄利克雷回歸遗座。使用 xgboost 時舀凛，建議通過將 nrounds 參數(shù)設(shè)置為 0 來自動調(diào)整它，盡管這是以增加運(yùn)行時間為代價的（實際上比預(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

方法真的不錯宽堆，最后腌紧，大家周末愉快

生活很好，有你更好