基于OmicVerse的Bulk RNA-seq数据细胞类型解卷积与模拟单细胞生成-优快云博客

本文链接：https://blog.youkuaiyun.com/Ljugg/article/details/146377270

🔬 基于OmicVerse的Bulk RNA-seq数据细胞类型解卷积与模拟单细胞生成

1. 数据准备与依赖库

import omicverse as ov
import scanpy as sc
import pandas as pd

核心库说明：

OmicVerse：专注于多组学数据整合，其tl.deconvolution模块支持基于单细胞参考的Bulk数据细胞类型分解
Scanpy：单细胞数据处理神器，用于数据标准化、质控和整合
Pandas：用于表格数据读取与处理

2. 数据读取与格式标准化

2.1 Bulk RNA-seq数据加载

bulk_data = pd.read_csv('bulk_rna_seq.csv', index_col=0)

数据格式要求：
CSV文件中行名为基因名，列名为样本名，数值为基因表达量（如原始计数或TPM）
注意事项：
若数据为基因-样本矩阵，需通过.T转置为AnnData要求的样本-基因格式

2.2 单细胞参考数据加载

single_cell_ref = sc.read_10x_mtx(
    'single_cell_data/filtered_gene_bc_matrices/hg19/',
    var_names='gene_symbols', cache=True
)
single_cell_ref.var_names_make_unique()  # 确保基因名唯一

数据来源：
10x Genomics标准输出的过滤后矩阵（filtered_gene_bc_matrices）
格式说明：
包含3个文件：matrix.mtx（表达矩阵）、barcodes.tsv（细胞 barcode）、features.tsv（基因信息）

3. 数据预处理

3.1 Bulk数据标准化

bulk_adata = sc.AnnData(bulk_data.T)  # 转置为样本×基因矩阵
sc.pp.normalize_total(bulk_adata, target_sum=1e4)  # 标准化到总计数10,000
sc.pp.log1p(bulk_adata)  # 对数转换（log(x+1)）

标准化逻辑：
常用的CPM（Counts Per Million）标准化变体，使不同样本间总表达量可比
对数转换作用：
压缩高表达基因动态范围，稳定方差，符合下游统计模型假设

3.2 单细胞参考数据预处理

# 1. 计算质控指标（如线粒体基因比例）
sc.pp.calculate_qc_metrics(
    single_cell_ref, 
    qc_vars=['mt'],  # 假设线粒体基因以'mt-'开头
    percent_top=None, 
    log1p=False, 
    inplace=True
)

# 2. 过滤低质量细胞（线粒体基因<20%）
single_cell_ref = single_cell_ref[single_cell_ref.obs['pct_counts_mt'] < 20, :]

# 3. 标准化与对数转换
sc.pp.normalize_total(single_cell_ref, target_sum=1e4)
sc.pp.log1p(single_cell_ref)

质控参数解析：
- pct_counts_mt < 20%：排除凋亡或受损细胞（高线粒体基因表达）
- 参考单细胞预处理流程，通常还需结合n_genes_by_counts（基因数）和total_counts（总UMI）过滤
标准化一致性：
与Bulk数据使用相同的标准化方法，确保量纲一致

4. 细胞类型解卷积

4.1 加载细胞类型注释

cell_type_annotation = pd.read_csv('cell_type_annotation.csv', index_col=0)
single_cell_ref.obs['cell_type'] = cell_type_annotation

注释文件格式：
需为两列CSV，第一列为细胞barcode，第二列为细胞类型（示例如下）：
```
barcode,cell_type
ATGCTG-1,Fibroblast
TGACGT-1,T_cell
...
```
注意事项：
确保barcode与单细胞矩阵完全匹配，可通过single_cell_ref.obs_names检查一致性

4.2 执行解卷积

deconvolution_result = ov.tl.deconvolution(
    bulk_adata, 
    single_cell_ref, 
    'cell_type'  # 单细胞数据中存储细胞类型的列名
)

算法原理：
基于非负矩阵分解（NMF）或线性回归，推断每个Bulk样本中各细胞类型的比例
输出说明：
DataFrame格式，行名为Bulk样本，列名为细胞类型，数值为比例（总和为1）

5. 模拟单细胞数据生成

5.1 基于比例抽样

simulated_single_cells = []
for sample in deconvolution_result.index:
    proportions = deconvolution_result.loc[sample]
    for cell_type, prop in proportions.items():
        num_cells = int(prop * 100)  # 假设总细胞数为100（可调整）
        # 提取对应细胞类型的单细胞数据
        cells = single_cell_ref[single_cell_ref.obs['cell_type'] == cell_type].copy()
        # 随机抽样（若细胞数不足则重复）
        selected_cells = cells[:num_cells] if len(cells) >= num_cells else cells.sample(num_cells, replace=True)
        simulated_single_cells.append(selected_cells)

参数调优：
- num_cells = prop * 总细胞数：总细胞数可根据实际需求调整（如1000）
- 抽样策略：若某细胞类型在参考数据中细胞数不足，可通过replace=True重复抽样
局限性：
未考虑细胞异质性，可通过添加高斯噪声模拟生物学变异

5.2 整合与保存

simulated_adata = sc.concat(simulated_single_cells)
simulated_adata.write('simulated_single_cell_data.h5ad')

文件格式：
.h5ad是Scanpy推荐的存储格式，支持高效读写和元数据存储
后续分析：
可直接用于单细胞聚类（sc.tl.leiden）、差异表达（sc.tl.rank_genes_groups）等分析

6. 扩展与优化建议

单细胞参考数据预处理增强：

按GitHub流程添加环境RNA校正（Soupx）和Doublet去除（Scrublet）

# 示例：使用Soupx去除环境RNA（需安装soupx）
import soupx
soupx.pp.soup_x(single_cell_ref, inplace=True)

解卷积算法对比：
尝试多种工具（如CIBERSORTx、MuSiC），评估结果一致性
模拟数据改进：
- 保留原始单细胞的UMI计数分布，而非直接使用标准化后数据
- 引入细胞特异性 dropout 噪声

7. 完整代码流程图

Bulk RNA-seq数据 → 标准化 → 解卷积 → 细胞类型比例  
       ↓                          ↑  
单细胞参考数据 → 质控/标准化 → 细胞类型注释  
       ↓  
   模拟单细胞数据（按比例抽样）

通过以上步骤，我们实现了从Bulk数据到模拟单细胞数据的全流程分析，为缺乏单细胞数据时的下游分析提供了有效解决方案。实际应用中需根据数据特点调整参数，并通过交叉验证确保结果可靠性。