hello，大家好劲绪，今天給大家分享10X空間轉(zhuǎn)錄組一個(gè)新的分析點(diǎn)男窟，圖卷積網(wǎng)絡(luò)，關(guān)于這個(gè)贾富，不知道大家對(duì)機(jī)器學(xué)習(xí)了解多少歉眷，希望認(rèn)真做科研的人，能夠好好學(xué)習(xí)一些數(shù)學(xué)方面的知識(shí)祷安，對(duì)自己大有裨益姥芥。

先放一張效果圖

圖片.png

今天分享的文章在Integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network，看題目就知道汇鞭，圖卷積網(wǎng)絡(luò)是把空間的基因表達(dá)凉唐，空間位置和組織學(xué)等多方面的信息結(jié)合在一起，從而實(shí)現(xiàn)對(duì)10X空間轉(zhuǎn)錄組新的角度的分析霍骄，好了台囱，我們提煉一下關(guān)鍵信息，最后分享示例代碼读整。

Abstract

軟件叫SpaGCN簿训，a graph convolutional network approach that integrates gene expression，spatial location and histology in spatial transcriptomics data analysis米间。

通過圖卷積强品，軟件SpaGCN從其相鄰spot聚合每個(gè)spot的基因表達(dá)，從而能夠識(shí)別具有一致表達(dá)和組織學(xué)的空間域（這也是我們做空間分析重要的一方面）屈糊。 隨后的domain引導(dǎo)差異表達(dá)分析然后檢測(cè)在已識(shí)別domain中具有豐富表達(dá)模式的基因的榛。（也就是說，識(shí)別相似的鄰域來劃分空間區(qū)域逻锐，基本是在替代普通的聚類）夫晌。

Introduction

1雕薪、了解組織中不同細(xì)胞的相對(duì)位置對(duì)于理解疾病病理學(xué)至關(guān)重要，因?yàn)榭臻g信息有助于了解細(xì)胞的基因表達(dá)如何受其周圍環(huán)境的影響以及相鄰區(qū)域如何在基因表達(dá)水平上相互作用 晓淀。（這也是空間轉(zhuǎn)錄組最大的意義）所袁。

2、獲得空間轉(zhuǎn)錄組信息主要有兩種策略 （1）single-molecule fluorescence in situ hybridization (smFISH) based techniques, such as MERFISH2 and seqFISH3, which measure expression level for hundreds of genes with 45 subcellular spatial resolution in a single cell 和（2）spatial barcoding followed by next generation sequencing based techniques, such as SLIDE-seq4 and 10X Genomics Visium, which measure the expression level for thousands of genes in captured locations, referred to as spots（當(dāng)然我們更加關(guān)心的是10X的空間轉(zhuǎn)錄組）凶掰。這些不同的空間轉(zhuǎn)錄組學(xué)技術(shù)使揭示異質(zhì)組織復(fù)雜的轉(zhuǎn)錄結(jié)構(gòu)成為可能燥爷，并增強(qiáng)了我們對(duì)疾病細(xì)胞機(jī)制的理解。

3锄俄、在空間轉(zhuǎn)錄組學(xué)研究中局劲，一個(gè)重要的步驟是確定空間domain，定義為在基因表達(dá)和組織學(xué)上具有空間一致性的區(qū)域奶赠。 （這個(gè)非常重要鱼填，是一切研究的基礎(chǔ)）。當(dāng)然毅戈，識(shí)別空間區(qū)域就會(huì)涉及到空間基因表達(dá)苹丸，空間位置和組織學(xué)方面的內(nèi)容。傳統(tǒng)的聚類方法（K-means and Louvain’s method苇经，在10X單細(xì)胞的分析中廣泛運(yùn)用）只分析基因表達(dá)的特征赘理，很顯然這在空間轉(zhuǎn)錄組的分析上是不夠的，To account for spatial dependency of gene expression, new methods have been developed.專門針對(duì)空間的聚類方法現(xiàn)在已經(jīng)有了stlearn,在聚類之前使用從組織學(xué)圖像中提取的特征以及相鄰點(diǎn)的表達(dá)到個(gè)空間平滑的基因表達(dá)數(shù)據(jù) ;還有BayesSpace扇单，運(yùn)用Bayesian approach通過對(duì)物理上靠近的點(diǎn)賦予更高權(quán)重的先驗(yàn)進(jìn)行聚類分析商模；還有隱馬爾可夫隨機(jī)場(chǎng)方法來模擬基因表達(dá)的空間依賴性等等，盡管這些方法已經(jīng)對(duì)空間劃分了區(qū)域蜘澜，但是缺少生物學(xué)解釋施流。

4、To link spatial domains with biological functions at the gene expression level, it is crucial to identify genes that show enriched expression in the identified domains.由于組織中細(xì)胞類型的空間變異鄙信，不同domain之間基因表達(dá)的差異主要由細(xì)胞類型組成的變化驅(qū)動(dòng)瞪醋。 另一方面，關(guān)于空間位置和相應(yīng)組織學(xué)的信息允許構(gòu)建基于解剖學(xué)的組織分類装诡，這為細(xì)胞類型組成提供了有用的視角银受。其中雖然stlearn已經(jīng)結(jié)合了空間基因表達(dá)，空間位置和組織學(xué)進(jìn)行聚類鸦采，但是細(xì)胞類型差異與組織組織結(jié)構(gòu)之間的推定對(duì)應(yīng)關(guān)系尚不清楚宾巍。而且許多空間區(qū)域在細(xì)胞類型方面高度混合。如果沒有進(jìn)一步的下游基因水平分析渔伯，stlearn檢測(cè)到的空間domain仍然缺乏可解釋性顶霞。 最近，Trendsceek咱旱、SpatialDE 和SPARK等新方法被開發(fā)用于檢測(cè)空間可變基因 (SVG)确丢。這些方法獨(dú)立檢查每個(gè)基因并返回一個(gè) p 值來表示基因的空間變異性。 然而吐限，由于缺乏對(duì)組織分類的考慮鲜侥，這些方法檢測(cè)到的基因沒有保證的空間表達(dá)模式，因此很難將這些基因用于進(jìn)一步的生物學(xué)研究(這也是空間高變基因的檢出結(jié)果很難運(yùn)用的原因)诸典。

5描函、不能將空間domian識(shí)別和 SVG 檢測(cè)視為單獨(dú)的問題，利用軟件SpaGCN狐粱，這是一種基于圖卷積網(wǎng)絡(luò)的方法舀寓，可以聯(lián)合考慮這兩個(gè)問題。 使用帶有附加迭代聚類層的圖卷積網(wǎng)絡(luò)肌蜻，SpaGCN 首先通過構(gòu)建表示數(shù)據(jù)空間依賴性的無向加權(quán)圖將基因表達(dá)互墓、空間位置和組織學(xué)整合在一起來識(shí)別空間域。 對(duì)于每個(gè)空間域蒋搜，SpaGCN 然后通過域信息引導(dǎo)的差異表達(dá)分析檢測(cè)在域中豐富的 SVG 與其周圍區(qū)域篡撵。 SpaGCN 還可以選擇檢測(cè)在給定域中唯一表達(dá)的meta基因。 空間domain以及為這些域檢測(cè)到的相應(yīng) SVG 和meta基因提供了關(guān)于組織中基因表達(dá)空間梯度的全面圖景豆挽。

我們來看看文章里面的分析結(jié)果

SpaGCN軟件的介紹

圖片.png

SpaGCN 首先構(gòu)建一個(gè)圖來表示所有樣本（基于測(cè)序的spot）的關(guān)系育谬，同時(shí)考慮空間位置和組織學(xué)信息；接下來帮哈，SpaGCN 利用圖卷積層聚合來自相鄰spot的基因表達(dá)信息膛檀；然后，SpaGCN 使用聚合的基因表達(dá)矩陣使用無監(jiān)督迭代聚類算法對(duì)spot進(jìn)行聚類娘侍；每個(gè)cluster被視為一個(gè)空間domain咖刃，然后 SpaGCN 從中檢測(cè)通過差異表達(dá)分析在域中豐富的 SVG。

圖片.png

當(dāng)單個(gè)基因不能標(biāo)記一個(gè)空間domain的表達(dá)模式時(shí)私蕾，SpaGCN會(huì)構(gòu)建一個(gè)由多個(gè)SVG組合而成的meta基因來表示該域的基因表達(dá)僵缺。 由于斑點(diǎn)/細(xì)胞的表達(dá)譜受到其局部微環(huán)境的嚴(yán)重影響，SpaGCN 還提供了每個(gè)空間域內(nèi)的子cluster檢測(cè)選項(xiàng)踩叭。 還可以檢測(cè) SVG 以幫助理解每個(gè)子空間domain的功能磕潮。

接下來是軟件的一些運(yùn)用

應(yīng)用于小鼠嗅球數(shù)據(jù)

圖片.png

K-means 只識(shí)別了 3 個(gè)主要空間域，只有很少的點(diǎn)被分配到了域 1 和 3容贝。 Louvain 的方法識(shí)別了 5 個(gè)主要空間域自脯。 然而，由于它沒有考慮空間和組織學(xué)信息斤富，域 2膏潮、3 和 4 的邊界模糊，并且比 SpaGCN 有更多的異常值满力。 相比之下焕参，SpaGCN 檢測(cè)到的域與生物已知的 MOB 5 層結(jié)構(gòu)更吻合轻纪。

然后對(duì)小鼠嗅球數(shù)據(jù)檢測(cè)高變基因

圖片.png

all genes show strong specificity for the corresponding domain（用SpaGCN檢測(cè)到的高變基因），當(dāng)然作者也比較了SpatialDE and SPARK檢測(cè)空間高變基因的結(jié)果叠纷，當(dāng)然刻帚，SpaGCN最具有解釋性。相比之下涩嚣，SpaGCN 檢測(cè)到的所有 SVG 都是特定于領(lǐng)域的崇众，根據(jù)對(duì)層結(jié)構(gòu)的了解，提供了合理的解釋航厚。 注意到顷歌，如果用戶將domian 2和domain 4 組合為 SVG 檢測(cè)中的目標(biāo)域，則 SpaGCN 也可以檢測(cè)到具有清晰但非domain特定空間模式的信息較少的 SVG幔睬，例如 PCP4眯漩。

在小鼠后腦數(shù)據(jù)中的應(yīng)用（還是SpaGCN的結(jié)果更加合理）

圖片.png

同樣的空間高變基因的檢測(cè)比較，The Moran’s I values of SpaGCN detected SVGs are much higher than those detected by SPARK and SpatialDE (median of 0.50 for SpaGCN against 0.21 for SPARK and 0.16 for SpatialDE).

圖片.png

對(duì) SPARK 和 SpatialDE 檢測(cè)到的 SVG 進(jìn)行仔細(xì)檢查后發(fā)現(xiàn)溪窒，大多數(shù) SVG 都存在先前在 MOB 數(shù)據(jù)集中觀察到的兩個(gè)問題之一：它們 (1) 僅在少數(shù)幾個(gè)點(diǎn)中表達(dá)或在大多數(shù)點(diǎn)中高度表達(dá)坤塞，表明 SPARK 和 SpatialDE 的高誤報(bào)率或 (2) 空間可變，但在多個(gè)相鄰空間域中表達(dá)澈蚌，使其難以解釋摹芙。這兩種方法的另一個(gè)限制是來自 SPARK 的 FDR 調(diào)整 p 值和來自 SpatialDE 的 q 值不提供信息。具有相似 p 值/q 值的基因不一定顯示相似的空間模式宛瞄，較小的 p 值/q 值并不能保證更好的空間模式（下圖）浮禾。

圖片.png

圖片.png

SPARK 和 SpatialDE 的 p 值和 q 值分布高度偏向 0（下圖）。

圖片.png

相比之下份汗，SpaGCN 檢測(cè)到的 SVG 在特定空間域中富集（下圖）盈电，

圖片.png

圖片.png

它們的表達(dá)模式可轉(zhuǎn)移到小鼠后腦中的相鄰組織切片（下圖）。

圖片.png

此外杯活，SpaGCN 中實(shí)現(xiàn)的多域自適應(yīng)過濾標(biāo)準(zhǔn)使其能夠消除誤報(bào) SVG匆帚，并確保所有檢測(cè)到的 SVG 具有清晰的空間表達(dá)模式。

為了說明適當(dāng)過濾的重要性旁钧，我們以domain 1吸重、5 和 8 為例。對(duì)于這些domain中的每一個(gè)歪今，SpaGCN 都檢測(cè)到在該區(qū)域中富含的單個(gè) SVG（下圖）嚎幸。

圖片.png

PVALB 在域 1 中富集，TRiM62 在域 8 中富集寄猩。雖然域 1 和域 8 彼此相鄰嫉晶，但這兩個(gè) SVG 仍然可以很好地標(biāo)記這些域。 NRGN 是 SpaGCN 針對(duì)域 5 和 7 檢測(cè)到的 SVG。域 5 和 7 中 NRGN 的高表達(dá)也表明這兩個(gè)域在神經(jīng)解剖學(xué)上相似——均由皮層和海馬的錐體層組成替废。皮層和海馬體都是位于大腦曲面上的區(qū)域箍铭。此后腦組織切片具有域 5 中曲面的頂部和域 7 中曲面的底部。域 5 和域 7 在完整的 3D 重建中是連續(xù)的椎镣，由于方式人為分離該部分被切斷坡疼。因此，除了 NRGN 之外衣陶，SpaGCN 還檢測(cè)到許多其他 SVG，如 APP闸氮、ATP6V1G2剪况、CALM2、CHN1蒲跨、CLSTN1译断、ARPP21、CYP46A1或悲、DCLK1孙咪、LINGO1 和 MARCKS，這并不奇怪巡语，它們?cè)趦蓚€(gè)域中都高度表達(dá) 5和 7翎蹈。SpaGCN 中獨(dú)特而強(qiáng)大的 SVG 檢測(cè)程序可確保不會(huì)遺漏此類基因。

（重點(diǎn)）meta gene 的識(shí)別

SpaGCN 沒有識(shí)別域 0 的任何 SVG男公。然而荤堪，我們推斷由多個(gè)基因組合形成的meta gene可能比任何單個(gè)基因更好地揭示空間模式。

圖片.png

首先枢赔，通過降低過濾閾值澄阳，SpaGCN 識(shí)別出在結(jié)構(gòu)域 0 下部高度表達(dá)的 KLK6。使用 KLK6 作為起始基因踏拜，SpaGCN 使用一種新方法找到了 KLK6碎赢、MBP 和 MBP 基因表達(dá)的對(duì)數(shù)線性組合。 ATP1B1速梗，準(zhǔn)確標(biāo)記了空間域 0肮塞。在該元基因中，KLK6 和 MBP 被認(rèn)為是陽性標(biāo)記镀琉，因?yàn)樗鼈冊(cè)谟?0 的某些點(diǎn)上高表達(dá)峦嗤，而 ATP1B1 被認(rèn)為是陰性標(biāo)記，因?yàn)樗饕谄渌麉^(qū)域表達(dá)比域 0屋摔。以前的研究表明 KLK6 和 MBP 的表達(dá)僅限于少突膠質(zhì)細(xì)胞烁设，而 ATP1B1 主要在神經(jīng)元和星形膠質(zhì)細(xì)胞中表達(dá)。這與域 0 代表由少突膠質(zhì)細(xì)胞主導(dǎo)且?guī)缀鯖]有神經(jīng)元細(xì)胞體的白質(zhì)的事實(shí)產(chǎn)生共鳴。因此装黑，構(gòu)成這個(gè)元基因的基因具有有意義的生物學(xué)解釋副瀑。使用這種元基因檢測(cè)程序，我們還檢測(cè)了域 2恋谭、7糠睡、8 和 9 的元基因，并發(fā)現(xiàn)這些元基因可轉(zhuǎn)移到相鄰的組織切片（下圖）疚颊。

圖片.png

一個(gè)spot的表達(dá)譜和生物學(xué)功能受其相鄰點(diǎn)的嚴(yán)重影響狈孔。周圍的spot可以觸發(fā)響應(yīng)路徑或向該點(diǎn)發(fā)出信號(hào)以執(zhí)行某些任務(wù)（其實(shí)就是我們經(jīng)常談到的臨近通訊）。盡管 SpaGCN 檢測(cè)到的一個(gè)空間域中的spot在空間上是一致的材义，并且具有相似的基因表達(dá)模式均抽，但由于它們周圍的spot不同，它們可能仍然具有不同的功能其掂。例如油挥，與位于空間域內(nèi)部的spot相比，位于空間域邊界附近的spot可能具有不同的功能款熬。為了更多地了解不同鄰域?qū)pot的影響深寥，我們進(jìn)行了子域檢測(cè)。例如贤牛，域 2 位于組織切片的中心并被多個(gè)其他空間域包圍惋鹅。因此，域 2 中的spot的相鄰環(huán)境會(huì)發(fā)生變化殉簸。如下圖所示负饲，域 2 被分成 5 個(gè)子域，它們位于域 2 的中心或不同邊界區(qū)域喂链，這表明spot鄰域的差異導(dǎo)致域內(nèi)異質(zhì)性返十。為每個(gè)子域檢測(cè)到的 SVG 可以幫助我們了解每個(gè)子域內(nèi)spot的基因表達(dá)變異性。

圖片.png

LIBD 人類背外側(cè)前額葉皮層數(shù)據(jù)的應(yīng)用（同樣的套路）

圖片.png

空間高變基因的檢測(cè)椭微，Patterns of SVGs for other domains are not very clear

圖片.png

SVGs detected by SPARK lack domain specificity（下圖）

圖片.png

SpatialDE detected 3,654 SVGs with 806 of them having q-values equal to 0, but these genes do not necessarily show better spatial pattern than genes with larger q-values

圖片.png

這兩種方法檢測(cè)到的基因缺乏區(qū)分表達(dá)中不同程度空間變異性的能力洞坑，因?yàn)樗鼈兊?p 值和 q 值分布高度偏向于 0

圖片.png

下圖顯示 SpaGCN 檢測(cè)到的 SVG 的 Moran's I 值顯著高于 SpatialDE 和 SPARK 檢測(cè)到的值。

圖片.png

有些區(qū)域可以識(shí)別單個(gè)specific gene（下圖）

圖片.png

但更合理的解釋是空間高變的meta gene

圖片.png

圖片.png

圖片.png

在人原發(fā)性胰腺癌組織中的應(yīng)用

圖片.png

上圖shows that domain 2 detected by SpaGCN has a better correspondence to the cancer region than clusters reported in the original study蝇率。

空間高變基因的檢測(cè)迟杂，SpaGCN的可解釋性最好

圖片.png

圖片.png

圖片.png

在 MERFISH 小鼠下丘腦數(shù)據(jù)中的應(yīng)用

圖片.png

Method，主要關(guān)注一下軟件的算法本慕，當(dāng)然非常難排拷，我們下一篇詳細(xì)分享。

示例代碼

安裝

pip3 install SpaGCN
#Note: you need to make sure that the pip is for python3锅尘，or we should install SpaGCN by
python3 -m pip install SpaGCN
pip3 install SpaGCN
#If you do not have permission (when you get a permission denied error), you should install SpaGCN by
pip3 install --user SpaGCN

加載,Import python modules

import os,csv,re
import pandas as pd
import numpy as np
import scanpy as sc
import math
import SpaGCN as spg
from scipy.sparse import issparse
import random, torch
import warnings
warnings.filterwarnings("ignore")
import matplotlib.colors as clr
import matplotlib.pyplot as plt
import SpaGCN as spg
#In order to read in image data, we need to install some package. Here we recommend package "opencv"
#inatll opencv in python
#!pip3 install opencv-python
import cv2

讀取數(shù)據(jù)

當(dāng)前版本的 SpaGCN 需要三個(gè)輸入數(shù)據(jù)：

• 基因表達(dá)矩陣（n by k）：expression_matrix.h5监氢；
• samplespositions.txt 的空間坐標(biāo)；
• 組織學(xué)圖像（可選）：histology.tif，可以是 tif 或 png 或 jepg浪腐。
• 基因表達(dá)數(shù)據(jù)可以存儲(chǔ)為 AnnData 對(duì)象纵揍。 AnnData 將數(shù)據(jù)矩陣 .X 與觀察注釋 .obs、變量 .var 和非結(jié)構(gòu)化注釋 .uns 一起存儲(chǔ)议街。

"""
#Read original 10x_h5 data and save it to h5ad
from scanpy import read_10x_h5
adata = read_10x_h5("../tutorial/data/151673/expression_matrix.h5")
spatial=pd.read_csv("../tutorial/data/151673/positions.txt",sep=",",header=None,na_filter=False,index_col=0) 
adata.obs["x1"]=spatial[1]
adata.obs["x2"]=spatial[2]
adata.obs["x3"]=spatial[3]
adata.obs["x4"]=spatial[4]
adata.obs["x5"]=spatial[5]
#Select captured samples
adata=adata[adata.obs["x1"]==1]
adata.var_names=[i.upper() for i in list(adata.var_names)]
adata.var["genename"]=adata.var.index.astype("str")
adata.write_h5ad("../tutorial/data/151673/sample_data.h5ad")
"""
#Read in gene expression and spatial location
adata=sc.read("../tutorial/data/151673/sample_data.h5ad")
#Read in hitology image
img=cv2.imread("../tutorial/data/151673/histology.tif")

Integrate gene expression and histology into a Graph

#Set coordinates
adata.obs["x_array"]=adata.obs["x2"]
adata.obs["y_array"]=adata.obs["x3"]
adata.obs["x_pixel"]=adata.obs["x4"]
adata.obs["y_pixel"]=adata.obs["x5"]
x_array=adata.obs["x_array"].tolist()
y_array=adata.obs["y_array"].tolist()
x_pixel=adata.obs["x_pixel"].tolist()
y_pixel=adata.obs["y_pixel"].tolist()

#Test coordinates on the image
img_new=img.copy()
for i in range(len(x_pixel)):
    x=x_pixel[i]
    y=y_pixel[i]
    img_new[int(x-20):int(x+20), int(y-20):int(y+20),:]=0

cv2.imwrite('./sample_results/151673_map.jpg', img_new)

• “s”參數(shù)確定在計(jì)算每?jī)蓚€(gè)點(diǎn)之間的歐幾里得距離時(shí)賦予組織學(xué)的權(quán)重泽谨。 's = 1' 表示組織學(xué)像素強(qiáng)度值與 (x,y) 坐標(biāo)具有相同的尺度方差，而 's' 值越高表示尺度方差越大特漩，因此吧雹，在計(jì)算歐幾里得距離時(shí)，組織學(xué)的權(quán)重越高 .
• “b”參數(shù)確定提取顏色強(qiáng)度時(shí)每個(gè)點(diǎn)的面積涂身。

#Calculate adjacent matrix
s=1
b=49
adj=spg.calculate_adj_matrix(x=x_pixel,y=y_pixel, x_pixel=x_pixel, y_pixel=y_pixel, image=img, beta=b, alpha=s, histology=True)
#If histlogy image is not available, SoaGCN can calculate the adjacent matrix using the fnction below
#adj=calculate_adj_matrix(x=x_pixel,y=y_pixel, histology=False)
np.savetxt('./data/151673/adj.csv', adj, delimiter=',')

Spatial domain detection using SpaGCN(空間區(qū)域檢測(cè))

adata=sc.read("./data/151673/sample_data.h5ad")
adj=np.loadtxt('./data/151673/adj.csv', delimiter=',')
adata.var_names_make_unique()
spg.prefilter_genes(adata,min_cells=3) # avoiding all genes are zeros
spg.prefilter_specialgenes(adata)
#Normalize and take log for UMI
sc.pp.normalize_per_cell(adata)
sc.pp.log1p(adata)

設(shè)置超參數(shù)

• p：社區(qū)貢獻(xiàn)的總表達(dá)百分比吮炕。
• l：控制p的參數(shù)。

p=0.5 
#Find the l value given p
l=spg.search_l(p, adj, start=0.01, end=1000, tol=0.01, max_run=100)

• n_clusters: Number of spatial domains wanted.
• res: Resolution in the initial Louvain's Clustering methods. If the number of clusters is known, we can use the spg.search_res() fnction to search for suitable resolution(optional)

#If the number of clusters known, we can use the spg.search_res() fnction to search for suitable resolution(optional)
#For this toy data, we set the number of clusters=7 since this tissue has 7 layers
n_clusters=7
#Set seed
r_seed=t_seed=n_seed=100
#Seaech for suitable resolution
res=spg.search_res(adata, adj, l, n_clusters, start=0.7, step=0.1, tol=5e-3, lr=0.05, max_epochs=20, r_seed=r_seed, t_seed=t_seed, n_seed=n_seed)

運(yùn)行SpaGCN

clf=spg.SpaGCN()
clf.set_l(l)
#Set seed
random.seed(r_seed)
torch.manual_seed(t_seed)
np.random.seed(n_seed)
#Run
clf.train(adata,adj,init_spa=True,init="louvain",res=res, tol=5e-3, lr=0.05, max_epochs=200)
y_pred, prob=clf.predict()
adata.obs["pred"]= y_pred
adata.obs["pred"]=adata.obs["pred"].astype('category')
#Do cluster refinement(optional)
#shape="hexagon" for Visium data, "square" for ST data.
adj_2d=spg.calculate_adj_matrix(x=x_array,y=y_array, histology=False)
refined_pred=spg.refine(sample_id=adata.obs.index.tolist(), pred=adata.obs["pred"].tolist(), dis=adj_2d, shape="hexagon")
adata.obs["refined_pred"]=refined_pred
adata.obs["refined_pred"]=adata.obs["refined_pred"].astype('category')
#Save results
adata.write_h5ad("./sample_results/results.h5ad")

Plot spatial domains

adata=sc.read("./sample_results/results.h5ad")
#Set colors used
plot_color=["#F56867","#FEB915","#C798EE","#59BE86","#7495D3","#D1D1D1","#6D1A9C","#15821E","#3A84E6","#997273","#787878","#DB4C6C","#9E7A7A","#554236","#AF5F3C","#93796C","#F9BD3F","#DAB370","#877F6C","#268785"]
#Plot spatial domains
domains="pred"
num_celltype=len(adata.obs[domains].unique())
adata.uns[domains+"_colors"]=list(plot_color[:num_celltype])
ax=sc.pl.scatter(adata,alpha=1,x="y_pixel",y="x_pixel",color=domains,title=domains,color_map=plot_color,show=False,size=100000/adata.shape[0])
ax.set_aspect('equal', 'box')
ax.axes.invert_yaxis()
plt.savefig("./sample_results/pred.png", dpi=600)
plt.close()

#Plot refined spatial domains
domains="refined_pred"
num_celltype=len(adata.obs[domains].unique())
adata.uns[domains+"_colors"]=list(plot_color[:num_celltype])
ax=sc.pl.scatter(adata,alpha=1,x="y_pixel",y="x_pixel",color=domains,title=domains,color_map=plot_color,show=False,size=100000/adata.shape[0])
ax.set_aspect('equal', 'box')
ax.axes.invert_yaxis()
plt.savefig("./sample_results/refined_pred.png", dpi=600)
plt.close()

圖片.png

Identify SVGs

#Read in raw data
raw=sc.read("../tutorial/data/151673/sample_data.h5ad")
raw.var_names_make_unique()
raw.obs["pred"]=adata.obs["pred"].astype('category')
raw.obs["x_array"]=raw.obs["x2"]
raw.obs["y_array"]=raw.obs["x3"]
raw.obs["x_pixel"]=raw.obs["x4"]
raw.obs["y_pixel"]=raw.obs["x5"]
#Convert sparse matrix to non-sparse
raw.X=(raw.X.A if issparse(raw.X) else raw.X)
raw.raw=raw
sc.pp.log1p(raw)

• target：用于識(shí)別 SVG 的目標(biāo)域访得。
• min_in_group_fraction：最小組內(nèi)表達(dá)分?jǐn)?shù)。
• min_in_out_group_ratio：最小值（組內(nèi)表達(dá)分?jǐn)?shù)）/（組外表達(dá)分?jǐn)?shù)）陕凹。
• min_fold_change：最泻芬帧（組內(nèi)表達(dá)式）/（組外表達(dá)式）。
• r：檢測(cè)一個(gè)點(diǎn)的相鄰點(diǎn)的半徑皿渗。

#Use domain 0 as an example
target=0
#Set filtering criterials
min_in_group_fraction=0.8
min_in_out_group_ratio=1
min_fold_change=1.5
#Search radius such that each spot in the target domain has approximately 10 neighbors on average
adj_2d=spg.calculate_adj_matrix(x=x_array, y=y_array, histology=False)
start, end= np.quantile(adj_2d[adj_2d!=0],q=0.001), np.quantile(adj_2d[adj_2d!=0],q=0.1)
r=spg.search_radius(target_cluster=target, cell_id=adata.obs.index.tolist(), x=x_array, y=y_array, pred=adata.obs["pred"].tolist(), start=start, end=end, num_min=10, num_max=14,  max_run=100)
#Detect neighboring domains
nbr_domians=spg.find_neighbor_clusters(target_cluster=target,
                                   cell_id=raw.obs.index.tolist(), 
                                   x=raw.obs["x_array"].tolist(), 
                                   y=raw.obs["y_array"].tolist(), 
                                   pred=raw.obs["pred"].tolist(),
                                   radius=r,
                                   ratio=1/2)

nbr_domians=nbr_domians[0:3]
de_genes_info=spg.rank_genes_groups(input_adata=raw,
                                target_cluster=target,
                                nbr_list=nbr_domians, 
                                label_col="pred", 
                                adj_nbr=True, 
                                log=True)
#Filter genes
de_genes_info=de_genes_info[(de_genes_info["pvals_adj"]<0.05)]
filtered_info=de_genes_info
filtered_info=filtered_info[(filtered_info["pvals_adj"]<0.05) &
                            (filtered_info["in_out_group_ratio"]>min_in_out_group_ratio) &
                            (filtered_info["in_group_fraction"]>min_in_group_fraction) &
                            (filtered_info["fold_change"]>min_fold_change)]
filtered_info=filtered_info.sort_values(by="in_group_fraction", ascending=False)
filtered_info["target_dmain"]=target
filtered_info["neighbors"]=str(nbr_domians)
print("SVGs for domain ", str(target),":", filtered_info["genes"].tolist())
filtered_info

#Plot refinedspatial domains
color_self = clr.LinearSegmentedColormap.from_list('pink_green', ['#3AB370',"#EAE7CC","#FD1593"], N=256)
for g in filtered_info["genes"].tolist():
    raw.obs["exp"]=raw.X[:,raw.var.index==g]
    ax=sc.pl.scatter(raw,alpha=1,x="y_pixel",y="x_pixel",color="exp",title=g,color_map=color_self,show=False,size=100000/raw.shape[0])
    ax.set_aspect('equal', 'box')
    ax.axes.invert_yaxis()
    plt.savefig("./sample_results/"+g+".png", dpi=600)
    plt.close()

圖片.png

最重要的部分障癌， Identify Meta Gene

#Use domain 2 as an example
target=2
meta_name, meta_exp=spg.find_meta_gene(input_adata=raw,
                    pred=raw.obs["pred"].tolist(),
                    target_domain=target,
                    start_gene="GFAP",
                    mean_diff=0,
                    early_stop=True,
                    max_iter=3,
                    use_raw=False)

raw.obs["meta"]=meta_exp
#Plot meta gene
g="GFAP"
raw.obs["exp"]=raw.X[:,raw.var.index==g]
ax=sc.pl.scatter(raw,alpha=1,x="y_pixel",y="x_pixel",color="exp",title=g,color_map=color_self,show=False,size=100000/raw.shape[0])
ax.set_aspect('equal', 'box')
ax.axes.invert_yaxis()
plt.savefig("./sample_results/"+g+".png", dpi=600)
plt.close()

raw.obs["exp"]=raw.obs["meta"]
ax=sc.pl.scatter(raw,alpha=1,x="y_pixel",y="x_pixel",color="exp",title=meta_name,color_map=color_self,show=False,size=100000/raw.shape[0])
ax.set_aspect('equal', 'box')
ax.axes.invert_yaxis()
plt.savefig("./sample_results/meta_gene.png", dpi=600)
plt.close()

圖片.png

Multiple tissue sections analysis(聯(lián)合分析)

圖片.png

讀取數(shù)據(jù)

adata1=sc.read("./data/Mouse_brain/MA1.h5ad")
adata2=sc.read("./data/Mouse_brain/MP1.h5ad")
img1=cv2.imread("./data/Mouse_brain/MA1_histology.tif")
img2=cv2.imread("./data/Mouse_brain/MP1_histology.tif")

Extract color intensity

b=49
s=1
x_pixel1=adata1.obs["x4"].tolist()
y_pixel1=adata1.obs["x5"].tolist()
adata1.obs["color"]=spg.extract_color(x_pixel=x_pixel1, y_pixel=y_pixel1, image=img1, beta=b)
z_scale=np.max([np.std(x_pixel1), np.std(y_pixel1)])*s
adata1.obs["z"]=(adata1.obs["color"]-np.mean(adata1.obs["color"]))/np.std(adata1.obs["color"])*z_scale

x_pixel2=adata2.obs["x4"].tolist()
y_pixel2=adata2.obs["x5"].tolist()
adata2.obs["color"]=spg.extract_color(x_pixel=x_pixel2, y_pixel=y_pixel2, image=img2, beta=b)
z_scale=np.max([np.std(x_pixel2), np.std(y_pixel2)])*s
adata2.obs["z"]=(adata2.obs["color"]-np.mean(adata2.obs["color"]))/np.std(adata2.obs["color"])*z_scale
del img1, img2

Modify coordinates to combine 2 sections

from anndata import AnnData
adata1.obs["x_pixel"]=x_pixel1
adata1.obs["y_pixel"]=y_pixel1
adata2.obs["x_pixel"]=x_pixel2-np.min(x_pixel2)+np.min(x_pixel1)
adata2.obs["y_pixel"]=y_pixel2-np.min(y_pixel2)+np.max(y_pixel1)
adata1.var_names_make_unique()
adata2.var_names_make_unique()
adata_all=AnnData.concatenate(adata1, adata2,join='inner',batch_key="dataset_batch",batch_categories=["0","1"])

Integrate gene expression and histology into a Graph

X=np.array([adata_all.obs["x_pixel"], adata_all.obs["y_pixel"], adata_all.obs["z"]]).T.astype(np.float32)
adj=spg.pairwise_distance(X)

Spatial domain detection using SpaGCN

sc.pp.normalize_per_cell(adata_all, min_counts=0)
sc.pp.log1p(adata_all)
p=0.5 
#Find the l value given p
l=spg.search_l(p, adj, start=0.01, end=1000, tol=0.01, max_run=100)
res=1.0
seed=100
random.seed(seed)
torch.manual_seed(seed)
np.random.seed(seed)
clf=spg.SpaGCN()
clf.set_l(l)
clf.train(adata_all,adj,init_spa=True,init="louvain",res=res, tol=5e-3, lr=0.05, max_epochs=200)
y_pred, prob=clf.predict()
adata_all.obs["pred"]= y_pred
adata_all.obs["pred"]=adata_all.obs["pred"].astype('category')


colors_use=['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728', '#9467bd', '#8c564b', '#e377c2', '#bcbd22', '#17becf', '#aec7e8', '#ffbb78', '#98df8a', '#ff9896', '#bec1d4', '#bb7784', '#0000ff', '#111010', '#FFFF00',   '#1f77b4', '#800080', '#959595', 
 '#7d87b9', '#bec1d4', '#d6bcc0', '#bb7784', '#8e063b', '#4a6fe3', '#8595e1', '#b5bbe3', '#e6afb9', '#e07b91', '#d33f6a', '#11c638', '#8dd593', '#c6dec7', '#ead3c6', '#f0b98d', '#ef9708', '#0fcfc0', '#9cded6', '#d5eae7', '#f3e1eb', '#f6c4e1', '#f79cd4']
num_celltype=len(adata_all.obs["pred"].unique())
adata_all.uns["pred_colors"]=list(colors_use[:num_celltype])
ax=sc.pl.scatter(adata_all,alpha=1,x="y_pixel",y="x_pixel",color="pred",show=False,size=150000/adata_all.shape[0])
ax.set_aspect('equal', 'box')
ax.axes.invert_yaxis()
plt.savefig("./sample_results/mouse_barin_muti_sections_domains.png", dpi=600)
plt.close()

圖片.png

非常棒的方法稿壁，大家一定要試一試

生活很好，有你更好