学习笔记：inferCNV

本来对4vCPUs | 32GB的配置没报太大期望，结果本地花了5个多小时才跑完的infercnv数据，服务器上只用了1个半小时就搞定了（果然是我的本太废柴了）。那就趁此机会把单细胞cnv分析捋捋清楚。
p.s. 不得不说kinesin老师的docker镜像真香，谁用谁知道
重点：本文仅为个人学习笔记，代码参考自健明老师，非原创！！！如有理解错误欢迎批评指正~参考文章如下：

首先，infercnv的计算方法来自于2014的science关于GBM文章（PMID: 24925914），该文章首次利用单细胞转录组数据推断CNV。在Supplementary Materials中有如下描述：

1.png

其实很简单，就是将所有基因根据染色体上的位置进行排序，第i个基因的cnv值就是其上游50个基因和下游50个基因的平均表达值（因此每条染色体开头和结尾的50个基因是无法计算的）。
注意这里的 𝐸𝑘(𝑜𝑗)是relative normalized expression，因此如果是seruat对象，需要GetAssayData(object, slot = "data")，如果直接拿到csv表达矩阵，那就edgeR::cpm一下。
文中的公式可以转换为下面的代码，更易于理解~

cnv <- lapply(51:(nrow(dat)-50), function(i){
    this_cnv <- unlist( lapply(1:ncol(dat), function(j){
      sum(dat[(i-50):(i+50),j])/101
    }))
    return(this_cnv)
  })

2.png

然后这里说他们的分析纳入了6000个基因，如果按照上述100基因为步长的计算方法，因为每条染色体的基因数分布不一样，所以最后拿到的是大概350个可计算的cnv。为了避免受到文库质量的影响造成cnv数值的偏移，他们对数据进行了中心化（就是scale了一下吧）。
文章中作者纳入~6000个基因进行分析，为啥不是所有基因？可能没必要为了变异度不高的基因降低运算速度？那是不是可以根据方差rank进行过滤呢？
其实我自己的分析中只考虑删除了在大多数细胞中都不表达的基因，一般来说会过滤到7-8000个基因左右。

subraw.data[apply(subraw.data,1, function(x) sum(x>1) > 50),]

理解了原理之后，其实走流程并不困难，毕竟我们可以走健明大神铺好的路~流程代码单细胞天地都有我就不贴了，不过下面这个基因位置参考文件的代码太牛了，给AnnoProbe打call，大家感受下：

library(AnnoProbe)
geneInfor=annoGene(rownames(dat),"SYMBOL",'human')
colnames(geneInfor)
geneInfor=geneInfor[with(geneInfor, order(chr, start)),c(1,4:6)]
geneInfor=geneInfor[!duplicated(geneInfor[,1]),]
length(unique(geneInfor[,1]))

运行代码很简单：

infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics
                             out_dir=my_dir,
                             cluster_by_groups=FALSE,
                             denoise=TRUE,
                             HMM=TRUE)

这里需要注意：

官方推荐smart-seq2的cutoff设置为1,10X的cutoff设置为0.1（其实10X这种测序深度用来估算cnv真的合适吗）；
HMM是隐马尔科夫模型，就到这吧，我是真没看懂HMM具体算法，反正你要是时间够就HMM=TRUE，不然的话可以先设置HMM=FALSE看一下大致结果。

写的比较简单，只是想从原理上去理解一下代码，主要是打个卡证明今天有学习。。。

最后编辑于：2020.11.02 17:32:25