支持中心

简体中文

English

简体中文

English

Appearance

ataCNV 分析教程：基于 scATAC-seq 的肿瘤拷贝数变异解析 ​

作者: SeekGene
时长: 13 分钟
字数: 3.0k 字
更新: 2026-02-27
阅读: 0 次
ATAC + RNA 双组学 Notebooks 分析指南 拷贝数变异分析

背景介绍 ​

染色体拷贝数变异（CNV）是肿瘤基因组不稳定性的典型标志，表现为染色体片段的大规模扩增或缺失。CNV 通过改变癌基因（如 MYC）或抑癌基因（如 TP53）的拷贝数，引发基因剂量效应，驱动肿瘤发生、进展与治疗抵抗。在单细胞分辨率下解析 CNV 对揭示肿瘤异质性尤为关键：依托恶性细胞特异的 CNV 图谱（如 7 号染色体扩增、10 号染色体缺失），可在复杂微环境中区分肿瘤细胞与非恶性细胞（如 T 细胞、成纤维细胞），进一步解析亚克隆结构及其表观调控差异。

相较于基于 scRNA-seq 的 CNV 推断，scATAC-seq 以基因组覆盖度为直接信号，不受转录程序与瞬时表达波动影响；在非转录区域也有更均匀的覆盖，因此对大片段与细粒度 CNV 的分辨率与稳健性更高。ataCNV 是专门面向 scATAC-seq 数据推断 CNV 的工具之一，利用基因组固定窗口的覆盖度计数作为 DNA 拷贝状态的代理信号，在单细胞层面识别染色体片段的扩增/缺失事件，从而实现恶性与非恶性细胞的鉴别与亚克隆精细刻画。

目的与意义 ​

  • 细胞级肿瘤鉴别：基于 CNV 图谱将恶性细胞与微环境细胞精准区分
  • 亚克隆异质性解析：重建拷贝数进化轨迹，理解克隆竞争与分化关系
  • 多组学联动解读：将 CNV 与染色质可及性、细胞注释/分群结果整合，阐明基因组结构改变对表观调控与细胞状态的影响
  • 临床转化价值：为预后分层、疗效评估与潜在靶点挖掘提供结构变异证据

ataCNV 算法原理 ​

ataCNV 主要包括两个核心步骤：

标准化步骤：

  • 在 100 万碱基对（Mbp）的基因组区间上生成单细胞读段计数矩阵
  • 根据基因组区间的可比对性和零值条目数量对细胞和基因组区间进行过滤
  • 通过一阶动态线性模型对计数矩阵进行平滑处理，减少极端噪声
  • 执行按细胞进行的局部回归，消除 GC 含量引起的潜在偏差
  • 利用正常细胞数据进行标准化，从染色质可及性等其他混杂因素中解卷积拷贝数信号

分段步骤：

  • 应用多样本 BIC-seq 算法对所有单细胞进行联合分段
  • 估算每个细胞中获取片段的拷贝比
  • 计算每个细胞群的负荷评分，将具有高 CNA 负荷的细胞群归类为恶性细胞

参数设置 ​

请先挂载将要分析的云平台相关流程数据，挂载方式请参照 jupyter 使用教程

具体参数填写：

  • rds：挂载的流程相关的 rds 文件，是流程相关 Seurat 对象；
  • meta：挂载的流程相关的 meta 文件，是流程相关 Seurat 对象里面的 meta.data；
  • outdir：输出目录路径，用于存放分析结果；
  • species：物种信息，仅支持人类（human）；
  • SAMplenames：样本名称，支持多个样本，用逗号分隔；
  • mode：分析模式，支持三种模式：
    • normalcells：使用指定的正常细胞类型进行标准化
    • allcells：使用所有细胞进行聚类后识别正常细胞
    • none：不进行标准化
  • celltype_col：细胞类型列名，选择要分析的细胞类型对应的标签；
  • cluster_col：细胞聚类列名，用于细胞分群信息；
  • Normal_celltype：正常细胞类型名称，当 mode 为"normalcells"时使用
R
# Parameters
#rds = "/PROJ2/STANDARD/scATAC/Signac/example_data/25020230_CC_T/result.rds"
#meta = "/PROJ2/STANDARD/scATAC/Signac/example_data/25020230_CC_T/meta.tsv"
#outdir = "/PROJ2/FLOAT/advanced-analysis/Seek_module_dev/ataCNV/demo/output/"
#species = "human"
#samplenames = "CC_T_A"
#celltype_col = "CellAnnotation"
#cluster_col = ""
#Normal_celltype = "T_cells"
#mode = "normalcells"
R
# Parameters
rds = "/PROJ2/FLOAT/advanced-analysis/test/2388-5504-6420-584654738725821050/input.rds"
meta = "/PROJ2/FLOAT/advanced-analysis/test/2388-5504-6420-584654738725821050/meta.txt"
outdir = "/PROJ2/FLOAT/advanced-analysis/test/2388-5504-6420-584654738725821050/output/"
species = "human"
samplenames = "joint14_A"
mode = "normalcells"
celltype_col = "CellAnnotation"
cluster_col = "wknn_res.0.5_d30_l2_50"
Normal_celltype = "T_cells"

环境准备 ​

R 包加载

请选择 common_r 这个环境进行该整合教程的学习

数据预处理 ​

加载并处理单细胞 ATAC-seq 数据，为 CNV 分析做准备。主要步骤包括：

  1. 数据加载：读取 Seurat 对象和元数据文件
  2. 细胞过滤：根据质控标准过滤低质量细胞
  3. 基因组设置：根据物种信息设置相应的基因组版本
  4. 计数矩阵提取：从 Seurat 对象中提取基因活性计数矩阵
R
options(warn = -1)
suppressWarnings({
  suppressMessages({
      library(data.table)
      library("AtaCNV")
      library(Matrix)
      library(Seurat,lib.loc="/PROJ2/FLOAT/shumeng/apps/miniconda3/envs/Geneactivity/lib/R/library")
      library(Signac,lib.loc="/PROJ2/FLOAT/shumeng/apps/miniconda3/envs/Geneactivity/lib/R/library")
      library(dplyr)
      library(png)
      library(ComplexHeatmap)
      library(RColorBrewer)
           })
})
options(warn = 0)

开始 CNV 分析 ​

开始执行 ataCNV 的 CNV 推断流程。该流程将 scATAC-seq 片段数据转换为基因组固定窗口的覆盖度矩阵，通过噪声控制与事件识别，在单细胞层面解析染色体拷贝数变异。

执行步骤概述:

  • 数据标准化：根据选择的模式对计数矩阵进行标准化处理
  • CNV 推断：调用 BIC-seq 算法进行分段分析，识别拷贝数变异事件
  • 结果输出：生成包含各窗口 CNV 状态的矩阵表格，支持后续可视化与细胞分群整合

关键参数说明:

  • 窗口大小：默认 1Mbp，平衡 CNV 分辨率与统计稳健性
  • 标准化模式：根据是否有正常细胞选择不同的标准化策略
  • 分段算法：使用 BIC-seq 算法进行联合分段分析
  • 染色体处理：自动处理所有常染色体和性染色体
R
# 创建输出目录
dir.create("../result", recursive = TRUE, showWarnings = FALSE)

# 数据加载和预处理
suppressWarnings({
  suppressMessages({
      samplenames <- unlist(strsplit(samplenames, split = ","))
      meta=read.table(meta, sep = "\t", header = T)
      rownames(meta)=meta$barcodes
  })
})
R
# 加载计数矩阵数据
suppressWarnings({
  suppressMessages({
      if(length(samplenames) >1){
          count=list()
          for(i in samplenames){
              count[[i]] <- readRDS(paste0(outdir, i, "_count.rds"))
              count[[i]] <- as.data.frame(count[[i]])
              count[[i]]$barcode=rownames(count[[i]])
              count=data.table::rbindlist(count)
              count <- as.data.frame(count)
              rownames(count)=count$barcode
              count$barcode=NULL
          }
      }else{
          count <- readRDS(paste0(outdir, samplenames, "_count.rds"))
          count <- as.data.frame(count)
      }
      meta=meta[rownames(meta) %in% rownames(count),]
  })
})
R
# 设置基因组版本
if (species == "human") {
  genome <- "hg38"
}else {
  stop(paste0("Error: Unsupported species '", species, "'. Supported values: human"))
}
R
# 数据标准化
options(warn = -1)
suppressWarnings({
  suppressMessages({
      if(mode=="normalcells"){
          norm_re <- AtaCNV::normalize(count,
                                       genome=genome, mode="normal cells",
                                       normal_cells=(meta[[celltype_col]]==Normal_celltype),
                                       output_dir=paste0(outdir,"dir"),
                                       output_name="norm_re.rds"
                                      )
      }else if(mode=="allcells"){
          norm_re <- AtaCNV::normalize(count,
                                       genome=genome, mode="all cells",
                                       cell_cluster=meta[[cluster_col]],
                                       output_dir=paste0(outdir,"dir"),
                                       output_name="norm_re.rds"
                                      )
      }else if(mode=="none"){
          norm_re <- AtaCNV::normalize(count,
                                       genome=genome, mode=mode,
                                       cell_cluster=meta[[cluster_col]],
                                       output_dir=paste0(outdir,"dir"),
                                       output_name="norm_re.rds"
                                      )
      }
          })
})
options(warn = 0)
output
[1] "Step 1: filtering bins and cells..."
[1] " number of bin before filter: 3102"
[1] " number of bin after filter: 2673"
[1] " number of cells before filter: 16653"
[1] " number of cells after filter: 16653"
[1] "Step 2: smoothing read count signals..."
[1] "Step 3: calculating baseline of normal cells..."
[1] " Baseline is calculated using median read count of normal cells"
[1] "Step 4: normalizing against baseline..."
[1] "99th percentile: 1.83815598812295"
[1] "Step 5: GC correction..."
[1] "Step 6: saving results..."
R
calculate_CNV <- function (norm_count, baseline, BICseq = "default", lambda = 5, 
    min_interval = 2, smooth = F, downsample = 0, output_dir = "./", 
    output_name = "CNV_result.rds") 
{
    if (BICseq == "default") {
        if (Sys.info()["sysname"] == "Windows") {
            BICseq <- system.file("BICseq2", "mbic-seq.exe", 
                package = "AtaCNV")
        }
        else if (Sys.info()["sysname"] == "Linux") {
            BICseq <- system.file("BICseq2", "MBICseq", package = "AtaCNV")
            Sys.chmod(BICseq, mode = "0755")
        }
    }
    if (!dir.exists(output_dir)) {
        dir.create(output_dir)
    }
    if (is.vector(baseline)) {
        baseline <- matrix(data = baseline, nrow = nrow(norm_count), 
            ncol = ncol(norm_count), byrow = T)
    }
    bins <- str2bin(colnames(norm_count))
    chr_bkp <- bin2chrbkp(bins)
    print("Running BICseq2 for segmentation...")
    bkp <- c()
    lambda <- lambda
    tmp_dir <- output_dir
    k <- 0
    if (downsample > 0) {
        idx <- sample(1:nrow(norm_count), downsample)
        norm_count_origin <- norm_count
        norm_count <- norm_count[idx, ]
        baseline_origin <- baseline
        baseline <- baseline[idx, ]
    }
    for (i in unique(bins$chr)) {
        print(i)
        if (sum(bins$chr == i) <= 2) {
            next
        }
        norm_count_ <- norm_count[, bins$chr == i]
        bic_input <- matrix(ncol = 2 * nrow(norm_count_) + 2, 
            nrow = ncol(norm_count_))
        bic_input[, 1:2] <- as.matrix(str2bin(colnames(norm_count_))[, 
            2:3])
        bic_input[, 2 * c(1:nrow(norm_count_)) + 1] <- t(norm_count_)
        bic_input[, 2 * c(1:nrow(norm_count_)) + 2] <- t(baseline[, 
            bins$chr == i])
        s <- c(1:ncol(bic_input))
        s[1:2] <- c("start", "end")
        s[2 * c(1:nrow(norm_count_)) + 1] <- rownames(norm_count_)
        s[2 * c(1:nrow(norm_count_)) + 2] <- paste0(rownames(norm_count_), 
            "_ref")
        colnames(bic_input) <- s
        write.table(x = bic_input, file = paste0(tmp_dir, i, 
            ".bin"), quote = F, sep = "\t", row.names = F, col.names = T, 
            append = F)
        system(paste0(BICseq, " -i ", paste0(tmp_dir, i, ".bin"), 
            " -l ", lambda))
        if(file.exists(paste0(tmp_dir, i, ".bin_seg"))){
            seg_re <- read.table(paste0(tmp_dir, i, ".bin_seg"), 
            header = T, sep = "\t")
            bkp_ <- cumsum(seg_re$binNum)
            bkp <- merge_bkp(bkp, bkp_ + k, min_interval = 2) 
        }else{
            warning(paste("Skipping chromosome", i, "because .bin_seg file is missing"))
            k <- k + sum(bins$chr == i)
        }
        
    }
    for (i in unique(bins$chr)) {
        file.remove(paste0(tmp_dir, i, ".bin"))
        file.remove(paste0(tmp_dir, i, ".bin_seg"))
    }
    bkp <- merge_bkp(chr_bkp, bkp, min_interval = 2)
    bkp <- sort(bkp)
    print(paste0("number of breakpoints: ", length(bkp)))
    if (downsample > 0) {
        norm_count <- norm_count_origin
        baseline <- baseline_origin
    }
    CNV_re <- matrix(0, nrow = nrow(norm_count), ncol = length(bkp) - 
        1)
    for (i in 1:(length(bkp) - 1)) {
        if (bkp[i] + 1 == bkp[i + 1]) {
            CNV_re[, i] <- norm_count[, bkp[i + 1]]/baseline[, 
                bkp[i + 1]]
            next
        }
        temp <- apply(norm_count[, (bkp[i] + 1):bkp[i + 1]], 
            1, mean)/apply(baseline[, (bkp[i] + 1):bkp[i + 1]], 
            1, mean)
        CNV_re[, i] <- temp
    }
    CNV <- norm_count
    for (i in 1:(length(bkp) - 1)) {
        CNV[, (bkp[i] + 1):bkp[i + 1]] <- matrix(CNV_re[, i], 
            nrow = nrow(CNV), ncol = bkp[i + 1] - bkp[i])
    }
    if (smooth) {
        print("smoothing...")
        CNV <- smooth_CNV(CNV)
    }
    print("Saving results...")
    CNV <- round(CNV, digit = 3)
    result <- list(copy_ratio = CNV, bkp = bkp, CNV_seg = CNV_re)
    saveRDS(result, file = paste0(output_dir, output_name))
    return(result)
}
R
# CNV分段分析
suppressWarnings({
  suppressMessages({
      seg_re <- calculate_CNV(norm_count=norm_re$norm_count,
                        baseline=norm_re$baseline,
                        output_dir=paste0(outdir,"dir/"),
                        output_name="seg_re.rds")
     })
})
output
[1] "Running BICseq2 for segmentation..."
[1] "chr1"
[1] "chr2"
[1] "chr3"
[1] "chr4"
[1] "chr5"
[1] "chr6"
[1] "chr7"
[1] "chr8"
[1] "chr9"
[1] "chr10"
[1] "chr11"
[1] "chr12"
[1] "chr13"
[1] "chr14"
[1] "chr15"
[1] "chr16"
[1] "chr17"
[1] "chr18"
[1] "chr19"
[1] "chr20"
[1] "chr21"
[1] "chr22"
[1] "chrX"
[1] "chrY"
[1] "number of breakpoints: 55"
[1] "Saving results..."

热图可视化 ​

用热图可视化细胞拷贝数变异情况，热图左侧色带分别表示将细胞按照细胞类型和细胞聚类情况进行排列，每行表示一个细胞，每列表示一个染色体臂，热图的颜色表示 copy ratio。copy ratio 值在 1 附近，表示正常的 2 倍体；copy ratio 在 0-1 之间，表示染色质存在缺失，越靠近 0（即热图颜色越蓝），表示缺失越严重；copy ratio 在 1-2 之间表示染色质存在扩增， 越靠近 2（即热图颜色越红），表示扩增越严重。

R
# 生成CNV热图
options(warn = -1)
suppressWarnings({
  suppressMessages({
      plot_heatmap(copy_ratio=seg_re$copy_ratio, 
             cell_cluster=meta[[cluster_col]],
             cell_annotation=meta[[celltype_col]],
             output_dir="../result/", 
             output_name="copy_ratio.png")
      #plot=readPNG("../result/copy_ratio.png")
      plot <- png::readPNG("../result/copy_ratio.png")
      options(repr.plot.width = 12, repr.plot.height = 8)
      grid::grid.raster(plot)
     })
})
options(warn = 0)
output
[1] "Step 1: clustering..."
[1] "Step 2: plotting..."

结果文件 ​

├── copy_ratio.png 拷贝数变异热图

结果目录：目录下包括此项分析所涉及的所有图片以及表格等结果文件

附录 ​

.txt ：结果数据表格文件，文件以制表符（Tab）分隔。unix/Linux/Mac 用户使用 less 或 more 命令查看；windows 用户使用高级文本编辑器 Notepad++ 等查看，也可以用 Microsoft Excel 打开。

.pdf ：结果图像文件，矢量图，可以放大和缩小而不失真，方便用户查看和编辑处理，可使用 Adobe Illustrator 进行图片编辑，用于文章发表等。

.rds : 包含基因活力 assay 的 Seurat 对象，需要在 R 环境中打开查看和用于进一步分析。

0 条评论·0 条回复