Nature | 5月14日发表：空间转录组文章复现（四）

天意生信云

发布于 2025-06-12 13:27:54

13600

代码可运行

运行总次数：0

代码可运行

在前三篇文章中，我们按照文章的结果顺序依次复现了Figure1、Figure2、Figure3。Figure1构建了人类胎儿大脑皮层的空间图谱，并构建了三个层级的聚类体系。Figure2识别出一组兴奋性神经元（EN）亚型，根据其分布模式分析出了大脑皮层的渐进形成过程。Figure3根据动态变化党的皮质层层级基因表达证明了神经亚型的区域特异性。基于以上的复现结果，我们今天继续完成Figure4的绘制，Figure4展示了V1与V2之间的明显转录边界。

V1与V2之间的明显转录边界

为了探究相邻皮层区域的精细分化，作者重点分析了初级视皮层（V1）和次级视皮层（V2）。

首先对GW20 样本的 数据集中特定神经元亚群（H3_annotation）进行筛选，并分析它们在不同皮层区域（Cortical area）的分布比例。

import scanpy as sc
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import openpyxl

adata = sc.read_h5ad("./gw20.h5ad")
data_want = {'H3_annotation': ['EN-IT-L3-c1', 'EN-IT-L3-c4', 'EN-IT-L3/4-c4', 
                               'EN-IT-L3/4-c1', 'EN-ET-L5/6-c5', 'EN-ET-L5/6-c4', 
                               'EN-ET-SP-2-c2', 'EN-ET-SP-2-c4']}
data_want = pd.DataFrame(data_want)
print(data_want)

filtered_df = adata.obs[(adata.obs['H2_annotation'] == 'EN-IT-L3') ]
print(filtered_df
      
# Create tuples for filtering
filtered_adata_obs = pd.merge(adata.obs, data_want, on=['H3_annotation'])
print(filtered_adata_obs)

# 1. Creating the 'sample_region' column
filtered_adata_obs['sample_region'] = filtered_adata_obs['sample'].astype(str) + '-' + filtered_adata_obs['region'].astype(str)
print(filtered_adata_obs)

meta_info = pd.read_excel("meta info.xlsx")
filtered_adata_obs = pd.merge(filtered_adata_obs, meta_info, on=['sample_region'])
print(filtered_adata_obs.head())

merged_encluster_ca_H3 = filtered_adata_obs.groupby(['Cortical area', 'H3_annotation']).size().reset_index(name='ca_count')
merged_encluster_ca_H3.to_csv("encluster_ca_H3.csv", index=False)

merged_encluster_H3 = filtered_adata_obs.groupby(['H3_annotation']).size().reset_index(name='H3_count')
merged_encluster_H3.to_csv("encluster_H3.csv", index=False)

count_ca_H3 = pd.merge(merged_encluster_ca_H3, merged_encluster_H3, on=["H3_annotation"])
count_ca_H3['percent'] = count_ca_H3['ca_count'] / count_ca_H3['H3_count']
print(count_ca_H3.head)

# Check the sum of the 'percent' column for each 'H2_annotation'
percent_sum_check = count_ca_H3.groupby(['H3_annotation'])['percent'].sum().reset_index(name='percent_sum')

# Display the sums to check if any do not equal 1
print(percent_sum_check)

count_ca_H3.to_csv("H3ca_histogram.csv")

根据输出表格，绘制堆叠柱状图，展示不同H3亚群在皮质区域（Cortical area）比例分布。

# 加载所需的R包
library(tidyverse)    # 用于数据处理和绘图（包含ggplot2）
library(readxl)       # 用于读取Excel文件（尽管本例未直接使用）
library(ggplot2)      # 用于绘图（已包含在tidyverse中，但显式加载以确保兼容）

# 读取CSV文件
count_ca_H3 <- read.csv("H3ca_histogram.csv")

# 移除第一列（假设为无关索引列）
count_ca_H3 <- count_ca_H3[, -1]

# 设置H3_annotation和Cortical.area为因子，并指定顺序
count_ca_H3$H3_annotation <- factor(
  count_ca_H3$H3_annotation,
  levels = c(
    'EN-IT-L3-c0', 'EN-IT-L3-c3', 'EN-IT-L3/4-c3', 'EN-IT-L3/4-c0',
    'EN-ET-L5/6-c4', 'EN-ET-L5/6-c3', 'EN-ET-SP-2-c1', 'EN-ET-SP-2-c3'
  )
)

count_ca_H3$Cortical.area <- factor(
  count_ca_H3$Cortical.area,
  levels = c("PFC", "PMC/M1", "Par", "Temp", "Occi")
)

# 绘制堆叠柱状图
stacked_barplot <- ggplot(count_ca_H3, aes(fill = Cortical.area, y = percent, x = H3_annotation)) +
  geom_bar(position = "fill", stat = "identity") +  # 堆叠柱状图，比例填充
  scale_fill_manual(values = c("red", "yellow", "cyan", "coral", "royalblue")) +  # 自定义颜色
  theme(
    axis.text.x = element_text(angle = 90, vjust = 0.5),  # X轴标签旋转90度
    axis.title.y = element_blank(),                       # 隐藏Y轴标题
    plot.margin = unit(c(1, 1, 1, 5), "cm"),             # 设置边距
    plot.title.position = "plot",                         # 标题位置
    plot.title = element_text(hjust = 0.4),               # 标题居中偏左
    legend.position = "right"                             # 图例置于右侧
  )

# 显示图表
print(stacked_barplot)

# 保存图表为PDF
ggsave("Cortical Area Prop.pdf", stacked_barplot, width = 5.4, height = 3.8, units = "in")

Figure4b

第20周（GW20）阶段的MERFISH分析显示，在V1和V2之间存在一个清晰的转录边界，该边界由四对相互对立的兴奋性神经元（EN）亚型组成。

import numpy as np
import scanpy as sc
import matplotlib.pyplot as plt
from multiprocessing import Pool
from matplotlib_scalebar.scalebar import ScaleBar


adata = sc.read('gw20.h5ad')

def make_plot(sample, region):
    adata1 = adata[(adata.obs['sample']==sample) & (adata.obs.region==region)].copy()
    types = ['EN-IT-L3-c1', 'EN-IT-L3/4-c4', 'EN-ET-L5/6-c5', 'EN-ET-SP-2-c2', 'EN-IT-L3-c4', 'EN-IT-L3/4-c1', 'EN-ET-L5/6-c4', 'EN-ET-SP-2-c4']
    colors = ['#FF97EE', '#FF00FF', '#BC00BC', '#A800A8', '#ACFC64', '#00FF00', '#00C800', '#007E00']
    color_dict = dict(zip(types, colors))
    unique_annotations = adata.obs['H3_annotation'].dropna().unique()
    remaining_types = np.setdiff1d(unique_annotations, types)
    color_dict.update(dict(zip(remaining_types, ['grey'] * len(remaining_types))))
    plot = sc.pl.embedding(adata1, basis="spatial", color = 'H3_annotation', groups = types, show = False, s=2, palette = color_dict, alpha=1); plt.axis('off');
    plot.set_aspect('equal');
    handles, labels = plt.gca().get_legend_handles_labels();
    order = [labels.index(i) for i in types];    
    plt.legend([handles[idx] for idx in order],[labels[idx] for idx in order], loc = 'center', fontsize=2, ncol = 2, bbox_to_anchor=(1.0,1.0), markerscale=0.25);
    plot.get_figure().gca().set_title('');
    scalebar = ScaleBar(1, "um", fixed_value=500, location = 'lower right');
    plot.add_artist(scalebar);
    plt.tight_layout();
    plt.savefig(sample + '_' + region + '.png', dpi=500); plt.clf()


sections = ['FB080-O1c', 'FB080-O1b', 'FB080-O1d', 'FB121-O1']


def main():
  with Pool(4) as pool:
    pool.starmap(make_plot, zip([i.split('-')[0] for i in sections], [i.split('-')[1] for i in sections]))

if __name__=="__main__":
    main()

Figure4c

import numpy as np
import scanpy as sc
import matplotlib.pyplot as plt
from multiprocessing import Pool
from matplotlib_scalebar.scalebar import ScaleBar
from itertools import repeat


adata = sc.read('./gw20.h5ad')

def make_plot(sample, region, i):
    adata1 = adata[(adata.obs['sample']==sample) & (adata.obs.region==region)].copy()
    types = ['EN-IT-L3-c1', 'EN-IT-L3/4-c4', 'EN-ET-L5/6-c5', 'EN-ET-SP-2-c2', 'EN-IT-L3-c4', 'EN-IT-L3/4-c1', 'EN-ET-L5/6-c4', 'EN-ET-SP-2-c4']
    colors = ['#FF97EE', '#FF00FF', '#BC00BC', '#A800A8', '#ACFC64', '#00FF00', '#00C800', '#007E00']
    types = [types[i], types[i+4]]
    colors = [colors[i], colors[i+4]]         
    color_dict = dict(zip(types, colors))
    unique_annotations = adata.obs['H3_annotation'].dropna().unique()
    remaining_types = np.setdiff1d(unique_annotations, types)
    color_dict.update(dict(zip(remaining_types, ['grey'] * len(remaining_types))))
    plot = sc.pl.embedding(adata1, basis="spatial", color = 'H3_annotation', groups = types, show = False, s=2, palette = color_dict, alpha=1); plt.axis('off');
    plot.set_aspect('equal');
    handles, labels = plt.gca().get_legend_handles_labels();
    order = [labels.index(i) for i in types];    
    plt.legend([handles[idx] for idx in order],[labels[idx] for idx in order], loc = 'center', fontsize=2, ncol = 2, bbox_to_anchor=(1.0,1.0), markerscale=0.25);
    plot.get_figure().gca().set_title('');
    scalebar = ScaleBar(1, "um", fixed_value=500, location = 'lower right');
    plot.add_artist(scalebar);
    plt.tight_layout();
    plt.savefig(sample + '_' + region + '_'+str(i)+ '.png', dpi=500); plt.clf()


samples = ['FB080-O1c', 'FB080-O1b', 'FB080-O1d', 'FB121-O1']


def main():
for sample in samples:
    with Pool(4) as pool:
      pool.starmap(make_plot, zip(repeat(sample.split('-')[0]), repeat(sample.split('-')[1]), range(4)))

if __name__=="__main__":
    main()

Figure4d

在所有皮层层级及下板层（SP）均表现出快速切换，尤以第3层和第4层的变化最为显著。

library(ggplot2)
library(reticulate)
anndata = import("anndata")

adata = anndata$read_h5ad("gw20.h5ad")
obs = adata$obs
obs = obs[obs$sample == "FB080" & obs$region == "O1c", ]
obs = obs[!is.na(obs$v1_v2_dist), ]
types = c('EN-IT-L3-c1', 'EN-IT-L3/4-c4', 'EN-ET-L5/6-c5', 'EN-ET-SP-2-c2', 'EN-IT-L3-c4', 'EN-IT-L3/4-c1', 'EN-ET-L5/6-c4', 'EN-ET-SP-2-c4')
colors = c('#FF97EE', '#FF00FF', '#BC00BC', '#A800A8', '#ACFC64', '#00FF00', '#00C800', '#007E00')
obs1 = obs[obs$H3_annotation %in% types,]
obs1$H3 = factor(obs1$H3_annotation, levels = types)
names(colors) = types
pdf(file = 'ridge.pdf', width=5.58, height=3.5)
print(ggplot(obs1, aes(x = v1_v2_dist, stat(count), color=H3_annotation, fill =H3_annotation )) + 
        geom_density(alpha=0.3)+  xlab('V1-V2 Depth') + ylab('Count') + 
        coord_flip() + scale_y_reverse() + scale_fill_manual(values=colors) + scale_color_manual(values=colors)+
        theme(legend.title = element_blank()))
dev.off()

Figure4e

对亚型中的两个簇来源于相同的H2细胞类型，且层级分布一致。

library(ggplot2)
library(reticulate)
anndata = import("anndata")

adata = anndata$read_h5ad("gw20.h5ad")
obs = adata$obs
obs = obs[obs$sample == "FB080" & obs$region == "O1c", ]
obs = obs[!is.na(obs$cortical_depth), ]
types = c('EN-IT-L3-c1', 'EN-IT-L3/4-c4', 'EN-ET-L5/6-c5', 'EN-ET-SP-2-c2', 'EN-IT-L3-c4', 'EN-IT-L3/4-c1', 'EN-ET-L5/6-c4', 'EN-ET-SP-2-c4')
obs1 = obs[obs$H3 %in% types,]
obs1$H3 = factor(obs1$H3_annotation, levels = types)
colors = c('#FF97EE', '#FF00FF', '#BC00BC', '#A800A8', '#ACFC64', '#00FF00', '#00C800', '#007E00')
names(colors) = types
pdf(file = 'ridge.pdf', width=5.58, height=3.5)
print(ggplot(obs1, aes(x = cortical_depth, stat(count), color=H3_annotation, fill =H3_annotation )) + 
        geom_density(alpha=0.3)+  xlab('Cortical Depth') + ylab('Count') + 
        coord_flip() + scale_y_reverse() + scale_fill_manual(values=colors) + scale_color_manual(values=colors)+
        geom_vline(xintercept=c(0.7354541900310434, 0.5103376408816647, 0.4110935664713694), linetype="dashed", color = "#1f77b4")+
        theme(legend.title = element_blank()))
dev.off()

Figure4f

V2富集的亚型在其他皮层区域也广泛分布，甚至在前部区域也有一定存在，这与之前描述的沿前后轴（AP轴）呈连续梯度分布的模式一致。相比之下，V1特异性的EN亚型仅存在于V1区域，其他区域均未检测到。这一结果表明，V1的区域特化可能依赖于一种高度局部化的调控程序，区别于一般的连续AP轴模式，使V1在皮层区域中特殊独立。

import scanpy as sc
import matplotlib.pyplot as plt
from multiprocessing import Pool
from matplotlib_scalebar.scalebar import ScaleBar
import matplotlib.colors
from itertools import repeat
import matplotlib


adata = sc.read('gw20.h5ad')

def make_plot(sample, region, gene, cutoff, cmap):
    adata1 = adata[(adata.obs['sample']==sample) & (adata.obs.region==region)].copy()
    sc.pp.scale(adata1, zero_center=True, max_value=cutoff)
    vmin = adata1.X.min()
    vmax = adata1.X.max()
    plot = sc.pl.embedding(adata1, basis="spatial", use_raw = False, color = gene, show = False, s=2, color_map = cmap, alpha=1, vmin=vmin, vmax=vmax, colorbar_loc = None);
    norm = matplotlib.colors.Normalize(vmin=vmin, vmax=vmax);
    sm = plt.cm.ScalarMappable(cmap=cmap, norm=norm); 
    sm.set_array([]);
    plt.colorbar(sm, location = 'top', orientation = 'horizontal', label = gene, shrink = 0.3);
    plt.axis('off');
    plot.set_aspect('equal');
    plot.get_figure().gca().set_title('');
    scalebar = ScaleBar(1, "um", fixed_value=500, location = 'lower right');
    plot.add_artist(scalebar);
    plt.tight_layout();
    plt.savefig(sample + '_' + region + '_'+gene+ '.png', dpi=500); plt.clf()


#sections = ['FB080-O1c', 'FB080-O1b', 'FB080-O1d', 'FB121-O1']

genes = ['NEFM', 'ETV1', 'UNC5C', 'PDE1A', 'NPY', 'IGFBPL1', 'TSHZ3', 'GABRA5']
cutoffs = [6,6,8,6,8,6,6,6]

cmap = ['YlGnBu']*8

def main():
    with Pool(8) as pool:
      pool.starmap(make_plot, zip(repeat('FB080'), repeat('O1c'), genes, cutoffs, cmap))

if __name__=="__main__":
    main()

Figure4h，i，j

作者还在与MERFISH分析相邻的连续切片上进行了Visium空间转录组测序，拓展空间分析至全转录组水平。结果验证了MERFISH所发现的清晰边界及基因表达模式，并揭示了更多差异表达基因。该边界在所有皮层层级及下板层形成一条直线将V1和V2分开，但内带区（IZ）、外侧和内侧次脑室区（oSVZ、iSVZ）及室管区（VZ）则表现均质。这一观察由Visium点的UMAP降维图支持，显示V1和V2形成两条平行轨迹，并在下板层及更深处分叉。

#Fig5_i-j.R
#Author : Monica D Manam
#         Walsh Lab-BCH
#         2024
#Package Credits : Seurat - Satija Lab
#Visium tutorial : https://satijalab.org/seurat/articles/spatial_vignette ####

#Import Libs
library(Seurat)
library(SeuratData)
library(ggplot2)
library(patchwork)
library(dplyr)
set.seed(1234)

#Read RDS from Zenodo Upload 
Visium_A1_brain_011124 <- readRDS("Visium_A1_brain_011124.rds")

A1_SpatialHeatMap <- SpatialFeaturePlot(Visium_A1_brain_011124, 
                                        max.cutoff = "1e+05",
                                        features = "nCount_Spatial",
                                        image.alpha = 0,
                                        pt.size.factor = 2.5) + theme(legend.position = "right")


#2.Spatial_slide
#Fig5_j
p1 <- DimPlot(Visium_A1_brain_011124, reduction = "umap", label = TRUE)
p2 <- SpatialDimPlot(Visium_A1_brain_011124, label = TRUE, label.size = 3)
A1_slide <- p1 + p2

#3.
#SpatilDimPlot
A1_spatialdimplot <- SpatialDimPlot(
  Visium_A1_brain_011124, 
  label = TRUE, 
  label.size = 6, 
  image.alpha = 0,
  pt.size.factor = 2.0,
  stroke = NA
)
#multiple
SpatialDimPlot(A1_brain, cells.highlight = CellsByIdentities(A1_brain),
               facet.highlight = TRUE,
               ncol = 5)

#Annotation : 
Visium_A1_brain_011124_annotated <- RenameIdents(Visium_A1_brain_011124, `0` = "iSVZ", `1` = "IZ", 
                                       `2` = "oSVZ-1", `3` = "VZ", `4` = "oSVZ-2",
                                       `5` = "SP-V1", `6` = "oSVZ-L1", 
                                       `7` = "oSVZ-L2", `8` = "L4-V2", 
                                       `9` = "L5/6-V1", `10` = "L4-V1", 
                                       `11` = "oSVZ-3", `12` = "L2/3-V2", 
                                       `13` = "SP-V2", `15` = "L5/6-V2", 
                                       `16` = "L2-V1", `17` = "L3-V1", 
                                       `18` = "iSVZ-L")

fig_5_i <- DimPlot(Visium_A1_brain_011124_annotated, label = TRUE)

#MARKER GENE ANALYSIS 
table(Visium_A1_brain_011124@active.ident)
Visium_A1_brain_011124_annotated.markers <- FindAllMarkers(Visium_A1_brain_011124_annotated, only.pos = TRUE)

#Filtering for markers grouped by clusters with an average log2 fold change greater than 1.
Visium_A1_brain_011124_annotated.markers %>%
  group_by(cluster) %>%
  dplyr::filter(avg_log2FC > 1) %>%
  slice_head(n = 10) %>%
  ungroup() -> top10
DoHeatmap(Visium_A1_brain_011124_annotated, features = top10$gene) + NoLegend()

table(Visium_A1_brain_011124_annotated@active.ident, Visium_A1_brain_011124_annotated@meta.data)
# Store cluster identities in object@meta.data$my.clusters
Visium_A1_brain_011124_annotated[["my.clusters"]] <- Idents(Visium_A1_brain_011124_annotated)

#CALL the function to flip the image
bb5=rotateSeuratImage(Visium_A1_brain_011124_annotated,rotation = "180")

#Plot Spatial 
fig_5_h <- SpatialFeaturePlot(object = bb5, 
                              features = c("AB13BP", "PDZRN4", "TAFA2", 
                                           "CCBE1", "THSD7B", "IL1RAP",
                                           "FLRT2"),
                              ncol = 3,
                              stroke = NA,
                              max.cutoff = "1e+05",
                              image.alpha = 0,
                              pt.size.factor = 2.5)


#Method to flip Image
# Code to Flipping Images, if the slide is flipped.
# see https://github.com/satijalab/seurat/issues/2702 for more.

# flip_angle %in% c(180, "R90", "L90", "Hf", "Vf")

rotimat=function(foo,rotation){
if(!is.matrix(foo)){
    cat("Input is not a matrix")
    return(foo)
  }
if(!(rotation %in% c("180","Hf","Vf", "R90", "L90"))){
    cat("Rotation should be either L90, R90, 180, Hf or Vf\n")
    return(foo)
  }
if(rotation == "180"){
    foo <- foo %>% 
      .[, dim(.)[2]:1] %>%
      .[dim(.)[1]:1, ]
  }
if(rotation == "Hf"){
    foo <- foo %>%
      .[, dim(.)[2]:1]
  }

if(rotation == "Vf"){
    foo <- foo %>%
      .[dim(.)[1]:1, ]
  }
if(rotation == "L90"){
    foo = t(foo)
    foo <- foo %>%
      .[dim(.)[1]:1, ]
  }
if(rotation == "R90"){
    foo = t(foo)
    foo <- foo %>%
      .[, dim(.)[2]:1]
  }
return(foo)
}

rotateSeuratImage = function(seuratVisumObject, slide = "slice1", rotation="Vf"){
if(!(rotation %in% c("180","Hf","Vf", "L90", "R90"))){
    cat("Rotation should be either 180, L90, R90, Hf or Vf\n")
    return(NULL)
  }else{
    seurat.visium = seuratVisumObject
    ori.array = (seurat.visium@images)[[slide]]@image
    img.dim = dim(ori.array)[1:2]/(seurat.visium@images)[[slide]]@scale.factors$lowres
    new.mx <- c()  
    # transform the image array
    for (rgb_idx in 1:3){
      each.mx <- ori.array[,,rgb_idx]
      each.mx.trans <- rotimat(each.mx, rotation)
      new.mx <- c(new.mx, list(each.mx.trans))
    }
    
    # construct new rgb image array
    new.X.dim <- dim(each.mx.trans)[1]
    new.Y.dim <- dim(each.mx.trans)[2]
    new.array <- array(c(new.mx[[1]],
                         new.mx[[2]],
                         new.mx[[3]]), 
                       dim = c(new.X.dim, new.Y.dim, 3))
    
    #swap old image with new image
    seurat.visium@images[[slide]]@image <- new.array
    
    ## step4: change the tissue pixel-spot index
    img.index <- (seurat.visium@images)[[slide]]@coordinates
    
    #swap index
    if(rotation == "Hf"){
      seurat.visium@images[[slide]]@coordinates$imagecol <- img.dim[2]-img.index$imagecol
    }
    
    if(rotation == "Vf"){
      seurat.visium@images[[slide]]@coordinates$imagerow <- img.dim[1]-img.index$imagerow
    }
    
    if(rotation == "180"){
      seurat.visium@images[[slide]]@coordinates$imagerow <- img.dim[1]-img.index$imagerow
      seurat.visium@images[[slide]]@coordinates$imagecol <- img.dim[2]-img.index$imagecol
    }
    
    if(rotation == "L90"){
      seurat.visium@images[[slide]]@coordinates$imagerow <- img.dim[2]-img.index$imagecol
      seurat.visium@images[[slide]]@coordinates$imagecol <- img.index$imagerow
    }
    
    if(rotation == "R90"){
      seurat.visium@images[[slide]]@coordinates$imagerow <- img.index$imagecol
      seurat.visium@images[[slide]]@coordinates$imagecol <- img.dim[1]-img.index$imagerow
    }
    
    return(seurat.visium)
  }  
}