单细胞肿瘤异质性分析¶

# === 单细胞肿瘤异质性分析（Scanpy + InferCNV） ===
import scanpy as sc                    # 单细胞分析
import numpy as np                     # 数值计算

# 1. 加载和预处理
adata = sc.read_10x_h5("filtered_feature_bc_matrix.h5")  # 10X数据
adata.var_names_make_unique()          # 基因名去重

sc.pp.filter_cells(adata, min_genes=200)    # 过滤低质量细胞
sc.pp.filter_genes(adata, min_cells=3)      # 过滤低表达基因

# 计算线粒体基因比例（肿瘤细胞可能偏高）
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], inplace=True)
adata = adata[adata.obs.pct_counts_mt < 20].copy()  # 过滤高线粒体

# 标准处理
sc.pp.normalize_total(adata, target_sum=1e4)  # 归一化
sc.pp.log1p(adata)                     # 对数变换
sc.pp.highly_variable_genes(adata, n_top_genes=3000)  # 高变异基因
sc.pp.scale(adata, max_value=10)       # 标准化
sc.tl.pca(adata, svd_solver="arpack")  # PCA降维
sc.pp.neighbors(adata, n_pcs=30)       # 构建KNN图
sc.tl.umap(adata)                      # UMAP降维
sc.tl.leiden(adata, resolution=0.8)    # Leiden聚类

# 2. 区分肿瘤细胞和非肿瘤细胞（用marker基因）
# 免疫细胞 marker
immune_markers = ["PTPRC", "CD3D", "CD8A", "CD4", "CD14", "CD68"]
# 上皮/肿瘤 marker
tumor_markers = ["EPCAM", "KRT18", "KRT19", "MUC1"]

sc.pl.dotplot(adata, immune_markers + tumor_markers,
              groupby="leiden")        # 点图查看marker表达

# 3. 手动注释细胞类型（根据marker表达）
cell_type_map = {
    "0": "Tumor", "1": "Tumor", "2": "T_cell",
    "3": "Macrophage", "4": "Tumor", "5": "Fibroblast",
    "6": "B_cell", "7": "Endothelial"  # 根据实际marker调整
}
adata.obs["cell_type"] = adata.obs["leiden"].map(cell_type_map)

# 4. 分析肿瘤细胞亚群
tumor_cells = adata[adata.obs["cell_type"] == "Tumor"].copy()
sc.tl.leiden(tumor_cells, resolution=0.5, key_added="tumor_subclone")
sc.pl.umap(tumor_cells, color="tumor_subclone",
           title="肿瘤细胞亚克隆")     # 可视化亚克隆

# === InferCNV 从 scRNA-seq 推断 CNV ===

# 1. 准备输入文件
# annotations.txt: 每个细胞的类型注释（肿瘤/正常）
# gene_order.txt: 基因在染色体上的位置

# 2. 运行 InferCNV（R 脚本）
cat > run_infercnv.R << 'RSCRIPT'
library(infercnv)                      # 加载InferCNV

# 创建 InferCNV 对象
infercnv_obj <- CreateInfercnvObject(
  raw_counts_matrix = "raw_counts.txt",    # 原始计数矩阵
  annotations_file = "annotations.txt",     # 细胞注释
  gene_order_file = "gene_order.txt",       # 基因位置
  ref_group_names = c("T_cell", "Macrophage")  # 正常参考细胞
)

# 运行分析
infercnv_obj <- infercnv::run(
  infercnv_obj,                        # InferCNV对象
  cutoff = 0.1,                        # 基因表达下限
  out_dir = "infercnv_output/",        # 输出目录
  cluster_by_groups = TRUE,            # 按组聚类
  denoise = TRUE,                      # 降噪
  HMM = TRUE,                         # 用HMM识别CNV状态
  num_threads = 8                      # 线程数
)
RSCRIPT

Rscript run_infercnv.R                 # 执行R脚本