InferCNV 肿瘤单细胞 CNV — 从单细胞转录组推断肿瘤基因组拷贝数变异¶

# 安装 InferCNV（通过 Bioconductor）
if (!require("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("infercnv")

# 安装 JAGS（贝叶斯推断依赖，必须先安装系统包）
# Ubuntu/Debian 系统：
# sudo apt-get install jags

# 安装 rjags（R 接口）
install.packages("rjags")

# 验证安装
library(infercnv)
packageVersion("infercnv")

# 安装辅助包
install.packages(c('ggplot2', 'dplyr', 'pheatmap'))
BiocManager::install("ComplexHeatmap")

library(infercnv)              # 加载 InferCNV

# ── 准备输入文件 ──────────────────────────────────────
# InferCNV 需要三个输入文件：
# 1. 表达矩阵（原始计数，行=基因，列=细胞）
# 2. 注释文件（每个细胞的类型：tumor 或 normal）
# 3. 基因顺序文件（每个基因的染色体位置）

# 假设从 Seurat 对象提取
library(Seurat)
seurat_obj <- readRDS('tumor_seurat.rds')

# 提取原始计数矩阵
raw_counts <- GetAssayData(
  seurat_obj,
  assay = 'RNA',
  slot = 'counts'               # 必须用原始计数（整数），不能用标准化后的
)

# 创建细胞注释文件（data.frame：行=细胞barcode，列=细胞类型）
cell_annotations <- data.frame(
  cell_type = seurat_obj@meta.data$cell_type,
  row.names = rownames(seurat_obj@meta.data)
)
write.table(
  cell_annotations,
  'cell_annotations.txt',
  sep = '\t',
  quote = FALSE,
  col.names = FALSE              # 不写列名
)

library(infercnv)

# ── 下载基因位置文件（hg38）─────────────────────────
# 从 InferCNV 官方仓库下载
# https://github.com/broadinstitute/inferCNV/tree/master/example/hg38_gencode_v27.txt
# 格式：基因名<tab>染色体<tab>起始位置<tab>终止位置

# 示例（hg38 基因位置文件片段）：
# WASH7P   chr1   14362   29370
# MIR1302-2HG  chr1  30366  30503
# ...

# ── 创建 InferCNV 对象 ───────────────────────────────
infercnv_obj <- CreateInfercnvObject(
  raw_counts_matrix = raw_counts,          # 原始计数矩阵（Seurat 中提取）
  annotations_file = 'cell_annotations.txt',  # 细胞注释文件（cell_type 列）
  delim = '\t',                            # 分隔符
  gene_order_file = 'hg38_gencode_v27.txt',   # 基因染色体位置文件
  ref_group_names = c('T_cell', 'B_cell', 'Fibroblast')  # 正常参考细胞类型
)

# ── 运行 InferCNV 分析 ───────────────────────────────
infercnv_obj_result <- infercnv::run(
  infercnv_obj,
  cutoff = 0.1,                  # 基因过滤：至少 10% 的参考细胞表达（10X 用 0.1）
  out_dir = 'infercnv_output/',  # 输出目录
  cluster_by_groups = TRUE,      # 按细胞类型分组聚类显示
  num_threads = 8,               # 并行线程数
  # 降噪设置
  denoise = TRUE,                # 去噪（去除不确定的小波动）
  HMM = TRUE,                   # 隐马尔科夫模型预测 CNV 状态（推荐，更准确）
  HMM_type = 'i6',              # HMM 类型（i6 = 6 状态模型，分辨更细）
  analysis_mode = 'subclusters', # 聚类模式（subclusters = 更精细）
  resume_mode = FALSE,           # 是否从上次断点继续
  output_format = 'pdf'          # 输出图片格式
)

# 输出文件说明（在 infercnv_output/ 目录中）：
# infercnv.png/pdf         →  主要热图（基因组位置 × 细胞）
# infercnv.observations.txt→  肿瘤细胞的 CNV 矩阵
# infercnv.references.txt  →  参考细胞的 CNV 矩阵
# HMM_CNV_predictions.HMMi6.*.pred_cnv_genes.dat  →  HMM 预测的 CNV 基因

# ── 解读结果 ──────────────────────────────────────────
# 热图解读：
# 红色（>1）= 拷贝数增加（扩增）
# 蓝色（<1）= 拷贝数减少（缺失）
# 白色（≈1）= 正常二倍体

# ── 在 Scanpy/Python 中使用 InferCNV 结果 ────────────

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

# 读取 InferCNV 的 CNV 矩阵（观测细胞）
cnv_obs = pd.read_csv(
    'infercnv_output/infercnv.observations.txt',
    sep='\t', index_col=0
)
# 行=基因，列=细胞，转置后行=细胞，列=基因
cnv_obs = cnv_obs.T

# 计算每个细胞的 CNV 分数（离均值的偏差越大，CNV 越多）
cnv_score = cnv_obs.var(axis=1)    # 方差作为 CNV 强度指标

# 将 CNV 分数加入 AnnData
adata = sc.read_h5ad('tumor_annotated.h5ad')
# 按 barcode 对齐
common_cells = list(set(adata.obs_names) & set(cnv_score.index))
adata_sub = adata[common_cells].copy()
adata_sub.obs['CNV_score'] = cnv_score.loc[common_cells].values

# 可视化 CNV 分数分布
sc.pl.umap(
    adata_sub,
    color='CNV_score',             # CNV 分数越高 = 越可能是肿瘤细胞
    color_map='Reds',              # 红色热图
    title='InferCNV Score (高=肿瘤细胞)'
)

# 用 CNV 分数区分肿瘤 vs 正常细胞
cnv_threshold = np.percentile(cnv_score.values, 75)  # top 25% 认定为肿瘤
adata_sub.obs['is_tumor'] = adata_sub.obs['CNV_score'] > cnv_threshold
sc.pl.umap(adata_sub, color='is_tumor', title='推断肿瘤细胞（红=肿瘤）')

# 输入要求
# 1. 原始整数计数矩阵（不能是标准化后的数据）
# 2. 细胞注释文件（包含肿瘤细胞和正常参考细胞的标注）
# 3. 基因染色体位置文件（从官方 GitHub 下载 hg38/mm10 版本）

# 关键参数
# cutoff=0.1   → 10X 数据
# cutoff=1.0   → Smart-seq2 数据
# HMM=TRUE     → 开启 HMM 模型（推荐，更准确）
# denoise=TRUE → 去噪（减少假阳性）

# 输出核心文件
# infercnv.png               → 主要热图
# infercnv.observations.txt  → 肿瘤细胞 CNV 矩阵