單細胞分析實錄(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="\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的想法。

因水平有限，有錯誤的地方，歡迎批評指正！