{ "cells": [ { "cell_type": "markdown", "id": "94a3feb5", "metadata": {}, "source": [ "---\n", "title: CopyscAT 分析：基于 scATAC-seq 的 CNV 推断与肿瘤识别\n", "author: SeekGene\n", "date: 2026-01-29\n", "tags:\n", " - ATAC + RNA 双组学\n", " - 分析指南\n", " - Notebooks\n", " - 拷贝数变异分析\n", "---\n" ] }, { "cell_type": "markdown", "id": "b0426a4a-48d1-4565-bf34-1227be2f6491", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "# CopyscAT 分析：基于 scATAC-seq 的 CNV 推断与肿瘤识别" ] }, { "cell_type": "markdown", "id": "35f271f8-2f87-4532-ab89-97ed5261a7d1", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## 背景介绍\n", "\n", "染色体拷贝数变异（CNV）是肿瘤基因组不稳定性的典型标志，表现为染色体片段的大规模扩增或缺失。CNV 通过改变癌基因（如 MYC）或抑癌基因（如 TP53）的拷贝数，引发基因剂量效应，驱动肿瘤发生、转移及治疗抵抗。在单细胞分辨率下解析 CNV 对揭示肿瘤异质性尤为关键：依托恶性细胞特异的 CNV 事件（如 7 号染色体扩增、10 号染色体缺失），可在复杂微环境中区分肿瘤细胞与非恶性细胞（如 T 细胞、成纤维细胞），进而解析亚克隆特异性表观调控网络及其与微环境的动态互作。\n", "\n", "Copy-scAT 面向 scATAC-seq 数据的 CNV 推断，通过整合染色质开放信号与基因组结构特征，在单细胞水平构建高精度 CNV 推断模型，实现恶性与非恶性细胞的鉴别与亚克隆精细刻画。\n", "\n", "### 核心原理\n", "\n", "- 局灶性 CNV 检测:\n", " - 利用 1 Mbp 分箱策略校正技术偏倚，识别局灶性扩增事件（如 EGFR 双微体、染色体外 DNA/ecDNA 区域）。\n", "- 大尺度 CNV 推断:\n", " - 基于归一化后的染色体臂水平信号强度，结合高斯混合模型聚类，判别分段性（如 EGFR 扩增）及全染色体臂（如 7 号染色体增益）的拷贝数变异。\n", "- 信号整合与校正:\n", " - 通过变点分析识别显著信号峰值，结合归一化覆盖矩阵，提升 CNV 检测的灵敏度与特异性。" ] }, { "cell_type": "markdown", "id": "6ece47dc", "metadata": {}, "source": [ "## 参数设置\n", "\n", "请先挂载将要分析的云平台相关流程数据，挂载方式请参照 jupyter 使用教程\n", "\n", "具体参数填写：\n", "\n", "- rds：挂载的流程相关的 rds 文件，RDS 文件是标准的 Seurat 对象文件，数据已经过预处理和多个样本合并\n", "- meta：挂载的流程相关的 meta 文件，是流程相关 rds 文件对应的 meta.data ；\n", "- species：物种信息,仅支持人；\n", "- SAMplenames:样本名,注意，为了兼容 SeekArc 输出的 fragment 文件命名方式，样本名需要包含 A，如下所示的\"joint14_A\",\"joint22_A\"，不能只写成\"joint14\",\"joint22\";\n", "- reduction：降维方式，在可视化环节，将每个细胞 CNV 情况映射到降维图上(umap/tsne/lsi/atacumap/atactsne/wnnumap/wnntsne)，方便直观比较不同细胞类型 CNV 的差异;\n", "- celltype_col：meta.data 中细胞类型列名，选择要分析的细胞类型或者聚类对应的标签。假如想对注释好的细胞类型进行分析，则选择其对应的标签，例如 CellAnnotation；\n", "- group_col:meta.data 中细胞分组的列名，样本分组情况，便于热图可视化 CNV 时，细胞的分组展示\n", "- resolution_cluster：meta.data 中细胞聚类列名；\n", "- minimumSegments:最小片段数​​：细胞过滤阈值，默认保留片段数 ≥ 40 的细胞，推荐直接使用默认参数\n", "- methods：auto 和 reference 两种选择，reference 表示以正常细胞作为参考来计算 SNV，auto 则表示不需要正常细胞做参考；\n", "- Normal_celltype:如果上面 methods=”reference”，这里则选择具体用什么细胞用作参考，一般建议用免疫 B 细胞或 T 细胞等作为参考。\n", "- file_path:待分析样本的 fragment 文件,支持多个样本；" ] }, { "cell_type": "code", "execution_count": null, "id": "a8cfbb99", "metadata": { "tags": [ "remove-cell" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "# Parameters\n", "rds = \"/PROJ2/FLOAT/advanced-analysis/test/2388-5504-6349-581745468808125050/input.rds\"\n", "meta = \"/PROJ2/FLOAT/advanced-analysis/test/2388-5504-6349-581745468808125050/meta.txt\"\n", "#outdir = \"/PROJ2/FLOAT/advanced-analysis/test/2388-5504-6349-581745468808125050/output/\"\n", "species = \"human\"\n", "samplenames = c(\"joint14_A\",\"joint22_A\")#为了兼容seekARC输出的fragment文件命名方式，样本名需要包含A，如下所示的\"joint14_A\",\"joint22_A\"，不能只写成\"joint14\",\"joint22\";\n", "reduction = \"wnnumap\"\n", "celltype_col = \"CellAnnotation\"\n", "group_col = \"\"\n", "resolution_cluster = \"wknn_res.0.5_d30_l2_50\"\n", "minimumSegments = \"40\"\n", "methods = \"reference\"\n", "Normal_celltype = \"T_cells\"\n", "file_path=c(\"/PROJ2/FLOAT/advanced-analysis/test/2388-5504-6345-581728456711101050/joint14_A_fragments.tsv.gz\",\"/PROJ2/FLOAT/advanced-analysis/test/2388-5504-6345-581728456711101050/joint22_A_fragments.tsv.gz\")" ] }, { "cell_type": "markdown", "id": "88215398", "metadata": {}, "source": [ "## 环境准备\n", "\n", "**R 包加载**\n", "\n", "请选择 common_r 这个环境进行该整合教程的学习" ] }, { "cell_type": "code", "execution_count": null, "id": "8ec39d76-7137-43ea-a27a-1b5bd78e0c9d", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "options(warn = -1)\n", "suppressWarnings({\n", " suppressMessages({\n", " invisible(capture.output(library(CopyscAT)))\n", " library(GenomicRanges)\n", " library(Seurat, lib.loc = \"/PROJ2/FLOAT/shumeng/apps/miniconda3/envs/Geneactivity/lib/R/library\",quietly = TRUE)\n", " library(tidyverse)\n", " library(ggplot2)\n", " library(ComplexHeatmap)\n", " library(circlize)\n", " library(RColorBrewer)\n", " library(base64enc, quietly = TRUE)\n", " library(DT, quietly = TRUE)\n", " library(htmltools)\n", " library(dplyr)\n", " library(data.table)\n", " source(\"/PROJ2/FLOAT/advanced-analysis/Seek_module_dev/CopyscAT/Module/bin/cnv_caller_functions_packaged.R\")\n", " })\n", "}) \n", "# 5. 根据物种加载基因组数据库\n", "suppressMessages({\n", " if(species==\"human\"){\n", " library(BSgenome.Hsapiens.UCSC.hg38, quietly = TRUE)\n", " } else {\n", " message(\"请安装相应物种的数据库\")\n", " }\n", "})\n", "options(warn = 0)" ] }, { "cell_type": "code", "execution_count": null, "id": "cba55682-5c97-467e-b1a4-dab935a6fa97", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "genomeFilepath=\"/PROJ2/FLOAT/advanced-analysis/Seek_module_dev/CopyscAT/Module/bin/\"\n", "dir.create(\"./result\", recursive = TRUE, showWarnings = FALSE)\n", "setOutputFile(\"./result\",\"result\")\n", "resolution.cluster=resolution_cluster" ] }, { "cell_type": "code", "execution_count": null, "id": "1574b959-444f-4f19-a7a6-9dcaec840d51", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "if (species == \"human\") {\n", " genomeFile <- paste0(genomeFilepath, \"hg38_chrom_sizes.tsv\")\n", " cytobandFile <- paste0(genomeFilepath, \"hg38_1e+06_cytoband_densities_granges.tsv\")\n", " cpgFile = paste0(genomeFilepath, \"hg38_1e+06_cpg_densities.tsv\")\n", "} else {\n", " stop(paste0(\"Error: Unsupported species '\", species, \"'. Supported species are human.\"))\n", "}" ] }, { "cell_type": "markdown", "id": "15e455a2", "metadata": {}, "source": [ "## Fragment 文件预处理\n", "\n", "SeekArctools 或 cellranger arc/ATAC 生成的 fragments 文件包含所有条形码（barcode）的片段信息，其中仅有一部分为真实细胞。为节省计算资源并与后续细胞注释保持一致，进行 CNV 分析前需先对 fragments 进行过滤：仅保留通过质控并用于下游分析的“细胞白名单”条形码。建议依据已确定的细胞元数据/注释表（如 QC 合格且被纳入分析的 barcodes 列表）对子集化 fragments，从而避免将低质量条形码或非细胞背景噪声引入 CNV 推断。" ] }, { "cell_type": "markdown", "id": "6369de94", "metadata": {}, "source": [ "### Fragment 文件过滤\n", "\n", "SeekArctools 或 cellranger arc/ATAC 生成的 fragments 文件包含所有条形码（barcode）的片段信息，其中仅有一部分为真实细胞。为节省计算资源并与后续细胞注释保持一致，进行 CNV 分析前需先对 fragments 进行过滤：仅保留通过质控并用于下游分析的“细胞白名单”条形码。建议依据已确定的细胞元数据/注释表（如 QC 合格且被纳入分析的 barcodes 列表）对子集化 fragments，从而避免将低质量条形码或非细胞背景噪声引入 CNV 推断。" ] }, { "cell_type": "code", "execution_count": null, "id": "fab04fb2", "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "meta <- read.table(meta,sep=\"\\t\",header=T) #读取挂载的meta文件，该文件包含挂载项目数据所有样本的barcode信息\n", "Barcodes <- meta$barcodes #提取全部Barcode\n", "\n", "# 修改循环，为每个样本分别保存文件\n", "for(sample in samplenames){\n", " file_path <- grep(pattern = sample, file_paths, value = TRUE) # 注意这里应该是 file_paths（复数）\n", " frag <- fread(file_path, header=F)\n", " \n", " # 修正这里的列名赋值语法错误\n", " colnames(frag) <- c(\"chr\", \"start\", \"end\", \"Barcode\", \"reads\") # 使用 <- 赋值\n", " \n", " frag_filter <- frag[Barcode %in% Barcodes]\n", " \n", " # 为每个样本分别保存文件，命名格式：{sample}_filtered_fragments.tsv\n", " output_file <- paste0(\"./\", sample, \"_filtered_fragments.tsv\")\n", " write.table(frag_filter, \n", " file = output_file, \n", " quote = FALSE, \n", " row.names = FALSE, \n", " col.names = TRUE, # 写入列名，因为每个文件都是独立的\n", " sep = \"\\t\")\n", "}" ] }, { "cell_type": "markdown", "id": "c41a74f6", "metadata": {}, "source": [ "### Fragment 文件 bin 处理\n", "\n", "### 说明\n", "\n", "CopyscAT 对 scATAC-seq 数据进行 CNV 分析之前，需要对过滤后发 fragment 进一步处理，将每条染色体分成 1000kb 长度的 bin，统计每个 bin 的转座酶切割事件，生成一个新的矩阵文件，行是细胞，列是染色体 bin 的编号,如 chr1_0 表示对第一天染色体将前 1000kb 碱基合并成一个 bin，“0”表示这个 bin 的起始碱基，值表示对应 bin 中该 barcode 中检测到的转座事件总数。如下所示：
\n", "Cell_id chr1_0 chr1_1000000 chr1_2000000 chr1_3000000 chr1_4000000
\n", "GTAACCACGACCGCTAG_1 321.0 6057.0 4540.0 1942.0 2437.0
\n", "AACACACCGACCAGCGT_1 368.0 998.0 490.0 0.0 0.0
\n", "ACGTAGTATGTGAGGTC_1 363.0 1123.0 831.0 56.0 0.0
\n", "AACACACTACATGATTG_1 400.0 2629.0 476.0 658.0 172.0
\n", "TTACGCCCGACCTGCGA_1 261.0 1960.0 2450.0 1916.0 0.0
\n", "CAGGTATGCTTGACTGT_1 429.0 670.0 286.0 1018.0 0.0
\n", "AGTCAACGCTACCGACC_1 2264.0 14874.0 17016.0 11043.0 4900.0
\n", "CTTGAGACGACCGTGGT_1 921.0 4590.0 4726.0 1870.0 1255.0
\n", "CAACCAATACATTTGCC_1 762.0 4810.0 6431.0 5360.0 827.0
\n", "\n", "### 处理方式\n", "\n", "去 CopyscAT 官方 github 复制 process_fragment_file.py 脚本和 hg38_chrom_sizes.tsv 文件，链接：\n", "https://github.com/spcdot/CopyscAT/blob/master/process_fragment_file.py \n", "https://github.com/spcdot/CopyscAT/blob/master/hg38_references/hg38_chrom_sizes.tsv\n", "\n", "然后按要求设定参数，例如：\n", "- python ./process_fragment_file.py -i ./${SAMple}_filtered_fragments.tsv -o ./${SAMple}_filtered_output_1000kb.tsv -b 1000000 -f 0 -g ./hg38_chrom_sizes.tsv" ] }, { "cell_type": "markdown", "id": "d3a4c607", "metadata": {}, "source": [ "## CNV 分析\n", "\n", "基于前述背景与参数设定，我们开始执行 Copy-scAT 的 CNV 推断流程。该流程将 scATAC-seq 片段数据转换为基因组固定窗口的覆盖度矩阵，通过多层次的分析策略，在单细胞层面解析染色体拷贝数变异。\n", "\n", "### 执行步骤概述\n", "- 数据输入：使用上部分处理结果文件${SAMple}_filtered_output_1000kb.tsv 作为输入，支持多个样本的批量处理。\n", "- 局灶性 CNV 检测：基于归一化覆盖矩阵中的显著信号峰值（如染色体外 DNA/ecDNA 区域），通过变点分析识别局灶性扩增事件。\n", "- 大尺度 CNV 推断：利用归一化后的染色体臂水平信号强度，结合高斯混合模型聚类，判别分段性（如 EGFR 扩增）及全染色体臂（如 7 号染色体增益）的拷贝数变异。\n", "- 结果整合：区分恶性与非恶性细胞群，输出单细胞 CNV 图谱及亚克隆表观特征关联分析。\n", "\n", "### 关键参数说明\n", "- 分箱策略：采用 1 Mbp 分箱策略，平衡 CNV 分辨率与统计稳健性\n", "- 信号校正：通过归一化处理校正技术偏倚，提升 CNV 检测的准确性\n", "- 聚类分析：结合高斯混合模型，支持局灶性与大尺度 CNV 的联合检测\n", "- 质量控制：自动过滤低质量信号，确保结果的可靠性" ] }, { "cell_type": "code", "execution_count": null, "id": "b345abf5", "metadata": { "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "suppressWarnings({\n", " suppressMessages({\n", " sc=readRDS(rds)\n", " meta=read.table(meta, sep = \"\\t\", header = T)\n", " rownames(meta)=meta$barcodes\n", " sc=AddMetaData(sc, metadata=meta)\n", " sc@meta.data$Sample <- gsub(\"_arc$\", \"_A\", sc@meta.data$Sample)\n", " sc@meta.data$num <- sub(\"^.*_(\\\\d+)$\", \"\\\\1\", sc@meta.data$barcodes)\n", " })\n", "})" ] }, { "cell_type": "code", "execution_count": null, "id": "f33c6055-241d-40b5-bf12-4979b0a99a49", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "if (length(samplenames) > 1) {\n", " inputdata <- list()\n", " # 根据 samplename 长度动态生成 cellSuffix\n", " suffixes <- paste0(\"_\", seq_along(samplenames)) # 生成如 \"_1\", \"_2\", ..., \"_n\"\n", " # 初始化环境（需根据样本数动态调整 cellSuffix）\n", " initialiseEnvironment(\n", " genomeFile = genomeFile,\n", " cytobandFile = cytobandFile,\n", " cpgFile = cpgFile,\n", " binSize = 1e6,\n", " minFrags = 1e4,\n", " cellSuffix = suffixes, # 动态传入后缀\n", " lowerTrim = 0.5,\n", " upperTrim = 0.8\n", " )\n", " for (i in seq_along(samplenames)) {\n", " inputdata[[i]] <- readInputTable(paste0(outdir, samplenames[[i]], \"_output_1000kb.tsv\"))\n", " inputdata[[i]]$Barcode=rownames(inputdata[[i]])\n", " }\n", " scData <- rbindlist(inputdata)\n", " scData <- as.data.frame(scData)\n", " rownames(scData)=scData$Barcode\n", " #setDT(scData, keep.rownames = TRUE) \n", " scData$Barcode=NULL\n", "} else {\n", " # 单个样本时直接使用 \"_1\"\n", " initialiseEnvironment(\n", " genomeFile = genomeFile,\n", " cytobandFile = cytobandFile,\n", " cpgFile = cpgFile,\n", " binSize = 1e6,\n", " minFrags = 1e4,\n", " cellSuffix = \"_1\", # 单个样本固定后缀\n", " lowerTrim = 0.5,\n", " upperTrim = 0.8\n", " )\n", " scData <- readInputTable(paste0(outdir, samplenames, \"_output_1000kb.tsv\"))\n", " rownames(scData) <- paste0(rownames(scData), \"_1\")\n", " sc@meta.data$barcodes=paste0(sc@meta.data$barcodes,\"_1\")\n", "}\n" ] }, { "cell_type": "code", "execution_count": null, "id": "458a3eaf-0cb4-4ac4-85d0-00963664fb32", "metadata": { "editable": true, "scrolled": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "#PART 1: INITIAL DATA NORMALIZATION\n", "#step 1 normalize the matrix\n", "suppressWarnings({\n", " suppressMessages({\n", " scData_k_norm <- normalizeMatrixN(scData,logNorm = FALSE,maxZero=2000,imputeZeros = FALSE,blacklistProp = 0.8,blacklistCutoff=125,dividingFactor=1,upperFilterQuantile = 0.95)\n", " #collapse into chromosome arm level\n", " summaryFunction<-cutAverage\n", " scData_collapse<-collapseChrom3N(scData_k_norm,summaryFunction=summaryFunction,binExpand = 1,minimumChromValue = 100,logTrans = FALSE,tssEnrich = 1,logBase=2,minCPG=300,powVal=0.73) \n", " })\n", "}) " ] }, { "cell_type": "code", "execution_count": null, "id": "d2bf0ceb-af6b-4c2b-bb69-b9fa361177db", "metadata": { "editable": true, "scrolled": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "options(repr.plot.width = 12, repr.plot.height = 6)\n", "suppressWarnings({\n", " suppressMessages({\n", " #apply additional filters\n", " scData_collapse<-filterCells(scData_collapse,minimumSegments = minimumSegments,minDensity =0.1)\n", " cells=rownames(sc@meta.data[which(sc@meta.data$barcodes %in% colnames(scData_collapse)),])\n", " sc_filter=subset(sc, cells=cells)\n", " #show unscaled chromosome list\n", " graphCNVDistribution(scData_collapse,outputSuffix = (\"_violinsn2\"))\n", " #compute centers\n", " median_iqr <- computeCenters(scData_collapse,summaryFunction=summaryFunction)\n", " })\n", "}) " ] }, { "cell_type": "code", "execution_count": null, "id": "7e56de42-5559-46f3-bc8b-2ce1d5ecf856", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "options(warn = -1)\n", "suppressWarnings({\n", " suppressMessages({\n", " nmf_results<-identifyNonNeoplastic(scData_collapse,methodHclust=\"ward.D\",cutHeight = 0.4)\n", " #?identifyNonNeoplastic\n", " write.table(x=rownames_to_column(data.frame(nmf_results$cellAssigns),var=\"Barcode\"),file=str_c(scCNVCaller$locPrefix,scCNVCaller$outPrefix,\"_nmf_clusters.csv\"),quote=FALSE,row.names = FALSE,sep=\",\")\n", " #print(paste(\"Normal cluster is: \",nmf_results$clusterNormal))\n", " })\n", "}) \n", "options(warn = 0)" ] }, { "cell_type": "code", "execution_count": null, "id": "261316d0-db5c-4dfc-bd22-cbca3ff976dd", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "options(warn = -1)\n", "suppressWarnings({\n", " suppressMessages({\n", " #PART 2: ASSESSMENT OF CHROMOSOME-LEVEL CNVs \n", " if(methods == \"auto\"){\n", " #print(\"methods is auto\")\n", " #OPTION 1: identify chromosome-level amplifications using all cells to generate 'normal' control\n", " #identify chromosome-level amplifications\n", " candidate_cnvs<-identifyCNVClusters(scData_collapse,median_iqr,useDummyCells = TRUE,propDummy=0.25,minMix=0.01,deltaMean = 0.03,deltaBIC2 = 0.25,bicMinimum = 0.1, subsetSize=600,fakeCellSD = 0.08, uncertaintyCutoff = 0.55,summaryFunction=summaryFunction,maxClust = 4,mergeCutoff = 3,IQRCutoff= 0.2,medianQuantileCutoff = 0.4)\n", " #cleanup step\n", " candidate_cnvs_clean<-clusterCNV(initialResultList = candidate_cnvs,medianIQR = candidate_cnvs[[3]],minDiff=1.5) #= 1.5)\n", " #final results and annotation\n", " final_cnv_list<-annotateCNV4(candidate_cnvs_clean, saveOutput=TRUE,outputSuffix = \"clean_cnv\",sdCNV = 0.5,filterResults=TRUE,filterRange=0.8)\n", " \n", " }else if(methods == \"reference\"){\n", " #print(\"methods is reference\")\n", " Idents(sc_filter)=sc_filter@meta.data[[celltype_col]]\n", " table(Idents(sc_filter))\n", " Normal_sc=subset(sc_filter, idents=Normal_celltype)\n", " Normal_cells=Normal_sc@meta.data$barcodes\n", " median_iqr <- computeCenters(scData_collapse %>% dplyr::select(chrom,Normal_cells),summaryFunction=summaryFunction)\n", " #setting medianQuantileCutoff to -1 and feeding non-neoplastic barcodes in as normalCells can improve accuracy of CNV calls\n", " candidate_cnvs<-identifyCNVClusters(scData_collapse,median_iqr,useDummyCells = TRUE,propDummy=0.25,minMix=0.01,deltaMean = 0.03,deltaBIC2 = 0.25,bicMinimum = 0.1, subsetSize=800,fakeCellSD = 0.09, uncertaintyCutoff = 0.65,summaryFunction=summaryFunction,maxClust = 4,mergeCutoff = 3,IQRCutoff = 0.25,medianQuantileCutoff = -1,normalCells=Normal_cells) \n", " candidate_cnvs_clean<-clusterCNV(initialResultList = candidate_cnvs,medianIQR = candidate_cnvs[[3]],minDiff=1.0) #= 1.5) \n", " #to save this data you can use annotateCNV4 as per usual\n", " #final_cnv_list<-annotateCNV4(candidate_cnvs_clean, saveOutput=TRUE,outputSuffix = \"clean_cnv\",sdCNV = 0.6,filterResults=TRUE,filterRange=0.4)\n", " final_cnv_list<-annotateCNV4B(candidate_cnvs_clean, Normal_cells, saveOutput=TRUE,outputSuffix = \"clean_cnv\",sdCNV = 0.6,filterResults=TRUE,filterRange=0.4,minAlteredCellProp = 0.5)\n", "\n", " }\n", " })\n", "}) \n", "options(warn = 0)" ] }, { "cell_type": "markdown", "id": "8cd1f09a-70c2-4647-8165-72f0e78dae01", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### 热图可视化" ] }, { "cell_type": "markdown", "id": "f82ba3a3-e248-4086-8889-bfd3fb309b5f", "metadata": {}, "source": [ "用热图可视化细胞拷贝数变异情况，copy number 值在 2 附近，表示正常的 2 倍体；copy number 在 1-2 之间，表示染色质存在缺失，越靠近 1，表示缺失越严重；copy number 在 2-3 之间表示染色质存在扩增， 越靠近 3 表示扩增越严重。热图左侧色带分别表示将细胞按照细胞类型和细胞聚类情况进行排列，每行表示一个细胞，每列表示一个染色体臂，热图的颜色表示 CNV 分数，颜色越红表示扩增越严重，颜色越蓝表示缺失越严重。" ] }, { "cell_type": "code", "execution_count": null, "id": "30bf7640-59ae-4d0b-9bc6-d03d89b7294d", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "options(warn = -1)\n", "ht_opt$message <- FALSE\n", "\n", "# 1. 读取数据\n", "cnv_data <- final_cnv_list[[3]]\n", "cnv_matrix <- cnv_data %>%\n", " tibble::column_to_rownames(\"rowname\") %>%\n", " as.matrix()\n", "\n", "# 2. 提取元数据\n", "if(nchar(group_col) > 0){\n", " if(nchar(celltype_col) > 0){\n", " metadata <- sc_filter@meta.data[,c(\"raw_Sample\", celltype_col, \"barcode\", resolution.cluster, group_col)]\n", " } else {\n", " metadata <- sc_filter@meta.data[,c(\"raw_Sample\", \"CellAnnotation\", \"barcode\", resolution.cluster, group_col)]\n", " }\n", "} else {\n", " if(nchar(celltype_col) > 0){\n", " metadata <- sc_filter@meta.data[,c(\"raw_Sample\", celltype_col, \"barcode\", resolution.cluster)]\n", " } else {\n", " metadata <- sc_filter@meta.data[,c(\"raw_Sample\", \"CellAnnotation\", \"barcode\", resolution.cluster)]\n", " }\n", "}\n", "\n", "# 3. 动态生成颜色 - 修正所有colorRampPalette函数\n", "# 用于Group的颜色\n", "if(nchar(group_col) > 0) {\n", " n_groups <- length(unique(metadata[[group_col]]))\n", " group_colors <- colorRampPalette(brewer.pal(min(9, 3), \"Set1\"))(n_groups)\n", "}\n", "\n", "# 用于Cluster的颜色\n", "n_clusters <- length(unique(metadata[[resolution.cluster]]))\n", "clusters_colors <- colorRampPalette(brewer.pal(min(8, 3), \"Set2\"))(n_clusters)\n", "\n", "# 用于Sample的颜色\n", "n_samples <- length(unique(metadata$Sample))\n", "sample_colors <- colorRampPalette(brewer.pal(min(12, 3), \"Paired\"))(n_samples)\n", "\n", "# 用于Cell type的颜色\n", "if(nchar(celltype_col) > 0) {\n", " n_celltypes <- length(unique(metadata[[celltype_col]]))\n", " cell_type_colors <- colorRampPalette(brewer.pal(min(12, 3), \"Set3\"))(n_celltypes)\n", "} else {\n", " n_celltypes <- length(unique(metadata$CellAnnotation))\n", " cell_type_colors <- colorRampPalette(brewer.pal(min(12, 3), \"Set3\"))(n_celltypes)\n", "}\n", "\n", "# 4. 准备热图颜色映射\n", "temp <- log(cnv_matrix + 0.01)\n", "d <- parallelDist::parallelDist(temp)\n", "hc_re <- hclust(d, method = \"ward.D2\")\n", "color_thre = c(min(cnv_matrix), max(cnv_matrix))\n", "c <- colorRampPalette(rev(RColorBrewer::brewer.pal(n = 7, name = \"RdYlBu\")))(100)\n", "colormap <- circlize::colorRamp2(c(1:100) * (color_thre[2] - color_thre[1]) * 0.01 + color_thre[1], c)\n", "\n", "# 5. 创建注释 - 隐藏细胞类型图例\n", "if(nchar(group_col) > 0){\n", " if(nchar(celltype_col) > 0){\n", " ha = rowAnnotation(\n", " Sample = metadata$Sample,\n", " Group = metadata[[group_col]],\n", " Cell_type = metadata[[celltype_col]],\n", " Cluster = metadata[[resolution.cluster]],\n", " col = list(\n", " Group = setNames(group_colors, unique(metadata[[group_col]])),\n", " Cluster = setNames(clusters_colors, unique(metadata[[resolution.cluster]])),\n", " Sample = setNames(sample_colors, unique(metadata$Sample)),\n", " Cell_type = setNames(cell_type_colors, unique(metadata[[celltype_col]]))\n", " ),\n", " show_annotation_name = TRUE,\n", " show_legend = c(TRUE, TRUE, TRUE, TRUE) # 隐藏Cell_type的图例\n", " )\n", " } else {\n", " ha = rowAnnotation(\n", " Sample = metadata$Sample,\n", " Group = metadata[[group_col]],\n", " Cell_type = metadata$CellAnnotation,\n", " Cluster = metadata[[resolution.cluster]],\n", " col = list(\n", " Group = setNames(group_colors, unique(metadata[[group_col]])),\n", " Cluster = setNames(clusters_colors, unique(metadata[[resolution.cluster]])),\n", " Sample = setNames(sample_colors, unique(metadata$Sample)),\n", " Cell_type = setNames(cell_type_colors, unique(metadata$CellAnnotation))\n", " ),\n", " show_annotation_name = TRUE,\n", " show_legend = c(TRUE, TRUE, TRUE, TRUE) # 隐藏Cell_type的图例\n", " )\n", " }\n", "} else {\n", " if(nchar(celltype_col) > 0){\n", " ha = rowAnnotation(\n", " Sample = metadata$Sample,\n", " Cell_type = metadata[[celltype_col]],\n", " Cluster = metadata[[resolution.cluster]],\n", " col = list(\n", " Cluster = setNames(clusters_colors, unique(metadata[[resolution.cluster]])),\n", " Sample = setNames(sample_colors, unique(metadata$Sample)),\n", " Cell_type = setNames(cell_type_colors, unique(metadata[[celltype_col]]))\n", " ),\n", " show_annotation_name = TRUE,\n", " show_legend = c(TRUE, TRUE, TRUE) # 隐藏Cell_type的图例\n", " )\n", " } else {\n", " ha = rowAnnotation(\n", " Sample = metadata$Sample,\n", " Cell_type = metadata$CellAnnotation,\n", " Cluster = metadata[[resolution.cluster]],\n", " col = list(\n", " Cluster = setNames(clusters_colors, unique(metadata[[resolution.cluster]])),\n", " Sample = setNames(sample_colors, unique(metadata$Sample)),\n", " Cell_type = setNames(cell_type_colors, unique(metadata$CellAnnotation))\n", " ),\n", " show_annotation_name = TRUE,\n", " show_legend = c(TRUE, TRUE, TRUE) # 隐藏Cell_type的图例\n", " )\n", " }\n", "}\n", "\n", "# 6. 绘制热图\n", "ht <- Heatmap(\n", " cnv_matrix,\n", " name = \"Copy Number\",\n", " col = colormap,\n", " left_annotation = ha,\n", " show_row_names = FALSE,\n", " cluster_rows = TRUE,\n", " cluster_columns = FALSE,\n", " row_split = metadata[,2],\n", " use_raster = TRUE,\n", " row_title = NULL,\n", " column_names_rot = 90,\n", " show_row_dend = FALSE,\n", " border = FALSE,\n", " column_names_gp = gpar(fontsize = 8),\n", " heatmap_legend_param = list(\n", " title = \"Copy number\",\n", " at = c(1, 1.5, 2, 2.5, 3),\n", " legend_height = unit(4, \"cm\"),\n", " grid_width = unit(0.5, \"cm\")\n", " )\n", ")\n", "# 8. 保存图片\n", "invisible(pdf(\"./result/cnv_heatmap_annotated.pdf\", width = 15, height = 10))\n", "invisible(draw(ht))\n", "invisible(dev.off())\n", "draw(ht)\n", "options(warn = 0)" ] }, { "cell_type": "markdown", "id": "00e5a4c7-ea65-4aab-a388-413d74782d4f", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "### 可视化每条染色体的 CNV score" ] }, { "cell_type": "markdown", "id": "ca7ca79f-81a2-4502-b6ae-805f85157286", "metadata": {}, "source": [ "对每个染色体区域进行 CNV 评分分级（<1.5 为缺失（1），1.5-2.5 为正常（2），>2.5 为扩增（3）），通过将每个细胞每条染色体上的 CNV 分数映射到 UMAP 图中，可以直观看出哪些细胞簇哪条染色体发生了明显的染色体拷贝数变异。" ] }, { "cell_type": "code", "execution_count": null, "id": "edd38e75-c5b6-4a41-8a1e-fbe567933781", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [ "remove-input" ], "vscode": { "languageId": "r" } }, "outputs": [], "source": [ "options(warn = -1)\n", "suppressWarnings({\n", " suppressMessages({\n", "\n", " celltype_colors <- c(\n", " \"1\" = \"#7A94F4\", # 蓝色\n", " \"2\" = \"#E4EAEA\", # 橙色\n", " \"3\" = \"#F769A3\"\n", " )\n", " library(RColorBrewer)\n", " library(scales) \n", " base_colors <- c(\"#EDD496\",\"#65D1BF\",\"#7ADAEA\",\"#F2AFAF\",\"#EA8F8F\",\n", " \"#C1A9D3\",\"#F2DA9E\",\"#D9EF78\",\"#EF9ED9\",\"#C9BCC6\",\"#8BD198\")\n", "\n", " generate_cols <- function(n) {\n", " if (n <= 0) return(character(0))\n", " if (n <= length(base_colors)) {\n", " base_colors[1:n]\n", " } else {\n", " colorRampPalette(base_colors)(n)\n", " }\n", " }\n", " # 创建目录保存结果\n", " output_dir <- file.path(\"./result\", \"cnv_plots\")\n", " dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)\n", "\n", " # 定义保存和编码图像的函数\n", " save_and_encode_plot <- function(plot, filename, width, height) {\n", " # 保存图片文件\n", " full_path_pdf <- paste0(filename, \".pdf\")\n", " temp_png <- tempfile(fileext = \".png\")\n", " \n", " # 保存PDF\n", " ggsave(full_path_pdf, plot, width = width, height = height)\n", " # 临时保存PNG用于base64编码\n", " ggsave(temp_png, plot, width = width, height = height, dpi = 300)\n", " \n", " # 将图片转换为base64\n", " png_data <- readBin(temp_png, \"raw\", file.info(temp_png)$size)\n", " base64_img <- base64encode(png_data)\n", " \n", " # 删除临时PNG文件\n", " unlink(temp_png)\n", " \n", " return(list(\n", " base64 = base64_img,\n", " pdf_path = full_path_pdf\n", " ))\n", " }\n", " # 生成细胞类型注释图颜色\n", " if (nchar(celltype_col) > 0) {\n", " unique_celltype <- sort(unique(sc_filter@meta.data[[celltype_col]]))\n", " clusters_cols <- generate_cols(length(unique_celltype))\n", " identity_colors <- setNames(clusters_cols, unique_celltype)\n", " celltype_umap <- DimPlot(sc_filter, \n", " group.by = celltype_col, \n", " cols = identity_colors,\n", " label = TRUE,\n", " reduction = reduction) + NoLegend()\n", " } else {\n", " unique_celltype <- sort(unique(sc_filter@meta.data[[\"CellAnnotation\"]]))\n", " clusters_cols <- generate_cols(length(unique_celltype))\n", " identity_colors <- setNames(clusters_cols, unique_celltype)\n", " celltype_umap <- DimPlot(sc_filter,\n", " group.by = \"CellAnnotation\",\n", " cols = identity_colors,\n", " label = TRUE,\n", " reduction = reduction) + NoLegend()\n", " }\n", " # 生成聚类分群图颜色\n", " unique_clusters <- sort(unique(sc_filter@meta.data[[resolution.cluster]]))\n", " clusters_cols <- generate_cols(length(unique_clusters))\n", " identity_colors <- setNames(clusters_cols, unique_clusters)\n", " cluster_umap <- DimPlot(sc_filter,\n", " group.by = resolution.cluster,\n", " cols = identity_colors,\n", " label = TRUE,\n", " reduction = reduction) + NoLegend()\n", " plot_filename2 <- file.path(output_dir, paste0(resolution.cluster, \"_umap\"))\n", " plot_filename3 <- file.path(output_dir, paste0(\"celltype\", \"_umap\"))\n", " result2 <- save_and_encode_plot(\n", " plot = cluster_umap,\n", " filename = plot_filename2,\n", " width = 5,\n", " height = 4\n", " )\n", " result3 <- save_and_encode_plot(\n", " plot = celltype_umap,\n", " filename = plot_filename3,\n", " width = 5,\n", " height = 4\n", " )\n", " # 创建数据框用于展示\n", " image_df <- data.frame(\n", " Chromosome = character(),\n", " CNV_UMAP = character(),\n", " cell_culster=character(),\n", " celltype_culster=character(),\n", " stringsAsFactors = FALSE\n", " )\n", " # 初始化存储列表\n", " char <- list()\n", " # 遍历每个染色体列\n", " for(i in colnames(cnv_matrix)) {\n", " # 处理单个染色\n", " char[[i]] <- cnv_data %>%\n", " dplyr::select(rowname, i) %>%\n", " mutate(\n", " copy_number = case_when(\n", " .data[[i]] < 1.5 ~ \"1\",\n", " .data[[i]] >= 1.5 & .data[[i]] <= 2.5 ~ \"2\",\n", " .data[[i]] > 2.5 ~ \"3\"\n", " )\n", " )\n", " colnames(char[[i]]) <- c(\"Barcodes\", paste0(i,\"_scores\"), i)\n", " rownames(char[[i]]) <- char[[i]]$Barcodes\n", " char[[i]] <- char[[i]][,-1]\n", " sc_filter <- AddMetaData(sc_filter, metadata=char[[i]])\n", " \n", " # 准备绘图数据\n", " max.col=max(cnv_matrix)\n", " #p=FeaturePlot(sc_filter, features=paste0(i,\"_scores\"),min.cutoff = 0,max.cutoff = max.col,raster=T)\n", " p=DimPlot(sc_filter, group.by=i,cols=celltype_colors,reduction = reduction)\n", " # 保存图片并获取base64编码\n", " plot_filename <- file.path(output_dir, paste0(i, \"_umap\"))\n", " result <- save_and_encode_plot(\n", " plot = p,\n", " filename = plot_filename,\n", " width = 5,\n", " height = 4\n", " )\n", "\n", " # 添加到数据框\n", " img_tag <- sprintf(\n", " '\n", " \n", " ',\n", " result$pdf_path,\n", " result$base64\n", " )\n", " \n", " img_tag2 <- sprintf(\n", " '\n", " \n", " ',\n", " result2$pdf_path,\n", " result2$base64\n", " )\n", " img_tag3 <- sprintf(\n", " '\n", " \n", " ',\n", " result3$pdf_path,\n", " result3$base64\n", " )\n", " image_df <- rbind(image_df, data.frame(\n", " Chromosome = i,\n", " CNV_UMAP = img_tag,\n", " cluster_UMAP = img_tag2,\n", " celltype_UMAP = img_tag3\n", " ))\n", " \n", " # 也保存普通PDF（与原代码相同）\n", " #ggsave(paste0(outdir, i, \"_umap.pdf\"), plot = p, width = 5, height = 4)\n", " }\n", " \n", " # 创建交互式表格展示结果\n", " # 6. 创建交互式表格\n", " cnv_table <- DT::datatable(\n", " image_df,\n", " escape = FALSE,\n", " options = list(\n", " pageLength = 1,\n", " autoWidth = TRUE,\n", " scrollX = TRUE,\n", " dom = 'Bfrtip',\n", " server = TRUE, # 强制服务器端模式\n", " deferRender = TRUE, # 延迟渲染\n", " scroller = TRUE, # 滚动加载\n", " columnDefs = list(\n", " # 仅针对前两列（索引0和1)\n", " list(targets = 0, width = \"4%\"),\n", " list(targets = 1, width = \"32%\"),\n", " list(targets = 2, width = \"32%\"),\n", " list(targets = 2, width = \"32%\")\n", " )\n", " ),\n", " rownames = FALSE,\n", " filter = 'top',\n", " caption = htmltools::tags$caption(\n", " style = 'caption-side: top; text-align: center; font-size: 1.2em; font-weight: bold;',\n", " '染色体区域CNV拷贝数UMAP可视化'\n", " )\n", " )\n", " \n", " # 显示表格\n", " cnv_table\n", " })\n", "})\n", "options(warn = 0)" ] }, { "cell_type": "markdown", "id": "60694d7d-d5d1-49b5-946a-e0571778f28a", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## 结果文件\n", "\n", "├── ***cnv_plots***
\n", "│ ├── chr*_umap.pdf 分别展示每条染色体臂的 cnv 变异情况
\n", "└── cnv_heatmap_annotated.pdf 拷贝数变异热图展示
\n", "\n", "结果目录：目录下包括此项分析所涉及的所有图片以及表格等结果文件\n" ] }, { "cell_type": "markdown", "id": "cb5c527e-89ad-4adb-bbf2-8ec24180df29", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## 文献案例解析\n", "\n", "* 《Single-cell multi-omics sequencing uncovers region-specific plasticity of glioblastoma for complementary therapeutic targeting》\n", " \n", "* 《Single-cell landscape of primary central nervous system diffuse large B-cell lymphoma》\n", "\n", " 在这项研究中，作者用 inferCNV 对 scRNA-seq 做拷贝数变异分析，用 Copy-scAT 对 scATAC-seq 做拷贝数变异分析，用于区分肿瘤细胞和非肿瘤细胞\n", "\n", "\n", "\n", "## 参考资料\n", "\n", "[1] Nikolic A, Singhal D, Ellestad K, et al.Copy-scAT: Deconvoluting single-cell chromatin accessibility of genetic subclones in cancer[J].Science advances, 2021, 7(42):eabg6045.DOI: 10.1126/sciadv.abg6045." ] }, { "cell_type": "markdown", "id": "059455e6-b826-4582-b83d-aa42d6a4594c", "metadata": { "editable": true, "slideshow": { "slide_type": "" }, "tags": [] }, "source": [ "## 附录\n", "\n", ".txt ：结果数据表格文件，文件以制表符（Tab）分隔。unix/Linux/Mac 用户使用 less 或 more 命令查看；windows 用户使用高级文本编辑器 Notepad++ 等查看，也可以用 Microsoft Excel 打开。\n", "\n", ".pdf ：结果图像文件，矢量图，可以放大和缩小而不失真，方便用户查看和编辑处理，可使用 Adobe Illustrator 进行图片编辑，用于文章发表等。\n", "\n", ".rds : 包含基因活力 assay 的 Seurat 对象，需要在 R 环境中打开查看和用于进一步分析。" ] } ], "metadata": { "kernelspec": { "display_name": "R.4.3.3", "language": "R", "name": "r.4.3.3" }, "language_info": { "codemirror_mode": "r", "file_extension": ".r", "mimetype": "text/x-r-source", "name": "R", "pygments_lexer": "r", "version": "4.3.3" } }, "nbformat": 4, "nbformat_minor": 5 }