hello，大家好，今天給大家?guī)硪粋€(gè)10X空間轉(zhuǎn)錄組去除污染的方法，具體怎么做，大家往下看，文章在SpotClean adjusts for spot swapping in spatial transcriptomics data其實(shí)關(guān)于去除污染的方法，很多都是針對(duì)單細(xì)胞數(shù)據(jù)游離污染的，大家可以參考我之前講的公開課。

企業(yè)微信截圖_16353858689413.png

好了，我們來看看10X空間轉(zhuǎn)錄組的污染以及如何去除。

summary

空間轉(zhuǎn)錄組學(xué) (ST) 是一種強(qiáng)大且廣泛使用的方法，用于分析組織中的全基因組基因表達(dá)，在分子醫(yī)學(xué)和腫瘤診斷中具有新興應(yīng)用。 最近的空間轉(zhuǎn)錄組學(xué)實(shí)驗(yàn)利用包含數(shù)千個(gè)點(diǎn)的載玻片，這些點(diǎn)具有結(jié)合 mRNA 的點(diǎn)特異性條形碼。 理想情況下，一個(gè)點(diǎn)的唯一分子標(biāo)識(shí)符測量點(diǎn)特定的表達(dá)，但由于附近點(diǎn)的bleed，這通常不是這種理想情況，可以將其稱為點(diǎn)交換。 所以建議使用SpotClean調(diào)整spot swapping ，并在此過程中提高進(jìn)行下游分析的靈敏度和精度。 （也就是說空間轉(zhuǎn)錄組的污染主要來自于周周圍spot mRNA的游離）。

簡介

在典型的 ST 實(shí)驗(yàn)中，將新鮮冷凍（或 FFPE）組織切片并放置在含有spot的載玻片上，每個(gè)spot包含數(shù)百萬個(gè)具有該spot獨(dú)特空間條形碼的捕獲寡核苷酸。組織成像，通常通過蘇木精和伊紅 (H&E) 染色。成像后，組織被透化以釋放 mRNA，然后與捕獲寡核苷酸結(jié)合，生成一個(gè) cDNA 文庫，該文庫由保留空間信息的條形碼結(jié)合的轉(zhuǎn)錄本組成。 ST 實(shí)驗(yàn)的數(shù)據(jù)由組織圖像和從每個(gè)點(diǎn)收集的 RNA 測序數(shù)據(jù)組成。處理 ST 數(shù)據(jù)的第一步是組織檢測，其中包含組織的載玻片上的點(diǎn)與沒有組織的背景點(diǎn)區(qū)分開來。然后在包含組織的每個(gè)點(diǎn)上的唯一分子標(biāo)識(shí)符 (UMI) 計(jì)數(shù)用于下游分析 (下圖)。

圖片.png

理想情況下，給定spot的基因特異性 UMI 將代表該基因在該點(diǎn)的表達(dá)，而沒有組織的spot將不顯示 UMI。實(shí)際情況并非如此。從附近spot流出的信使 RNA 會(huì)導(dǎo)致 UMI 計(jì)數(shù)的大量污染，這里將這種現(xiàn)象稱為點(diǎn)交換。下圖顯示了死后人類大腦組織樣本中的點(diǎn)交換證據(jù)，該組織樣本被描述為空間 LIBD 的一部分，該項(xiàng)目旨在定義六層人類背外側(cè)前額葉皮層 (DLPFC) 中基因表達(dá)的。

圖片.png

• 注： Data from the human dorsolateral prefrontal cortex profiled in the spatialLIBD experiment,sample LIBD_151507. (a) UMI count densities for tissue and background spots show relatively high counts in the background. (b) UMI total counts in the background decrease with increasing distance from the tissue;the perimeter delineating tissue and background is shown in white. (c) Counts of the top 50 genes from a select tissue region (upper), from a nearby background region (middle), and from a distant background region (bottom) show the similarity between expression in tissue spots and nearby background spots due to spot swapping from tissue to background, an effect that decreases as distance from the tissue increases. (d) Tissue and background spots are not distinguished visually via UMAP. (e) Graph-based clustering of all spots identifies 9 clusters. (f) Spots on the slide are colored by their cluster membership shown in (e). Black arrows highlight areas of spot swapping of signal from tissue to background. Spots on the perimeter (shown in white) have been removed from the summaries shown here to ensure that the effects shown are not due to spots on the tissue-background boundary.

具體來說，上圖a 顯示背景點(diǎn)的 UMI 計(jì)數(shù)（在沒有污染的情況下為零）與組織點(diǎn)的計(jì)數(shù)相比較高；并且計(jì)數(shù)隨著與組織距離的增加而減少（上圖b）。上圖c 顯示了組織區(qū)域、附近背景區(qū)域和遠(yuǎn)處背景區(qū)域中 50 個(gè)基因的 UMI 計(jì)數(shù)分布。由于組織和附近背景之間的表達(dá)相似性，組織和背景斑點(diǎn)不容易區(qū)分（上圖d）。這在上圖f 中再次得到強(qiáng)調(diào)，where spots on the slide are 46 colored by membership in the graph-based clusters shown in 上圖e。其他數(shù)據(jù)的測試顯示了類似結(jié)果；下表顯示，在大多數(shù)數(shù)據(jù)集中，背景點(diǎn)中 UMI 計(jì)數(shù)的比例在 5% 到 20% 之間。

圖片.png

下面又是一些證明的例子

圖片.png

圖片.png

圖片.png

圖片.png

上面的結(jié)果表明從組織到背景發(fā)生點(diǎn)交換，但評(píng)估從組織spot到組織spot的點(diǎn)交換程度更具挑戰(zhàn)性。雖然 SpotClean 模型提供了一個(gè)估計(jì)值（下表），

圖片.png

但也考慮了在 spatialLIBD 項(xiàng)目中識(shí)別的組織特異性標(biāo)記基因。在沒有點(diǎn)交換的情況下，層特定標(biāo)記的表達(dá)在該層中應(yīng)該很高，而在其他層中應(yīng)該很低（或關(guān)閉）。當(dāng)點(diǎn)交換發(fā)生時(shí)，附近層中的標(biāo)記表達(dá)相對(duì)較高。這在 GFAP 中很明顯，例如，一種已知在白質(zhì) (WM) 和 DLPFC 的第一個(gè)注釋層 (Layer1) 中上調(diào)的標(biāo)記。下圖顯示了 GFAP 在 WM 和 Layer1 斑點(diǎn)中的高表達(dá)，

圖片.png

正如預(yù)期的那樣，但在與 WM 和 Layer1 相鄰的組織spot中也有相對(duì)高的表達(dá)，GFAP 表達(dá)隨著與 WM（或 Layer1）的距離增加而降低。雖然相鄰組織spot中標(biāo)記物表達(dá)的一些增加可能是由于這些斑點(diǎn)處存在 WM（或 Layer1）細(xì)胞，但應(yīng)該注意到表達(dá)速率衰減到背景斑點(diǎn)（其中不存在細(xì)胞）類似于衰減到相鄰組織區(qū)域的速率。因此，相鄰組織斑點(diǎn)中可能存在的 WM（或 Layer1）細(xì)胞不足以完全解釋觀察到的表達(dá)模式。 WM 標(biāo)記、MOBP（上圖）以及 13 個(gè)附加標(biāo)記（下圖）顯示了類似的結(jié)果。

圖片.png

為了更直接地量化點(diǎn)交換的程度，這里進(jìn)行了嵌合實(shí)驗(yàn)，其中人類和小鼠組織在樣品制備過程中連續(xù)放置。對(duì)于每個(gè)實(shí)驗(yàn)，我們對(duì) H&E 圖像進(jìn)行注釋以識(shí)別物種特異性區(qū)域，并計(jì)算點(diǎn)交換讀數(shù)的比例（人類spot中的小鼠特異性讀數(shù)、小鼠spot中的人類特異性讀數(shù)和背景spot中的讀數(shù)）。這是點(diǎn)交換讀數(shù) (LPSS) 比例的下限，因?yàn)樗豢紤]物種內(nèi)的點(diǎn)交換（例如，來自人類點(diǎn) t' 處的探針結(jié)合的人類點(diǎn) t 的讀數(shù)）；在這些實(shí)驗(yàn)中，LPSS 的范圍在 26-37% 之間。總之，組織和背景表達(dá)的比較、標(biāo)記基因分析和嵌合實(shí)驗(yàn)（下圖）表明點(diǎn)交換影響 ST 實(shí)驗(yàn)中的 UMI 計(jì)數(shù)。

圖片.png

這種令人討厭的變異性降低了下游分析的能力和精度（下圖）。

圖片.png

圖片.png

圖片.png

為調(diào)整 RNA-seq 實(shí)驗(yàn)中已知污染源而開發(fā)的統(tǒng)計(jì)方法不適應(yīng)點(diǎn)交換中固有的空間依賴性，因此，在此設(shè)置中是不夠的。為了調(diào)整 ST 實(shí)驗(yàn)中點(diǎn)交換的影響，作者開發(fā)了 SpotClean。該方法在 R 包 R/spotClean 中實(shí)現(xiàn)。 SpotClean 是根據(jù)模擬和案例研究數(shù)據(jù)進(jìn)行評(píng)估的。在 SimI 中，假設(shè)局部污染遵循Gaussian kernel，則生成污染計(jì)數(shù)； SimII-IV 放寬了高斯假設(shè)。在 SimV 中，對(duì)平均表達(dá)在載玻片上有系統(tǒng)變化的基因的污染計(jì)數(shù)進(jìn)行了模擬。下表顯示了模擬數(shù)據(jù)集中真實(shí)基因表達(dá)和去污染基因表達(dá)之間的均方誤差 (MSE)，

圖片.png

圖片.png

圖片.png

圖片.png

這些結(jié)果表明 SpotClean 提供了更好的表達(dá)估計(jì)；和下圖表明，SpotClean 表達(dá)估計(jì)提高了識(shí)別空間變化基因的精度。

圖片.png

SpotClean 對(duì)下游分析的好處也在案例研究數(shù)據(jù)中得到了說明。 具體來說，SpotClean 增加了標(biāo)記基因表達(dá)的特異性，增加了識(shí)別 DE 基因的能力，并提高了斑點(diǎn)注釋的準(zhǔn)確性。 下圖表明SpotClean 通過維持 WM 和 Layer1 中的表達(dá)水平并減少其他層中的虛假表達(dá)，提高了 SpatialLIBD 數(shù)據(jù)中 GFAP 的特異性。

圖片.png

圖片.png

上圖b 和下圖考慮了已知在原始和 SpotClean 凈化數(shù)據(jù)中 WM 和 Layer6 之間差異表達(dá)（DE）的基因；

圖片.png

SpotClean 導(dǎo)致已知 DE 基因的倍數(shù)變化增加和 p 值減小。 嵌合數(shù)據(jù)集提供了額外的例子。 特別是，SpotClean 減少了嵌合數(shù)據(jù)集中點(diǎn)交換 UMI 計(jì)數(shù)的比例。 類似的結(jié)果顯示在上圖中，考慮到人類特異性和小鼠特異性基因在人類特異性和小鼠特異性位點(diǎn)的表達(dá)。 通過 SpotClean 凈化的數(shù)據(jù)顯示小鼠組織中人類基因的表達(dá)減少，而人體組織沒有減少，反之亦然。

將空間轉(zhuǎn)錄組學(xué)應(yīng)用于個(gè)性化醫(yī)療引起了相當(dāng)大的興趣，例如患者腫瘤活檢的分子分析以指導(dǎo)診斷和精準(zhǔn)治療。 SpotClean 在精確點(diǎn)注釋至關(guān)重要的此類應(yīng)用中展示了巨大的優(yōu)勢(shì)。下圖顯示了人類乳腺癌樣本（導(dǎo)管癌），其中腫瘤的診斷、范圍和侵襲性通常是通過病理學(xué)家對(duì) H&E 圖像的評(píng)估來估計(jì)的。

圖片.png

空間轉(zhuǎn)錄組學(xué)可以提供額外的信息，包括識(shí)別惡性細(xì)胞的細(xì)微集合，但需要準(zhǔn)確的spot注釋才能使這些信息在臨床實(shí)踐中有用，尤其是不要過度強(qiáng)調(diào)腫瘤負(fù)荷。上圖顯示了使用 SpotClean 數(shù)據(jù)注釋的spot與使用未通過 SpotClean 凈化的數(shù)據(jù)注釋的spot。未經(jīng)凈化的數(shù)據(jù)將許多spot誤認(rèn)為是惡性的，包括那些在腫瘤周圍含有良性炎癥細(xì)胞的spot，而 SpotClean 凈化的數(shù)據(jù)更接近于在 H&E 圖像上識(shí)別惡性細(xì)胞。下圖顯示，在沒有 SpotClean 的情況下，原始數(shù)據(jù)中超過 13% 的標(biāo)記為惡性的spot可能是由于spot交換造成的誤判。

圖片.png

空間轉(zhuǎn)錄組學(xué)為解決生物學(xué)問題和加強(qiáng)患者護(hù)理提供了前所未有的機(jī)會(huì)，但必須調(diào)整點(diǎn)交換引起的偽影，以確保從這些強(qiáng)大的實(shí)驗(yàn)中獲得最大的信息。 SpotClean 提供更準(zhǔn)確的表達(dá)估計(jì)，從而提高下游分析的能力和精度。

最后看看示例代碼

安裝和加載

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("SpotClean")
library(SpotClean)

Short Demo

在這里，在捆綁的示例數(shù)據(jù)上快速演示了一般的 SpotClean 工作流程。 然后在下面進(jìn)行分步說明。

# Not run

# Load example data
data(mbrain_raw)
data(mbrain_slide_info)

# Visualize raw data
mbrain_obj <- CreateSlide(count_mat = mbrain_raw, 
                          slide_info = mbrain_slide_info)
VisualizeSlide(slide_obj = mbrain_obj)
VisualizeHeatmap(mbrain_obj,rownames(mbrain_raw)[1])

# Decontaminate raw data
decont_obj <- SpotClean(mbrain_obj)

# Visualize decontaminated gene
VisualizeHeatmap(decont_obj,rownames(mbrain_raw)[1])

# Visualize the estimated per-spot contamination rate
VisualizeHeatmap(decont_obj,decont_obj@metadata$contamination_rate, 
                 logged = FALSE, legend_title = "contamination rate",
                 legend_range = c(0,1))

# (Optionally) Transform to Seurat object for downstream analyses
seurat_obj <- ConvertToSeurat(decont_obj,image_dir = "/path/to/spatial/folder")

Running Speed

計(jì)算速度與輸入數(shù)據(jù)集的大小和結(jié)構(gòu)有關(guān)，主要由組織spot的數(shù)量驅(qū)動(dòng)。 SpotClean 不需要并行計(jì)算，因此不會(huì)占用過多 CPU 或內(nèi)存資源。 作為參考，SpotClean 在默認(rèn)基因過濾下在中等大小的數(shù)據(jù)集（大約 30,000 個(gè)基因和 2,000 個(gè)組織spot）上運(yùn)行所需的時(shí)間不到 15 分鐘。

Situations you should think twice about before applying SpotClean

• $組織斑點(diǎn)過多（背景斑點(diǎn)不夠)：雖然觀察到的數(shù)據(jù)是具有固定列數(shù)（spot）的單個(gè)矩陣，但未知參數(shù)的數(shù)量與組織spot的數(shù)量成正比。 在所有spot都被組織覆蓋的極端情況下，有比觀察到的數(shù)據(jù)值更多的未知參數(shù)。 在這種情況下，受污染的表達(dá)式與真實(shí)表達(dá)式混淆，SpotClean 估計(jì)變得不可靠。 因此建議輸入數(shù)據(jù)至少有 25% 的斑點(diǎn)未被組織占據(jù)，以便 SpotClean 從背景斑點(diǎn)中獲得足夠的信息來估計(jì)污染。
• Lowly-expressed genes： Lowly-expressed genes typically contain relatively less information and relatively more noise than highly-expressed genes. SpotClean by default only keeps highly-expressed and highly-variable genes for decontamination. It can be forced to run on manually-specified lowly-expressed genes. However, even in this case, expression for the lowly-expressed genes is typically not changed very much. Given the high sparsity in most lowly expressed genes, there is not enough information available to confidently reassign UMIs in most cases. However, we do not filter genes by sparsity because there can be interesting genes highly concentrated in a small tissue region. In cases like this, SpotClean is effective at adjusting for spot swapping in these regions. If the defaults are not appropriate, users can either adjust the expression cutoffs or manually specify genes to decontaminate.
• Inference based on sequencing depth： SpotClean reassigns bled-out UMIs to their tissue spots of origin which changes the estimated sequencing depth of tissue spots after decontamination, since most estimations of sequencing depth rely on total expressions at every spot. As a result, decontamination can be considered as another type of normalization and might conflict with existing sequencing depth normalization methods.

Recommended applications

SpotClean 通過校正點(diǎn)交換來改進(jìn)表達(dá)的估計(jì)。 換句話說，SpotClean 通過增強(qiáng)高表達(dá)區(qū)域的信號(hào)并減少未表達(dá)區(qū)域的測量信號(hào)來降低噪聲。 因此，SpotClean 將提高基于標(biāo)記基因的推斷的準(zhǔn)確性，包括組織類型注釋、細(xì)胞類型分解、與單細(xì)胞 RNA-seq 數(shù)據(jù)的集成以及相關(guān)的下游分析。 SpotClean 還改進(jìn)了空間可變和差異表達(dá) (DE) 基因的識(shí)別。 注意到，在某些情況下，由于spot組內(nèi)的變異性增加，與已知 DE 基因相關(guān)的 p 值可能會(huì)在 SpotClean 后略微增加（即真正表達(dá)的區(qū)域變得更高表達(dá)；并且從未表達(dá)的區(qū)域中去除信號(hào)）。

鑒于cluster主要由相對(duì)較少的高表達(dá)基因決定，SpotClean 不會(huì)在大多數(shù)數(shù)據(jù)集中顯著改變cluster。 雖然cluster的定義可能會(huì)稍微好一些，但在大多數(shù)情況下，在應(yīng)用 SpotClean 后，看不到cluster數(shù)量和/或cluster之間關(guān)系的差異。

step-to-step analyse

加載10X空間數(shù)據(jù)

# Not run

raw_mat <- Read10xRaw(count_dir = "/path/to/raw_feature_bc_matrix/")
slide_info <- Read10xSlide(tissue_csv_file="/path/to/tissue_positions_list.csv", 
                           tissue_img_file="/path/to/tissue_lowres_image.png",
                           scale_factor_file="/path/to/scalefactors_json.json")

# Compare with bundled example data
data(mbrain_raw)
data(mbrain_slide_info)
slide_info$slide$total_counts <- colSums(
    raw_mat[rownames(mbrain_raw),mbrain_slide_info$slide$barcode]
)

identical(raw_mat[rownames(mbrain_raw),], mbrain_raw)
identical(slide_info$slide, mbrain_slide_info$slide)

Create the slide object

data(mbrain_raw)
data(mbrain_slide_info)

slide_obj <- CreateSlide(mbrain_raw, mbrain_slide_info)
slide_obj

Visualize the slide object

VisualizeSlide(slide_obj)

圖片.png

Function VisualizeLabel() shows the spot labels. You can specify the column name of character labels in the slide metadata, or manually provide a vector of character labels corresponding to each spot. For example, we can plot their tissue/background labels, which has been pre-stored in the input slide information:

VisualizeLabel(slide_obj,"tissue")

圖片.png

Function VisualizeHeatmap() draws a heatmap of values at every spot in the 2-D slide space. Similar to VisualizeLabel(), you can specify the column name of numerical values in the slide metadata, or manually provide a vector of numerical values corresponding to each spot. For example, we can plot the total UMI counts in every spot, which again has been pre-stored in the input slide information:

VisualizeHeatmap(slide_obj,"total_counts")

圖片.png

You can also provide a certain gene name appearing in the raw count matrix in input slide object to VisualizeHeatmap(). For example, the expression of the Mbp gene can be visualized:

VisualizeHeatmap(slide_obj,"Mbp")

圖片.png

VisualizeLabel() and VisualizeHeatmap() both support manual label/value inputs, subsetting spots to plot, title and legend name modification. VisualizeHeatmap() also supports different color palettes (rainbow vs. viridis) and log-scaling options. These visualization functions return ggplot2 objects, which can be further modified by users.

Decontaminate the data

SpotClean() is the main function for performing decontamination. It takes the slide object of raw data as input together with some parameters for controlling optimization and convergence, and returns a slide object with decontaminated gene expressions and other model-related parameters and statistics appending to the slide information. Detailed parameter explanations can be found by running ?SpotClean. Here we set maxit=10 and candidate_radius=20 to save computation time. In practice, SpotClean() by default evaluates a series of candidate radii and automatically chooses the best one. The default maximum number of iterations is 30, which can be extended if convergence has not been reached.

decont_obj <- SpotClean(slide_obj, maxit=10, candidate_radius = 20)

meta數(shù)據(jù)現(xiàn)在包含更多信息，包括來自 SpotClean 模型的參數(shù)估計(jì)和污染水平的測量。

We can visualize the Mbp gene expressions after 10 iterations of decontamination:

VisualizeHeatmap(decont_obj,"Mbp")

圖片.png

Estimate contamination levels in observed data

Our model is able to estimate the proportion of contaminated expression at each tissue spot (i.e. expression at a tissue spot that orginated from a different spot due to spot swapping):

summary(decont_obj@metadata$contamination_rate)

This indicates around 30% of UMIs at a tissue spot in the observed data came from spot swapping contamination, averaging across all tissue spots.

ARC score

還提供了另一種污染水平的主觀估計(jì)，稱為環(huán)境 RNA 污染 (ARC) 分?jǐn)?shù)。 它可以使用函數(shù) ARCScore() 計(jì)算，也是去污輸出的一部分。 直觀地，ARC 分?jǐn)?shù)是觀察到的組織斑點(diǎn)中污染比例的保守下限。 ARC 分?jǐn)?shù)還可以應(yīng)用于基于液滴的單細(xì)胞數(shù)據(jù)，以在用空液滴替換背景點(diǎn)時(shí)估計(jì)環(huán)境 RNA 污染。 可以通過運(yùn)行 ?ARCScore 找到詳細(xì)信息。

ARCScore(slide_obj)
## [1] 0.05160659

This indicates at least 5% expressions in observed data came from spot swapping contamination.

Convert to Seurat object for downstream analyses

ConvertToSeurat() can be used to convert our slide object to a Seurat spatial object. Note that Seurat requires full input of the spatial folder. In the above example, this is the spatial folder.

seurat_obj <- ConvertToSeurat(decont_obj,image_dir = "/path/to/spatial/folder")

生活很好，有你更好