这一个分析，我目前只在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="\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
输出结果有几个文件很重要，后面画进化树会用到：
-
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的想法。
因水平有限，有错误的地方，欢迎批评指正！
网友评论