应用空间统计学分析空间表达数据
空间信息在空间转录组中的运用
Giotto|| 空间表达数据分析工具箱
SPOTlight || 用NMF解卷积空间表达数据
空间转录组教程|| stLearn ：空间轨迹推断
Seurat 新版教程：分析空间转录组数据
单细胞转录组数据分析|| scanpy教程：空间转录组数据分析
10X Visium:空间转录组样本制备到数据分析
定量免疫浸润在单细胞研究中的应用
在之前的文章中，我们提出地理学三大定律是完全适用于空间表达数据的。分析空间表达数据，如果离开空间信息，只用其表达矩阵那么单细胞的所有分析点当然是完全能跑得通的，但是有两点我们需要追问：
- 这样做的生物学意义是什么
- 既然你忽视了空间数据，为什么要做空间表达，而不是只做表达
这两个问题值得我们思考，也值得我们做一些探索：把空间信息纳入到分析中来。这里我们再一次思考空间信息所带来的新的可能。首先，我们来熟悉一下空间表达数据中包含的数据类型：
我们看到图象/空间/表达这三种数据类型都可以对应到矩阵的形式上，也就是在这里我们面对的是三个矩阵。既可以对他们分别做分析，也可以将他们联系在一起分析。结合空间数据当然是我们喜闻乐见的了，但是我们先来看看图象数据的分析。这一块，我推荐阅读公号【单细胞组学】的一篇文章：H&E染色图像分割与空间转录组共定位揭示肿瘤组织异质性，在这篇文章中作者演示了如何利用H&E染色图象数据来辅助聚类。为了致敬原作，文章部分内容我们将在下面复现
- H&E图像辅助细胞聚类
首先，我们配置好演示用的R包和数据。
library(Seurat)
library(SeuratData)
library(cowplot)
library(ggplot2)
library(sp)
library(grid)
library(raster)
library(spatstat)
library('SPARK')
library(spdep)
library(GWmodel)
AvailableData()
# InstallData('stxBrain')
stxBrain.SeuratData::anterior1 -> brain
brain
An object of class Seurat
31053 features across 2696 samples within 1 assay
Active assay: Spatial (31053 features, 0 variable features)
head(brain@meta.data)
orig.ident nCount_Spatial nFeature_Spatial slice region
AAACAAGTATCTCCCA-1 anterior1 13069 4242 1 anterior
AAACACCAATAACTGC-1 anterior1 37448 7860 1 anterior
AAACAGAGCGACTCCT-1 anterior1 28475 6332 1 anterior
AAACAGCTTTCAGAAG-1 anterior1 39718 7957 1 anterior
AAACAGGGTCTATATT-1 anterior1 33392 7791 1 anterior
AAACATGGTGAGAGGA-1 anterior1 20955 6291 1 anterior
查看Seurat图象数据的结构：
其中image中是图象RGB三原色信息，这个三维的矩阵。
str(brain@images$anterior1@image)
num [1:599, 1:600, 1:3] 0.722 0.722 0.722 0.718 0.718 ...
我们可以直接将其画出来：
grid.raster(brain@images$anterior1@image)
我们看到染色是有颜色深浅的，那么就可以用RGB三原色信息的矩阵数据来做聚类，也就是识别出着色的不同模式。
# 提取矩阵信息
img <- brain@images$anterior1@image
imgDm <- dim(img)
# Assign RGB channels to data frame
imgRGB <- data.frame(
x = rep(1:imgDm[2], each = imgDm[1]),
y = rep(imgDm[1]:1, imgDm[2]),
R = as.vector(img[,,1]),
G = as.vector(img[,,2]),
B = as.vector(img[,,3])
)
#kmeans 聚类 ，当然你可以用其他的聚类算法
kClusters <- 6
kMeans <- kmeans(imgRGB[, c("R", "G", "B")], centers = kClusters)
#col2hex(brain@images$anterior1@image)
kColours <- rgb(kMeans$centers[kMeans$cluster,])
unicolor <- unique(kColours)
map(1:length(unicolor),function(i){
binary_col <- kColours
#set other clusters as background color.
binary_col[kColours %in% unicolor[-i]] <- "#FAF9F9"
p <- ggplot(data = imgRGB, aes(x = x, y = y)) +
geom_point(colour = binary_col) +
labs(title = paste("k-Means Clus: ",i)) +
xlab("x") +
ylab("y")
p
}) %>% cowplot::plot_grid(plotlist = .)
仅用图象数据也可以完成聚类，只是这里的聚类没有对应到spot的barcode上，但是不同部位的细胞对苏木精和伊红的感知不同，那么也会聚成不同的类。
所以有了染色数据也是可以用来做数据探索的。
- 空间变异基因
至于空间信息，我们不妨理解为采样点，不同分辨率的技术，理解为一个采样点采到了多少个细胞。
coordinates <- GetTissueCoordinates(object = brain)
ggplot(coordinates,aes(imagecol,imagerow)) + geom_point()+ theme_bw()
所谓空间变异基因，是指在空间中分布是由显著差异的基因，举个例子就明白了：
head(brain@images$anterior1@coordinates[,c(2:3)])
row col
AAACAAGTATCTCCCA-1 50 102
AAACACCAATAACTGC-1 59 19
AAACAGAGCGACTCCT-1 14 94
AAACAGCTTTCAGAAG-1 43 9
AAACAGGGTCTATATT-1 47 13
AAACATGGTGAGAGGA-1 62 0
> brain
An object of class Seurat
31053 features across 2696 samples within 1 assay
Active assay: Spatial (31053 features, 0 variable features)
indf <-brain@images$anterior1@coordinates[,c(2:3)]
brain_sparkX <- sparkx(brain@assays$Spatial@counts,as.matrix(indf),
numCores=1,option="mixture")
brain_sparkX$stats %>% as.data.frame() %>% arrange(projection) %>% head(4)%>% rownames() -> mdown
brain_sparkX$stats %>% as.data.frame() %>% arrange(desc(projection)) %>% head(3) %>% rownames() -> mtop
dim(brain_sparkX$stats)
SpatialFeaturePlot(brain, features = c(mdown[2:4],mtop))
我们看到上面一行的三个基因就没有多少空间特异性，而下面的就很明显。显然这对我们寻在某位置的基因表达模式是很有帮助的，加以扩展，如何看一个基因集的空间模式，这一模式对应的生物学意义是什么？ 这里我们用的方法是广义空间线性模型（generalized spatial linear models,GSLM），这一方法在单细胞转录组中的应用被封装在R包SPARK中，文章见：
Shiquan Sun, Jiaqiang Zhu and Xiang Zhou#. Statistical analysis of spatial expression pattern for spatially resolved transcriptomic studies, 2020, Nature Methods, in press.
Jiaqiang Zhu, Shiquan Sun and Xiang Zhou#. Scalable and robust non-parametric detection of spatial expression patterns for large spatial transcriptomic studies., 2020
- 地理加权回归
地理加权回归（Geographically weighted regression, GWR)是一种空间分析技术，广泛应用于地理学及涉及空间模式分析的相关学科。GWR通过建立空间范围内每个点处的局部回归方程，来探索研究对象在某一尺度下的空间变化及相关驱动因素，并可用于对未来结果的预测。由于它考虑到了空间对象的局部效应，因此其优势是具有更高的准确性。【来自百科】
在上一步中，我们检测到了空间特异的基因，那么这些基因之间或他们与空间中某一事件的关系是怎样的，这时候我们往往会想到用回归的方法。在空间分析中，观测数据一般按照给定的地理位置作为采样单元进行采样，随着地理位置的变化，变量间的关系或者结构会发生改变，即GIS中所说的“空间非平稳性”。这种空间非平稳性普遍存在于空间数据中，如不同省份的AIDS发病率、湖泊不同深度TN含量、城市工业区与非工业区PM2.5浓度等等。如果采用传统的线性回归模型来分析空间数据，一般很难得到令人满意的结果，因为全局模型在分析前就假定了变量间的关系具有“各向同性”，所得结果只是研究区域内的某种“平均”。因此，有必要采用一种新的局部回归方法，来应对空间数据自身的这种属性。GWR模型便顺势被研究者提出并加以大量实践和验证【同样，来自百科】。总而言之，空间信息很重要。
我们知道单细胞数据分析过程中，落脚点往往是在某个基因集上面，这里我们也选一个基因集来做地理加权回归。
brain <- FindVariableFeatures(brain)
InDevar <- VariableFeatures(brain)[1:10]
Devar= "Hpca"
InDevar <- setdiff(InDevar,c(Devar,"1500015O10Rik","2900040C04Rik","Hba-a1","Hba-a2","Hbb-bs"))
featherdf <- Seurat::FetchData(brain, vars = c("imagerow", "imagecol",Devar, InDevar))
head(featherdf)
imagerow imagecol Hpca Ttr S100a5 Enpp2 Ptgds Doc2g Kl Sst Cdhr1 Fabp7
AAACAAGTATCTCCCA-1 386 439.0 0 2 0 3 16 1 0 0 0 0
AAACACCAATAACTGC-1 442 144.0 1 3 27 3 13 25 0 2 40 24
AAACAGAGCGACTCCT-1 163 410.5 19 1 0 11 25 0 0 3 0 0
AAACAGCTTTCAGAAG-1 343 108.4 3 2 0 2 19 3 0 1 4 17
AAACAGGGTCTATATT-1 367 122.6 1 2 0 3 8 2 0 0 2 7
AAACATGGTGAGAGGA-1 460 76.4 0 4 53 3 5 6 0 0 4 71
构建空间数据对象。这里用的是R包sp, 在sp中定义了一个空间对象基础类Spatial，由两个solt 构成：bbox和proj4string, 在Spatial类的基础上，分别扩展为点线面和栅格4种空间数据类型，分别为SpatialPoints/ SpatialLines/SpatialPolygons/SpatialGirds。这里我们主要用的是SpatialPoints数据类型，当然如果分完群之后，形成了基本格局和分区是可以构建后面的几种数据类型的。
Sptm <- SpatialPoints(featherdf)
Sptm<- SpatialPointsDataFrame(Sptm ,featherdf)
Sptm
class : SpatialPointsDataFrame
features : 2696
extent : 139, 541, 76.4, 492 (xmin, xmax, ymin, ymax)
crs : NA
variables : 12
names : imagerow, imagecol, Hpca, Ttr, S100a5, Enpp2, Ptgds, Doc2g, Kl, Sst, Cdhr1, Fabp7
min values : 138.640278405, 76.41996724, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
max values : 540.567997997, 492.237532229, 269, 12442, 88, 602, 1707, 85, 41, 305, 57, 117
plot(Sptm)
大头儿子思考者
?model.selection.gwr
model.sel <- model.selection.gwr(Devar,InDevar,data =Sptm,kernel = "gaussian" ,adaptive = T,bw=1000)
# str(model.sel)
sorted.m <- model.sort.gwr(model.sel,numVars = length(InDevar),ruler.vector = model.sel[[2]][,2])
# sorted.m
model.l <- sorted.m[[1]]
?model.view.gwr
model.view.gwr(Devar,InDevar,model.list =model.l )
以上执行了地理加权回归的模型选择（也可以叫做变量选择）过程，就是哪些变量是值得纳入模型中的，同时给出了一个判断依据。
和
plot(sorted.m[[2]][,2],col="black",pch=20,lty=5,type="b")
可以看出曲线大概在40的时候是平稳的，对应上图的一个组合： Hpca~Ttr+Ptgds+Fabp7+Enpp2+Cdhr1+Sst+Kl+S100a5。这个回归公式，我们在教科书上就见过吧。
选定模型后，选择带宽（Bandwidth）。
Sptm@coords <- Sptm@coords[,1:2] # 前面的建立SpatialPointsDataFrame 需要调整
bw.gwr <- bw.gwr(Hpca ~ Ttr+ S100a5+ Enpp2 +Ptgds+ Doc2g +Kl +Sst+ Cdhr1 ,
data = Sptm,approach = "CV",kernel = "gaussian" ,adaptive = T )
下面执行地理加权回归
gwr.res1 <- gwr.basic(Hpca ~ Ttr+Ptgds+Fabp7+Enpp2+Cdhr1+Sst+Kl+S100a5 , data = Sptm,bw = bw.gwr,kernel = "gaussian" ,adaptive = T)
gwr.res1
***********************************************************************
* Package GWmodel *
***********************************************************************
Program starts at: 2020-12-23 00:01:48
Call:
gwr.basic(formula = Hpca ~ Ttr + Ptgds + Fabp7 + Enpp2 + Cdhr1 +
Sst + Kl + S100a5, data = Sptm, bw = bw.gwr, kernel = "gaussian",
adaptive = T)
Dependent (y) variable: Hpca
Independent variables: Ttr Ptgds Fabp7 Enpp2 Cdhr1 Sst Kl S100a5
Number of data points: 2696
***********************************************************************
* Results of Global Regression *
***********************************************************************
Call:
lm(formula = formula, data = data)
Residuals:
Min 1Q Median 3Q Max
-19.29 -9.60 -3.70 4.45 254.43
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 14.75829 0.47249 31.24 < 2e-16 ***
Ttr -0.00128 0.00402 -0.32 0.750
Ptgds -0.01903 0.00423 -4.50 7.1e-06 ***
Fabp7 -0.26754 0.05811 -4.60 4.3e-06 ***
Enpp2 -0.01390 0.07748 -0.18 0.858
Cdhr1 -0.47453 0.07557 -6.28 4.0e-10 ***
Sst 0.03153 0.01326 2.38 0.018 *
Kl 0.44907 0.58133 0.77 0.440
S100a5 0.00305 0.06787 0.04 0.964
---Significance stars
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 16.1 on 2687 degrees of freedom
Multiple R-squared: 0.0584
Adjusted R-squared: 0.0556
F-statistic: 20.8 on 8 and 2687 DF, p-value: <2e-16
***Extra Diagnostic information
Residual sum of squares: 696982
Sigma(hat): 16.1
AIC: 22647
AICc: 22647
***********************************************************************
* Results of Geographically Weighted Regression *
***********************************************************************
*********************Model calibration information*********************
Kernel function: gaussian
Adaptive bandwidth: 17 (number of nearest neighbours)
Regression points: the same locations as observations are used.
Distance metric: Euclidean distance metric is used.
****************Summary of GWR coefficient estimates:******************
Min. 1st Qu. Median 3rd Qu. Max.
Intercept -2.60e-01 2.97e+00 9.86e+00 1.55e+01 62.74
Ttr -9.08e-01 2.16e-03 1.73e-01 7.72e-01 8.40
Ptgds -8.63e-01 -8.71e-02 -2.53e-02 -7.34e-05 1.19
Fabp7 -3.71e+00 -1.04e-01 1.35e-01 6.52e-01 13.28
Enpp2 -7.75e+00 -6.78e-02 1.77e-01 6.30e-01 2.39
Cdhr1 -7.90e+00 -5.44e-01 8.22e-02 2.55e+00 90.57
Sst -1.57e+00 -4.05e-02 1.10e-02 4.70e-02 1.83
Kl -2.94e+01 -2.08e+00 1.85e-01 1.93e+00 65.92
S100a5 -4.55e+01 -1.52e+00 -1.29e-03 2.75e+00 140.35
************************Diagnostic information*************************
Number of data points: 2696
Effective number of parameters (2trace(S) - trace(S'S)): 766
Effective degrees of freedom (n-2trace(S) + trace(S'S)): 1930
AICc (GWR book, Fotheringham, et al. 2002, p. 61, eq 2.33): 21251
AIC (GWR book, Fotheringham, et al. 2002,GWR p. 96, eq. 4.22): 20395
Residual sum of squares: 247514
R-square value: 0.666
Adjusted R-square value: 0.533
***********************************************************************
Program stops at: 2020-12-23 00:01:51
地理统计的基本信息，我们看到有几个基因在模型中是比较显著的，画一下看看。
genes <- c("Hpca" ,"Cdhr1", "Fabp7", "Ptgds")
SpatialFeaturePlot(brain, features = genes)
我们还可以在空间位置上细化每一个变量的参数估计和诊断信息。
spplot(gwr.res1$SDF,key.space="right")
借助地理加权回归，我们可以感受到基因的表达是有空间异质性的。那么，之前直接拿一堆基因做的相关性网络就显得略失稳妥了。
空间统计已经形成一个独立的学科门类，空间表达数据如能利用空间统计的基本概念与模型，一定会带来新的角度。本文只是做了一些简单探索，甚至空间统计的许多基础概念都只是一笔带过，算是抛砖引玉吧。当然，生搬硬套模型也会贻笑大方。应用模型的标准不是代码跑不跑得通，而是该模型能给我们带来怎样的神奇体验。我相信每一门学科都是人类的一双眼睛，让我们得以看见这平凡世界的离奇的美。
https://xzhoulab.github.io/SPARK/01_about/
https://www.biorxiv.org/content/10.1101/810903v1.full
http://spatialgiotto.rc.fas.harvard.edu/giotto_spatial_genes.html
https://zhuanlan.zhihu.com/p/90627938
R语言空间数据处理与分析实践教程，卢宾宾编著
