hello，今天給大家分享一個檢測空間高變基因的方法笑跛，分享這個內(nèi)容的主要原因是下一篇我們要利用空間高變基因運(yùn)用到我們的空間轉(zhuǎn)錄組分析付魔，文獻(xiàn)在Statistical analysis of spatial expression patterns for spatially resolved transcriptomic studies聊品，nature method，IF27几苍，軟件在SPARK翻屈，不過最近作者推出了SPARK加強(qiáng)版，SPARK-X妻坝，我們后面也分享一下代碼伸眶，其實(shí)關(guān)于10X空間轉(zhuǎn)錄組檢測空間高變基因，大家可以參考之前的文章10X空間轉(zhuǎn)錄組-----空間高變基因檢測之SpatialDE刽宪、10X空間轉(zhuǎn)錄組之基因的空間表達(dá)模式赚抡、trendsceek、HMRF纠屋，每種算法都有其優(yōu)勢和適用場景涂臣。對于10X空間轉(zhuǎn)錄組的空間高變基因該如何運(yùn)用，我們下一篇分享售担，這一篇來看看這個軟件赁遗。
SPARK軟件算法通過廣義線性空間模型直接對空間的count數(shù)據(jù)進(jìn)行建模，用于識別從各種空間轉(zhuǎn)錄組學(xué)技術(shù)生成的數(shù)據(jù)中基因的空間表達(dá)模式族铆。通訊作者是來自于美國密西根大學(xué)生物統(tǒng)計學(xué)的周翔教授（中國人）岩四。
圖片.png
Abstract
鑒定出空間轉(zhuǎn)錄組數(shù)據(jù)中呈現(xiàn)具有空間表達(dá)模式的基因是表征復(fù)雜組織的空間轉(zhuǎn)錄整體概況的重要第一步。在本文中作者提出了一種統(tǒng)計方法——SPARK哥攘，用于識別從空間轉(zhuǎn)錄組數(shù)據(jù)中基因的空間表達(dá)模式剖煌。SPARK通過廣義線性空間模型直接對空間轉(zhuǎn)錄組數(shù)據(jù)進(jìn)行建模。它依靠最新開發(fā)的統(tǒng)計公式進(jìn)行假設(shè)檢驗(yàn)逝淹，可有效控制I型錯誤并具有較高的統(tǒng)計能力耕姊。利用基于懲罰擬似然的高效計算算法，SPARK的應(yīng)用可以擴(kuò)展到包含數(shù)萬個樣本的數(shù)萬個基因的數(shù)據(jù)集栅葡。本文中作者使用SPARK分析了四個已發(fā)表的空間轉(zhuǎn)錄組數(shù)據(jù)集茉兰，發(fā)現(xiàn)它的功效可能是現(xiàn)有方法的十倍之多，并揭示了新的生物學(xué)發(fā)現(xiàn)欣簇。
introduction
有效的基因空間表達(dá)模式分析（簡寫為SE）面臨著重要的統(tǒng)計和計算挑戰(zhàn)规脸。當(dāng)前僅有的兩種SE分析方法：SpatialDE和Trendsceek都存在能效不足的問題，并且在應(yīng)用于具有數(shù)百個基因且跨數(shù)百個空間位置的大型數(shù)據(jù)集時熊咽，在計算方面都存在一定的挑戰(zhàn)莫鸭。SE分析的一般步驟包括：空間表達(dá)的模型構(gòu)建，對鑒定出的具有空間表達(dá)模式基因的進(jìn)行統(tǒng)計檢驗(yàn)以及可視化横殴。作者提出了一種基于核函數(shù)（kernels）進(jìn)行空間模式識別的SPARK算法被因，是一種旨在解決上述這些問題并更好地控制I類錯誤的新方法，可以從空間轉(zhuǎn)錄組數(shù)據(jù)中提取具有空間表達(dá)模式的基因。具體地說（圖1）氏身，SPARK首先以原始計數(shù)的形式對空間數(shù)據(jù)進(jìn)行建模巍棱，該數(shù)據(jù)呈現(xiàn)過度分散的Poisson分布并包含一個未知參數(shù)，該參數(shù)表示了每個位置的潛在空間隨機(jī)過程蛋欣，包括三個組成部分：1.觀察到的解釋變量（例如批效應(yīng)或細(xì)胞周期信息）航徙；2.由于隨機(jī)噪聲引起的殘余誤差；3.用于描述數(shù)據(jù)空間自相關(guān)模式的廣義線性空間模型（generalized linear spatial model陷虎，GLSM）到踏。GLSM模型是基于核函數(shù)展開的，而鑒定具有空間表達(dá)模式基因的統(tǒng)計能力將取決于空間核函數(shù)K(s)的選擇尚猿。具體地說窝稿，不同的核函數(shù)適用于具有不同空間表達(dá)模式的基因，核函數(shù)與目標(biāo)基因所展現(xiàn)的潛在空間模式的匹配程度將決定該方法鑒定空間可變基因的能力凿掂。為了確卑槔疲可靠地識別各種空間表達(dá)模式基因，作者考慮使用總共十個不同的空間核函數(shù)庄萎，包括五個具有不同周期性參數(shù)的周期核函數(shù)和五個具有不同平滑度參數(shù)的高斯核函數(shù)踪少。為了克服解高維積分的難題，使利用GLSM實(shí)現(xiàn)可擴(kuò)展的估計和推斷成為可能糠涛，作者基于懲罰擬似然算法（penalized quasi-likelihood algorithm援奢，PQL）來進(jìn)行參數(shù)估計，通過柯西組合規(guī)則（Cauchy combination rule）得到一個最終校正后的P值
圖片.png
- 注：SPARK方法的示意圖
Result
模擬數(shù)據(jù)：與Trendsceek方法不同忍捡，SpatialDE和SPARK具有基礎(chǔ)的數(shù)據(jù)生成模型集漾。但SPARK與SpatialDE不同的是，SPARK模型直接對數(shù)據(jù)進(jìn)行計數(shù)砸脊，并依靠適當(dāng)?shù)慕y(tǒng)計過程來獲得校準(zhǔn)的P值具篇。作者使用模擬數(shù)據(jù)來評估SPARK的性能，并將其與兩種現(xiàn)有方法SpatialDE和Trendsceek進(jìn)行了比較脓规。在這里作者主要考察方法兩個方面的性能：第一個是其控制I類錯誤的能力栽连，第二是SPARK方法鑒定空間可變基因的能力险领。從圖2b和c可以看出侨舆，在原假設(shè)條件下，SPARK在整個轉(zhuǎn)錄組范圍內(nèi)得到了最為準(zhǔn)確的P值绢陌。Trendsceek方法中的Markvario和Vmark檢驗(yàn)統(tǒng)計得到了較為合理的P值挨下，而其他Emark統(tǒng)計和markcorr統(tǒng)計得出的P值則稍微偏保守。相比之下脐湾，SpatialDE產(chǎn)生的P值過于保守臭笆。圖2d和e表征的是不同方法鑒別不同空間表達(dá)模式的空間可變基因的能力。結(jié)果證明，在不同F(xiàn)DR下愁铺， SPARK相比其他方法都更為有效鹰霍。由于Trendsceek的性能相對較差。作者對第二組模擬數(shù)據(jù)進(jìn)行了分析茵乱，這組數(shù)據(jù)模擬完全基于原始Trendsceek文獻(xiàn)茂洒，在這組模擬數(shù)據(jù)中也得到了類似的結(jié)論∑拷撸總體而言督勺，SPARK在檢測空間可變基因時比其他兩種方法有更好控制I類錯誤的能力，并且檢出能力更佳斤贰。
圖片.png
真實(shí)數(shù)據(jù)1-小鼠嗅球數(shù)據(jù)：它由在260個位點(diǎn)上測得的11,274個基因組成智哀。從圖3a中可以看出，與模擬數(shù)據(jù)的結(jié)果一致荧恍，SPARK和Trendsceek均產(chǎn)生了較準(zhǔn)確的P值瓷叫，而SpatialDE則沒有。與SpatialDE和Trendsceek相比送巡，SPARK鑒定出更多的具有空間可變基因（圖3b）赞辩。并且，僅由SPARK授艰、SpatialDE和兩者共同檢測出的空間可變基因在整個空間位點(diǎn)上的表達(dá)量辨嗽。作者發(fā)現(xiàn)大多數(shù)僅由SpatialDE檢測到的空間可變基因的表達(dá)水平趨于零，這些基因在空間表達(dá)譜只在一個或兩個位點(diǎn)上有表達(dá)淮腾，表明這些可能是潛在的假陽信號糟需。相反，僅由SPARK檢測到的空間可變基因通常具有與被兩種方法同時檢測到的空間可變基因相當(dāng)?shù)谋磉_(dá)水平谷朝。
圖片.png
為了評估由SPARK鑒定的空間可變基因的質(zhì)量洲押，作者又對由SPARK鑒定的772 基因進(jìn)行了聚類分析，獲得了三種主要的空間表達(dá)模式（下圖）：模式一代表二尖瓣細(xì)胞層圆凰；模式二代表腎小球?qū)予菊剩荒Ｊ饺砹罴?xì)胞層。通過對三個層的先前已知的標(biāo)記基因的表達(dá)譜圖進(jìn)行展示专钉，可以清晰地看出三種空間表達(dá)模式（圖3c和d）挑童。進(jìn)一步，在僅由SPARK檢測到的空間可變基因中跃须，三個模式每個模式隨機(jī)挑選20個基因進(jìn)行表達(dá)譜展示站叼。結(jié)果顯示，幾乎所有這些基因都顯示出清晰的空間表達(dá)模式菇民，從而確認(rèn)了SPARK的更高功效尽楔。
作者提供了另外三條證據(jù)來驗(yàn)證SPARK檢測空間可變基因的可靠性（圖3g）投储。首先，檢查了該小鼠嗅球數(shù)據(jù)原文章中提出的嗅覺系統(tǒng)中比較顯著的標(biāo)記基因阔馋。作者對分析結(jié)果和此基因列表取交集玛荞，結(jié)果顯示此基因列表中的十個的基因中SPARK能檢測到八個，SpatialDE僅檢測到三個呕寝，而Trendsceek未檢測到冲泥。再者，作者獲得了最近的有關(guān)小鼠嗅球的單細(xì)胞RNA測序研究中鑒定的2,030種細(xì)胞類型特異性標(biāo)記基因的列表壁涎。結(jié)果顯示凡恍，僅由SPARK鑒定出基因中有55％位于列表中，而僅由SpatialDE鑒定只有20％怔球。第三嚼酝，作者在Harmonizome數(shù)據(jù)庫中獲得了3,262個與嗅覺系統(tǒng)相關(guān)的基因。同樣竟坛，僅由SPARK鑒定有26％闽巩，而僅由SpatialDE鑒定只有20％。這三個附加的驗(yàn)證分析證明了SPARK在鑒定SE基因方面的確更有優(yōu)勢担汤。
最后涎跨，作者進(jìn)行了對由SPARK和SpatialDE方法鑒定的空間可變基因的功能富集分析。對于SPARK崭歧，富集得到了1,023個GO注釋和79個KEGG通路（圖3h）隅很；而對SpatialDE鑒定的基因分析，得到了87個GO注釋和兩個KEGG通路率碾。僅由SPARK識別的基因富集得到的許多GO注釋或KEGG通路與突觸組織和嗅球發(fā)育直接相關(guān)叔营，包括嗅葉發(fā)育和催產(chǎn)素信號傳導(dǎo)途徑。分別單獨(dú)選取模式I–III中的基因進(jìn)行的深入富集分析提供了其他生物學(xué)見解所宰，比如其揭示了突觸組織在二尖瓣細(xì)胞層中的關(guān)鍵作用绒尊、細(xì)胞連接和突觸對神經(jīng)層連接的重要性以及樹突形態(tài)發(fā)生和突觸-樹突可塑性對于顆粒層的重要性∽兄啵總的來說婴谱，基因的功能富集分析也彰顯出SPARK的優(yōu)勢。
圖片.png
和第一個數(shù)據(jù)的驗(yàn)證方式一樣躯泰，作者提供了另外三個證據(jù)來驗(yàn)證SPARK檢測到的SE基因（圖4c）谭羔。首先，作者將數(shù)據(jù)原文中得到的14個顯著與癌癥相關(guān)的基因與三種方法得到的結(jié)果進(jìn)行交叉驗(yàn)證：SPARK檢測到其中的十個基因斟冕，SpatialDE檢測到七個基因口糕，Trendsceek檢測到兩個基因。SpatialDE和Trendsceek都沒有檢測到三個與癌癥相關(guān)的著名基因：SCGB2A2磕蛇、KRT17和MMP14景描。圖4e展示了五個僅由SPARK識別的腫瘤基因（HLA-B、EEF1A1秀撇、ERBB2超棺、MMP14和CD44）的空間表達(dá)譜。圖4e左上角的圖展示了對樣本進(jìn)行染色后的樣本圖呵燕，深色表示潛在的腫瘤區(qū)域棠绘。其顯示出腫瘤區(qū)域和這5種腫瘤基因的空間表達(dá)譜的模式相似。其次再扭，作者在CancerMine數(shù)據(jù)庫中收集了1144個與乳腺癌相關(guān)的基因列表氧苍。僅由SPARK鑒定的基因有14％在列表中，而僅由SpatialDE鑒定的基因有10％在列表中泛范。其中让虐，由大量的先前文獻(xiàn)所證實(shí)的著名原癌基因ERBB2基因僅被SPARK鑒定。第三罢荡，作者從Harmonizome數(shù)據(jù)庫中收集了3538個與乳腺癌相關(guān)的基因列表赡突。僅由SPARK鑒定的基因的44％在列表中，而僅由SpatialDE鑒定的基因的37％在列表中区赵〔宴郑總體而言，這三個附加證據(jù)為SPARK的強(qiáng)大功能提供了支持笼才。
最后漱受，作者進(jìn)行了功能富集分析（圖4f）。在FDR值為5％時骡送，SPARK識別出的空間可變基因功能富集后得到542個GO注釋和20條KEGG通路拜效，而SpatialDE識別了266個GO注釋和三個KEGG通路。功能富集分析結(jié)果顯示各谚，僅由SPARK發(fā)現(xiàn)的許多富集的基因集與細(xì)胞外基質(zhì)的組織和免疫應(yīng)答有關(guān)紧憾，突顯了它們在癌癥發(fā)展和轉(zhuǎn)移中的重要性。
討論
在本文中昌渤，作者提出了一種統(tǒng)計方法SPARK赴穗，用于識別從各種空間轉(zhuǎn)錄組數(shù)據(jù)中具有空間表達(dá)模式的基因。SPARK通過廣義線性空間模型直接對空間數(shù)據(jù)進(jìn)行建模膀息。它依靠最新開發(fā)的統(tǒng)計公式進(jìn)行假設(shè)檢驗(yàn)般眉，可有效控制I類錯誤并具有較高的統(tǒng)計能力。利用基于懲罰擬似然的高效計算算法潜支，使SPARK的應(yīng)用可以擴(kuò)展到包含成千上萬個樣本中成千上萬個基因的數(shù)據(jù)集甸赃。使用SPARK分析已發(fā)布的空間轉(zhuǎn)錄組數(shù)據(jù)集，發(fā)現(xiàn)它的功效可能是現(xiàn)有方法的十倍之多冗酿，并揭示了現(xiàn)有方法無法揭示的生物學(xué)發(fā)現(xiàn)埠对。SPARK直接對原始計數(shù)進(jìn)行建模络断，以說明在空間數(shù)據(jù)中觀察到的均值-方差相關(guān)性，從而使該方法檢測空間可變基因的效能增加项玛。作者綜合了不同核函數(shù)的P值貌笨，以確保SPARK在各種可能的空間表達(dá)模式下的具有穩(wěn)定性能，可靠地識別各種空間表達(dá)模式的基因襟沮，核函數(shù)的類型和數(shù)量的選擇可以根據(jù)不同的研究需求來自行調(diào)整锥惋。總的來說开伏，SPARK是當(dāng)前用于識別空間轉(zhuǎn)錄組數(shù)據(jù)中具有空間表達(dá)模式基因的一種較好的統(tǒng)計分析方法膀跌。
最后來看看示例代碼(SPARK)
SPARK 通過廣義空間線性模型直接對從各種空間解析轉(zhuǎn)錄組學(xué)技術(shù)生成的計數(shù)數(shù)據(jù)進(jìn)行建模。 它依賴于用于可擴(kuò)展計算的懲罰準(zhǔn)似然算法和最近開發(fā)的用于假設(shè)檢驗(yàn)的統(tǒng)計公式固灵，提供對 I 類錯誤的有效控制并產(chǎn)生高統(tǒng)計能力穷蛹。
圖片.png
加載
library('SPARK')
load("./Layer2_BC_Count.rds")
rawcount[1:5,1:5]
17.907x4.967 18.965x5.003 18.954x5.995 17.846x5.993 20.016x6.019
GAPDH 1 7 5 1 2
USP4 1 0 0 0 0
MAPKAPK2 1 1 0 0 1
CPEB1 0 0 0 0 0
LANCL2 0 0 0 0 0
提取每個樣本的注釋信息萧锉，即位置或坐標(biāo)
## extract the coordinates from the rawdata
info <- cbind.data.frame(x=as.numeric(sapply(strsplit(colnames(rawcount),split="x"),"[",1)),
y=as.numeric(sapply(strsplit(colnames(rawcount),split="x"),"[",2)),
total_counts=apply(rawcount,2,sum))
rownames(info) <- colnames(rawcount)
## filter genes and cells/spots and
spark <- CreateSPARKObject(counts=rawcount,
location=info[,1:2],
percentage = 0.1,
min_total_counts = 10)
## total counts for each cell/spot
spark@lib_size <- apply(spark@counts, 2, sum)
## Take the first ten genes as an example
#spark@counts <- spark@counts[1:10,]
Fit the statistical model under the null hypothesis.
## Estimating Parameter Under Null
spark <- spark.vc(spark,
covariates = NULL,
lib_size = spark@lib_size,
num_core = 5,
verbose = F)
測試空間表達(dá)的模式基因揖铜。 默認(rèn)情況下膳音，核矩陣通過坐標(biāo)自動計算，并檢查核矩陣的正定義大审。 還有一個選項(xiàng)可以由用戶提供內(nèi)核矩陣蘸际。
## Calculating pval
spark <- spark.test(spark,
check_positive = T,
verbose = F)
head(spark@res_mtest[,c("combined_pvalue","adjusted_pvalue")])
combined_pvalue adjusted_pvalue
GAPDH 7.477171e-09 4.233403e-06
MAPKAPK2 1.016078e-01 1.000000e+00
MCL1 1.149519e-08 6.078909e-06
TMEM109 4.303998e-01 1.000000e+00
TMEM189 6.189064e-01 1.000000e+00
ITPK1 7.213287e-01 1.000000e+00
再來看看示例代碼SPARK-X
SPARK-X 建立在強(qiáng)大的協(xié)方差測試框架之上，以對通過不同技術(shù)收集的各種空間轉(zhuǎn)錄組學(xué)數(shù)據(jù)進(jìn)行建模徒扶。 它依賴于可擴(kuò)展計算的代數(shù)創(chuàng)新以及新開發(fā)的用于假設(shè)檢驗(yàn)的統(tǒng)計公式粮彤，產(chǎn)生校準(zhǔn)良好的 p 值并產(chǎn)生高統(tǒng)計功效。 SPARK-X 具有很高的計算效率姜骡，并且是唯一可用于 HDST 數(shù)據(jù)的 SE 方法导坟。
圖片.png
library('SPARK')
load("./CN24_D1_unmodgtf_filtered_red_ut_HDST_final_clean.rds")
sp_count[26:30,1:5]
1000x100 1000x103 1000x113 1000x114 1000x116
Gm42418 1 . 2 . .
Gm10925 . . . . .
Gm7135 . . . . .
Atrx . . . . .
Celf2 . . . . .
View the coordinates information location, each row denotes a gene and each column represents a spot.
## extract the coordinates from the rawdata
info <- cbind.data.frame(x=as.numeric(sapply(strsplit(colnames(sp_count),split="x"),"[",1)),
y=as.numeric(sapply(strsplit(colnames(sp_count),split="x"),"[",2)))
rownames(info) <- colnames(sp_count)
location <- as.matrix(info)
mt_idx <- grep("mt-",rownames(sp_count)) ###去除線粒體基因
if(length(mt_idx)!=0){
sp_count <- sp_count[-mt_idx,]
}
####Analyze the data with SPARK-X
sparkX <- sparkx(sp_count,location,numCores=1,option="mixture")
## ===== SPARK-X INPUT INFORMATION ====
## number of total samples: 177455
## number of total genes: 19913
## Running with single core, may take some time
## Testing With Projection Kernel
## Testing With Gaussian Kernel 1
## Testing With Gaussian Kernel 2
## Testing With Gaussian Kernel 3
## Testing With Gaussian Kernel 4
## Testing With Gaussian Kernel 5
## Testing With Cosine Kernel 1
## Testing With Cosine Kernel 2
## Testing With Cosine Kernel 3
## Testing With Cosine Kernel 4
## Testing With Cosine Kernel 5
輸出結(jié)果
head(sparkX$res_mtest)
combinedPval adjustedPval
Rcn2 0.11886766 1
Mycbp2 0.03704549 1
Mprip 0.10150035 1
Mroh1 0.92949857 1
Zfp560 0.28341540 1
Cacna1e 0.08444910 1
好了，下一篇我們來對空間高變基因進(jìn)行實(shí)際運(yùn)用
生活很好圈澈，有你更好
