ixxmu
commented
11 months ago infercnv 三回首：深入理解infercnv为何能发nature by 生信小博士
大家好，不知你是否还记得，前两次关于infercnv的介绍。请看这里：
如果下载了示例数据，并且你已经跑了上述代码，不难得到这张图：
高考作文问你对这张图的理解。
01
关于我对infercnv的理解：其实结果是哪种图真的不重要，重点是对图形的解读。比如本文中出现的第一张图，里面有个绿色圈圈，可以看到柱子里深灰和浅灰部分（就是我标注的A 和B部分）的突变图景肉眼几乎不可分辨，但是却被分属于两群不同的细胞群。
这就说明：通过肉眼观察，从A和B的突变图景上来看，A和B的突变图景是一致的，那么我们就可以把A亚群细胞与B亚群细胞合并，把他们当作恶性细胞。
但是 柱子里所有的深灰和浅灰 所对应的组别别是seurat分群里面的2群和4群
这就提示我们本文中的seurat分群结果并不一定能把恶性细胞单纯分出来
（也许我们在seurat里可以把所有分辨率都试一下，理论上应该可以分出恶性细胞）。所以我们要对seurat的分群结果进行重新命名，2群和4群中我画出来的A和B是恶性细胞，但是2群和4群中剩下的其他部分细胞不是。
其实，你仔细去看这张图的话，最左边还有一根柱子，这个柱子不同的颜色对应着infercnv突变图景的聚类结果。你会发现最左边的柱子上存在A1 和B1部分（用红字标出来的部分）正好对应着A1和B1部分的细胞。这里的A1 和B1的聚类结果就是我们想要定义的恶性细胞。我们可以通过代码获取A1和B1部分细胞id（这是一个思考题，如何获取细胞id见文末），进而命名为恶性细胞。
就本数据集而言，通过这样的分析，大家能明白为啥我们要做cnvscore了吗：因为单纯通过seruat分群，只能找出2群和4群为恶性细胞。但实际上只有2群和4群中的A部分和B部分为恶性细胞。 所以说infercnv备受CNS青睐，是因为它精准的预测恶性细胞的能力，而这种能力是普通的降维聚类分群所达不到的(感兴趣的同学可以在seurat中把分辨率调大一点再试试)。
致敬！
Anoop P. Patel, Itay Tirosh, et al. Single-cell RNA-seq highlights intratumoral heterogeneity in primary glioblastoma. Science. 2014 Jun 20: 1396-1401
Tirosh I et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016 Apr 8;352(6282):189-96
Tirosh I et al. Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. Nature. 2016 Nov 10;539(7628):309-313. PubMed PMID: 27806376; PubMed Central PMCID: PMC5465819.
Venteicher AS, Tirosh I, et al. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science. 2017 Mar 31;355(6332).PubMed PMID: 28360267; PubMed Central PMCID: PMC5519096.
Puram SV, Tirosh I, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell. 2017 Dec 14;171(7):1611-1624.e24. PubMed PMID: 29198524; PubMed Central PMCID: PMC5878932.
02
ok，我这里只是稍带了一下infercnv的厉害之处，下面我们进入第二个问题，你知道为什么为什么有的文章会使用另外一种图吗？：
大家可以看到最左边的柱子，本文第一个图是多色的，而这个只有靛青色。你知道这分别代表什么意思吗？欢迎来稿，说出你的见解，下期见分晓~。
2.5
有的时候还会出现下面这种图，看着就很诡异：
但是真不用害怕，这些图都是正常的，这是由于参数的原因！你要学会找到你到底需要哪一张图。
出现各种奇怪的图，基本就两个原因：1.参数原因（cluster_by_groups和leiden_resolution参数特别重要） 2.细胞数量原因 。详情见官方链接https://www.youtube.com/watch?v=-qOcHAavZT8
03
下面列举七种不同参数出的不同结果
当然你可以在每个代码里面都加上 denoise=T参数 进行降噪处理
接下来就展示一下不同的代码参数，所获得的突变图景，稍有差别。
0
如果没有cluster_by_groups=T参数的话,图形丑且不能用。
# cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsinfercnv_obj2 = infercnv::run(infercnv_obj,cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsout_dir= "infercnv_output_noclusterbygroup", # dir is auto-created for storing outputs# cluster_by_groups=T , # clusterhclust_method="ward.D2", plot_steps=F)
后面都加上这个参数：
1
# cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsinfercnv_obj2 = infercnv::run(infercnv_obj,cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsout_dir= "infercnv_output", # dir is auto-created for storing outputscluster_by_groups=T , # clusterhclust_method="ward.D2", plot_steps=F)
2
infercnv_obj2 = infercnv::run(infercnv_obj,cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsout_dir= "infercnv_output2", # dir is auto-created for storing outputscluster_by_groups=T , # clusterleiden_resolution = 0.01,hclust_method="ward.D2", plot_steps=F)
3
infercnv_obj2 = infercnv::run(infercnv_obj,cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsout_dir= "infercnv_outputres0.6", # dir is auto-created for storing outputscluster_by_groups=T , # clusterleiden_resolution = 0.6,hclust_method="ward.D2", plot_steps=F)
4
infercnv_obj2 = infercnv::run(infercnv_obj,cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsout_dir= "infercnv_outputres0.01", # dir is auto-created for storing outputscluster_by_groups=T , # clusterleiden_resolution = 0.01,HMM = TRUE,plot_steps=F)#subcluster
5
infercnv_obj2 = infercnv::run(infercnv_obj,cutoff=0.1, # use 1 for smart-seq, 0.1 for 10x-genomicsout_dir="infercnv_output_mycodes", # dir is auto-created for storing outputsanalysis_mode="subclusters", # 更好方法，按照 cnv 亚群进行预测，而不是样本。能更好的区分肿瘤异质性。cluster_by_groups=T, # 根据细胞类型对肿瘤细胞执行单独的聚类，如 cell annotations 文件中定义的那样denoise=T, # 去噪处理num_threads=18, # 设置线程数, 多线程运行，加快计算速度HMM=T,# hclust_method="ward.D2",plot_steps=F)
6
infercnv_obj2 = infercnv::run(infercnv_obj,cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsout_dir= "infercnv_outputres0.1", # dir is auto-created for storing outputscluster_by_groups=T , # clusterleiden_resolution = 0.1,#hclust_method="HMM",plot_steps=F)
7
infercnv_obj2 = infercnv::run(infercnv_obj,cutoff=0.1, # cutoff=1 works well for Smart-seq2, and cutoff=0.1 works well for 10x Genomicsout_dir= "infercnv_outputres0.2", # dir is auto-created for storing outputscluster_by_groups=T , # clusterleiden_resolution = 0.2,HMM = TRUE, plot_steps=F)
04
题外话：如果超过5k个细胞，不建议用r跑程序。如果超过1w个恶性细胞，普通电脑（我的OMEN游戏本）根本带不起来哈。2w个恶性细胞的话，普通服务器带不起来哈，得用高性能计算。当然你可以设置HMM=F，这样内存消耗会小很多，然后再试试。
如果使用python的话，应该会快很多，可能内存要求不那么高，我目前还没有试验~。哪天来个超大数据集，说不定会逼得我用python。
本文中的思考题，提示：如果想获取A1和B1部分细胞id，你可以检查以下圈出来的这两个文件，可以根据里面的行名获取A1和B1部分细胞id。欢迎来稿，使用详细代码写出如何获取这些恶性细胞的id。
https://mp.weixin.qq.com/s/IujcPhKM1SHKgPLTwGyZfg