内容更新----CNV推断之inferCNVpy

#!/usr/bin/env python
# -*- coding: utf-8 -*-
####zhaoyunfei
####20240321

import argparse
import os
import infercnvpy as cnv
import scanpy as sc
from matplotlib import pyplot as plt

# 创建参数解析器
parser = argparse.ArgumentParser(
    description="运行inferCNVpy分析单细胞CNV",
    formatter_class=argparse.ArgumentDefaultsHelpFormatter
)

# 添加参数
parser.add_argument("-i", "--input_file", required=True, help="输入h5ad文件路径")
parser.add_argument("-c", "--cell_type_col", required=True, help="细胞类型所在列名")
parser.add_argument("-r", "--reference_types", required=True, nargs="+", help="参考细胞类型列表")
parser.add_argument("-o", "--output_dir", default="infercnv_results", help="输出目录")
parser.add_argument("-w", "--window_size", type=int, default=100, help="基因组窗口大小")
parser.add_argument("--resolution", type=float, default=1.0, help="聚类分辨率(值越大聚类越多)")

# 解析参数
args = parser.parse_args()

# 创建输出目录
os.makedirs(args.output_dir, exist_ok=True)

# 读取数据
print(f"读取输入文件: {args.input_file}")
adata = sc.read_h5ad(args.input_file)

# 验证细胞类型列
if args.cell_type_col not in adata.obs.columns:
    raise ValueError(f"列 '{args.cell_type_col}' 不存在于 adata.obs")

# 验证参考细胞类型
missing_refs = [rt for rt in args.reference_types if rt not in adata.obs[args.cell_type_col].unique()]
if missing_refs:
    raise ValueError(f"参考类型 {missing_refs} 不存在于列 '{args.cell_type_col}'")

# 运行inferCNV
print("正在运行inferCNV分析...")
cnv.tl.infercnv(
    adata,
    reference_key=args.cell_type_col,
    reference_cat=args.reference_types,
    window_size=args.window_size,
    step=args.window_size,
    out_dir=args.output_dir
)

# 计算CNV分数
cnv.tl.cnv_score(adata)

# CNV聚类
print(f"使用分辨率 {args.resolution} 进行CNV聚类...")
cnv.tl.pca(adata)
cnv.pp.neighbors(adata)
cnv.tl.leiden(adata, key_added="cnv_leiden", resolution=args.resolution)

# 输出聚类信息
n_clusters = adata.obs["cnv_leiden"].nunique()
print(f"在分辨率 {args.resolution} 下获得 {n_clusters} 个CNV聚类")

# 可视化
print("生成可视化图表...")
sc.pl.umap(adata, color=["cnv_leiden", args.cell_type_col], 
           save="_cnv_clusters.png", show=False, wspace=0.5)
cnv.pl.chromosome_heatmap(adata, groupby="cnv_leiden", 
                         save="_cnv_heatmap.png", show=False)

# 保存聚类结果
cluster_file = os.path.join(args.output_dir, "cnv_clusters.tsv")
adata.obs[["cnv_leiden"]].to_csv(cluster_file, sep="\t")
print(f"聚类结果已保存至 {cluster_file}")

# 保存完整结果
output_file = os.path.join(args.output_dir, "infercnv_results.h5ad")
adata.write(output_file)
print(f"完整结果已保存至 {output_file}")

print("CNV分析完成!")