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

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

1. 先来看一下例子

Single-cell analysis reveals new evolutionary complexity in uveal melanoma

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

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

这两篇文章在做完这个分析之后，并没有说太多的东西，只是简单说明了一下哪些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="	",
                                    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

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

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格式)，和一张表格，如下图所示：

先来看表格，这里面有很多编码，先从四位编码来看，比如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，我的脚本大概可以得到以下的信息

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

基因层面

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

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

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

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