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

空间转录组数据入口：从Space Ranger到三个核心文件

做空间转录组分析，第一步就是搞清楚数据长什么样。说实话，我刚接触10x Visium数据时，也被一堆文件搞晕过——h5、mtx、tsv，到底该用哪个？

这一章，我就带你把这几个核心文件彻底吃透。你搞懂了它们，后面80%的数据处理问题都能自己解决。

Space Ranger到底输出了什么

Space Ranger做完比对和计数后，会在输出目录下生成一套标准文件。我个人习惯先扫一眼这个目录结构：

outs/
├── filtered_feature_bc_matrix/   # 过滤后的表达矩阵
│   ├── barcodes.tsv
│   ├── features.tsv
│   └── matrix.mtx
├── raw_feature_bc_matrix/         # 未过滤的原始矩阵
│   ├── barcodes.tsv
│   ├── features.tsv
│   └── matrix.mtx
├── molecule_info.h5               # 分子级别的H5文件
├── spatial/
│   ├── tissue_positions_list.csv
│   ├── scalefactors_json.json
│   └── tissue_hires_image.png
├── analysis/
│   ├── clustering/
│   └── pca/
└── web_summary.html

嗯，这里要注意——我们真正要用的核心数据，只有三个文件格式：H5、barcodes、features、matrix。其他那些图像和聚类结果，说实话是给预览用的。

为什么有filtered和raw两套？

raw保存了所有检测到的barcode，包括那些可能是空液滴的背景噪音。filtered则是Space Ranger根据算法判断的“有效组织区域”。我在项目中通常直接用filtered版本，除非你要做自己的空液滴过滤。

核心原则：先检查filtered_feature_bc_matrix，如果发现细胞数明显偏少（比如芯片上明明有组织，却只检出几百个spot），再回去看raw_feature_bc_matrix。

H5文件：一个容器装所有

我个人非常喜欢H5格式——一个文件就把基因表达矩阵、barcode列表、feature信息全打包了。你想想看，如果是mtx格式需要三个文件，万一传输时丢了一个，数据就废了。

拿Python读取H5文件很简单：

import h5py

file_path = "filtered_feature_bc_matrix.h5"
with h5py.File(file_path, "r") as f:
    # 看一眼里面有什么
    print(list(f.keys()))
    
    # 矩阵数据
    matrix = f["matrix"]
    data = matrix["data"][:]
    indices = matrix["indices"][:]
    indptr = matrix["indptr"][:]
    shape = matrix["shape"][:]
    
    # 条码和特征
    barcodes = [b.decode() for b in f["matrix/barcodes"][:]]
    features = f["matrix/features"]
    gene_ids = [g.decode() for g in features["id"][:]]
    gene_names = [n.decode() for n in features["name"][:]]
    
    print(f"矩阵形状: {shape}")
    print(f"条码数: {len(barcodes)}")
    print(f"基因数: {len(gene_names)}")

输出大概长这样：

['matrix']
矩阵形状: [36601, 4992]
条码数: 4992
基因数: 36601

你看，一个H5文件里包含了3万多个基因和近5000个spot的信息。我建议你在拿到H5文件后，第一件事就是打印出这些基本信息，确认数据量是否合理。

小技巧：如果你用的是Scanpy，直接用sc.read_10x_h5()就能读取。但我还是推荐你先用h5py看一眼原始结构，这样出问题时更容易定位。

H5文件内部采用了一种叫CSR（Compressed Sparse Row）的压缩格式来存储矩阵。为什么要这么做？因为表达矩阵里绝大多数位置的值都是0——这不表达了，没检测到的基因占99%以上。稀疏存储能把文件体积缩小几十倍。

barcodes.tsv：每个spot的身份ID

barcodes.tsv，说白了就是每个捕获区域的条码序列。一行一个barcode，没有表头，没有额外信息。

AAACCAACTCAGAGAT-1
AAACCAACTCGGGCAT-1
AAACCAAGTCGACTCT-1
AAACCAAGTCGGCATC-1
...

它有两种情况：

filtered版本：只有被判定为“有效组织”的barcode，数量一般在1000~5000之间
raw版本：包含芯片上所有可能捕获到的barcode，Visium芯片是4992个spot，raw版本会包含这些全部（以及一些空液滴噪音）

这里有个坑：barcode的命名规则是{实际条码}-{样本编号}。我印象最深的一次，是帮同事排查数据整合问题，发现两个样本的barcode竟然有重复——原来是样本编号没加上去。你以后合并多个样本时，一定记得在barcode后面标记来源。

注意：不要直接修改barcode的内容！你后续的定位坐标、图像配准全部依赖这个条码去关联。改错了，数据就全乱了。

features.tsv：基因目录和feature类型

features.tsv在旧版Space Ranger里叫genes.tsv，v3之后改成features.tsv——因为现在不光能测基因表达，还能加抗体标签（ADT）、CRISPR扰动等信息。

这个文件有三列，用Tab分隔：

ENSG00000243485  MIR1302-2HG      Gene Expression
ENSG00000237613  FAM138A          Gene Expression
ENSG00000186092  OR4F5            Gene Expression
...

三列分别是：

列号	名称	含义
1	id	Ensembl基因ID（或feature的唯一标识符）
2	name	基因符号（如GAPDH、ACTB）
3	feature_type	通常是"Gene Expression"，也可能是"Antibody Capture"等

我个人习惯读取时直接转成DataFrame：

import pandas as pd

features = pd.read_csv(
    "filtered_feature_bc_matrix/features.tsv.gz",
    sep="\t",
    header=None,
    names=["gene_id", "gene_name", "feature_type"]
)
print(features.head())

输出：

            gene_id    gene_name     feature_type
0  ENSG00000243485  MIR1302-2HG  Gene Expression
1  ENSG00000237613      FAM138A  Gene Expression
2  ENSG00000186092        OR4F5  Gene Expression

为什么要保留gene_id？我在项目中遇到过基因名重复的情况——比如同一个symbol对应多个Ensembl ID。这时候用gene_id做去重才是准确的，不要偷懒只用gene_name。

matrix.mtx：真正的表达矩阵

这是三个文件里最“硬核”的一个。它采用Market Exchange格式，专门存储稀疏矩阵。

打开看一眼：

%%MatrixMarket matrix coordinate integer general
%metadata_json: {"software_version": "2.0.0", "format_version": 3}
36601 4992 29457423
1 1 1
1 2 1
1 3 2
2 1 1
3 1 1
...

解读一下这个文件的结构：

第一行：固定标识，声明这是稀疏矩阵
第二行（%开头）：元数据注释，记录软件版本
第三行：三个数字——基因数条码数非零元素个数
从第四行开始：每行三列——基因索引条码索引表达计数

注意，索引是从1开始的！这个跟Python习惯不一样。我曾经因为没注意这一点，直接用索引去取基因，结果第一个基因怎么也匹配不上——嗯，后来才发现差了1。

读取建议：推荐用scipy.io.mmread()直接读取mtx文件，它自动帮你处理索引偏移。不推荐自己手动解析，除非你特别熟悉稀疏矩阵格式。

非零元素个数是多少呢？拿上面的例子来说：36601个基因、4992个spot，如果全满应该有1.8亿个值。但实际只有2945万个非零元素——稀疏度只有约1.6%。你看，这也是为什么我们必须用稀疏格式来存储。

三个文件的关联关系

现在我们把它们串起来讲。matrix.mtx里的第一列对应features.tsv的第几行，第二列对应barcodes.tsv的第几行。所以：

matrix.mtx[行索引i, 列索引j] = 表达的UMI计数
这个计数对应的是barcodes.tsv中第j行的barcode
对应的基因是features.tsv中第i行的基因

我画了一张图，帮你理清这三者的关系：

梳理一下整个数据价值链：

H5文件：把三个文件打包成一个，方便传输和管理
barcodes.tsv：标识每个空间spot
features.tsv：定义每个基因/feature
matrix.mtx：存储稀疏表达矩阵

在实际项目中，你直接读取H5是最省事的。但当你需要和空间坐标、组织图像对齐时，就不可避免地要操作barcodes.tsv——因为tissue_positions_list.csv里用的正是这些barcode。

数据质量检查：拿到数据后先做什么

我记得刚入行时，导师跟我说过一句话：“不要信任任何你还没检查过的数据。”后来在项目中，这句话帮我避开了很多坑。

我建议你拿到数据后，先做以下三件事：

# 1. 检查barcode数量是否匹配
import scanpy as sc

adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")
print(f"barcode数: {adata.n_obs}")
print(f"基因数: {adata.n_vars}")

# 2. 查看UMI分布
adata.obs["n_counts"] = adata.X.sum(axis=1).A1
print(f"UMI中位数: {adata.obs['n_counts'].median():.0f}")
print(f"UMI范围: {adata.obs['n_counts'].min()} - {adata.obs['n_counts'].max()}")

# 3. 检查基因检出比例
print(f"检测到基因中位数: {adata.obs['n_genes'].median():.0f}")

输出示例：

barcode数: 4992
基因数: 36601
UMI中位数: 1245
UMI范围: 0 - 15843
检测到基因中位数: 892

这里有几个数值需要你警惕：

UMI中位数 < 500：可能是组织质量差，或者切片时丢失了组织。
检测到基因中位数 < 300：同样提示数据质量偏低。
大量barcode的UMI为0：检查是否包含了太多空液滴，考虑是否应该使用filtered版本。

注意事项：如果用的是raw_feature_bc_matrix，里面可能包含大量UMI=0的barcode。直接聚类会把这些空液滴也算进去，结果会严重失真。

我个人习惯，在做任何下游分析前，先用H5文件跑一遍这个质量检查。如果数据出问题，越早发现修复成本越低。你要知道，等到聚完类才发现数据不对，那回头重做的代价就太大了。

总结一下这一章的核心内容：

H5文件：一站式容器，推荐作为主要数据入口
barcodes.tsv：spot的身份ID，后续所有空间关联的桥梁
features.tsv：基因目录，解决基因名对应问题
matrix.mtx：稀疏表达的数值核心（千万级别的零值不用存）