脚本更新----空间转录组的Radial Distance分析

作者，Evil Genius

今天我们需要更新的分析内容，空间距离分析。

空间距离分析包括空间的基因距离与细胞距离。

Radial Distance分析的基本思路：

定义中心点： 通常，空间转录组数据会选择组织切片的某个特定区域作为参考中心，可能是组织的几何中心、特定的解剖学结构位置或某个感兴趣区域（例如肿瘤中心），核心在于中心区域的无规则结构的定义。

计算各个位置的辐射距离： 一旦确定了参考点，Radial Distance分析通过计算每个空间位置到这个中心点的距离来进行。这些位置可以是Visium平台上单个探针捕获到的空间位置。

基因表达的空间分布： 结合每个空间位置的基因表达数据，可以分析特定基因的表达强度是否随着距离中心点的增加或减少而发生变化。这种分析可以帮助揭示一些潜在的空间模式，例如某些基因在远离组织中心的区域表达增高或降低。

Radial Distance分析的步骤：

数据预处理：

需要先对Visium平台生成的空间转录组数据进行标准化处理（如去噪声、批次效应校正等），确保分析的准确性。

确定参考点和距离计算：

确定一个感兴趣的参考点（例如，肿瘤组织的中心或其他解剖标志），然后计算其他每个空间位置到该点的辐射距离。

基因表达与距离的关系：

对比不同辐射距离的区域内基因的表达水平。这可能是单个基因的分析，也可以是某一组基因（例如免疫反应基因、代谢通路相关基因等）。

可视化分析：

通过热图、散点图、箱型图等可视化方式，展示基因表达与空间距离的关系。还可以进行聚类分析，揭示哪些区域基因表达模式有显著差异。

应用实例：

肿瘤微环境：研究肿瘤组织中心与边缘不同区域的基因表达变化。例如，某些免疫细胞标志基因可能在离肿瘤中心较远的区域上调，而其他类型的基因（如肿瘤相关基因）可能在中心区域表达更高。
组织发育：在胚胎发育或组织再生过程中，Radial Distance分析可以帮助识别在组织不同半径区域内的基因表达差异，揭示发育过程中不同区域的生物学差异。

首先我们用python代码实现这个内容，大家做好基础分析，下面是一个示例。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scanpy as sc
from scipy.spatial.distance import cdist
from shapely.geometry import Polygon, Point

# 1. 载入空间转录组数据
adata = sc.read('your_visium_data.h5ad')

# 假设adata.obs中包含坐标信息（例如，x和y坐标）
coordinates = adata.obs[['x', 'y']].values  # 提取x和y坐标

# 2. 定义参考区域
# 假设参考区域是一个多边形，定义一个多边形的顶点（示例区域）
polygon_points = np.array([[100, 100], [200, 100], [200, 200], [100, 200]])  # 自定义区域顶点

# 创建一个Shapely Polygon对象
polygon = Polygon(polygon_points)

# 3. 计算每个点到参考区域的最小距离
# 使用Shapely的distance方法计算每个点到多边形区域的最短距离
min_distances = []
for coord in coordinates:
    point = Point(coord)  # 创建一个Shapely Point对象
    min_distance = point.distance(polygon)  # 计算到多边形的最短距离
    min_distances.append(min_distance)

# 将最短距离加入adata.obs中
adata.obs['min_distance_to_region'] = min_distances

# 4. 选择感兴趣的基因
genes_of_interest = ['GeneA', 'GeneB', 'GeneC']  # 根据需要调整基因名
gene_expression = adata.raw[:, genes_of_interest].to_df()

# 5. 按最小距离分组并计算每组的基因表达平均值
# 假设将距离分为若干组（例如，每50μm一组）
distance_bins = np.arange(0, np.max(min_distances), 50)  # 按最小距离分组
binned_distances = np.digitize(min_distances, distance_bins)

# 将基因表达按最小距离分组，计算每组的平均表达
expression_by_distance = []
for i in range(1, len(distance_bins)):
    group_indices = np.where(binned_distances == i)[0]
    group_expression = gene_expression.iloc[group_indices].mean(axis=0)
    expression_by_distance.append(group_expression)

# 转换为DataFrame方便查看
expression_by_distance_df = pd.DataFrame(expression_by_distance, columns=genes_of_interest)
expression_by_distance_df['distance_bin'] = distance_bins[1:]

# 6. 可视化
# 绘制基因表达随最小距离变化的曲线
plt.figure(figsize=(10, 6))
for gene in genes_of_interest:
    plt.plot(expression_by_distance_df['distance_bin'], expression_by_distance_df[gene], label=gene)

plt.xlabel('Min Distance to Region (μm)')
plt.ylabel('Gene Expression (Mean)')
plt.title('Gene Expression vs Min Distance to Region')
plt.legend()
plt.grid(True)
plt.show()

解析：1、定义参考区域：

polygon_points：我们用一个简单的矩形多边形作为示例。你可以根据实际的参考区域调整这个多边形的坐标，或者使用不同的形状（例如圆形、复杂多边形等）。

使用shapely.geometry.Polygon来定义这个区域，shapely.geometry.Point用来定义每个空间位置，并计算每个点到这个区域的最短距离。

2、计算每个点到区域的最短距离：

通过point.distance(polygon)计算空间点到参考区域的最短距离。

将每个点的最短距离存储在adata.obs['min_distance_to_region']中。

3、按最小距离分组：

使用np.digitize函数按最小距离将空间位置分组（例如每50μm为一组）。

计算每组的基因表达均值。

4、可视化：

使用Matplotlib绘制基因表达随最小距离变化的曲线。每个基因的表达值会随着距离变化显示在图中。

注意事项：

参考区域形状：在实际应用中，你可以根据需要选择不同形状的区域（例如用圆形、椭圆形、多边形等）。Shapely库支持这些操作。

数据格式：adata对象假设包含坐标和基因表达数据，确保数据结构适合该代码。如果需要处理不同格式的数据，可以调整加载方式。

分组方式：distance_bins分组方式可以根据数据的特性调整，选择合适的间隔大小。

基于R的示例代码

# 安装并加载必要的R包
# install.packages(c("ggplot2", "sf", "dplyr", "Seurat"))
library(ggplot2)
library(sf)
library(dplyr)
library(Seurat)

# 1. 从RDS文件加载空间转录组数据
# 假设RDS文件包含一个Seurat对象，使用readRDS()读取数据
seurat_data <- readRDS("spatial_data.rds")

# 2. 获取坐标和基因表达数据
# 假设Seurat对象中包含空间坐标（'spatial'），并且数据存储在'RNA' assay中
coordinates <- seurat_data@images[[1]]@coordinates  # 获取坐标信息
coordinates <- as.data.frame(coordinates)  # 转换为data.frame
coordinates$x <- coordinates$`1`  # 假设x坐标在'1'列
coordinates$y <- coordinates$`2`  # 假设y坐标在'2'列

# 获取基因表达数据（假设我们关心的是原始表达数据）
gene_expression <- seurat_data@assays$RNA@counts  # 获取基因表达数据

# 将基因表达数据转换为data.frame形式
gene_expression_df <- as.data.frame(as.matrix(gene_expression))

# 3. 定义参考区域（多边形）
# 假设参考区域是一个矩形多边形（你可以根据实际情况修改为不同形状）
polygon_points <- matrix(c(100, 100, 200, 100, 200, 200, 100, 200), ncol = 2, byrow = TRUE)  # 定义一个矩形区域
polygon <- st_sfc(st_polygon(list(polygon_points)), crs = 4326)  # 创建一个SF多边形对象

# 4. 计算每个点到参考区域的最小距离
# 将每个点转换为sf对象
coordinates_sf <- st_as_sf(coordinates, coords = c("x", "y"), crs = 4326)

# 计算每个点到多边形区域的最小距离
coordinates$min_distance_to_region <- st_distance(coordinates_sf, polygon) %>% as.numeric()  # 最短距离

# 5. 按最小距离分组，并计算每组的基因表达均值
# 将最小距离分为若干组（例如，每50单位距离一组）
coordinates$distance_bin <- cut(coordinates$min_distance_to_region, breaks = seq(0, max(coordinates$min_distance_to_region), by = 50), include.lowest = TRUE)

# 计算每个距离组的平均基因表达值
expression_by_distance <- coordinates %>%
  cbind(gene_expression_df) %>%
  group_by(distance_bin) %>%
  summarise(across(everything(), mean, na.rm = TRUE))

# 6. 可视化
# 绘制基因表达随最小距离变化的趋势
ggplot(expression_by_distance, aes(x = distance_bin)) +
  geom_line(aes(y = GeneA, color = "GeneA"), size = 1) +
  geom_line(aes(y = GeneB, color = "GeneB"), size = 1) +
  geom_line(aes(y = GeneC, color = "GeneC"), size = 1) +
  labs(title = "Gene Expression vs Min Distance to Region",
       x = "Min Distance to Region (units)",
       y = "Gene Expression (Mean)") +
  scale_color_manual(values = c("GeneA" = "blue", "GeneB" = "red", "GeneC" = "green")) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))  # 旋转x轴标签

生活很好，有你更好