如何在Python脚本中获取用户输入阈值处理h5ad单细胞数据

需求说明 #

我目前通过手动硬编码参数的方式对h5ad格式的单细胞数据进行过滤，希望实现一个自动化函数，只需传入以下参数即可完成过滤并保存结果：

• h5ad文件的输入路径
• 每个细胞需保留的最小基因数（仅保留表达基因数≥该值的细胞）
• 每个基因需保留的最小细胞数（仅保留在≥该值的细胞中表达的基因）

现有手动处理代码 #

导入依赖模块 #

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

sc.settings.verbosity = 3  # 日志级别：0=错误, 1=警告, 2=信息, 3=提示
sc.settings.set_figure_params(dpi=80, color_map='Greys')
sc.logging.print_versions()

读取h5ad数据 #

adata = sc.read_h5ad('/Single_cell/cerebellar_development/508cf892-174c-45ab-b2dc-05f54f1ee7ed/aldinger20.processed.h5ad')
adata

数据结构信息：

obs: 'orig.ident', 'nCount_RNA', 'nFeature_RNA', 'sample_id', 'percent.mito', 'S.Score', 'G2M.Score', 'Phase', 'CC.Difference', 'nCount_SCT', 'nFeature_SCT', 'age', 'figure_clusters', 'sex', 'type', 'experiment', 'fig_cell_type', 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt'
var: 'name', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts'
obsm: 'X_pca', 'X_umap', 'X_tsne'

手动执行过滤 #

sc.pp.filter_cells(adata, min_genes=500)
sc.pp.filter_genes(adata, min_cells=10)
print(f"过滤后细胞数：{adata.n_obs}，基因数：{adata.n_vars}")

自动化函数实现 #

以下是封装好的自动化函数，包含数据读取、过滤、结果保存功能，同时输出过滤前后的统计信息：

import scanpy as sc

def filter_and_save_h5ad(input_path, output_path, min_genes_per_cell, min_cells_per_gene):
    # 读取数据
    adata = sc.read_h5ad(input_path)
    print(f"原始数据：细胞数={adata.n_obs}，基因数={adata.n_vars}")
    
    # 执行过滤
    sc.pp.filter_cells(adata, min_genes=min_genes_per_cell)
    sc.pp.filter_genes(adata, min_cells=min_cells_per_gene)
    print(f"过滤后数据：细胞数={adata.n_obs}，基因数={adata.n_vars}")
    
    # 保存过滤后的结果
    adata.write_h5ad(output_path)
    print(f"过滤后数据已保存至：{output_path}")

函数使用示例 #

# 示例参数：输入路径、输出路径、每个细胞最小基因数500、每个基因最小细胞数10
filter_and_save_h5ad(
    input_path='/Single_cell/cerebellar_development/508cf892-174c-45ab-b2dc-05f54f1ee7ed/aldinger20.processed.h5ad',
    output_path='/Single_cell/cerebellar_development/aldinger20.filtered.h5ad',
    min_genes_per_cell=500,
    min_cells_per_gene=10
)

内容的提问来源于stack exchange，提问作者PesKchan