infercnv的cluster_by_groups参数影响后续层次聚类文件读取函数

生信技能树

发布于 2022-03-03 13:27:31

1.4K0

发布于 2022-03-03 13:27:31

文章被收录于专栏：生信技能树

因为教程跨越了不同时间周期，软件更新，数据集的特异性，导致很多小伙伴follow不同系统的教程会得到不一样的报错。

一般来说，我会使用500个两种不同的正常血液细胞作为inferCNV算法的对照，然后在被计算拷贝数的上皮细胞里面混入300个两种不同的正常血液细胞作为控制条件，代码如下所示：

dat=cbind(epiMat,TcellsMat,monoMat)
# 其中 sce 是除去了两种不同的正常血液细胞 后的其它细胞，被预测拷贝数情况
groupinfo=data.frame(v1=colnames(dat),
                     v2=c( sce$celltype ,
                          rep('spike-Tcells',300),
                          rep('ref-Tcells',500),
                          rep('spike-mono',300),
                          rep('ref-mono',500)))

因为被预测的上皮细胞有可能是分群了的，其实是除去了两种不同的正常血液细胞后的所有的其它细胞，有分群很正常。

所以我们在运行了infercnv流程的时候，会设置 cluster_by_groups=T ：

library(infercnv)
infercnv_obj = CreateInfercnvObject(raw_counts_matrix=expFile,
annotations_file=groupFiles,
delim="\t",
gene_order_file= geneFile,
ref_group_names=c('ref-Tcells',
'ref-mono'))  ## 这个取决于自己的分组信息里面的


# cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomics
infercnv_obj2 = 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= "infercnv_output",  # dir is auto-created for storing outputs
cluster_by_groups=T ,   # cluster
hclust_method="ward.D2", 
plot_steps=F)

以前我们的教程里面，其实并没有设置 cluster_by_groups=T ，因为我们是对全部的上皮细胞进行预测，这个时候并不会对它进行分群。所以也没必要去 cluster_by_groups=T ，但是最近我发现很多其它非上皮细胞也是有可能有拷贝数的，所以就预测了除去了两种不同的正常血液细胞后的所有的其它细胞，所以就设置 cluster_by_groups=T 。

这样得到的 inferCNV 的 dendrogram文件就不能使用之前的代码读取：

 infercnv.dend <- read.dendrogram(file = "inferCNV_output/infercnv.observations_dendrogram.txt")

而是需要使用ape包的read.tree函数，如下所示：

myTree <- ape::read.tree(file = "inferCNV_output/infercnv.observations_dendrogram.txt")
u

我们可以看到读入的 inferCNV 的 dendrogram文件其实是9个内容：

> length(myTree)
[1] 9
> unlist(lapply(myTree, function(x){length(x$tip.label)}))
[1]  23 858  57  19  71  10  14 300 300

但是如果我们看前面的细胞分组信息，如下所示：

> as.data.frame(table(gp$V2))
           Var1 Freq
1        Bcells   14
2      CD24-epi   71
3       cyc-epi   57
4          endo   10
5      ESR1-epi    1
6         fibro   23
7       lum-epi  858
8        plasma   19
9      ref-mono  500
10   ref-Tcells  500
11          SMC    2
12   spike-mono  300
13 spike-Tcells  300

应该是13个细胞类型，其中两个 'ref-Tcells' 和 'ref-mono'是另外的文件，不会在infercnv.observations_dendrogram.txt里面，而ESR1-epi 和 SMC 这两个细胞亚群因为细胞数量太少，直接被过滤了。所以就是读入的 inferCNV 的 dendrogram文件的9个内容。

如果你确实觉得我的教程对你的科研课题有帮助，让你茅塞顿开，或者说你的课题大量使用我的技能，烦请日后在发表自己的成果的时候，加上一个简短的致谢，如下所示：

We thank Dr.Jianming Zeng(University of Macau), and all the members of his bioinformatics team, biotrainee, for generously sharing their experience and codes.

十年后我环游世界各地的高校以及科研院所（当然包括中国大陆）的时候，如果有这样的情谊，我会优先见你。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-02-15，如有侵权请联系 cloudcommunity@tencent.com 删除