📚 课程目录第3章数据获取与格式

第三章：数据获取与格式——从公共数据库到本地加载

做空间转录组实战，第一步就是搞到数据。这事听起来简单，但我在辅导学员时发现，光是下载和读取就能卡住一大半人。就因为它不是单一格式，不同平台、不同版本、不同组织来源，文件长相千奇百怪。这一章我就带你把它彻底捋清楚。

核心要点：空间转录组数据 = 表达矩阵 + 空间坐标 + 图像（可选）。无论从哪个数据库下载，最后读进 Python 时你都要看清这三样东西有没有齐全。

3.1 公共数据库怎么挑？

目前主流的空间转录组数据仓库有两个：GEO（Gene Expression Omnibus）和 Single Cell Portal。GEO 是老牌选手，大部分10x Visium、Slide-seq、MERFISH 的原始数据都往这存。Single Cell Portal 则是 Broad 研究所的亲儿子，专注单细胞和空间数据，交互做得更友好。

我个人习惯先查 Single Cell Portal，因为它会把数据整理成统一的 H5AD 格式，下载后几乎不用改代码就能用。但如果项目很新，比如今年刚发的空间数据，往往只有 GEO 才有。这时候你得做好手动整理元数据的准备。

数据库	网址	文件特点	适用场景
GEO	www.ncbi.nlm.nih.gov/geo/	通常提供原始 FASTQ、处理后的 MTX 或 H5 文件；图像可能以 supplemental 形式提供	经典数据集、新发表数据
Single Cell Portal	singlecell.broadinstitute.org	统一 H5AD，附带 tsne/umap 坐标和 metadata	快速验证、教学数据集

小技巧：在 GEO 搜索时，关键词加上 "spatial transcriptomics"、"Visium" 或捕获芯片型号，比搜 "sequencing" 更精准。我踩过多次坑，搜出来的结果一半是单细胞。

3.2 下载数据——别急着点“Download”

很多新手一进 GEO 页面就直接点 Download 超链接，结果下回来一堆 SRA 文件，还得用 fastq-dump 转换，非常麻烦。我建议你优先找 Supplementary file 区域的 .h5 或 .mtx.gz。这些才是分析可以直接用的计数矩阵。

举个例子，10x Visium 的标准数据包里通常包含以下文件：

sample/
├── filtered_feature_bc_matrix.h5      # 官方推荐，包含基因-细胞矩阵
├── filtered_feature_bc_matrix/
│   ├── barcodes.tsv.gz                # 细胞条形码列表
│   ├── features.tsv.gz                # 基因特征列表
│   └── matrix.mtx.gz                  # 稀疏矩阵（三列格式）
├── spatial/
│   ├── tissue_positions_list.csv       # 空间坐标（barcode, x, y）
│   ├── scalefactors_json.json          # 图像缩放参数
│   └── tissue_hires_image.png          # 组织切片高分辨率图
└── raw_feature_bc_matrix.h5            # 未过滤的原始矩阵

我最常用的是 filtered_feature_bc_matrix.h5，它包含了经过空点过滤后的所有 barcode 和基因计数。但注意，有些课题组上传的是 raw_feature_bc_matrix.h5，你需要自己再过滤一遍空点——我在项目中遇到过这种“好心坑”，后来我每次下完数据都会先看一眼 H5 文件里的 cell barcode 数量来确认。

避坑指南：我曾经在 GEO 下了一个 MERFISH 数据集，满心欢喜地加载后发现矩阵是 2D 连续坐标但图像文件只有 thumbnail 分辨率。后来才读到 supplementary 说明里那句“high-res image available upon request”。所以下载前一定通读 README 或 metadata 表格，确认图像分辨率满足你的分析需求。

3.3 数据格式解析——三件套：H5AD、h5、mtx

空间转录组最常用的三种格式，我把它们拆开讲清楚。

3.4.1 H5AD（H5 封装的 AnnData）

这是 Scanpy 生态的默认格式。.h5ad 本质上是一个 HDF5 容器，里面装了四样东西：X（表达矩阵）、obs（细胞注释）、var（基因注释）、obsm（降维结果和空间坐标）。读起来只需一行代码：

import anndata
adata = anndata.read_h5ad("sample.h5ad")
print(adata)

输出示例：

AnnData object with n_obs × n_vars = 4720 × 18000
    obs: 'library_id', 'n_counts', 'n_genes', 'cluster'
    var: 'gene_ids', 'feature_types', 'genome'
    obsm: 'spatial', 'X_umap'

我建议你拿到 H5AD 后先检查 adata.obsm['spatial'] 的形状——它应该是 (n_obs, 2) 的二维坐标。如果只有一个维度，那很可能是被错误地存成了其他 key。

3.4.2 H5 格式（10x 官方的 HDF5）

这是 10x Genomics 原生的 .h5 文件，不同于 H5AD。它内部结构更接近“基因-细胞对”的压缩存储。用 Scanpy 的 read_10x_h5 函数可以直接读：

import scanpy as sc
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
adata.var_names_make_unique()
print(adata)

但这里有个坑：很多 H5 文件里的基因名不保证唯一（比如多个 ENST 对应同一 gene symbol）。我习惯在 read_10x_h5 后立刻跑 adata.var_names_make_unique()，否则后续合并多个样本时一定会报错。

3.4.3 MTX 格式（Market Exchange Format）

如果你下载的是 matrix.mtx.gz + barcodes.tsv.gz + features.tsv.gz 三件套，那就是 MTX 格式。它把表达矩阵写成稀疏格式，存储效率极高。读取方式：

import scanpy as sc
adata = sc.read_10x_mtx("filtered_feature_bc_matrix/")
print(adata)

MTX 的优点是非常轻量，适合在不稳定网络下断点续传。缺点是你得手动把 spatial 文件夹里的坐标数据拼接到 adata 上。比如：

import pandas as pd
pos = pd.read_csv("spatial/tissue_positions_list.csv",
                  header=None,
                  names=["barcode", "in_tissue", "row", "col", "pxl_row", "pxl_col"])
pos = pos[pos["in_tissue"] == 1]
pos = pos.set_index("barcode")
adata.obs = adata.obs.join(pos[["pxl_col", "pxl_row"]])
adata.obsm["spatial"] = adata.obs[["pxl_col", "pxl_row"]].values

嗯，这是比较繁琐的一步，但我曾经因为这个拼接粗心导致坐标轴顺序搞反，整个下游分析和可视化全反了。现在我会在拼接后先画一个 sc.pl.spatial(adata, color="n_counts") 来验证切片方向。

3.5 图像文件读取

空间转录组离不开高分辨率组织图像。10x Visium 提供 tissue_hires_image.png 和 tissue_lowres_image.png。加载时你需要把它接到 adata.uns["spatial"] 里。Scanpy 有封装好的函数帮你做：

import scanpy as sc
from scanpy import read
import matplotlib.pyplot as plt

adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
img = plt.imread("spatial/tissue_hires_image.png")

# 手动设置 spatial 库的 uns 结构
adata.uns["spatial"] = {
    "library_id": {
        "scalefactors": {
            "tissue_hires_scalef": 1.0,  # 通常从 scalefactors_json.json 读取
            "spot_diameter_fullres": 50.0
        },
        "images": {
            "hires": img
        }
    }
}
sc.pl.spatial(adata, img_key="hires", color="n_counts")

这里要注意 tissue_hires_scalef 的值必须跟 scalefactors_json.json 匹配，否则坐标会漂移。我每次都会打印 scalefactors_json.json 的内容检查一遍：

import json
with open("spatial/scalefactors_json.json") as f:
    sf = json.load(f)
print(sf)
# 输出类似：
# {'tissue_hires_scalef': 0.17012642, 'spot_diameter_fullres': 53.946408, ...}

3.6 知识体系一览

我把本章的核心内容提炼成一张图，方便你记忆整个流程。

3.7 实战建议：三步拿到可分析数据

确定数据库和数据集 ID：如果是重复我讲过的案例，直接用 sc.datasets.visium_sge(sample_id) 一行搞定。但实际操作中多数是自定义数据集。
下载并验证文件完整性：用 md5sum 检查 H5/MTX 文件的哈希值，跟网站提供的对比。我之前有次因为断网导致 H5 文件损坏，加载时直接 segmentation fault。
读入并检查 adata：打印一下 adata.obs.shape 和 adata.obsm 的 keys，确保 spatial 存在且长度匹配。不匹配的话回头查坐标文件的行数。

我的主力代码片段：

def load_visium_sample(path):
    adata = sc.read_10x_h5(path + "/filtered_feature_bc_matrix.h5")
    adata.var_names_make_unique()
    import json
    with open(path + "/spatial/scalefactors_json.json") as f:
        sf = json.load(f)
    img = plt.imread(path + "/spatial/tissue_hires_image.png")
    adata.uns["spatial"] = {
        "sample": {
            "scalefactors": sf,
            "images": {"hires": img}
        }
    }
    pos = pd.read_csv(path + "/spatial/tissue_positions_list.csv", header=None,
                      names=["barcode","in_tissue","row","col","pxl_row","pxl_col"])
    pos = pos[pos["in_tissue"]==1].set_index("barcode")
    adata.obs = adata.obs.join(pos[["pxl_col","pxl_row"]])
    adata.obsm["spatial"] = adata.obs[["pxl_col","pxl_row"]].values
    return adata

这个函数我用了两年，兼容大部分10x Visium GEO下载包。注意它假设了空间坐标列名是 pxl_col 和 pxl_row，但有些新版本改成了 x 和 y。遇到报错时先看看 CSV 的列数。

好了，数据到你手里了。接下来就可以开始质控、归一化、批次矫正……不过那是后面章节的事。你在实际动手时如果遇到某个格式读不进去，可以先试试用 h5py 或 pandas 直接 peek 一下文件结构，避免盲目依赖上层封装。我当年被一个不规范的 H5 文件卡了两天，最后发现是 matrix.mtx 的排序跟 barcodes.tsv 不一致——没办法，手动排序解决的。