单细胞转录组分析之CNV推断

[ixxmu]

https://mp.weixin.qq.com/s/O5xL-QQXROSlWfkVeioqRw

[ixxmu]

单细胞转录组分析之CNV推断 by 医学僧的科研日记

单细胞转录组公共数据库已经很多了，但单细胞多组学的数据大多数还没公开，如何利用已有工具让我们的研究更出彩成为了大家亟需解决的问题。在这里，继上篇单细胞转录组分析之细胞通讯，小王子给大家带来新的模块-CNV推断。

1.什么是CNV？

拷贝数变异（Copy number variation, CNV）：是由基因组发生重排而导致的,一般指长度1KB以上的基因组大片段的拷贝数增加或者减少, 主要表现为亚显微水平的缺失和重复。

CNV的发生机制就是非等位基因重组，第一次是在减数第一次分裂前期，一对同源染色体染色体上的非姐妹染色单体交叉互换，第二次是在减数第一次分裂后期，同源色体分离，非同源染色体自由组合。基因组上非等位的两个高度同源的DNA序列在减数分裂或者有丝分裂的过程中发生错误的配对，并发生序列交换，从而导致缺失、重复的出现。

PS：对背景知识不理解的小伙伴建议去知乎搜索专栏：焦老师讲遗传系列之12：拷贝数变异CNV

2.文献中怎么用CNV推断？

一篇21年的揭示转移性胰腺导管腺癌瘤内异质性的文章：

在进行细胞注释后，作者用inferCNV在所测的5个样本中依次进行了不同胰腺导管细胞亚群的CNV推断，发现不同样本具有显著的异质性。

3.该如何进行CNV推断？

由于篇幅有限，这里着重介绍下常用的inferCNV和copyKAT。

1.inferCNV

inferCNV前后一共有3篇CNS文章进行开发和改进，在这里，小王子分别给大家简述其演进过程以方便大家理解inferCNV进行CNV推断的原理：

(1)提出：Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014 Jun. PMID: 24925914a

第一版整体来讲的原理就是：使用101个基因的平均表达值表示中间基因的染色体CNV值，流程如下：

图1：原始数据标准化，去除测序深度不同造成的差异；

图2：根据公式把基因表达值转换为CNV值；

图3：对每个细胞的CNV值进行中心化处理，使细胞之间相同染色体区域的CNV值具有可比性；

图4：用肿瘤细胞的CNV值减去对应区域正常细胞的CNV；

蓝框左图：如果infercnv::run函数中的参数denoise=TRUE，则使用算法进一步去除背景噪音凸显CNV区域；

蓝框右图：如果infercnv::run函数中的参数HMM=TRUE，则使用隐马尔可夫模型（Hidden Markov Model, HMM）预测CNV区域，并用贝叶斯潜在混合模型（Bayesian Network Latent Mixture Model）对结果进行校正。

(2)改进：Dissectingthe multicellular ecosystem of metastatic melanoma by single-cell RNA-seq.Science. 2016 Apr 8.PMID: 27124452

第二版与第一版相比，有了以下两处改进:A. 为避免极端值（-3，-3）；B.5种CNV稳定的正常细胞作为基线值进行CNV矫正，用整体CNV和与同一肿瘤top10%细胞的平均CNV的相关性估计细胞的良恶性。

到第二版，inferCNV已经能够判断良恶性细胞了

(3)细化+升华：Single-CellTranscriptomic Analysis of Primary and Metastatic Head and Neck Cancer. Cell.2017 Dec 14.PMID: 29198524

与第二版使用同一肿瘤top10%细胞的平均CNV的相关性相比，第三版有了如下改进：

A.用整体CNV和与同一肿瘤所有细胞的平均CNV的相关性估计细胞的良恶性。

B.给定了以下判断良恶性的相关性参考阈值：

恶性细胞：整体CNV>0.05&CNV相关性>0.5

非恶性细胞：整体CNV<0.05&CNV相关性<0.5

未定义细胞：整体CNV>0.05&CNV相关性<0.5或者整体CNV>0.05&CNV相关性<0.5

#代码实战：比较简单，只需准备好count矩阵表达文件、#单个细胞类群的注释文件和基因在染色体上的位置文件即可，#运行速度看电脑配置和细胞数量，我示例是3317个细胞跑了20分钟。rm(list = ls())#安装infercnv包BiocManager::install("infercnv")#使用infercnv自带数据-----------------------------------------------------------------------library(infercnv)library(tidyverse)dir.create("inferCNV")dir.create("inferCNV/test")##读取示例数据目录exprMatrix = system.file("extdata", "oligodendroglioma_expression_downsampled.counts.matrix.gz", package="infercnv")cellAnnota = system.file("extdata", "oligodendroglioma_annotations_downsampled.txt", package="infercnv")geneLocate = system.file("extdata", "gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt", package="infercnv")#创建inferCNV对象，直接给相应的文件路径即可infercnv_obj = CreateInfercnvObject(delim = '\t',                                    raw_counts_matrix = exprMatrix,                                    annotations_file = cellAnnota,                                    gene_order_file = geneLocate,                                    ref_group_names = c("Microglia/Macrophage","Oligodendrocytes (non-malignant)"))
##分析细胞CNV#cutoff阈值，官网建议Smart-seq2数据选“1”, 10x Genomics数据选“0.1”infercnv_obj = infercnv::run(infercnv_obj,                             cutoff=1,                              out_dir='inferCNV/test',                              cluster_by_groups=TRUE,                             #analysis_mode="subclusters", #默认是"samples"                             denoise=TRUE,                             HMM=TRUE)library(RColorBrewer)library(ggsci)infercnv::plot_cnv(infercnv_obj,                   plot_chr_scale = T,                   output_filename = "inferCNV/test/better_plot",output_format = "pdf",                   custom_color_pal =  color.palette(c(pal_nejm()(8)[1],"white",pal_nejm()(8)[2]), c(2,2)))#使用TISCH中的Seurat(LIHC的GSE125449)数据-----------------------------------------------------------------------rm(list = ls())library(infercnv)library(tidyverse)library(Seurat)library(tidyverse)dir.create("inferCNV/valudation")#获得Count矩阵data <- Read10X_h5('LIHC_GSE125449_aPDL1aCTLA4_expression.h5')data <- CreateSeuratObject(data, project = "", min.cells = 3, min.features = 200)data #3835个细胞，17356个基因meta <- read.delim("LIHC_GSE125449_aPDL1aCTLA4_CellMetainfo_table.tsv", header=F)colnames(meta) <- meta[1,]meta <- meta[-1,]rownames(data)colnames(meta)colnames(meta) <- c("Cell","UMAP_1","UMAP_2","Celltype_malignancy",                      "Celltype_major_lineage","Celltype_minor_lineage",                    "Cluster","Celltype","Treatment","Patient","Sample",                    "Source","Age","Gender","Stage")data@assays$RNA@meta.features <- as.data.frame(rownames(data))identical(meta$Cell,colnames(data))meta <- meta[meta$Cell%in%colnames(data),]a <- data@meta.dataa <- cbind(a,meta)data@meta.data <- atable(data@meta.data$Celltype_malignancy)#数据质控+去除批次效应 ---------------------------------------------------------data[["percent.mt"]] <- PercentageFeatureSet(data, pattern = "^MT-") # 注意基因名称变化VlnPlot(data, group.by = "Patient", features = c("nFeature_RNA", "nCount_RNA", "percent.mt"), pt.size = 0, ncol = 3) +   NoLegend()##筛选合格细胞cat("Before filter:",nrow(data@meta.data),"cells\n")data <- subset(data, subset =                  nFeature_RNA > 200 &                  nFeature_RNA < 3000 &                  nCount_RNA > 500 &                  nCount_RNA < 4000 &                 percent.mt < 5)cat("After filter :",nrow(data@meta.data),"cells\n")data #保留17356个基因，3317个细胞。save(data,file = "GSE125449_count.rda")
dd <- data[["RNA"]]@countsdd2 <- data@meta.data  dd <- as.data.frame(dd)dd2 <- dd2[,c(4,7)]table(dd2$Celltype_malignancy)dd2$Celltype_malignancy <- factor(dd2$Celltype_malignancy,levels = c("Immune cells","Stromal cells","Malignant cells"))
load('GRCh38_v39.gtf.rda') #参考基因组，可自行下载dd3 <- gtf_data[,c(12,1:3)]dd3 <- distinct(dd3,gene_name,.keep_all = T)
gene <- intersect(rownames(dd),dd3$gene_name)dd <- dd[rownames(dd)%in%gene,]dd3 <- dd3[dd3$gene_name%in%gene,]
write.table(as.matrix.data.frame(dd), file ="inferCNV/valudation/counts.matrix",quote = F,sep = "\t",row.names = T,col.names = T)write.table(as.data.frame(dd3),file = "inferCNV/valudation/gencode.txt", quote = F,sep = "\t",row.names = F,col.names = F)write.table(as.data.frame(dd2),file = "inferCNV/valudation/annotation.txt", quote = F,sep = "\t",row.names = F,col.names = F)#运行infercnv-----------------------------------------------------------------------rm(list = ls())library(tidyverse)library(infercnv)memory.limit(10000000)#使用inferCNV之前，最好过滤掉低质量的细胞infercnv_obj = CreateInfercnvObject(raw_counts_matrix="inferCNV/valudation/counts.matrix",                                    annotations_file="inferCNV/valudation/annotation.txt",                                    delim="\t",                                    gene_order_file="inferCNV/valudation/gencode.txt",                                    ref_group_names=c("Immune cells") )#ref_group_names参数根据细胞注释文件填写，在本数据，我们假设用“Immune cells”作为参照，看“Stromal cells”和“Malignant cells”的恶性；#ref_group_names=NULL，则会选用所有细胞的平均表达作为参照，这种模式下，最好能确保细胞内在有足够的差异infercnv_obj = infercnv::run(infercnv_obj,                             cutoff=1,                              out_dir="try2",                             cluster_by_groups=TRUE,                              #analysis_mode="subclusters", #默认是"samples"                             denoise=TRUE,                             HMM=TRUE,                             num_threads=4)#去噪，HMM预测CNV这两项我一般都选上

2.copyKAT

与inferCNV相比，copyKAT是专门为单细胞CNV开发的，不需要给定正常参考细胞，总体效果好于inferCNV。

copyKAT主要原理：贝叶斯方法+层次聚类

a:首先输入ScRNA-seq的UMI表达矩阵，通过基因组坐标对它们进行排序，之后用Freeman-Tukey变换和多项式动态线性建模矫正矩阵中的异常值；

b:之后，将所有单细胞集中到几个小的亚群分类中，并使用高斯混合模型估算每个分类的方差，具有最小估计方差的聚类被定义为“标准的二倍体细胞”。

c:为了检测染色体断点，他们整合泊松-伽玛模型和马尔可夫链蒙特卡罗迭代生成每个基因窗口的后验均值，然后应用Kolmogorov-Smirnov检验对均值无显著差异的相邻窗口进行合并，然后计算每个窗口的最终拷贝数值，以此作为跨越每个细胞中相邻染色体断点的的所有基因的后验平均值。

#copyKAT-----------------------------------------------------------------------rm(list = ls())
#install_github("navinlabcode/copykat")library(copykat)load("GSE125449_count.rda")counts <- as.matrix(data@assays$RNA@counts)cnv <- copykat(rawmat = counts,ngene.chr = 5,               sam.name = "GSE125449_copykat",n.cores = 1)save(cnv,file="GSE125449_cnv_copyKAT.rda"

示例结果展示：可以直接识别出输入的细胞是非整倍体（恶性）、整倍体（良性）或者未识别

需要指出的是：单细胞转录组推断CNV用的是Count数据。不过很多数据库中的数据都已经经过log化，有文献提示利用log后的数据进行inferCNV和copyKAT分析对结果并没影响，不过还是建议用Count数据，因为它会自动进行log标准化。

小结

总的来说，单细胞转录组的CNV推断就是利用转录组的count表达矩阵，根据前后一些基因的表达推断该基因的CNV情况，基本原理就是CNV是以片段为单位的扩增和缺失，一旦该基因发生改变，周围的相邻基因很可能“同步”发生改变，论文中经常使用CNV推断判断细胞的良恶性，感兴趣的同学可以现学现用，立刻就能增加到你的研究中来！

本文主要参考：

1.知乎-焦老师讲遗传系列之12：拷贝数变异CNV

2.腾讯云-单细胞转录组高级分析四：scRNA数据推断CNV

3.简书-单细胞转录组鉴定肿瘤细胞：CopyKAT和InferCNV

4.简书-单细胞CNV分析-copyKAT

PS:我们在B站的同名账号“医学僧的科研日记”同步有公众号的讲解视频~目前已有10w余粉丝，感兴趣的小伙伴可以多多支持幺。

医学僧的心愿是为医学僧的支持者和科研爱好者营造一个勤学好问、互帮互助、共同进步的优良环境。但是人多了难免鱼龙混杂，良莠不齐，群里不乏有真心实意、脚踏实地的科研人，我们很尊重并欢迎他们，因此不希望被少数来发广告的、划水的、看热闹的、占个坑的所打扰，最终我们的心意付之东流，大家的热情消磨殆尽。为了解决这个问题，我们决定提高入群的门槛——收费进群，这样，愿意进群的同志相信也是个圈内活跃用户，圈内人相互交流，岂不快哉！

因此，为了提升大家的效率，调动大家的积极性，吸收真正的科研爱好者，形成一个良好的学习交流氛围，我们决定增加群的数量，并细化群的功能，暂时增加以下几个主题群：

• 文献下载-医学僧

• 机器学习（分类）-医学僧

• ggplot可视化-医学僧

• 热图可视化-医学僧

• TCGA挖掘-医学僧

• GEO挖掘-医学僧

• 肿瘤预后模型-医学僧

• 肿瘤分子亚型-医学僧

欢迎小伙伴们进群👏👏👏，每个群收费15元。请联系小编微信yixuesengkyrj或扫描下面二维码：