当前位置: 首页>大数据>正文

单细胞分析实录(13)- inferCNV结合UPhyloplot2分析肿瘤进化

这一个分析,我目前只在Nature Communications上面看到过两篇文章,都是2020年发表的。肿瘤进化是单细胞里面做肿瘤内部异质性一个比较创新的角度,我参与的课题中也用到了这个分析。公众号后台回复20210420获取今天的数据和部分代码。

1. 先来看一下例子

Single-cell analysis reveals new evolutionary complexity in uveal melanoma

单细胞分析实录(13)- inferCNV结合UPhyloplot2分析肿瘤进化,第1张

这个图并不复杂,单独看每一个样本的进化树,上方主干对应的CNV事件,是早期就发生的;下方分枝的CNV事件是后期发生的。其中蕴含的逻辑就是含有某个CNV的细胞占比越多,那么这个CNV就发生得越早;含有某个CNV的细胞占比越少,这个CNV就发生得越晚。

Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma

单细胞分析实录(13)- inferCNV结合UPhyloplot2分析肿瘤进化,第2张

这两篇文章在做完这个分析之后,并没有说太多的东西,只是简单说明了一下哪些CNV在主干,哪些在分枝(与前人研究对比),而且还是长短臂水平,没有细化到基因。

2. inferCNV预测CNV,并依据CNV聚类

这次教程用到的测试数据同上上篇一样,肿瘤细胞来源于4个病人。但是为了演示,我假设这些细胞都来源于一个病人,是不同的亚克隆。

library(infercnv)
#1
infercnv_obj = CreateInfercnvObject(raw_counts_matrix="oligodendroglioma_expression_downsampled.counts.matrix",
                                    annotations_file="oligodendroglioma_annotations_downsampled.txt",
                                    delim="\t",
                                    gene_order_file="gencode_downsampled.EXAMPLE_ONLY_DONT_REUSE.txt",
                                    ref_group_names=c("Microglia/Macrophage","Oligodendrocytes (non-malignant)")
                                    )
#2
infercnv_obj = infercnv::run(infercnv_obj,
                             cutoff=1, 
                             out_dir="try2",
                             cluster_by_groups=F, 
                             analysis_mode="subclusters",
                             denoise=TRUE,
                             HMM=TRUE,
                             num_threads=1)

运行的结果,保存在try2文件夹中。
注意两个参数cluster_by_groups=F,以及analysis_mode="subclusters",这个参数最终会将肿瘤细胞分为8个cluster(少数情况是7类,如果实在找不出进一步的差别),每个cluster有各自的CNV模式,如果analysis_mode="samples",则一个样本不同细胞最终预测的CNV模式是唯一的。另外需要注意的是,一般文章放的热图是去噪后的热图,那张图两种模式没什么区别,因为去噪和预测CNV在inferCNV里面是分开的两步。

subclusters模式的CNV预测如下图:infercnv.19_HMM_predHMMi6.rand_trees.hmm_mode-subclusters.Pnorm_0.5.repr_intensities.png

单细胞分析实录(13)- inferCNV结合UPhyloplot2分析肿瘤进化,第3张

输出结果有几个文件很重要,后面画进化树会用到:

  • 17_HMM_predHMMi6.rand_trees.hmm_mode-subclusters.cell_groupings
    包含了根据CNV分类的结果,一共两列,一列是类别名称(1.1.1.1, 1.1.1.2, 1.1.2.1, 1.1.2.2, 1.2.1.1, 1.2.1.2, 1.2.2.1, 1.2.2.2这8类),另一列是细胞编号。这个文件不止包含观测,还有参照,参照对应的行要去掉

  • HMM_CNV_predictions.HMMi6.rand_trees.hmm_mode-subclusters.Pnorm_0.5.pred_cnv_regions.dat

# cell_group_name cnv_name        state   chr     start   end
# all_observations.all_observations.1.1.1.1       chr1-region_1   2       chr1    14363   145116922
# all_observations.all_observations.1.1.1.1       chr1-region_3   3       chr1    151264273       156182587

第二列是CNV的name,唯一;第一列是CNV所属的group,示例在"subclusters"模式下有7个group;4 5 6列包含CNV的坐标;第三列表示状态:

# State 1: 0x: complete loss
# State 2: 0.5x: loss of one copy
# State 3: 1x: neutral
# State 4: 1.5x: addition of one copy
# State 5: 2x: addition of two copies
# State 6: 3x: essentially a placeholder for >2x copies but modeled as 3x
  • HMM_CNV_predictions.HMMi6.rand_trees.hmm_mode-subclusters.Pnorm_0.5.pred_cnv_genes.dat
# cell_group_name gene_region_name        state   gene    chr     start   end
# all_observations.all_observations.1.1.1.1       chr1-region_1   2       WASH7P  chr1    14363   29806
# all_observations.all_observations.1.1.1.1       chr1-region_1   2       LINC00115       chr1    14363   29806

每一个group(第一列), 每一个CNV片段(第二列)上面每一个基因(第四列)的CNV状态(第三列),文件中基因这一列是唯一的。相当于上一个文件细化到基因层面。

需要说明的是,上面三个文件只有第一个文件是画进化树需要的,后面两个文件是为了注释进化树的分枝。

3. UPhyloplot2画图

这个小软件其实很简单,两百行代码,只有一个功能就是画图,里面唯一的分析就是计算一下各种CNV cluster的比例,小于5%的cluster不画。
链接在:https://github.com/harbourlab/uphyloplot2

使用的时候,将主程序uphyloplot2.py和文件夹Inputs放在一起,上面提到.cell_groupings文件放到Inputs文件夹里面。

# #命令行运行

# less -S ./try2/17_HMM_predHMMi6.rand_trees.hmm_mode-subclusters.cell_groupings | grep "references" -v | less -S > ./uphyloplot2-2.3/Inputs/test.cell_groupings
# 
# cd uphyloplot2-2.3
# 
# #如果用python3,可能会遇到下面的报错,换成python2就行了
# /d/hsy/software/anaconda/python.exe uphyloplot2.py
# uphyloplot2 version 2.3
# Traceback (most recent call last):
#   File "uphyloplot2.py", line 237, in <module>
#     main()
#   File "uphyloplot2.py", line 166, in main
#     if len(data_row[0].split(".")) > longest_tree:
# IndexError: list index out of range
# 
# #python2
# /d/hsy/software/anaconda/envs/python2/python.exe uphyloplot2.py

结果包含一个图(svg格式),和一张表格,如下图所示:

单细胞分析实录(13)- inferCNV结合UPhyloplot2分析肿瘤进化,第4张

先来看表格,这里面有很多编码,先从四位编码来看,比如1.2.2.1和1.2.2.2这两个前三位是一样的,说明这两个CNV group比较像,而1.2.2就是推测的它俩的上一个状态,对应字母就是K和L的上一个状态是J。表格的第二列是占比,这两个group的占比都是11%,反映在图上,就是K和L的枝长一样(the length of tree branches is proportional to the number of cells in each subclone)。剩下的group以此类推。

四位编码最多是8类,所以树的末端最多8个分枝(本示例7类,1.1.2实在不能往下分),同时占比小于5%的不会输出到表格和图中(这个例子中的1.1.1.2,所以是6个终端分枝,7-1)。对于那些推测的状态,不计算占比,所以都是0,也不反映在枝长上,所以图中推测状态的枝长都是没有意义的。原图需要添加注释,调整分枝角度使其不那么紧凑,才算是合格的图。

延伸一下,这种方法做出的进化树只涉及细胞占比信息,和其他利用突变数量推测molecular time的算法不一样。

4. 进化树分枝注释

这就要用到另外两个.dat文件了

染色体长、短臂层面

第一种注释就是跟之前的文章一样,只注释出长短臂层面的CNV,我的脚本大概可以得到以下的信息

单细胞分析实录(13)- inferCNV结合UPhyloplot2分析肿瘤进化,第5张

再依次从上往下注释进化树的分枝即可,需要对照着上文那个表格(哪一个cluster对应哪一个字母),半成品是这样的

单细胞分析实录(13)- inferCNV结合UPhyloplot2分析肿瘤进化,第6张
基因层面

如果需要注释得更精细,可以把基因也加上,这里我就加到热图上吧。热图显示每个CNV subgroup的CNV情况,进化树显示这些CNV的发生先后,组合起来就是一张CNS级别的配图了

单细胞分析实录(13)- inferCNV结合UPhyloplot2分析肿瘤进化,第7张

本文CNV cluster注释的代码不无偿提供,有需要的朋友可以后台联系小编。

好啦,这一节就到这里,本文也是inferCNV系列的第三篇,后面暂时没有更新inferCNV的想法。

因水平有限,有错误的地方,欢迎批评指正!


https://www.xamrdz.com/bigdata/7ak1877218.html

相关文章: