作者，追風(fēng)少年i
hello，大家好，新的一周，新的開始，腦海里浮現(xiàn)周星馳的一句臺(tái)詞，只有越來越強(qiáng)大，才能越來越童話。
今天我們簡(jiǎn)單來學(xué)習(xí)一個(gè)簡(jiǎn)單的內(nèi)容，腫瘤進(jìn)化樹，先看示例。不過之前大家要先補(bǔ)充一下基礎(chǔ)知識(shí)，我列在下面，供大家參考。
10X單細(xì)胞（10X空間轉(zhuǎn)錄組）CNV分析之inferCNVpy
10X單細(xì)胞個(gè)性化分析之CNV篇
10X單細(xì)胞（10X空間轉(zhuǎn)錄組）CNV分析回顧之CopyKAT
10X單細(xì)胞（10X空間轉(zhuǎn)錄組）分析之將單細(xì)胞轉(zhuǎn)錄組映射到拷貝數(shù)進(jìn)化樹（scatrex）
10X空間轉(zhuǎn)錄組CNV和速率分析的空間圖譜
10X單細(xì)胞（10X空間轉(zhuǎn)錄組）數(shù)據(jù)分析之識(shí)別腫瘤細(xì)胞的CNV分析原理
其中文章10X單細(xì)胞個(gè)性化分析之CNV篇和10X單細(xì)胞（10X空間轉(zhuǎn)錄組）數(shù)據(jù)分析之識(shí)別腫瘤細(xì)胞的CNV分析原理主要講分析原理，建議大家好好看看。
先來看看范例
文章Single-cell analysis reveals new evolutionary complexity in uveal melanoma，2020年的NC文章。
圖片要說明的問題并不復(fù)雜，a 和 b圖就是簡(jiǎn)單的CNV事件熱圖，c圖是每一個(gè)樣本的進(jìn)化樹，上方主干對(duì)應(yīng)的CNV事件，是早期就發(fā)生的；下方分枝的CNV事件是后期發(fā)生的。其中蘊(yùn)含的邏輯就是含有某個(gè)CNV的細(xì)胞占比越多，那么這個(gè)CNV就發(fā)生得越早；含有某個(gè)CNV的細(xì)胞占比越少，這個(gè)CNV就發(fā)生得越晚。
文章Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma
和上一張圖的形式差不多，解讀也是一樣的，我們今天的目標(biāo)就是來實(shí)現(xiàn)進(jìn)化樹的結(jié)果
其中我們需要注意2點(diǎn)
- 進(jìn)化樹的繪制
- 長短臂注釋的由來
我們先來看看文章中的方法，在文章Single-cell analysis reveals new evolutionary complexity in uveal melanoma中的方法是這樣的：
inferCNV and clonality analysis:
“使用 10x 數(shù)據(jù)”部分（inferCNV，https://github.com/broadinstitute/inferCNV）中的建議，從 Seurat 對(duì)象中提取原始基因表達(dá)數(shù)據(jù)。對(duì)于每位患者，通過 CD3E 的表達(dá)高于平均表達(dá) 2 個(gè)標(biāo)準(zhǔn)差并且沒有 PRAME 和 HTR2B 的表達(dá)來選擇正常參考細(xì)胞(選擇reference)。對(duì)于 inferCNV 分析，使用了以下參數(shù)：“denoise”、默認(rèn)隱藏馬爾可夫模型 (HMM) 設(shè)置和“cutoff”值 0.1。為了減少假陽性 CNV 調(diào)用的可能性，實(shí)施了默認(rèn)的貝葉斯?jié)撛诨旌夏Ｐ蛠碜R(shí)別每個(gè)細(xì)胞中改變的后驗(yàn)概率。使用閾值的默認(rèn)值“0.5”過濾低概率 CNV。為了確定每個(gè)腫瘤中的克隆 CNV 變化，對(duì) HMM 生成的 CNV 使用了“subcluster”方法。 GRCh38 cytoband information用于將每個(gè) CNV 轉(zhuǎn)換為 p 或 q 臂水平變化，以便根據(jù)其位置進(jìn)行簡(jiǎn)化。每個(gè) CNV 都被注釋為gain或loss。數(shù)據(jù)轉(zhuǎn)換后，包含相同臂級(jí) CNV 的亞克隆被折疊，樹被重組以準(zhǔn)確表示亞克隆 CNV 架構(gòu)。該分析排除了線粒體CNV。對(duì)于數(shù)據(jù)可視化，開發(fā)了 UPhyloplot2 (https://github.com/harbourlab/UPhyloplot2) 繪圖算法以自動(dòng)生成腫瘤內(nèi)進(jìn)化樹。從 inferCNV HMM 子集群 CNV 預(yù)測(cè)算法策劃的臂級(jí) CNV 調(diào)用和每個(gè)子克隆中的細(xì)胞百分比用作輸入。為每個(gè)樣本生成可視化系統(tǒng)發(fā)育樹的可縮放矢量圖形 (.svg) 文件。臂長與細(xì)胞百分比加上間隔（圓直徑 + 5 像素）成正比。
在文章Single-cell RNA landscape of intratumoral heterogeneity and immunosuppressive microenvironment in advanced osteosarcoma中：
Single-cell copy-number variation (CNV) and clonality analysis:
使用 R 的 inferCNV 包（版本 1.2.2；https://github.com/broadinstitute/inferCNV/wiki）估計(jì)成骨細(xì)胞和軟骨母細(xì)胞腫瘤細(xì)胞中每個(gè)細(xì)胞的初始 CNV。計(jì)算成骨細(xì)胞和軟骨細(xì)胞腫瘤細(xì)胞的CNVs，并以免疫細(xì)胞作為參考。在過濾具有 < 2000 個(gè) UMI 的不合格細(xì)胞后，使用包括“denoise”、默認(rèn)隱藏馬爾可夫模型 (HMM) 設(shè)置和“cutoff”值 0.1 在內(nèi)的參數(shù)執(zhí)行 inferCNV 分析。為了減少假陽性 CNV 調(diào)用，實(shí)現(xiàn)了默認(rèn)的貝葉斯?jié)撛诨旌夏Ｐ蛠碜R(shí)別每個(gè)細(xì)胞中 CNV 變化的后驗(yàn)概率，默認(rèn)值為 0.5 作為閾值。為了推斷克隆單細(xì)胞CNV的變化，應(yīng)用“subcluster”方法根據(jù)HMM生成的CNV值推斷子簇細(xì)胞。注釋gene cytoband information信息，每個(gè) p 或 q 臂水平變化都根據(jù)其位置簡(jiǎn)單地轉(zhuǎn)換為等效的 CNV。每個(gè) CNV 都被注釋為gain或loss。數(shù)據(jù)轉(zhuǎn)換后，包含相同臂級(jí) CNV 的亞克隆被折疊，樹被重組以表示亞克隆 CNV 架構(gòu)。對(duì)于數(shù)據(jù)可視化，遵循了 Durante 等人開發(fā)的 UPhyloplot2 算法。 （https://github.com/harbourlab/UPhyloplot2）自動(dòng)生成腫瘤內(nèi)進(jìn)化樹。從 inferCNV HMM 子集群 CNV 預(yù)測(cè)算法策劃的臂級(jí) CNV 調(diào)用和每個(gè)子克隆中的細(xì)胞百分比用作輸入。
我們希望拿到的結(jié)果
信息量有點(diǎn)龐大，我們先來一步一步的解析：
首先來看基礎(chǔ)的CNV分析參數(shù)：
- “denoise”
- (HMM) 設(shè)置和“cutoff”值 0.1
- 貝葉斯?jié)撛诨旌夏Ｐ蛠碜R(shí)別每個(gè)細(xì)胞中 CNV 變化的后驗(yàn)概率，默認(rèn)值為 0.5 作為閾值
- 注意其中參考細(xì)胞的選擇，可以選擇免疫細(xì)胞
運(yùn)行實(shí)例
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)
注意兩個(gè)參數(shù)cluster_by_groups=F，以及analysis_mode="subclusters"，這個(gè)參數(shù)最終會(huì)將腫瘤細(xì)胞分為8個(gè)cluster(少數(shù)情況是7類，如果實(shí)在找不出進(jìn)一步的差別)，每個(gè)cluster有各自的CNV模式，如果analysis_mode="samples"，則一個(gè)樣本不同細(xì)胞最終預(yù)測(cè)的CNV模式是唯一的。另外需要注意的是，一般文章放的熱圖是去噪后的熱圖，那張圖兩種模式?jīng)]什么區(qū)別，因?yàn)槿ピ牒皖A(yù)測(cè)CNV在inferCNV里面是分開的兩步。
inferCNV分析完之后，我們一般會(huì)得到如下的三個(gè)文件：
17_HMM_predHMMi6.rand_trees.hmm_mode-subclusters.cell_groupings包含了根據(jù)CNV分類的結(jié)果，一共兩列，一列是類別名稱(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類)，另一列是細(xì)胞編號(hào)。這個(gè)文件不止包含觀測(cè)，還有參照，參照對(duì)應(yīng)的行要去掉。
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個(gè)group；4 5 6列包含CNV的坐標(biāo)；第三列表示狀態(tài)：
# 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
每一個(gè)group（第一列）, 每一個(gè)CNV片段（第二列）上面每一個(gè)基因（第四列）的CNV狀態(tài)（第三列），文件中基因這一列是唯一的。相當(dāng)于上一個(gè)文件細(xì)化到基因?qū)用?/strong>。
需要說明的是，上面三個(gè)文件只有第一個(gè)文件是畫進(jìn)化樹需要的，后面兩個(gè)文件是為了注釋進(jìn)化樹的分枝。
