单细胞 RNA-seq 质量控制：常见问题

import scanpy as sc
import anndata as ad

# 1. .h5ad 文件（AnnData 格式，推荐）
adata = sc.read_h5ad('data.h5ad')

# 2. .h5 文件（10X Genomics 输出）
adata = sc.read_10x_h5('filtered_feature_bc_matrix.h5')

# 3. 10X 目录（包含 barcodes.tsv, genes.tsv, matrix.mtx）
adata = sc.read_10x_mtx('path/to/10X/output/')

# 4. .csv 文件（基因 × 细胞矩阵）
import pandas as pd
counts = pd.read_csv('counts.csv', index_col=0)
adata = ad.AnnData(X=counts.T)  # 转置为细胞 × 基因

# 5. .txt 文件（类似 csv）
counts = pd.read_csv('counts.txt', sep='\t', index_col=0)
adata = ad.AnnData(X=counts.T)

# 验证加载结果
print(f"数据形状：{adata.shape}")
print(f"细胞数：{adata.n_obs}")
print(f"基因数：{adata.n_vars}")

import scanpy as sc

# 1. 加载多个数据集
adata1 = sc.read_h5ad('sample1.h5ad')
adata2 = sc.read_h5ad('sample2.h5ad')
adata3 = sc.read_h5ad('sample3.h5ad')

# 2. 添加批次标识
adata1.obs['batch'] = 'sample1'
adata2.obs['batch'] = 'sample2'
adata3.obs['batch'] = 'sample3'

# 3. 合并（推荐方法）
adata_combined = adata1.concatenate(adata2, adata3, index_unique=None)

print(f"合并后：{adata_combined.n_obs} 细胞")
print(f"批次分布:\n{adata_combined.obs['batch'].value_counts()}")

# 4. 执行 QC（一起过滤所有批次）
# ... QC 代码 ...

# 5. （可选）分回各批次
adata1_qc = adata_combined[adata_combined.obs['batch'] == 'sample1']
adata2_qc = adata_combined[adata_combined.obs['batch'] == 'sample2']

# 取基因交集
common_genes = adata1.var_names.intersection(adata2.var_names)
adata1 = adata1[:, common_genes]
adata2 = adata2[:, common_genes]

# 方案 1：检查是否保存了原始层
print("可用的层：")
print(list(adata.layers.keys()))

# 如果有 'counts' 或 'raw' 层
if 'counts' in adata.layers:
    adata.X = adata.layers['counts'].copy()

# 方案 2：从 log1p 转换回去（近似，不推荐）
# 如果数据是 log1p(counts + 1) 转换的
adata.X = np.expm1(adata.X) - 1  # 反向变换

# 方案 3：重新加载原始数据（最佳）
# 如果你有原始 .h5 文件
adata_raw = sc.read_10x_h5('raw_feature_bc_matrix.h5')

# 1. 检查基因名格式
print("前 20 个基因名：")
print(adata.var_names[:20])

# 2. 可能的情况：
# 情况 A: 基因名是 Ensembl ID（如 ENSG00000244734）
# 解决：需要转换基因 ID 为基因符号
# 使用 mygene 包或参考文件

# 情况 B: 基因名大小写不同
# 小鼠（小写 mt-）
adata.var['mt'] = adata.var_names.str.lower().str.startswith('mt-')

# 人类（大写 MT-）
adata.var['mt'] = adata.var_names.str.upper().str.startswith('MT-')

# 情况 C: 特殊格式（如 "mt-Co1"）
# 使用更灵活的匹配
adata.var['mt'] = adata.var_names.str.contains('mt-', case=False, regex=True)

# 3. 验证
mt_genes = adata.var_names[adata.var['mt']]
print(f"检测到 {len(mt_genes)} 个线粒体基因：")
print(mt_genes[:10])

# 错误 1: KeyError: 'mt'
# 原因：未先定义 mt 列
# 解决：
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)

# 错误 2: 内存不足
# 原因：数据太大
# 解决：分批计算或使用稀疏矩阵
from scipy import sparse
if not sparse.issparse(adata.X):
    adata.X = sparse.csr_matrix(adata.X)

# 错误 3: 数据类型不对
# 原因：X 不是计数矩阵
# 解决：检查数据
print(f"数据类型：{adata.X.dtype}")
print(f"最小值：{adata.X.min()}")
print(f"最大值：{adata.X.max()}")
# 如果是负数或小数，可能已归一化

# 1. 定义基因列表
hemoglobin_genes = ['HBA1', 'HBA2', 'HBB', 'HBD', 'HBE1', 'HBG1', 'HBG2']
cell_cycle_genes = ['MKI67', 'PCNA', 'MCM5', 'STMN1', 'TOP2A']

# 2. 计算自定义指标
def calculate_gene_percentage(adata, gene_list, var_name='custom'):
    """计算特定基因列表的 UMI 百分比"""
    # 找到存在的基因
    genes_present = [g for g in gene_list if g in adata.var_names]

    if not genes_present:
        print(f"警告：{var_name} 基因未找到")
        adata.obs[f'pct_counts_{var_name}'] = 0
        return

    # 计算这些基因的总 UMI
    gene_mask = adata.var_names.isin(genes_present)
    custom_counts = adata[:, gene_mask].X.sum(axis=1).A1
    total_counts = adata.X.sum(axis=1).A1

    # 计算百分比
    percentage = (custom_counts / total_counts) * 100
    adata.obs[f'pct_counts_{var_name}'] = percentage

    print(f"{var_name} 基因数：{len(genes_present)}")
    print(f"平均占比：{percentage.mean():.2f}%")

# 3. 使用
calculate_gene_percentage(adata, hemoglobin_genes, 'hb')
calculate_gene_percentage(adata, cell_cycle_genes, 'cellcycle')

# 4. 查看
print(adata.obs[['pct_counts_hb', 'pct_counts_cellcycle']].describe())

# 1. 检查数据分布
print("=== QC 指标统计 ===")
print(adata.obs[['n_genes_by_counts', 'n_counts', 'pct_counts_mt']].describe())

# 2. 使用更宽松的阈值（5 MAD）
N_MADS = 5  # 增加到 5

# 3. 分别检查各指标的过滤率
def mad_filter(data, n_mads, direction='both'):
    median = np.median(data)
    mad = np.median(np.abs(data - median))
    if direction == 'lower':
        return data > (median - n_mads * mad)
    elif direction == 'upper':
        return data < (median + n_mads * mad)
    else:
        lower = median - n_mads * mad
        upper = median + n_mads * mad
        return (data > lower) & (data < upper)

for metric in ['n_genes_by_counts', 'n_counts', 'pct_counts_mt']:
    if metric == 'pct_counts_mt':
        filter_mask = mad_filter(adata.obs[metric], 3, 'upper')
    else:
        filter_mask = mad_filter(adata.obs[metric], 3, 'both')

    retention = filter_mask.sum() / len(filter_mask) * 100
    print(f"{metric}: {retention:.1f}% 保留")

# 4. 检查不同批次的保留率
if 'batch' in adata.obs.columns:
    for batch in adata.obs['batch'].unique():
        batch_cells = adata.obs['batch'] == batch
        retention = (filter_mask[batch_cells].sum() / batch_cells.sum() * 100)
        print(f"  {batch}: {retention:.1f}%")

# 5. 如果某个批次保留率特别低，考虑单独 QC 或丢弃

# 可视化
import matplotlib.pyplot as plt
plt.hist(adata.obs['n_genes_by_counts'], bins=100)
plt.show()

# 如果出现双峰，可能原因：
# 1. 混入不同细胞类型
# 2. 技术批次效应
# 3. 部分样本质量差

# 方案 1: 分组过滤
if 'batch' in adata.obs.columns:
    # 对每个批次分别过滤
    adata_list = []
    for batch in adata.obs['batch'].unique():
        batch_adata = adata[adata.obs['batch'] == batch]
        # 对该批次应用 MAD 过滤
        # ...
        adata_list.append(batch_adata_filtered)

    # 合并
    adata_filtered = adata_list[0].concatenate(adata_list[1:], index_unique=None)

# 方案 2: 使用百分位数而非 MAD
lower = np.percentile(adata.obs['n_genes_by_counts'], 5)
upper = np.percentile(adata.obs['n_genes_by_counts'], 95)
filter_mask = (adata.obs['n_genes_by_counts'] > lower) & (adata.obs['n_genes_by_counts'] < upper)

# 方案 3: 基于建模的方法（如 GMM）
from sklearn.mixture import GaussianMixture

# 选择"高质量"峰
gmm = GaussianMixture(n_components=2, random_state=42)
gmm.fit(adata.obs[['n_genes_by_counts']])
labels = gmm.predict(adata.obs[['n_genes_by_counts']])

# 选择基因数较高的簇
high_quality_cluster = np.argmax(gmm.means_)
filter_mask = labels == high_quality_cluster

# 1. 检查过滤前该细胞群的 QC 特征
# 假设有细胞类型注释
if 'cell_type' in adata.obs.columns:
    for ct in adata.obs['cell_type'].unique():
        ct_cells = adata.obs['cell_type'] == ct
        print(f"\n{ct}:")
        print(f"  n_genes: {adata.obs[ct_cells]['n_genes_by_counts'].median():.0f}")
        print(f"  n_counts: {adata.obs[ct_cells]['n_counts'].median():.0f}")
        print(f"  mt%: {adata.obs[ct_cells]['pct_counts_mt'].median():.1f}%")

# 2. 调整过滤策略
# 如果某个细胞类型天生基因数低，不要使用固定的 n_genes 下限
# 如果某个细胞类型线粒体含量高，使用更高的 mt% 阈值

# 3. 基于细胞类型的分群过滤（如果有参考）
# 对不同细胞类型使用不同阈值
filter_mask = np.ones(adata.n_obs, dtype=bool)

for ct in adata.obs['cell_type'].unique():
    ct_cells = adata.obs['cell_type'] == ct
    ct_adata = adata[ct_cells]

    # 对该细胞类型计算 MAD
    median = np.median(ct_adata.obs['n_genes_by_counts'])
    mad = np.median(np.abs(ct_adata.obs['n_genes_by_counts'] - median))
    lower = median - 5 * mad  # 使用 5 MAD（宽松）

    ct_filter = ct_adata.obs['n_genes_by_counts'] > lower
    filter_mask[ct_cells] = ct_filter

import scrublet as scr

# 错误 1: 数据类型问题
counts_matrix = adata.X.toarray() if hasattr(adata.X, 'toarray') else adata.X
scrub = scr.Scrublet(counts_matrix)

# 错误 2: 双细胞率异常高
# 原因：可能是阈值设置不当
scrub = scr.Scrublet(counts_matrix)
doublet_scores, predicted_doublets = scrub.scrub_doublets()

# 可视化分布
scrub.plot_doublet_score_histogram()

# 手动调整阈值
# 在双峰之间的谷处选择阈值
threshold = 0.25  # 根据直方图调整
predicted_doublets_manual = doublet_scores > threshold

print(f"自动阈值双细胞率：{predicted_doublets.mean():.1%}")
print(f"手动阈值双细胞率：{predicted_doublets_manual.mean():.1%}")

# 错误 3: 期望双细胞率设置不当
expected_doublet_rate = 0.06  # 根据实际加载调整
scrub = scr.Scrublet(counts_matrix, expected_doublet_rate=expected_doublet_rate)

# 方法 1: 检查双细胞的基因表达特征
doublets = adata.obs['predicted_doublet']
adata_doublets = adata[doublets]

# 检查是否表达两个不同细胞类型的标志基因
# 例如：T 细胞（CD3D）+ B 细胞（MS4A1）
if 'CD3D' in adata.var_names and 'MS4A1' in adata.var_names:
    t_cell = adata_doublets[:, 'CD3D'].X > 0
    b_cell = adata_doublets[:, 'MS4A1'].X > 0
    both = (t_cell & b_cell).sum()
    print(f"同时表达 CD3D 和 MS4A1 的细胞：{both}")

# 方法 2: 检查 QC 指标
print("双细胞 QC 特征：")
print(adata_doublets.obs[['n_genes_by_counts', 'n_counts']].describe())
print("\n单细胞 QC 特征：")
print(adata[~doublets].obs[['n_genes_by_counts', 'n_counts']].describe())
# 双细胞应该有更高的 n_genes 和 n_counts

# 1. 基本统计
print("=== 过滤前后对比 ===")
print(f"细胞数：{adata.n_obs} → {adata_filtered.n_obs}")
print(f"保留率：{adata_filtered.n_obs / adata.n_obs * 100:.1f}%")
print(f"基因数：{adata.n_vars} → {adata_filtered.n_vars}")

# 2. QC 指标改善
print("\n=== MT% 改善 ===")
print(f"过滤前中位数：{adata.obs['pct_counts_mt'].median():.1f}%")
print(f"过滤后中位数：{adata_filtered.obs['pct_counts_mt'].median():.1f}%")

# 3. 可视化对比
import matplotlib.pyplot as plt

fig, axes = plt.subplots(2, 3, figsize=(15, 10))

# 过滤前
axes[0, 0].violinplot([adata.obs['n_genes_by_counts']])
axes[0, 0].set_title('n_genes (过滤前)')
axes[0, 1].violinplot([adata.obs['n_counts']])
axes[0, 1].set_title('n_counts (过滤前)')
axes[0, 2].violinplot([adata.obs['pct_counts_mt']])
axes[0, 2].set_title('MT% (过滤前)')

# 过滤后
axes[1, 0].violinplot([adata_filtered.obs['n_genes_by_counts']])
axes[1, 0].set_title('n_genes (过滤后)')
axes[1, 1].violinplot([adata_filtered.obs['n_counts']])
axes[1, 1].set_title('n_counts (过滤后)')
axes[1, 2].violinplot([adata_filtered.obs['pct_counts_mt']])
axes[1, 2].set_title('MT% (过滤后)')

plt.tight_layout()
plt.show()

# 4. 快速聚类检查（可选）
sc.pp.pca(adata_filtered)
sc.pp.neighbors(adata_filtered)
sc.tl.leiden(adata_filtered, resolution=0.5)
sc.tl.umap(adata_filtered)
sc.pl.umap(adata_filtered, color='leiden')

# 期望：聚类清晰，无明显技术伪影

# 如果保留率过低，检查
retention = adata_filtered.n_obs / adata.n_obs * 100
if retention < 70:
    print("警告：保留率过低！")
    print("建议：")
    print("1. 增加 MAD 倍数到 5")
    print("2. 检查是否有技术批次")
    print("3. 考虑组织特异性调整")
    print("4. 验证数据质量")

# 1. 检查版本
python -c "import scanpy as sc; print(sc.__version__)"

# 2. 更新到最新稳定版
pip install --upgrade scanpy

# 3. 或安装特定版本（如果需要）
pip install scanpy==1.10.0

# 4. 完整的环境重建（推荐）
conda create -n scrna-qc python=3.10 -y
conda activate scrna-qc
pip install scanpy==1.10.0 anndata==0.10.0

# 1. 使用稀疏矩阵
from scipy import sparse
if not sparse.issparse(adata.X):
    adata.X = sparse.csr_matrix(adata.X)

# 2. 分批处理
def batch_qc(adata, batch_size=10000):
    """分批计算 QC 指标"""
    results = []
    for i in range(0, adata.n_obs, batch_size):
        batch = adata[i:i+batch_size]
        sc.pp.calculate_qc_metrics(batch, qc_vars=['mt'], inplace=True)
        results.append(batch.obs.copy())

    import pandas as pd
    adata.obs = pd.concat(results)

# 3. 删除中间变量
import gc
del adata_large  # 删除大对象
gc.collect()  # 强制垃圾回收

# 4. 使用 Dask（适合超大数据）
# import dask.array as da
# adata.X = da.from_array(adata.X, chunks=(10000, 1000))

# CITE-seq 数据有两个矩阵
# RNA: adata.X
# Protein: adata.obsm['protein_counts']

# 1. RNA QC（标准流程）
adata.var['mt'] = adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)

# 2. 蛋白质 QC
protein_counts = adata.obsm['protein_counts']
adata.obs['n_proteins'] = (protein_counts > 0).sum(axis=1)
adata.obs['total_protein_counts'] = protein_counts.sum(axis=1)

# 3. 联合过滤
# 基于 RNA 和蛋白质指标
filter_mask = (
    mad_filter(adata.obs['n_genes_by_counts'], 3, 'lower') &
    mad_filter(adata.obs['n_proteins'], 3, 'lower')  # 蛋白质指标
)

adata_filtered = adata[filter_mask].copy()

# 4. 可视化
import matplotlib.pyplot as plt
plt.scatter(
    adata.obs['n_genes_by_counts'],
    adata.obs['n_proteins'],
    alpha=0.5
)
plt.xlabel('n_genes')
plt.ylabel('n_proteins')
plt.show()

#!/usr/bin/env python3
"""自动化单细胞 QC 脚本"""

import scanpy as sc
import numpy as np
import sys

def mad_filter(data, n_mads=3, direction='both'):
    """MAD 过滤"""
    median = np.median(data)
    mad = np.median(np.abs(data - median))
    if direction == 'lower':
        return data > (median - n_mads * mad)
    elif direction == 'upper':
        return data < (median + n_mads * mad)
    else:
        lower = median - n_mads * mad
        upper = median + n_mads * mad
        return (data > lower) & (data < upper)

def auto_qc(input_file, output_file, species='human', n_mads=3):
    """自动化 QC"""
    print(f"=== QC: {input_file} ===")

    # 1. 加载
    adata = sc.read_h5ad(input_file)
    print(f"加载：{adata.n_obs} 细胞 × {adata.n_vars} 基因")

    # 2. 识别线粒体基因
    if species == 'human':
        adata.var['mt'] = adata.var_names.str.startswith('MT-')
    elif species == 'mouse':
        adata.var['mt'] = adata.var_names.str.startswith('mt-')
    else:
        adata.var['mt'] = adata.var_names.str.contains('mt-', case=False)

    # 3. 计算 QC 指标
    sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], inplace=True)

    # 4. 过滤
    filter_mask = (
        mad_filter(adata.obs['n_genes_by_counts'], n_mads, 'lower') &
        mad_filter(adata.obs['n_genes_by_counts'], n_mads, 'upper') &
        mad_filter(adata.obs['n_counts'], n_mads, 'lower') &
        mad_filter(adata.obs['n_counts'], n_mads, 'upper') &
        mad_filter(adata.obs['pct_counts_mt'], n_mads, 'upper')
    )

    # 5. 应用过滤
    adata_filtered = adata[filter_mask].copy()
    sc.pp.filter_genes(adata_filtered, min_cells=10)

    # 6. 保存
    adata_filtered.write(output_file)

    retention = filter_mask.sum() / len(filter_mask) * 100
    print(f"保存：{adata_filtered.n_obs} 细胞 ({retention:.1f}%)")
    print(f"输出：{output_file}\n")

if __name__ == '__main__':
    # 命令行使用
    if len(sys.argv) < 3:
        print("用法：python auto_qc.py input.h5ad output.h5ad [species] [n_mads]")
        sys.exit(1)

    input_file = sys.argv[1]
    output_file = sys.argv[2]
    species = sys.argv[3] if len(sys.argv) > 3 else 'human'
    n_mads = int(sys.argv[4]) if len(sys.argv) > 4 else 3

    auto_qc(input_file, output_file, species, n_mads)

# 使用：
# python auto_qc.py data_raw.h5ad data_qc.h5ad human 3