阅读 22 SEO

学习笔记:inferCNV

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

  1. https://www.jianshu.com/p/05b728c39eee
  2. https://mp.weixin.qq.com/s/Qns9TCSgNg_CQuwQxQbVnw
  3. https://mp.weixin.qq.com/s/aWM5P244HPjWAmPwqacq-Q
  4. https://mp.weixin.qq.com/s/3zM7sfhX-v1FmKcxC8Ae_w

首先,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)

这里需要注意:

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

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

作者:cnwxy

原文链接:https://www.jianshu.com/p/80ca4d59e25f

文章分类
后端
文章标签
版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌抄袭侵权/违法违规的内容, 请发送邮件至 gxwowoo@163.com 举报,一经查实,本站将立刻删除。
相关推荐