hello，大家好，又是周五，一周又將過去了，時(shí)間總是匆匆，片刻不停留，我們唯一能做的，就是在時(shí)間流逝的過程中，留下點(diǎn)印記，雖然我們窮，找不下對(duì)象，但是過的開心也是成功~~~~
今天我們來學(xué)習(xí)一個(gè)很好的分析點(diǎn)，空間基因表達(dá)模式與共定位（共定位可以是細(xì)胞類型的共定位，也可以是配受體的共定位），參考的文章是Scalable and model-free detection of spatial patterns and colocalization,方法非常棒，推薦給大家。
空間組學(xué)分析的一個(gè)重要步驟是表征空間表達(dá)模式和共定位。已經(jīng)開發(fā)了幾種方法來識(shí)別空間可變基因。Trendsceek 使用permutation test來檢測點(diǎn)的空間分布與其基于標(biāo)記點(diǎn)過程的表達(dá)水平之間的顯著依賴性。Sepal 通過擴(kuò)散時(shí)間對(duì)空間可變基因進(jìn)行排序，理由是具有空間模式的基因比具有隨機(jī)空間分布的基因需要更多時(shí)間才能達(dá)到同質(zhì)狀態(tài)。SpatialDE 和 SPARK 都利用高斯過程回歸作為空間協(xié)方差結(jié)構(gòu)的基礎(chǔ)數(shù)據(jù)生成模型。SpatialDE 將表達(dá)變異性分解為空間方差和噪聲，并通過比較有和沒有空間分量的可能性來估計(jì)統(tǒng)計(jì)顯著性。SPARK 通過廣義線性空間誤差模型擴(kuò)展 SpatialDE，能夠直接模擬原始計(jì)數(shù)并調(diào)整協(xié)變量。SPARK-X 檢查表達(dá)式協(xié)方差矩陣和距離協(xié)方差矩陣的相似性，并測試它們是否比偶然預(yù)期的更相似。這些方法的統(tǒng)計(jì)能力高度依賴于空間協(xié)方差模型，即它們與真實(shí)的潛在表達(dá)模式的匹配程度。盡管考慮了多個(gè)內(nèi)核，包括具有不同平滑度參數(shù)的高斯、線性和周期性內(nèi)核，以確保識(shí)別各種空間模式，但統(tǒng)計(jì)能力將大大降低，以識(shí)別由那些預(yù)定義的內(nèi)核函數(shù)建模不佳的空間模式。此外，空間協(xié)方差模型是建立在細(xì)胞距離上的，這會(huì)將真實(shí)的表達(dá)差異與由細(xì)胞密度差異驅(qū)動(dòng)的表達(dá)差異混淆?？紤]到細(xì)胞密度不均勻，MERINGUE 基于空間鄰域圖計(jì)算空間自相關(guān)和互相關(guān)，以識(shí)別空間可變基因和基因相互作用。
章節(jié)1：空間基因表達(dá)模式
library(devtools)
install_github("liuqivandy/SpaGene")
Load the library and the data
library(SpaGene)
load("BreastCancer/bc_raw.rds")###Seurat分析出來的rds
Find spatially variable genes and patterns
brain10x_sv<-SpaGene(count,location)
# the most significant spt
head(brain10x_sv$spagene_res[order(brain10x_sv$spagene_res$adjp),])
## score zval pval adjp
## COL12A1 8.48 -16.17570 3.741996e-59 5.534038e-55
## FN1 8.52 -16.06389 2.284824e-58 1.689513e-54
## COL3A1 8.72 -15.50484 1.608692e-54 7.930316e-51
## COL1A2 8.88 -15.05759 1.538985e-51 5.690014e-48
## B2M 9.16 -14.27491 1.568364e-46 3.865755e-43
## COL1A1 9.16 -14.27491 1.568364e-46 3.865755e-43
空間表達(dá)模式繪圖
pattern<-FindPattern(bc_spagene)
PlotPattern(pattern,location)
Top five genes falling into each pattern
Top five genes falling into each pattern
top5<-apply(pattern$genepattern,2,function(x){names(x)[order(x,decreasing=T)][1:5]})
library(pheatmap)
pheatmap(pattern$genepattern[rownames(pattern$genepattern)%in%top5,])
章節(jié)2 ：Identification of spatially colocalized ligand-receptor pairs
SpaGene 以識(shí)別由共定位配體和受體對(duì)介導(dǎo)的細(xì)胞間通訊。SpaGene 通過空間轉(zhuǎn)錄組學(xué)從示例數(shù)據(jù)中發(fā)現(xiàn)了 35 種配體-受體相互作用。兩個(gè)最重要的配體-受體對(duì)是 Igfbp5-Cav1 (adjp=3e-31) 和 Apoe-Lrp6 (adjp=1e-18)，兩者都發(fā)生在 ONL 和 GL 之間。Apoe 已知富含 ONL 和 GL，并且被 SpaGene 鑒定為非常重要（adjp=1e-50）。大多數(shù)具有高 Apoe 表達(dá)的斑點(diǎn)都被具有高 Lrp6 表達(dá)的斑點(diǎn)包圍，表明它們之間存在潛在的相互作用。Apoe-Lrp6 介導(dǎo) Wnt 信號(hào)傳導(dǎo)，這對(duì)于調(diào)節(jié)突觸完整性和認(rèn)知很重要。在 ONL 和 GL 層之間鑒定 Apoe-Lrp6 可能暗示 Wnt 信號(hào)傳導(dǎo)在建立外圍-CNS 嗅覺連接中的潛在調(diào)節(jié)。
Identify colocalized ligand-receptor pairs
load("LRpair_human.rds")
bc_lr<-SpaGene_LR(count,location,LRpair=LRpair)
# the most signficant colocalized LR pairs
head(bc_lr[order(bc_lr$adj),])
## score comm zval pval adjp
## COL1A2_ITGA11 28.36 18 -10.176989 1.256075e-24 2.348861e-21
## FN1_SDC2 28.68 19 -9.034053 8.271455e-20 7.733810e-17
## COL1A1_ITGA11 28.88 16 -8.319718 4.408658e-17 2.748063e-14
## FN1_CD44 28.92 18 -8.176851 1.456795e-16 6.810516e-14
## COL1A2_ITGB1 29.12 20 -7.462516 4.244301e-14 1.587369e-11
## COL1A1_ITGB1 29.16 15 -7.319648 1.243107e-13 3.320871e-11
Plot the ligand-receptor pair FN1-CD44
plotLR(count,location,LRpair=c("FN1","CD44"),alpha.min=0.5)
最后，我們以小鼠腦的數(shù)據(jù)為例分析一下空間模式與共定位
library(SpaGene)
load("brain10X/brain10x.rds")
brain10x_sv<-SpaGene(count,location)
# the most significant spt
head(brain10x_sv$spagene_res[order(brain10x_sv$spagene_res$adjp),])
## score zval pval adjp
## 3110035E14Rik 7.838590 -58.63093 0 0
## Igfbp2 9.747681 -41.03227 0 0
## Scg2 8.992579 -47.99306 0 0
## Ptma 9.705009 -41.42563 0 0
## Ngef 8.725417 -50.45585 0 0
## Dbi 7.745826 -59.48606 0 0
pattern<-FindPattern(brain10x_sv,nPattern = 15)
PlotPattern(pattern,location,pt.size = 0.5)
top5<-apply(pattern$genepattern,2,function(x){names(x)[order(x,decreasing=T)][1:5]})
library(pheatmap)
pheatmap(pattern$genepattern[rownames(pattern$genepattern)%in%top5,],fontsize_row = 6)
Identify colocalized ligand-receptor pairs
load("LRpair.rds")
brain10x_lr<-SpaGene_LR(count,location,LRpair=LRpair)
# the most signficant colocalized LR pairs
head(brain10x_lr[order(brain10x_lr$adj),]
## score comm zval pval adjp
## Cck_Cckbr 28.63822 214 -32.06833 6.095309e-226 1.075212e-222
## Apoe_Sdc4 28.74954 242 -30.67785 5.620267e-207 4.957076e-204
## Pdyn_Oprk1 28.85343 201 -29.38006 4.936498e-190 2.902661e-187
## Apoe_Abca1 28.95733 210 -28.08228 8.062429e-174 3.555531e-171
## Tac1_Tacr1 29.13173 204 -25.90386 3.012735e-148 1.062893e-145
## Nptx1_Nptxr 29.15028 204 -25.67211 1.197629e-145 3.521030e-143
plotLR(count,location,LRpair=c("Cck","Cckbr"),alpha.min=0.2,pt.size = 1)
生活很好，有你更好
