hello，大家好剑按，又到周五了疾就，今天給大家?guī)硪粋€(gè)新的分析內(nèi)容，空間高變基因聯(lián)合算法推斷區(qū)域識(shí)別的方法艺蝴，SpatialDE2猬腰，文章在SpatialDE2: Fast and localized variance component analysis of spatial transcriptomics,看到這個(gè)，大家首先聯(lián)想到的是軟件SpatialDE猜敢，其實(shí)這本身就是SpatialDE的升級(jí)版姑荷，關(guān)于SpatialDE大家可以參考我的文章10X空間轉(zhuǎn)錄組-----空間高變基因檢測(cè)之SpatialDE，好了缩擂，我們來看看SpatialDE2有什么改進(jìn)之處鼠冕。

放一張效果圖，劃分區(qū)域相當(dāng)不錯(cuò)胯盯，優(yōu)于Leiden

圖片.png

Abstract

空間轉(zhuǎn)錄組學(xué)現(xiàn)在是一項(xiàng)成熟的技術(shù)懈费，可以在復(fù)雜組織的組織學(xué)背景下分析基因表達(dá)的變化。典型分析工作流程從識(shí)別具有相似表達(dá)譜的組織區(qū)域開始博脑，然后檢測(cè)高度可變或空間可變的基因憎乙。空間轉(zhuǎn)錄組數(shù)據(jù)集規(guī)模和復(fù)雜性的快速增加要求以一致和集成的方式進(jìn)行這些分析步驟叉趣，當(dāng)前很多方法無法滿足這一要求泞边。為了解決這個(gè)問題，作者開發(fā)了SpatialDE2疗杉，它將組織區(qū)域的映射和空間高變基因檢測(cè)統(tǒng)一為集成軟件框架阵谚，同時(shí)推進(jìn)這兩個(gè)步驟的當(dāng)前算法。該模型在貝葉斯框架中制定，考慮了泊松計(jì)數(shù)噪聲梢什，同時(shí)與以前的方法相比提供了卓越的計(jì)算速度闻牡。使用模擬數(shù)據(jù)驗(yàn)證了 SpatialDE2，并證明了其在實(shí)際中的運(yùn)用價(jià)值绳矩。

Introduction

空間轉(zhuǎn)錄組學(xué)技術(shù)目前獲得了極大的關(guān)注，因?yàn)槠淠軌蛱剿骷?xì)胞局部環(huán)境中的細(xì)胞身份和功能玖翅。 在快速技術(shù)發(fā)展的推動(dòng)下翼馆，空間轉(zhuǎn)錄組現(xiàn)在允許并行研究成百上千個(gè)基因。 基于多重成像的技術(shù)金度，如 SeqFISH+ 或 MERFISH应媚，允許同時(shí)檢測(cè)數(shù)百個(gè)基因參數(shù)，提供亞細(xì)胞分辨率猜极。 基于 RNA 測(cè)序 (RNAseq) 的方法使用空間條形碼引物中姜，并允許對(duì)少數(shù)細(xì)胞進(jìn)行分辨率并接近單細(xì)胞分辨率。 這些方法包括空間轉(zhuǎn)錄組學(xué)跟伏、HDST 和 Slide-seq丢胚，分辨率為少數(shù)細(xì)胞，接近單細(xì)胞分辨率受扳。 特別是携龟，空間轉(zhuǎn)錄組學(xué)的改進(jìn)版本 10x Visium 由于其商業(yè)可用性和易用性而在社區(qū)中得到廣泛采用（大家用到的空間轉(zhuǎn)錄組大部分都是10X公司的）。

雖然最初的研究集中在相對(duì)均勻的小組織切片中的概念驗(yàn)證應(yīng)用勘高，但這些技術(shù)越來越多地應(yīng)用于具有不同結(jié)構(gòu)和區(qū)域的復(fù)雜組織或器官峡蟋，包括人腦、子宮內(nèi)膜和腫瘤微環(huán)境华望。這些數(shù)據(jù)允許解決一系列不同的問題蕊蝗，空間組學(xué)分析工作流程的典型起點(diǎn)是識(shí)別空間可變基因。此步驟產(chǎn)生與下游分析最相關(guān)的基因赖舟，例如組織中局部生態(tài)位的定義支持細(xì)胞分化和功能（這也是我特別強(qiáng)調(diào)大家重視空間高變基因的原因）蓬戚，但關(guān)于空間變異性的知識(shí)本身通常可以提供生物學(xué)見解建蹄，例如變成癌癥碌更。然而，現(xiàn)有的空間可變基因檢測(cè)方法通扯瓷鳎考慮空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)集中的整個(gè)視野痛单。隨著空間轉(zhuǎn)錄組學(xué)的視野不斷擴(kuò)大，并且這些方法被應(yīng)用于由不同細(xì)胞類型組成的區(qū)域組成的日益復(fù)雜的組織劲腿，單純地應(yīng)用空間可變基因檢測(cè)不再能產(chǎn)生相關(guān)的見解甘凭，因?yàn)橐炎R(shí)別的空間可變基因組主要是包含并非固有空間可變的細(xì)胞類型標(biāo)記弃鸦。因此调限，空間方差分析需要結(jié)合合適的計(jì)算方法來識(shí)別組織區(qū)域（說的很有道理）匠楚。

空間可變基因檢測(cè)和組織區(qū)域識(shí)別都在該領(lǐng)域受到關(guān)注。例如僻弹，SpatialDE1是最早檢測(cè)空間可變基因的計(jì)算解決方案之一，最近已通過Leiden 模型對(duì)其進(jìn)行了改進(jìn)。兩種模型都基于非參數(shù)高斯過程 (GP) 回歸房午，SPARK 還提供基于計(jì)數(shù)的可能性和更強(qiáng)大的統(tǒng)計(jì)測(cè)試。同樣基于 GP 回歸的 SVCA 和基于距離加權(quán)非參數(shù)回歸的 scHOT 擴(kuò)展了空間變量基因檢測(cè)的原理丹允，通過考慮相鄰細(xì)胞或voxels之間的相互作用郭厌，提供精度更高的空間基因表達(dá)變異分解。然而雕蔽，所有這些方法的共同點(diǎn)是它們不能擴(kuò)展到大型數(shù)據(jù)集折柠，部分原因是受 CPU 限制的實(shí)現(xiàn)無法利用現(xiàn)代高度并行化的 GPU 架構(gòu)，部分原因是算法效率低下批狐。同樣扇售，存在一系列聚類方法，但它們不太適合識(shí)別空間基因表達(dá)數(shù)據(jù)中的組織區(qū)域嚣艇。 Scanpy 和 Seurat 是兩個(gè)廣泛使用的 scRNA-seq 分析框架承冰，推薦使用 Leiden 聚類從空間組學(xué)中識(shí)別組織區(qū)域，這是一種不知道空間關(guān)系的算法髓废。 Giotto 為此任務(wù)提出了一種基于隱馬爾可夫隨機(jī)場(chǎng)的方法巷懈，從而強(qiáng)加了空間平滑性約束。然而慌洪，該模型要求用戶預(yù)先指定組織區(qū)域的數(shù)量顶燕，并且它采用了對(duì)計(jì)數(shù)數(shù)據(jù)來說不是最佳的高斯似然模型。最后冈爹，注意到缺乏能夠結(jié)合空間變量基因選擇和組織區(qū)域識(shí)別的集成軟件和工作流程涌攻。

基于這方面的確實(shí)，作者開發(fā)了SpatialDE2频伤，一個(gè)用于建目一眩空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)的靈活框架。 SpatialDE2 實(shí)現(xiàn)了兩個(gè)主要模塊憋肖，這兩個(gè)模塊共同提供了用于分析空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)的end-to-end工作流程：組織區(qū)域分割模塊和用于檢測(cè)空間可變基因的模塊因痛。與以前的方法相比，組織區(qū)域分割速度快岸更，并提供了改進(jìn)的可用性鸵膏。特別地，該模塊能夠在采用適當(dāng)?shù)幕谟?jì)數(shù)的可能性的同時(shí)自動(dòng)確定組織區(qū)域的數(shù)量怎炊。用于檢測(cè)空間可變基因的模塊通過提供技術(shù)創(chuàng)新和計(jì)算加速擴(kuò)展了先前的方法谭企，例如 SpatialDE 和 SVCA廓译。文章中使用模擬數(shù)據(jù)和對(duì)兩個(gè)真實(shí)世界數(shù)據(jù)集的應(yīng)用來驗(yàn)證 SpatialDE2 的分割和空間變量基因檢測(cè)模塊，其中包含小鼠大腦和人類子宮內(nèi)膜的 10x Visium 數(shù)據(jù)债查。在這些應(yīng)用中非区，軟件展示了與以前的方法相比，兩個(gè)模塊的速度和魯棒性都有所提高盹廷。 SpatialDE2 為評(píng)估組織區(qū)域內(nèi)的空間表達(dá)異質(zhì)性提供了一個(gè)集成的解決方案征绸，從而為處理復(fù)雜組織和大樣本提供了原則性的策略。

SpatialDE2算法原理

SpatialDE2 通過實(shí)現(xiàn)兩個(gè)無縫集成的分析模塊：組織區(qū)域分割模塊和空間可變基因檢測(cè)模塊俄占，實(shí)現(xiàn)了用于表征子區(qū)域空間異質(zhì)性的end-to-end工作流程（下圖）歹垫。

圖片.png

• 注：Top: SpatialDE2 accepts input spatial transcriptome profiles from a tissue sample, either in the form of a raw gene count matrix, or using cell type counts as provided from a computational deconvolution step (e.g. cell2location). Bottom: input data processing workflow.

這兩個(gè)模塊都直接對(duì)原始 mRNA 計(jì)數(shù)進(jìn)行建模，這些計(jì)數(shù)是從多路分解的空間轉(zhuǎn)錄組學(xué)工作流程中獲得的颠放，或作為輸入的成像技術(shù)。 可選地吭敢，SpatialDE2 還可以對(duì)從附加解卷積步驟獲得的細(xì)胞計(jì)數(shù)估計(jì)值進(jìn)行操作（這里用到的是cell2location）碰凶。

簡(jiǎn)而言之，空間組織區(qū)域分割模塊基于貝葉斯隱馬爾可夫隨機(jī)場(chǎng)鹿驼，將組織分割成不同的組織學(xué)區(qū)域欲低，同時(shí)明確考慮相鄰位置之間的空間平滑度（下圖）

圖片.png

• 注：SpatialDE2 segments the input tissue into connected and transcriptionally similar regions. Top: Schematic output of the segmentation with colours denoting identified tissue regions. Bottom: implementation of the spatial segmentation based on a Poisson Hidden Markov Random field to encode the assumption of spatial smoothness. Nodes correspond to locations with colour denoting the region label. The number of regions is determined by the model.

與當(dāng)前最廣泛使用的基于Leiden聚類的算法相比，SpatialDE2方法的主要優(yōu)勢(shì)是雙重的畜晰。 首先砾莱，與以前的方法不同，SpatialDE2 在分割步驟中明確說明了空間信息凄鼻。 其次腊瑟，SpatialDE2 為使用原則性貝葉斯方法實(shí)現(xiàn)的組織分割提供了一個(gè)連貫?zāi)Ｐ停撃Ｐ托枰粋€(gè)用戶定義的參數(shù)块蚌，該參數(shù)對(duì)空間平滑度的先驗(yàn)假設(shè)進(jìn)行編碼闰非。 相比之下， Leiden 聚類方法的結(jié)果會(huì)受到手動(dòng)定義處理步驟相互作用的非直觀影響峭范，每個(gè)處理步驟都由多個(gè)參數(shù)決定财松。 Leiden聚類工作流包括原始數(shù)據(jù)歸一化、降維纱控、k-最近鄰搜索和 Leiden 聚類算法本身的應(yīng)用辆毡。

The spatially variable gene detection module models variance components of individual genes within identified regions using an appropriate count-based likelihood.。 雖然 SpatialDE2 建立在高斯過程回歸的基礎(chǔ)上甜害，也被 SpatialDE舶掖、SVCA 和 SPARK 等方法使用，但該模型通過提供新穎的特征唾那、改進(jìn)的可擴(kuò)展性访锻、支持多個(gè)方差分量和 GPU 計(jì)算來概括這些方法褪尝。 核心部分，該模型將表達(dá)變異分解為不同的結(jié)構(gòu)化成分和一個(gè)模擬隨機(jī)變異性的噪聲項(xiàng)(下圖期犬，我還是特別推薦有能力的童鞋多多研究數(shù)學(xué)算法河哑，不然最根本的原理完全不明白)。

圖片.png

• 注：SpatialDE2 models spatial variance components of individual genes within tissue regions. Top: Schematic for the identification of spatial variance components of individual genes in specific tissue regions. Bottom: SpatialDE2 models expression variation within tissue regions by partitioning gene expression variation into one or multiple functional components (U1,..,Un). Each component is characterized by a covariance matrix that is parametrized by spatial or non-spatial covariates. The special case of spatially variable gene selection corresponds to a functional component parametrized by distance between locations.

根據(jù)一個(gè)或多個(gè)協(xié)方差矩陣的設(shè)計(jì)龟虎，可以將不同的測(cè)試形式化璃谨，以量化空間可變基因，或識(shí)別受細(xì)胞-細(xì)胞相互作用調(diào)節(jié)的基因鲤妥。

最后佳吞，注意到 SpatialDE2 提供了下游分析工具來幫助解釋。 這包括用于識(shí)別空間共變的基因棉安。 AEH 將特定基因的表達(dá)建模為來自定義數(shù)量的平滑空間模式之一底扳，并嘗試使用貝葉斯框架估計(jì)模式和基因?qū)δＪ降姆峙洹?/h5>

圖片.png

• 注：(D) Run time of SpatialDE2’s tissue region segmentation module and a Leiden clustering workflow for semi-synthetic dataset of increasing size and when using alternative compute environments. Clustering/segmentation was based on 2,000 genes. Leiden denotes the scanpy workflow. Only SpatialDE2 supports GPU computations. (E) Run time of SpatialDE2’s spatially variable gene selection module versus alternative methods. Considered is a dataset consisting of 200 genes for increasing numbers of locations and alternative computing environments. Only SpatialDE2 supports GPU computations.

Model validation using simulated data

首先，在無空間變量基因表達(dá)的零假設(shè)下使用模擬數(shù)據(jù)來確認(rèn)空間變量基因檢測(cè)模塊的統(tǒng)計(jì)校準(zhǔn)贡耽。 簡(jiǎn)而言之衷模，與SpatialDE1和SPARK 類似，SpatialDE2 increases the sensitivity of its spatially variable gene detection by testing multiple kernel matrices for each gene. SpatialDE2 implements two strategies to estimate statistical significance: Each kernel matrix is tested separately and the p-values are combined using the Cauchy combination蒲赂，or all kernel matrices are tested simultaneously using an omnibus test阱冶。 后一種選擇更快，因?yàn)橹贿M(jìn)行了一次測(cè)試滥嘴。 分析確認(rèn)這兩種策略都產(chǎn)生了校準(zhǔn)結(jié)果木蹬。

接下來，模擬了真實(shí)空間可變基因若皱，調(diào)整了來自 10x Visium 小鼠大腦數(shù)據(jù)集的經(jīng)驗(yàn)參數(shù)镊叁。 我們?cè)u(píng)估了 SpatialDE2、SpatialDE和SPARK 的靈敏度（統(tǒng)計(jì)性能）以檢測(cè)真正的空間可變基因走触。SpatialDE的靈敏度最低意系，而 SpatialDE2 和SPARK 的結(jié)果相當(dāng)。 默認(rèn)情況下饺汹，SPARK 通過執(zhí)行內(nèi)核矩陣的特征分解并將負(fù)特征值設(shè)置為零來強(qiáng)制內(nèi)核的正定性蛔添。 由于其默認(rèn)的周期內(nèi)核不是正定的，這會(huì)導(dǎo)致內(nèi)核定義不明確且無法解釋兜辞。 因此迎瞧，還在基準(zhǔn)測(cè)試中包含了沒有特征值裁剪的SPARK ，分析得到一致的結(jié)果逸吵。

Tissue region segmentation recovers known histological features along a continuous resolution gradient

接下來凶硅，測(cè)試來自小鼠大腦的 10x Visium 數(shù)據(jù)，以評(píng)估 SpatialDE2 和組織分割的替代方法扫皱。小鼠大腦非常適合此類基準(zhǔn)測(cè)試目的足绅，因?yàn)樗哂忻鞔_定義和注釋良好的不同區(qū)域捷绑。 使用空間平滑度的默認(rèn)設(shè)置，SpatialDE2 正確解析了主要的大腦區(qū)域（下圖）氢妈，特別是將海馬錐體/顆粒細(xì)胞與周圍的海馬結(jié)構(gòu)（區(qū)域 0）分開并識(shí)別了側(cè)腦室（區(qū)域 3）粹污。 為了進(jìn)行比較，還考慮了應(yīng)用于相同輸入數(shù)據(jù)的Leiden聚類工作流程（方法）首量。 從Leiden聚類獲得的分割解決方案未能解決一些預(yù)期的大腦區(qū)域壮吩，特別是海馬錐體/顆粒細(xì)胞和側(cè)腦室都沒有被識(shí)別。

圖片.png

• 注：(A) Clustering of mouse brain Visium data with SpatialDE2 using default parameters (left), SpatialDE2 with adjusted parameters and clustering results obtained using the Leiden workflow (right).Identified tissue regions are annotated in colour. To avoid clutter, only selected regions that are referred to in the main text are labeled in the legend (out of 18 in total) . (B) Corresponding reference annotation of the mouse brain regions obtained from the Allen Brain Atlas

為了提高空間分辨率加缘，調(diào)整了 SpatialDE2 的空間平滑度參數(shù)和Leiden聚類的分辨率參數(shù)鸭叙。 值得注意的是，SpatialDE2 在廣泛的平滑參數(shù)設(shè)置中產(chǎn)生了一致的分割拣宏，在低精度和高精度解決方案之間提供了有意義的interpolation沈贝。 正如預(yù)期的那樣，具有較低平滑度參數(shù)的分割解決了等皮質(zhì)層以及胼胝體（上圖A）勋乾，并進(jìn)一步將海馬錐體和顆粒細(xì)胞層分成不同的區(qū)域缀程。 另一方面，當(dāng)改變Leiden工作流程的分辨率參數(shù)時(shí)市俊，結(jié)果變化很大并且缺乏直觀的過渡。 此外滤奈，不管分辨率參數(shù)的設(shè)置如何摆昧，Leiden算法既不能分辨胼胝體，也不能分辨海馬錐體/顆粒細(xì)胞蜒程。

對(duì) SpatialDE2 在轉(zhuǎn)錄組距離背景下識(shí)別的區(qū)域的分析表明绅你，該模型獲得的解決方案在空間平滑度和轉(zhuǎn)錄相似性之間取得了適當(dāng)?shù)钠胶狻?特別是，如果這些區(qū)域在轉(zhuǎn)錄上非常相似昭躺，則模型將相同的聚類標(biāo)簽分配給空間不同的區(qū)域忌锯。 這在海馬結(jié)構(gòu)和等皮質(zhì)的第 1 層的情況下尤為明顯，它們被 SpatialDE2 分配到同一cluster领炫。

還通過在各個(gè)位置引導(dǎo) RNA-seq 讀數(shù)來評(píng)估 SpatialDE2 和Leiden聚類的穩(wěn)健性（下圖C）偶垮，并且評(píng)估了兩種方法對(duì)測(cè)序深度變化的敏感性（下圖D）。 總的來說帝洪，這些實(shí)驗(yàn)表明 SpatialDE2 至少與Leiden 工作流程一樣強(qiáng)大似舵。

圖片.png

• 注：(C) Assessment of the robustness of tissue region segmentations obtained using SpatialDE2 or Leiden clustering. Shown is the concordance across 100 bootstrap experiments, sampling sequencing reads with replacement from the full dataset (N=4.5·10⁷ reads total). Bar height denotes the mean split/join distance compared to the results obtained on the full dataset; error bars denote plus or minus one standard deviation. (D) Assessment of the robustness of the SpatialDE2 tissue region segmentation to downsampling of sequencing reads (out of N=45x10^6 reads total). The total number of sequencing reads in the mouse brain dataset was downsampled to the indicated amount. Shown are the segmentation results.

最后，使用 SpatialDE2 空間可變基因檢測(cè)模塊來測(cè)試已識(shí)別組織區(qū)域內(nèi)的空間可變基因葱峡。 這鑒定了多達(dá) 867 個(gè)空間可變基因（FDR<0.1%；下圖E），這些基因豐富了與神經(jīng)元和神經(jīng)組織相關(guān)的 GO terms琼稻，例如“神經(jīng)系統(tǒng)發(fā)育”或“突觸”。 幾個(gè)單獨(dú)的空間可變基因在特定組織區(qū)域中具有已知的功能提鸟，例如突觸結(jié)合蛋白、突觸蛋白仅淑、促甲狀腺激素釋放激素和丘腦中的血管緊張素原（下圖F）称勋。 總的來說，這些結(jié)果表明漓糙，即使在看似同質(zhì)的組織區(qū)域內(nèi)铣缠，在個(gè)體基因和通路水平上也存在表達(dá)變異的空間模式，這可能有助于在更精細(xì)的空間尺度上理解組織功能昆禽。

圖片.png

• 注：(E) Results from the identification of spatially variable genes within the tissue regions identified from the segmentation step as in (A). Bar heights denote the number of spatially variable genes within each region (FDR<0.1%; out of 12,682 genes assessed). (F) Selected examples of spatially variable genes identified in the thalamus. Shown is normalized relative expression of the indicated genes.

SpatialDE2 enables fine-grained analyses in complex human tissue

接下來蝗蛙，為了在更復(fù)雜的數(shù)據(jù)中測(cè)試 SpatialDE2，我們將該模型應(yīng)用于來自人類子宮內(nèi)膜的 Visium 載玻片醉鳖，這是一種高度分隔的組織捡硅，其中細(xì)胞身份取決于剛剛開始全面表征的空間環(huán)境。 最初盗棵，再次應(yīng)用 SpatialDE2 和Leiden聚類來分割組織區(qū)域（下圖）

圖片.png

• 注：Comparison of tissue segmentation obtained using SpatialDE2 (left) and Leiden clustering (right). SpatialDE2 was used with default settings, and the Leiden resolution parameter was selected to yield a biologically meaningful clustering for results using alternative parameters). Myometrium LQ denotes a region of the myometrium with low data quality壮韭。

為了注釋使用 SpatialDE2 識(shí)別的區(qū)域，結(jié)合了來自組織組織學(xué)的線索和從參考細(xì)胞類型的計(jì)算分配中獲得的細(xì)胞類型注釋（使用 cell2location）纹因；這確定了高精度組織亞結(jié)構(gòu)的生理學(xué)上有意義的注釋喷屋。例如，SpatialDE2 識(shí)別了區(qū)域?qū)?yīng)于腺體及其周圍區(qū)域瞭恰。在這兩個(gè)區(qū)域中屯曹，腺上皮細(xì)胞和成纖維細(xì)胞都存在，但比例不同惊畏。這些觀察到的細(xì)胞類型組成的逐漸變化是預(yù)期的恶耽，反映了 Visium 平臺(tái)的低分辨率（10 到 50每個(gè)位置的細(xì)胞。手動(dòng)檢查腺細(xì)胞的典型標(biāo)志物颜启，包括 EPCAM（一種上皮標(biāo)志物）和 PAEP（一種腺體細(xì)胞的典型分泌標(biāo)志物）偷俭，證實(shí)了這些細(xì)胞類型注釋。同樣缰盏，肌肉細(xì)胞標(biāo)志物 ACTA2 和 MYLK 主要是在子宮肌層中表達(dá)涌萤，從而提供額外的信心 SpatiaIDE2 鑒定了具有生理意義的組織區(qū)域。

與 SpatialDE2 相比口猜，Leiden聚類工作流程產(chǎn)生的結(jié)果與生物學(xué)預(yù)期不一致形葬。 例如，識(shí)別出的子宮肌層區(qū)域遠(yuǎn)不如 SpatialDE2 的平滑暮的，并且功能不同的子宮區(qū)域被分配到同一簇：子宮內(nèi)膜基底層的一部分與子宮肌層的一部分聚集在一起笙以，基底層的另一部分聚集在一起 部分功能層，腺上皮與管腔上皮聚集在一起(上圖),至關(guān)重要的是冻辩，雖然調(diào)整Leiden分辨率參數(shù)導(dǎo)致某些區(qū)域的部分改善 - 提高分辨率會(huì)產(chǎn)生更細(xì)粒度的基底層和功能性子宮內(nèi)膜層聚集猖腕，并分離管腔和腺上皮 - 這種調(diào)整同時(shí)導(dǎo)致子宮肌層過度聚集.

接下來拆祈，著手研究單個(gè)組織區(qū)域內(nèi)的空間基因表達(dá)。 SpatialDE2 測(cè)試確定了每個(gè)區(qū)域 53 到 409 個(gè)空間可變基因（FDR<0.001倘感，下圖）

圖片.png

• 注：Number of spatially variable genes within each region (FDR<0.01%).

GO 富集分析確定了幾個(gè)與金屬離子相關(guān)的 GO terms放坏，這些terms用于上皮中的空間可變基因

圖片.png

• 注：GO enrichment analysis of spatially variable genes in the named regions of C). GO terms containing metal-binding genes are highlighted. To avoid clutter, only the GO IDs are shown for the other terms.

這些terms中包含的八個(gè)空間可變基因中有七個(gè)是金屬硫蛋白，它們幾乎只在一小群管腔上皮spot中表達(dá)

圖片.png

• 注：Normalized relative expression of metallothioneins in the functional layer.

管腔上皮富含纖毛細(xì)胞老玛，其特性取決于卵巢激素雌激素和孕激素淤年。 金屬硫蛋白是多種組織中管腔上皮的一個(gè)特征，先前已被證明受孕酮調(diào)節(jié)蜡豹。

然后麸粮，嘗試使用 SpatialDE2 的 AEH 模塊將基因分組為表達(dá)模式。 這種方法是對(duì)將組織分割成離散組織區(qū)域的補(bǔ)充镜廉，特別是允許識(shí)別平滑過渡弄诲。 將 AEH 應(yīng)用于除被歸類為子宮肌層的那些點(diǎn)之外的所有點(diǎn)。 排除子宮肌層有兩個(gè)原因：它是一個(gè)相對(duì)均質(zhì)的組織娇唯，低質(zhì)量區(qū)域可能會(huì)干擾分析齐遵。 這確定了子宮內(nèi)膜中的五種空間模式。 其中四種模式似乎概括了功能層中腺體的分布塔插，但一種模式（模式 1）明顯不同梗摇，并且包含在功能層最內(nèi)腔部分高表達(dá)的基因，就在上皮下方

圖片.png

• 注：Normalized relative expression of the progesterone receptor PGRMC1 and other genes belonging to the same spatial pattern (top right)

這些基因中有孕酮受體 PGRMC1想许，這與子宮內(nèi)膜的功能層響應(yīng)激素而分化（蛻膜化）的事實(shí)一致伶授。

SpatialDE2 設(shè)計(jì)將空間組學(xué)概況作為輸入，從而使模型適用于一系列不同的技術(shù)伸刃，這些技術(shù)產(chǎn)生的表達(dá)估計(jì)可以通過基于計(jì)數(shù)的可能性來解釋。這既包括基于測(cè)序的分析逢倍，也包括成像技術(shù)捧颅。然而，特別是對(duì)于 Visium 或其他無法實(shí)現(xiàn)真正單細(xì)胞分辨率的技術(shù)较雕，相對(duì)粗略的分割可能會(huì)導(dǎo)致將特定組織分配給集群的限制碉哑，如具有“腺體（邊界）”集群的子宮內(nèi)膜數(shù)據(jù)集中所見.受此限制的啟發(fā)，因此著手探索如何將 SpatialDE2 與從利用參考 scRNA-seq 數(shù)據(jù)集估計(jì)細(xì)胞類型豐度的計(jì)算反卷積工作流程中獲得的細(xì)胞類型計(jì)數(shù)估計(jì)相結(jié)合亮蒋。具體來說扣典，將 SpatialDE2 應(yīng)用于通過 cell2location 在子宮內(nèi)膜數(shù)據(jù)集上獲得的細(xì)胞類型豐度估計(jì)。與將 SpatialDE2 直接應(yīng)用于基因計(jì)數(shù)相比慎玖，這導(dǎo)致了更高精度的分割贮尖。

圖片.png

• 注：Comparison of tissue segmentation obtained using SpatialDE2 (left) and Leiden clustering (right) based on cell2location cell type abundance estimates. Number of clusters was determined automatically for SpatialDE2, whereas the Leiden resolution parameter was tuned by hand to yield a biologically meaningful clustering

例如，SpatialDE2 現(xiàn)在將基底層分成三個(gè)不同的區(qū)域（圖中編號(hào)為 1,5 和 8）趁怔。 屬于這些區(qū)域的spot也以基因表達(dá)空間中的不同cluster為特征

圖片.png

• 注：UMAP embedding of individual locations based on their expression profiles. Color denotes the cluster labels assigned to one of three basal layer clusters by SpatialDE2.

表明這些區(qū)域之間存在相關(guān)的表型差異湿硝。 一致地薪前，觀察到這些區(qū)域之間細(xì)胞類型組成的顯著差異。 例如关斜，區(qū)域 8 富含上皮腺細(xì)胞示括，而區(qū)域 1 具有更多的血管周圍 STEAP4 細(xì)胞

圖片.png

• 注：Cell type composition of the three basal layer clusters as in (B). Shown are estimates of cell type fractions from cell2location across individual locations assigned to either cluster 1, 5 or 8

一般來說，大多數(shù)區(qū)域顯示出獨(dú)特的細(xì)胞類型組成痢畜，允許對(duì)組織組織學(xué)進(jìn)行細(xì)粒度解剖垛膝。 最后注意到，如果應(yīng)用于 cell2location 輸出并且錯(cuò)誤地將上皮和腺體以及部分功能層和基底層組合在一起丁稀，Leiden 聚類工作流程再次表現(xiàn)不佳吼拥，同時(shí)嚴(yán)重過度聚類子宮肌層。

Discussion

在這里二驰，展示了 SpatialDE2扔罪，一個(gè)用于分析空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)的多功能工具箱。 該方法的核心是兩個(gè)分析模塊——區(qū)域分割和空間可變基因識(shí)別桶雀。 與現(xiàn)有方法相比矿酵，這兩個(gè)模塊都具有重要的優(yōu)勢(shì)。 此外矗积，這些模塊是兼容的全肮，可以在一個(gè)連貫的軟件工作流程中有效組合。

區(qū)域分割模塊獨(dú)特地考慮了組織的空間結(jié)構(gòu)棘捣，直接對(duì) RNA 計(jì)數(shù)數(shù)據(jù)進(jìn)行操作辜腺，并且可以通過單個(gè)直觀參數(shù)進(jìn)行調(diào)整。 同時(shí)乍恐，該模塊與忽略空間坐標(biāo)的傳統(tǒng) Leiden 聚類一樣有效评疗。 此外，SpatialDE2 提供 GPU 加速計(jì)算選項(xiàng)茵烈。 類似地百匆，空間可變基因方向模塊提供比以前的方法（如 SpatialDE 或 SPARK）快幾個(gè)數(shù)量級(jí)的計(jì)算，同時(shí)提供可比較或更好的統(tǒng)計(jì)能力呜投。 SpatialDE2 還保留了這些基于內(nèi)核的方法的顯著特征加匈，例如通過設(shè)計(jì)相應(yīng)的內(nèi)核來測(cè)試特定類型的表達(dá)模式的能力

這兩個(gè)模塊同時(shí)進(jìn)行，可以有效實(shí)施工作流程仑荐，該工作流程包括首先分割組織切片雕拼，然后檢測(cè)每個(gè)區(qū)域內(nèi)的空間可變基因。 在小鼠大腦上展示了這個(gè)工作流程粘招，這是一種適合評(píng)估性能的特征明確的組織啥寇。 在這個(gè)數(shù)據(jù)集上，組織區(qū)域分割的表現(xiàn)與Leiden 聚類工作流程相當(dāng)，能夠恢復(fù)已知的生物學(xué)示姿，例如主要的大腦區(qū)域甜橱。 還將 SpatialDE2 工作流程應(yīng)用于人類子宮內(nèi)膜數(shù)據(jù)集，這是一種表征較少的組織栈戳。 在這里岂傲，組織區(qū)域分割證明優(yōu)于Leiden 聚類。 除了概括已知的生物學(xué)之外子檀，SpatialDE2 還對(duì)這種復(fù)雜組織產(chǎn)生了新的見解镊掖，例如金屬硫蛋白的表達(dá)模式和與孕酮受體在空間上共表達(dá)的幾個(gè)基因。

最后褂痰，展示了如何將 SpatialDE2 應(yīng)用于基于參考表達(dá)譜的空間轉(zhuǎn)錄組學(xué)計(jì)算反卷積的下游亩进，從而產(chǎn)生具有更細(xì)粒度結(jié)果的組織區(qū)域分割，從而實(shí)現(xiàn)組織組織學(xué)的詳細(xì)解剖缩歪。

最后看一下示例代碼

import SpatialDE

import numpy as np
import scipy
import pandas as pd
import scanpy as sc
import anndata as ad

from tqdm.auto import trange, tqdm

import matplotlib.pyplot as plt
from matplotlib.pyplot import rc_context
from IPython.display import set_matplotlib_formats
%matplotlib inline
set_matplotlib_formats('png')
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 200
mpl.rcParams['svg.fonttype'] = 'none'
mpl.rc('font', **{'family':'sans-serif','sans-serif':['Helvetica']})
from plotnine import *
import mizani

theme_set(theme_bw())
theme_update(legend_key=element_blank(), text=element_text(family="Helvetica"))
data = sc.read_visium("mouse_brain_visium_wo_cloupe_data/rawdata/ST8059048")
data.var_names_make_unique()
data.var["mt"] = data.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(data, qc_vars=["mt"], inplace=True)
sc.pp.filter_cells(data, min_counts=4000)
data = data[data.obs["pct_counts_mt"] < 20]
sc.pp.filter_genes(data, min_cells=100)
svg_full, _ = SpatialDE.test(data, omnibus=True)
svg_full["total_counts"] = np.asarray(data.X.sum(axis=0)).squeeze()
svg_full.to_pickle("ST8059048_svg_full.pkl")
svg_full = pd.read_pickle("ST8059048_svg_full.pkl")
vargenes = svg_full[svg_full.padj < 0.001].sort_values("total_counts", ascending=False).gene[:2000]
segm_default, _ = SpatialDE.tissue_segmentation(data, genes=vargenes, rng=np.random.default_rng(seed=42))
data.obs["segmentation_labels_default"] = data.obs.segmentation_labels
segm_adj, _ = SpatialDE.tissue_segmentation(data, genes=vargenes, rng=np.random.default_rng(seed=42), params=SpatialDE.TissueSegmentationParameters(smoothness_factor=5))
with rc_context({"figure.figsize": (2, 2)}):
    sc.pl.spatial(data, color="segmentation_labels_default", title="SpatialDE2", return_fig=True).savefig("figures/3a_default.svg")
    sc.pl.spatial(data, color="segmentation_labels", title="SpatialDE2 (adjusted)", return_fig=True).savefig("figures/3a_adjusted.svg")

圖片.png

data_normalized = data.copy()
sc.pp.normalize_total(data_normalized, target_sum=1e4, key_added="scaling_factor")
data_norm_reg = data_normalized.copy()
sc.pp.log1p(data_norm_reg)
variable = sc.pp.highly_variable_genes(data, flavor="seurat_v3", inplace=False, n_top_genes=2000)
data_norm_reg = data_normalized.copy()
sc.pp.log1p(data_norm_reg)
data_norm_reg = data_norm_reg[:, vargenes]
sc.pp.regress_out(data_norm_reg, ['total_counts'])
sc.pp.scale(data_norm_reg, max_value=10)
sc.tl.pca(data_norm_reg, n_comps=100)
sc.pp.neighbors(data_norm_reg, n_neighbors=20, n_pcs=80)
sc.tl.leiden(data_norm_reg, random_state=42)
sc.tl.umap(data_norm_reg)

with rc_context({"figure.figsize": (2, 2)}):
    sc.pl.spatial(data_norm_reg, color="leiden", title="Leiden", return_fig=True).savefig("figures/3a_leiden.svg")
    sc.pl.umap(data_norm_reg, size=3, color=["segmentation_labels_default", "segmentation_labels", "leiden"], title=["SpatialDE2", "SpatialDE2 (adjusted)", "Leiden"], return_fig=True, ncols=3).savefig("figures/s3c.svg")

圖片.png

resolutions = (0.1, 0.5, 1., 1.5, 2., 3., 4.)
smoothnesses = (0.5, 1., 2., 3., 5., 8., 10.)
data_norm_reg.obs["leiden_1.0"] = data_norm_reg.obs.leiden
with rc_context({"figure.figsize": (2, 2)}):
    for res in resolutions:
        sc.tl.leiden(data_norm_reg, random_state=42, resolution=res, key_added=f"leiden_{res}")
        sc.pl.spatial(data_norm_reg, color=f"leiden_{res}", return_fig=True, title=f"Leiden ({res})").savefig(f"figures/s3b_{res}.svg")

圖片.png

testdata = data.copy()
with rc_context({"figure.figsize": (2, 2)}):
    for smoothness in smoothnesses:
        segm_adj, _ = SpatialDE.tissue_segmentation(testdata, genes=vargenes, rng=np.random.default_rng(seed=42), params=SpatialDE.TissueSegmentationParameters(smoothness_factor=smoothness))
        sc.pl.spatial(testdata, color=f"segmentation_labels", return_fig=True, title=f"SpatialDE2 ({smoothness})").savefig(f"figures/s3a_{smoothness}.svg")

圖片.png

圖片.png

nresamples = 100
totals = (20e6, 10e6, 5e6, 1e6)
obs = data.obs.copy()
for i in trange(nresamples):
    bdata = data.copy()
    bdata.X.data = np.random.default_rng(seed=i).multinomial(bdata.X.data.sum(), bdata.X.data.astype(np.uint16) / bdata.X.data.astype(np.uint16).sum())
    bdata.obs.total_counts = bdata.X.sum(axis=1)
    
    segm, _ = SpatialDE.tissue_segmentation(bdata, genes=vargenes, rng=np.random.default_rng(seed=42))
    obs[f"segmentation_labels_bootstrap_{i}"] = bdata.obs.segmentation_labels
    
    sc.pp.normalize_total(bdata, target_sum=1e4, key_added="scaling_factor")
    bdata = bdata[:, vargenes]
    sc.pp.regress_out(bdata, ['total_counts'])
    sc.pp.scale(bdata, max_value=10)
    sc.tl.pca(bdata, n_comps=100)
    sc.pp.neighbors(bdata, n_neighbors=20, n_pcs=80)
    sc.tl.leiden(bdata, random_state=42)
    obs[f"leiden_bootstrap_{i}"] = bdata.obs.leiden
    
    for total in totals:
        cdata = sc.pp.downsample_counts(data, total_counts=total, copy=True, random_state=i)
        cdata.obs.total_counts = cdata.X.sum(axis=1)
        
        for smoothness in smoothnesses:
            segm, _ = SpatialDE.tissue_segmentation(cdata, genes=vargenes, rng=np.random.default_rng(seed=42), params=SpatialDE.TissueSegmentationParameters(smoothness_factor=smoothness))
            obs[f"segmentation_labels_{total:.0e}_{i}_{smoothness}"] = cdata.obs.segmentation_labels

        sc.pp.normalize_total(cdata, target_sum=1e4, key_added="scaling_factor")
        cdata = cdata[:, vargenes]
        sc.pp.regress_out(cdata, ['total_counts'])
        sc.pp.scale(cdata, max_value=10)
        sc.tl.pca(cdata, n_comps=100)
        sc.pp.neighbors(cdata, n_neighbors=20, n_pcs=80)
        for res in resolutions:
            sc.tl.leiden(cdata, random_state=42, resolution=res)
            obs[f"leiden_{total:.0e}_{i}_{res}"] = cdata.obs.leiden
obs.to_pickle("data_obs.pkl")
data.write_h5ad("data.h5ad")
data = ad.read_h5ad("data.h5ad")
obs = pd.read_pickle("data_obs.pkl")
data.obs = obs
def split_join_metric(obs1, obs2):
    cl1 = np.unique(obs1)
    cl2 = np.unique(obs2)
    
    cp1 = {c: 0 for c in cl1}
    cp2 = {c: 0 for c in cl2}
    for ccl1 in cl1:
        cobs1 = obs1 == ccl1
        for ccl2 in cl2:
            intersect = np.sum(cobs1 & (obs2 == ccl2))
            cp1[ccl1] = max(cp1[ccl1], intersect)
            cp2[ccl2] = max(cp2[ccl2], intersect)
    p1 = sum(cp1.values())
    p2 = sum(cp2.values())
    return 2 * len(obs1) - p1 - p2
data.obs["leiden"] = data_norm_reg.obs.leiden
dists = (pd.concat([
            pd.DataFrame.from_records([(i, split_join_metric(data.obs.segmentation_labels_default, data.obs[f"segmentation_labels_bootstrap_{i}"]), "SpatialDE") for i in range(nresamples)], columns=["resample", "dist", "method"]),
            pd.DataFrame.from_records([(i, split_join_metric(data.obs.leiden, data.obs[f"leiden_bootstrap_{i}"]), "Leiden") for i in range(nresamples)], columns=["resample", "dist", "method"])
        ], axis=0).groupby("method", as_index=False).
            agg(y=("dist", np.mean),
                ymin=("dist", lambda x: x.mean() - x.std()),
                ymax=("dist", lambda x: x.mean() + x.std()))
       )
(
    ggplot(
        dists,
        aes(x="method", y="y", ymin="ymin", ymax="ymax")) +
    geom_bar(stat='identity') +
    geom_errorbar() +
    labs(x=None, y="split / join distance") +
    theme(figure_size=(2,2))
).save("figures/3c.svg")
dists

圖片.png

dists_spatialde = pd.DataFrame.from_records([(total, i, sm, split_join_metric(data.obs.segmentation_labels_default, data.obs[f"segmentation_labels_{total:.0e}_{i}_{sm}"])) for total in totals for i in range(nresamples) for sm in smoothnesses], columns=["depth", "resample", "smoothness", "dist"])
mean_dists = dists_spatialde.groupby(["depth", "smoothness"], as_index=False).agg(mean_dist=("dist", np.mean))
spatialde_optimal = mean_dists.loc[mean_dists.groupby("depth")["mean_dist"].idxmin()]

dists_leiden = pd.DataFrame.from_records([(total, i, sm, split_join_metric(data.obs.leiden, data.obs[f"leiden_{total:.0e}_{i}_{sm}"])) for total in totals for i in range(nresamples) for sm in resolutions], columns=["depth", "resample", "res", "dist"])
mean_dists = dists_leiden.groupby(["depth", "res"], as_index=False).agg(mean_dist=("dist", np.mean))
leiden_optimal = mean_dists.loc[mean_dists.groupby("depth")["mean_dist"].idxmin()]
def label(x):
    if x == 0:
        return "0.0"
    exp = np.floor(np.log10(x))
    factor = x / (10 ** exp)
    return f"${factor:.1f}\\! \\times\\! 10^{{{int(exp)}}}$"
dists_spatialde = dists_spatialde.set_index(["depth", "smoothness"]).loc[spatialde_optimal.set_index(["depth", "smoothness"]).index].reset_index()
dists_leiden = dists_leiden.set_index(["depth", "res"]).loc[leiden_optimal.set_index(["depth", "res"]).index].reset_index()
dists_spatialde["method"] = "SpatialDE2"
dists_leiden["method"] = "Leiden"
(ggplot(
    pd.concat((dists_spatialde, dists_leiden), axis=0).groupby(["depth", "method"], as_index=False).agg(y=("dist", np.mean),
                                ymin=("dist", lambda x: x.mean() - x.std()),
                                ymax=("dist", lambda x: x.mean() + x.std())),
    aes("depth", "y", ymin="ymin", ymax="ymax", color="method")) + geom_line() + geom_errorbar(width=5e5) +
    labs(x="sequencing depth", y="split/join distance") +
    scale_x_continuous(labels=lambda x: [label(i) for i in x]) +
    theme(figure_size=(2,2), axis_text_x=element_text(rotation=20), legend_title=element_blank())
).save("figures/s3e.svg")
spatialde_optimal

圖片.png

leiden_optimal

圖片.png

with rc_context({"figure.figsize": (2, 2)}):
    for total, optimal in zip(totals, (3., 3., 3., 0.5)):
        sc.pl.spatial(data, color=f"segmentation_labels_{total:.0e}_0_{optimal}", return_fig=True, title=f"SpatialDE2 ({total:.0e} reads)").savefig(f"figures/3d_{total:.0e}.svg")

圖片.png

with rc_context({"figure.figsize": (2, 2)}):
    for total, optimal in zip(totals, (1., 1., 1., 0.5)):
        sc.pl.spatial(data, color=f"leiden_{total:.0e}_0_{optimal}", return_fig=True, title=f"Leiden ({total:.0e} reads)").savefig(f"figures/s3d_{total:.0e}.svg")

圖片.png

ulabels, lcounts = np.unique(data.obs.segmentation_labels_default, return_counts=True)
svg_default = []
for l in ulabels[lcounts > 10]:
    csvg, _ = SpatialDE.test(data[data.obs.segmentation_labels_default == l, :], omnibus=True)
    csvg["label"] = l
    svg_default.append(csvg)
svg_default = pd.concat(svg_default, axis=0, ignore_index=True)
svg_default.to_pickle("ST8059048_svg_default.pkl")

ulabels, lcounts = np.unique(data.obs.segmentation_labels, return_counts=True)
svg = []
for l in ulabels[lcounts > 10]:
    csvg, _ = SpatialDE.test(data[data.obs.segmentation_labels == l, :], omnibus=True)
    csvg["label"] = l
    svg.append(csvg)
svg = pd.concat(svg, axis=0, ignore_index=True)
svg.to_pickle("ST8059048_svg.pkl")
svg = pd.read_pickle("ST8059048_svg.pkl")
(
    ggplot(svg.assign(region=lambda x: pd.Categorical(x.label))[svg.padj < 0.001], aes("region")) +
        geom_bar() +
        ylab("spatially variable genes") +
        theme(axis_text_x=element_text(rotation=20, ha="right"), figure_size=(4, 2))
).save("figures/3e.svg")
from bioservices import BioMart
import goatools
from io import StringIO
s = BioMart()
s.add_dataset_to_xml("mmusculus_gene_ensembl")
s.add_attribute_to_xml("ensembl_gene_id")
s.add_attribute_to_xml("external_gene_name")
s.add_attribute_to_xml("entrezgene_accession")
s.add_attribute_to_xml("entrezgene_id")

s.add_filter_to_xml("ensembl_gene_id", ",".join(data.var.gene_ids.to_numpy()))
res = s.query(s.get_xml())
res = pd.read_table(StringIO(res), names=["ensembl_gene_id", "gene_name", "entrezgene_accession", "entrezgene_id"], dtype={"entrezgene_id": pd.Int32Dtype()})

def filter_entrez(group):
    group = group.drop("ensembl_gene_id", axis=1)
    if group.shape[0] == 1:
        ret = group
    else:
        matched = group[group.gene_name == group.entrezgene_accession]
        if matched.shape[0] > 0:
            ret = matched.iloc[[0],:] # no clear way to choose if there are multiple matches, so I just choose the first one
        else:
            ret = matched
    return ret.drop("entrezgene_accession", axis=1)
res = res.dropna().reset_index(drop=True).groupby("ensembl_gene_id").apply(filter_entrez)

data.var["entrez_id"] = (
    data.var.
    reset_index().
    rename(columns={"index":"gene_name", "gene_ids": "ensembl_gene_id"}).
    merge(
        res.assign(entrezgene_id=res.entrezgene_id.astype(pd.Int32Dtype())),
        on=["gene_name", "ensembl_gene_id"]).
    set_index("gene_name").
    rename_axis(index=None)["entrezgene_id"]
from goatools.base import download_go_basic_obo, download_ncbi_associations
from goatools.obo_parser import GODag
from goatools.anno.genetogo_reader import Gene2GoReader
obo_fname = download_go_basic_obo()
fin_gene2go = download_ncbi_associations()
obodag = GODag(obo_fname)
objanno = Gene2GoReader(fin_gene2go, taxids=[10090])
ns2assoc = objanno.get_ns2assc()
background_genes = data.var.entrez_id.to_numpy()
from goatools.goea.go_enrichment_ns import GOEnrichmentStudyNS
goeaobj = GOEnrichmentStudyNS(background_genes, ns2assoc, obodag, propagate_counts=False, alpha=0.01, methods=['fdr_by'])
go_enrichment = []
for label, g in svg.groupby("label"):
    enrichment = goeaobj.run_study(data.var.entrez_id[g.gene[g.padj < 0.001]].to_numpy().tolist(), prt=None)
    res = []
    for term in enrichment:
        res.append({
            "ont": term.NS,
            "go": term.GO,
            "desc": term.name,
            "study_count": term.study_count,
            "population_count": term.pop_count,
            "study_ratio": np.divide(*term.ratio_in_study),
            "population_ratio": np.divide(*term.ratio_in_pop),
            "pval": term.p_uncorrected,
            "padj": term.p_fdr_by,
            "label": label,
            "genes": data.var.query("entrez_id in @term.study_items").index.to_numpy()
        })
    go_enrichment.append(pd.DataFrame(res))
go_enrichment = pd.concat(go_enrichment, axis=0, ignore_index=True).assign(label=lambda x: pd.Categorical(x.label))
(
    ggplot(go_enrichment[(go_enrichment.padj < 1e-3) & (go_enrichment.ont == "CC")], aes("label", "desc", color="-np.log10(pval)", size="study_ratio")) +
    geom_point() +
    facet_wrap("ont", ncol=1, scales="free") +
    scale_size_area(breaks=[0.1, 0.2, 0.3, 0.4, 0.5], limits=(0.1, 0.5), oob=mizani.bounds.squish, name="fraction of sample") +
    scale_color_distiller(palette="YlGnBu", direction=1, limits=(0, 25), oob=mizani.bounds.squish, name=r"$-\log_{10}(p)$") +
    theme(figure_size=(2,7), subplots_adjust={'hspace': 0.1})
).save("figures/s3f_1.svg")
(
    ggplot(go_enrichment[(go_enrichment.padj < 1e-3) & (go_enrichment.ont != "CC")], aes("label", "desc", color="-np.log10(pval)", size="study_ratio")) +
    geom_point() +
    facet_wrap("ont", ncol=1, scales="free") +
    scale_size_area(breaks=[0.1, 0.2, 0.3, 0.4, 0.5], limits=(0.1, 0.5), oob=mizani.bounds.squish, name="fraction of sample") +
    scale_color_distiller(palette="YlGnBu", direction=1, limits=(0, 25), oob=mizani.bounds.squish, name=r"$-\log_{10}(p)$") +
    theme(figure_size=(2,7), subplots_adjust={'hspace': 0.2})
).save("figures/s3f_2.svg")
with rc_context({"figure.figsize": (2, 2)}):
    for gene in ("Syt7", "Syngr3", "Trh", "Agt"):
        sc.pl.spatial(data_normalized[data.obs.segmentation_labels == 5, :], color=gene, return_fig=True).savefig(f"figures/3f_{gene}.svg")

圖片.png

with pd.ExcelWriter("figures/table_s_mouse_brain.xlsx") as ew:
    svg_full[["gene", "pval", "padj"]].to_excel(ew, sheet_name="SV genes whole slice", index=False)
    svg[["gene", "pval", "padj", "label"]].rename(columns={"label": "region"}).to_excel(ew, sheet_name="SV genes per region", index=False)
    go_enrichment.rename(columns={"region": "region_name", "label": "region"}).to_excel(ew, sheet_name="SV genes GO enrichment", index=False)
    data.obs[["segmentation_labels_default", "segmentation_labels"]].reset_index().rename(columns={"index": "spatial_barcode", "segmentation_labels": "region", "segmentation_labels_names": "region_name"}).to_excel(ew, sheet_name="segmentation (counts)", index=False)

import SpatialDE

import numpy as np
import scipy
import pandas as pd
import scanpy as sc
import anndata as ad

from tqdm.auto import trange, tqdm

import matplotlib.pyplot as plt
from matplotlib.pyplot import rc_context
from IPython.display import set_matplotlib_formats
%matplotlib inline
set_matplotlib_formats('png')
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 200
mpl.rcParams['svg.fonttype'] = 'none'
mpl.rc('font', **{'family':'sans-serif','sans-serif':['Helvetica']})
from plotnine import *
import mizani

theme_set(theme_bw())
theme_update(legend_key=element_blank(),
             text=element_text(family="Helvetica"),
             axis_text=element_text(color="black"),
             panel_border=element_rect(color="black"),
             strip_background=element_rect(color="black")
            )

from io import StringIO
from itertools import chain
from bioservices import BioMart, UniProt
import goatools
a30_0 = sc.read_visium("152807")
a30_0.var_names_make_unique()
a30_0.var["mt"] = a30_0.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(a30_0, qc_vars=["mt"], inplace=True)
sc.pp.filter_cells(a30_0, min_counts=800)
a30_0 = a30_0[a30_0.obs["pct_counts_mt"] < 20]
sc.pp.filter_genes(a30_0, min_cells=100)
s = BioMart()
s.add_dataset_to_xml("hsapiens_gene_ensembl")
s.add_attribute_to_xml("ensembl_gene_id")
s.add_attribute_to_xml("external_gene_name")
s.add_attribute_to_xml("entrezgene_accession")
s.add_attribute_to_xml("entrezgene_id")

s.add_filter_to_xml("ensembl_gene_id", ",".join(a30_0.var.gene_ids.to_numpy()))
res = s.query(s.get_xml())
res = pd.read_table(StringIO(res), names=["ensembl_gene_id", "gene_name", "entrezgene_accession", "entrezgene_id"], dtype={"entrezgene_id": pd.Int32Dtype()})

def filter_entrez(group):
    group = group.drop("ensembl_gene_id", axis=1)
    if group.shape[0] == 1:
        ret = group.iloc[0,:]
    else:
        matched = group[group.gene_name == group.entrezgene_accession]
        ret = matched.iloc[0,:] # no clear way to choose if there are multiple matches, so I just choose the first one
    return ret.drop("entrezgene_accession")
res = res.dropna().reset_index(drop=True).groupby("ensembl_gene_id").apply(filter_entrez)

a30_0.var["entrez_id"] = (
    a30_0.var.
    reset_index().
    rename(columns={"index":"gene_name", "gene_ids": "ensembl_gene_id"}).
    merge(
        res.assign(entrezgene_id=res.entrezgene_id.astype(pd.Int32Dtype())),
        on=["gene_name", "ensembl_gene_id"]).
    set_index("gene_name").
    rename_axis(index=None)["entrezgene_id"]
svg_full, _ = SpatialDE.test(a30_0, omnibus=True)
svg_full["total_counts"] = np.asarray(a30_0.X.sum(axis=0)).squeeze()
svg_full.to_pickle("152807_svg_full.pkl")
svg_full = pd.read_pickle("152807_svg_full.pkl")


vargenes = svg_full[svg_full.padj < 0.001].sort_values("total_counts", ascending=False).gene[:2000]
segm, _ = SpatialDE.tissue_segmentation(a30_0, genes=vargenes, rng=np.random.default_rng(seed=42))
region_mapping = {0:"basal layer", 1: "myometrium LQ", 2: "myometrium", 3: "epithelium", 4: "functional layer", 7: "glands (border)", 11: "glands"}
a30_0.obs["segmentation_labels_names"] = a30_0.obs.segmentation_labels.cat.rename_categories(region_mapping)

with rc_context({"figure.figsize": (2,2)}):
    sc.pl.spatial(a30_0, color="segmentation_labels_names", title="SpatialDE2", return_fig=True).savefig("figures/4a_spatialde.svg", dpi=1000)

圖片.png

a30_0_normalized = a30_0.copy()
sc.pp.normalize_total(a30_0_normalized, target_sum=1e4, key_added="scaling_factor")
a30_0_norm_reg = a30_0_normalized.copy()
sc.pp.log1p(a30_0_norm_reg)
a30_0_norm_reg = a30_0_norm_reg[:, vargenes]
sc.pp.regress_out(a30_0_norm_reg, ['total_counts'])
sc.pp.scale(a30_0_norm_reg, max_value=10)
sc.tl.pca(a30_0_norm_reg, svd_solver='arpack')
sc.pp.neighbors(a30_0_norm_reg, n_neighbors=10, n_pcs=50)
sc.tl.leiden(a30_0_norm_reg, random_state=42, resolution=0.6)
with rc_context({"figure.figsize": (2,2)}):
    sc.pl.spatial(a30_0_norm_reg, color="leiden", title="Leiden", return_fig=True).savefig("figures/4a_leiden.svg", dpi=1000)

圖片.png

resolutions = (0.1, 0.3, 0.6, 0.8, 1., 1.5, 2.)
smoothnesses = (0.5, 1., 2., 3., 5., 8., 10.)


a30_0_norm_reg.obs["leiden_0.6"] = a30_0_norm_reg.obs.leiden
with rc_context({"figure.figsize": (2, 2)}):
    for res in resolutions:
        sc.tl.leiden(a30_0_norm_reg, random_state=42, resolution=res, key_added=f"leiden_{res}")
        sc.pl.spatial(a30_0_norm_reg, color=f"leiden_{res}", return_fig=True, title=f"Leiden ({res})").savefig(f"figures/s4d_{res}.svg", dpi=1000)

圖片.png

圖片.png

testdata = a30_0.copy()
with rc_context({"figure.figsize": (2, 2)}):
    for smoothness in smoothnesses:
        segm_adj, _ = SpatialDE.tissue_segmentation(testdata, genes=vargenes, rng=np.random.default_rng(seed=42), params=SpatialDE.TissueSegmentationParameters(smoothness_factor=smoothness))
        sc.pl.spatial(testdata, color=f"segmentation_labels", return_fig=True, title=f"SpatialDE2 ({smoothness})").savefig(f"figures/s4c_{smoothness}.svg", dpi=1000)

圖片.png

圖片.png

ulabels, lcounts = np.unique(segm.labels, return_counts=True)
svg_byregion = []
for l in ulabels[lcounts > 100]:
    res, _ = SpatialDE.test(a30_0[segm.labels == l, :], omnibus=True)
    res["label"] = l
    svg_byregion.append(res)
svg_byregion = pd.concat(svg_byregion, axis=0, ignore_index=True)
svg_byregion.to_pickle("15807_svg_byregion.pkl")
svg_byregion = pd.read_pickle("15807_svg_byregion.pkl")
(
    ggplot(svg_byregion.
               assign(region=lambda x: pd.Categorical(x.label).rename_categories(region_mapping))[svg_byregion.padj < 0.001],
           aes("region")) +
        geom_bar() +
        ylab("spatially variable genes") +
        theme(axis_text_x=element_text(rotation=20, ha="right"), figure_size=(3, 2))
).save("figures/4b.svg"
from goatools.base import download_go_basic_obo, download_ncbi_associations
from goatools.obo_parser import GODag
from goatools.anno.genetogo_reader import Gene2GoReader
obo_fname = download_go_basic_obo()
fin_gene2go = download_ncbi_associations()
obodag = GODag(obo_fname)
objanno = Gene2GoReader(fin_gene2go, taxids=[9606])
ns2assoc = objanno.get_ns2assc()
background_genes = a30_0.var.entrez_id.to_numpy()
from goatools.goea.go_enrichment_ns import GOEnrichmentStudyNS
goeaobj = GOEnrichmentStudyNS(background_genes, ns2assoc, obodag, propagate_counts=False, alpha=0.01, methods=['fdr_by'])
go_enrichment = []
for label, g in svg_byregion.groupby("label"):
    enrichment = goeaobj.run_study(a30_0.var.entrez_id[g.gene[g.padj < 0.001]].to_numpy().tolist(), prt=None)
    res = []
    for term in enrichment:
        res.append({
            "ont": term.NS,
            "go": term.GO,
            "desc": term.name,
            "study_count": term.study_count,
            "population_count": term.pop_count,
            "study_ratio": np.divide(*term.ratio_in_study),
            "population_ratio": np.divide(*term.ratio_in_pop),
            "pval": term.p_uncorrected,
            "padj": term.p_fdr_by,
            "label": label,
            "genes": a30_0.var.query("entrez_id in @term.study_items").index.to_numpy()
        })
    go_enrichment.append(pd.DataFrame(res))
go_enrichment = pd.concat(go_enrichment, axis=0, ignore_index=True).assign(label=lambda x: pd.Categorical(x.label))
go_enrichment["region"] = go_enrichment.label.cat.rename_categories(region_mapping)
(
    ggplot(go_enrichment[(go_enrichment.padj < 1e-2) & (go_enrichment.ont == "BP")].assign(
            label=lambda x: np.where(x.go.isin(("GO:0010273", "GO:0006882", "GO:0071294", "GO:0071280", "GO:0071276")), x.desc, x.go)),
        aes("region", "label", color="-np.log10(pval)", size="study_ratio")) +
    geom_point() +
    facet_wrap("ont", ncol=1, scales="free") +
    labs(size="fraction of sample", color=r"$-\log_{10}(p)$", y=None, x=None) + 
    theme(figure_size=(2.5,5), subplots_adjust={'hspace': 0.2}, axis_text_x=element_text(rotation=30, ha="right"))
).save("figures/4c.svg")
with rc_context({"figure.figsize": (2, 1)}):
    for gene in ("MT2A", "MT1E", "MT1M", "MT1F", "MT1G", "MT1H", "MT1X"):
        sc.pl.spatial(a30_0_normalized, color=gene, size=1.5, crop_coord=[0, 1000, 1500, 1200], return_fig=True).savefig(f"figures/4d_{gene}.svg", dpi=1000)

圖片.png

allsignifgenes = svg_full[svg_full.padj < 0.001].gene.to_numpy()
upper_patterns, upper_a30_0 = SpatialDE.spatial_patterns(a30_0[np.isin(segm.labels, (0, 3,4,7,8,9,11)),:], genes=allsignifgenes, rng=np.random.default_rng(seed=42), copy=True)
with rc_context({"figure.figsize": (2,1)}):
    for i in range(upper_patterns.patterns.shape[1]):
        sc.pl.spatial(upper_a30_0, color=f"spatial_pattern_{i}", title=f"spatial pattern {i}", crop_coord=[0, 1000, 1500, 1000], size=1.5, return_fig=True).savefig(f"figures/4e_pattern_{i}.svg", dpi=1000)

圖片.png

upper_a30_0_normalized = a30_0_normalized[np.isin(segm.labels, (0, 3,4,7,8,9,11)),:]
with rc_context({"figure.figsize": (2,1)}):
    for gene in allsignifgenes[upper_patterns.labels == 1]:
        sc.pl.spatial(upper_a30_0_normalized, color=gene, size=1.5, crop_coord=[0, 1000, 1500, 1000], return_fig=True).savefig(f"figures/4e_{gene}.svg", dpi=1000)

圖片.png

with rc_context({"figure.figsize":(2,2)}):
    for gene, vmax in zip(["EPCAM", "SCGB2A2", "PAEP", "ACTA2", "MYLK"], (15, 15, None, None, None)):
        sc.pl.spatial(a30_0_normalized, return_fig=True, color=gene, size=1.5, vmax=vmax).savefig(f"figures/s4b_{gene}.svg", dpi=1000)

圖片.png

rnakeys = ("q05_nUMI_factors", "mean_nUMI_factors", "q95_nUMI_factors")
weightkeys = ("q05_spot_factor", "mean_spot_factor", "q95_spot_factor")
cell2locationdata = ad.read_h5ad("20201207_LocationModelLinearDependentWMultiExperiment_19clusters_20952locations_19980genes/sp.h5ad")
cell2locationdata = cell2locationdata[cell2locationdata.obs["sample"] == "152807"]
a30_0_unfiltered = sc.read_visium("152807")
for key in chain(rnakeys, weightkeys):
    a30_0_unfiltered.obs = a30_0_unfiltered.obs.join(cell2locationdata.obs.set_index("spot_id").filter(like=key), sort=False)
a30_0_celltypernas = tuple(a30_0_unfiltered.obs.filter(like=k) for k in rnakeys)
a30_0_celltypeads = tuple(ad.AnnData(X=ctr.to_numpy().round().astype(np.int32), obs=pd.DataFrame(index=a30_0_unfiltered.obs_names), var=pd.DataFrame(index=ctr.columns), uns=a30_0_unfiltered.uns, obsm=a30_0_unfiltered.obsm) for ctr in a30_0_celltypernas)
slabs = a30_0.obs.segmentation_labels_names[a30_0.obs.segmentation_labels.isin((0, 4,7, 11))]
obs = a30_0_celltypeads[1].obs.loc[slabs.index,:]
obs["segmentation_labels_names"] = slabs
(
    ggplot(obs.
               join(pd.DataFrame(a30_0_celltypeads[1].X, columns=a30_0_celltypeads[1].var_names.str[17:], index=a30_0_celltypeads[1].obs_names)).
               reset_index().
               melt(id_vars=["segmentation_labels_names", "index"], var_name="celltype", value_name="mean_spot_factor").
               assign(segmentation_labels_names=lambda df: df.segmentation_labels_names.cat.remove_unused_categories()).
               groupby("index").
               apply(lambda g: g.assign(rel_mean_spot_factor=lambda g: g.mean_spot_factor / g.mean_spot_factor.sum())),
           aes("celltype", "rel_mean_spot_factor", color="segmentation_labels_names")) +
        geom_sina(scale="width", shape=".", alpha=0.6, size=0.01) +
        labs(x=None, y="fraction of transcripts") +
        guides(color=guide_legend(title="", override_aes={"alpha": 1, "shape": 'o', "size": 2})) +
      theme(figure_size=(10, 2), axis_text_x=element_text(rotation=30, ha="right"))
).save("figures/s4a.svg")
c2l_segm, _ = SpatialDE.tissue_segmentation(a30_0_celltypeads[1], rng=np.random.default_rng(seed=42))


with rc_context({"figure.figsize": (2,2)}):
    sc.pl.spatial(a30_0_celltypeads[1], color="segmentation_labels", title="SpatialDE2", return_fig=True).savefig("figures/5a_spatialde.svg", dpi=1000)

圖片.png

a30_0_celltypeweights = tuple(a30_0_unfiltered.obs.filter(like=k) for k in weightkeys)
a30_0_celltypeweightads = tuple(ad.AnnData(X=ctr.to_numpy(), obs=a30_0_unfiltered.obs, var=pd.DataFrame(index=ctr.columns), uns=a30_0_unfiltered.uns, obsm=a30_0_unfiltered.obsm) for ctr in a30_0_celltypeweights)

sc.pp.neighbors(a30_0_celltypeweightads[1], n_neighbors=10, n_pcs=0, random_state=42, metric="cosine")
sc.tl.leiden(a30_0_celltypeweightads[1], random_state=42, resolution=0.4)
with rc_context({"figure.figsize": (2,2)}):
    sc.pl.spatial(a30_0_celltypeweightads[1], color="leiden", title="Leiden", return_fig=True).savefig("figures/5a_leiden.svg", dpi=1000)

圖片.png

resolutions = (0.1, 0.2, 0.3, 0.4, 0.5, 0.8, 1, 1.2)
with rc_context({"figure.figsize": (2,2)}):
    for r in resolutions:
        sc.tl.leiden(a30_0_celltypeweightads[1], random_state=42, resolution=r, key_added=f"leiden_{r}")
        sc.pl.spatial(a30_0_celltypeweightads[1], color=f"leiden_{r}", title=f"Leiden ({r:.1f})", return_fig=True).savefig(f"figures/s5_leiden_resolution_{r:.1f}.svg", dpi=1000)

圖片.png

a30_0_umap = a30_0.copy()
a30_0_umap.obs["segmentation_labels_c2l"] = a30_0_celltypeads[1].obs.segmentation_labels
sc.pp.highly_variable_genes(a30_0_umap, flavor="seurat_v3", n_top_genes=2000)
sc.pp.normalize_total(a30_0_umap, target_sum=1e4, inplace=True)
sc.pp.log1p(a30_0_umap)
sc.pp.regress_out(a30_0_umap, ['total_counts'])
sc.pp.scale(a30_0_umap, max_value=10)
sc.pp.pca(a30_0_umap, n_comps=10)
sc.pp.neighbors(a30_0_umap, n_neighbors=15)
sc.tl.umap(a30_0_umap, min_dist=0.5)
with rc_context({"figure.figsize": (2,2)}):
    sc.pl.umap(a30_0_umap, color="segmentation_labels_c2l", groups=(1,5,8), size=10, title="", return_fig=True).savefig("figures/5b.svg")
(
    ggplot(a30_0_celltypeads[1].obs.
               query("segmentation_labels in (1,5,8)").
               join(pd.DataFrame(a30_0_celltypeads[1].X, columns=a30_0_celltypeads[1].var_names.str[17:], index=a30_0_celltypeads[1].obs_names)).
               reset_index().
               melt(id_vars=["segmentation_labels", "index"], var_name="celltype", value_name="mean_spot_factor").
               assign(segmentation_labels=lambda df: df.segmentation_labels.cat.remove_unused_categories()).
               groupby("index").
               apply(lambda g: g.assign(rel_mean_spot_factor=lambda g: g.mean_spot_factor / g.mean_spot_factor.sum())),
           aes("celltype", "rel_mean_spot_factor", color="segmentation_labels")) +
        geom_sina(scale="width", shape=".", alpha=0.6, size=0.01) +
        labs(x=None, y="fraction of transcripts") +
        guides(color=guide_legend(title="", override_aes={"alpha": 1, "shape": 'o', "size": 2})) +
      theme(figure_size=(6, 2), axis_text_x=element_text(rotation=30, ha="right"))
).save("figures/5c.svg")
(
    ggplot(a30_0_celltypeads[1].obs.
               join(pd.DataFrame(a30_0_celltypeads[1].X, columns=a30_0_celltypeads[1].var_names.str[17:], index=a30_0_celltypeads[1].obs_names)).
               reset_index().
               melt(id_vars=["segmentation_labels", "index"], var_name="celltype", value_name="mean_spot_factor").
               assign(segmentation_labels=lambda df: df.segmentation_labels.cat.remove_unused_categories()).
               groupby("index").
               apply(lambda g: g.assign(rel_mean_spot_factor=lambda g: g.mean_spot_factor / g.mean_spot_factor.sum())),
           aes("celltype", "rel_mean_spot_factor", color="segmentation_labels")) +
        geom_sina(scale="width", shape=".", alpha=0.6, size=0.01) +
        labs(x=None, y="fraction of transcripts") +
        guides(color=guide_legend(title="", override_aes={"alpha": 1, "shape": 'o', "size": 2}, nrow=1)) +
      theme(figure_size=(15, 2), axis_text_x=element_text(rotation=30, ha="right"), legend_position="top")
).save("figures/5d.svg")


with pd.ExcelWriter("figures/table_s_endometrium.xlsx") as ew:
    svg_full[["gene", "pval", "padj"]].to_excel(ew, sheet_name="SV genes whole slice", index=False)
    svg_byregion[["gene", "pval", "padj", "label"]].rename(columns={"label": "region"}).to_excel(ew, sheet_name="SV genes per region", index=False)
    go_enrichment.rename(columns={"region": "region_name", "label": "region"}).to_excel(ew, sheet_name="SV genes GO enrichment", index=False)
    a30_0.obs[["segmentation_labels", "segmentation_labels_names"]].reset_index().rename(columns={"index": "spatial_barcode", "segmentation_labels": "region", "segmentation_labels_names": "region_name"}).to_excel(ew, sheet_name="segmentation (counts)", index=False)
    a30_0_celltypeads[1].obs.reset_index().rename({"index": "spatial_barcode", "segmentation_labels": "region"}).to_excel(ew, sheet_name="segmentation (cell2location)", index=False)

有點(diǎn)難归薛，生活很好，有你更好