单细胞数据分析第一步：用Python scanpy正确读取10x数据，并保存为.h5ad文件

单细胞数据分析第一步：用Python scanpy正确读取10x数据，并保存为.h5ad文件

刚接触单细胞数据分析的Python用户，面对10x Genomics提供的原始数据文件（barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz）时，往往会感到无从下手。本文将带你从零开始，一步步完成从数据导入到持久化保存的全流程操作，让你在Jupyter Notebook或Python脚本中轻松构建scanpy AnnData对象。

1. 环境准备与数据理解

在开始之前，确保你已经安装了必要的Python库。如果你使用的是conda环境，可以通过以下命令安装：

conda install -c conda-forge scanpy pandas scipy

10x Genomics的单细胞数据通常包含三个核心文件：

• barcodes.tsv.gz：包含每个细胞的唯一标识符
• features.tsv.gz（或genes.tsv.gz）：包含基因的Ensembl ID和Gene Symbol
• matrix.mtx.gz：稀疏矩阵格式的基因表达数据

理解这些文件的含义至关重要，因为后续的每一步操作都建立在对数据结构的清晰认识上。

2. 数据导入与初步处理

2.1 文件路径设置与检查

首先，我们需要设置数据存放的路径，并检查文件是否存在：

import os from pathlib import Path data_dir = Path("/path/to/your/10x_data") required_files = ["barcodes.tsv.gz", "features.tsv.gz", "matrix.mtx.gz"] for file in required_files: if not (data_dir / file).exists(): raise FileNotFoundError(f"Missing required file: {file}")

2.2 使用scanpy直接读取10x数据

scanpy提供了便捷的read_10x_mtx函数来读取10x数据：

import scanpy as sc adata = sc.read_10x_mtx( data_dir, # 包含10x文件的目录 var_names="gene_symbols", # 使用Gene Symbol作为变量名 make_unique=True # 确保基因名唯一 )

注意：如果features.tsv.gz文件中的基因名有重复，make_unique=True会为重复的基因名添加后缀使其唯一。

3. 数据结构检查与调整

3.1 理解AnnData对象结构

读取后的数据存储在AnnData对象中，这是scanpy的核心数据结构。让我们查看一下它的基本属性：

print(f"Number of cells: {adata.n_obs}") print(f"Number of genes: {adata.n_var}") print(f"Matrix shape: {adata.X.shape}")

3.2 处理Ensembl ID和Gene Symbol

10x数据通常包含基因的Ensembl ID和Gene Symbol。我们可以将这些信息存储在AnnData对象的var属性中：

import pandas as pd # 读取features文件获取完整信息 features = pd.read_csv(data_dir / "features.tsv.gz", sep="\t", header=None) features.columns = ["gene_ids", "gene_symbols", "feature_type"] # 将信息添加到AnnData对象 adata.var["gene_ids"] = features["gene_ids"].values adata.var["feature_type"] = features["feature_type"].values

4. 数据转置与索引设置

4.1 理解矩阵方向

10x的matrix.mtx.gz文件通常是"基因×细胞"的格式，而scanpy的AnnData对象期望的是"细胞×基因"的格式。如果矩阵方向不正确，我们需要进行转置：

if adata.X.shape[0] == adata.n_var: # 检查是否需要转置 adata = adata.T print("Matrix has been transposed")

4.2 设置obs和var的索引

确保细胞和基因的索引正确设置：

# 设置细胞(barcode)索引 barcodes = pd.read_csv(data_dir / "barcodes.tsv.gz", header=None) adata.obs_names = barcodes[0].values # 设置基因索引 adata.var_names = features["gene_symbols"].values

5. 数据保存与后续使用

5.1 稀疏矩阵与密集矩阵转换

默认情况下，AnnData使用稀疏矩阵存储表达数据以节省内存。但在某些情况下，你可能需要转换为密集矩阵：

# 转换为密集矩阵（谨慎使用，会显著增加内存使用） adata.X = adata.X.toarray()

提示：除非必要，否则建议保持数据为稀疏格式，特别是对于大型单细胞数据集。

5.2 保存为.h5ad文件

最后，我们将处理好的数据保存为.h5ad格式：

output_path = data_dir / "processed_data.h5ad" adata.write(output_path) print(f"Data saved to {output_path}")

.h5ad文件是AnnData的专用格式，保存了所有数据、注释和元数据。后续使用时，只需一行代码即可重新加载：

adata = sc.read_h5ad(output_path)

6. 常见问题与解决方案

6.1 文件命名问题

如果scanpy无法自动识别你的10x文件，可能是因为文件名不符合预期。标准的10x文件命名应该是：

6.2 内存不足问题

处理大型单细胞数据集时，可能会遇到内存不足的情况。以下是一些优化建议：

• 始终使用稀疏矩阵格式
• 过滤掉低质量的细胞和基因
• 考虑使用Dask等工具进行分布式计算

6.3 基因名重复问题

当features.tsv.gz中存在重复的基因名时，可以采取以下策略：

# 方法1：使用make_unique参数（推荐） adata = sc.read_10x_mtx(data_dir, var_names="gene_symbols", make_unique=True) # 方法2：手动处理重复基因名 from collections import defaultdict count = defaultdict(int) new_names = [] for name in features["gene_symbols"]: count[name] += 1 new_names.append(f"{name}_{count[name]}" if count[name] > 1 else name) adata.var_names = new_names

7. 进阶技巧与最佳实践

7.1 添加样本元数据

在实际分析中，你可能需要为每个细胞添加样本信息：

# 假设我们有样本信息 sample_info = pd.DataFrame({ "sample_id": ["sample1"] * len(adata), "batch": ["batch1"] * len(adata) }, index=adata.obs_names) # 将信息添加到AnnData对象 adata.obs = pd.concat([adata.obs, sample_info], axis=1)

7.2 质量控制标记

在读取数据阶段就可以添加一些质量控制标记：

# 计算每个细胞的总计数 adata.obs["total_counts"] = adata.X.sum(axis=1).A1 # 计算检测到的基因数 adata.obs["n_genes_by_counts"] = (adata.X > 0).sum(axis=1).A1 # 计算线粒体基因比例 mito_genes = adata.var_names.str.startswith("MT-") adata.obs["pct_mito"] = ( adata[:, mito_genes].X.sum(axis=1).A1 / adata.obs["total_counts"] * 100 )

7.3 并行处理多个样本

如果你有多个10x样本需要合并分析，可以这样做：

import glob # 获取所有样本路径 sample_dirs = glob.glob("/path/to/samples/*") # 读取并合并所有样本 adatas = [sc.read_10x_mtx(d) for d in sample_dirs] adata = adatas[0].concatenate(adatas[1:])

8. 性能优化与调试技巧

8.1 加速数据读取

对于非常大的数据集，可以尝试以下优化：

# 使用更快的矩阵读取方法 from scipy.io import mmread mtx = mmread(data_dir / "matrix.mtx.gz").tocsr() # 然后手动构建AnnData import anndata adata = anndata.AnnData(X=mtx)

8.2 内存使用监控

在处理大数据时，监控内存使用很有帮助：

import psutil import humanize def print_memory_usage(): process = psutil.Process() mem = process.memory_info().rss print(f"Memory usage: {humanize.naturalsize(mem)}")

8.3 处理特殊字符

有时基因名可能包含特殊字符，需要清理：

import re # 清理基因名中的特殊字符 clean_names = [re.sub(r"[^a-zA-Z0-9_.-]", "_", name) for name in adata.var_names] adata.var_names = clean_names

掌握这些技巧后，你将能够高效地处理各种10x Genomics单细胞数据，为后续的分析打下坚实基础。记住，良好的数据导入和处理习惯是成功分析的一半。