隔離的第七日秃殉，孤獨(dú)如影隨形坝初，必須寫一些內(nèi)容來(lái)排解心中的無(wú)聊，上一篇我們?cè)敿?xì)回顧了copyCAT鳄袍，文章在10X單細(xì)胞（10X空間轉(zhuǎn)錄組）CNV分析回顧之CopyKAT绢要，還有分享的文章copyKAT推斷單細(xì)胞轉(zhuǎn)錄組腫瘤細(xì)胞CNV（自動(dòng)識(shí)別腫瘤normal和tumor）拗小，但是copyCAT的引用率還是比較低重罪，相比于inferCNV 2000多的引用率十籍，還是差很多蛆封，所以，inferCNV還是主流惨篱，我們這一篇就來(lái)回顧一下。
其實(shí)inferCNV已經(jīng)是非常老的一個(gè)軟件了围俘，14年就發(fā)表了文章，所以我們來(lái)分享其最新的python版本界牡。
加載,與scanpy無(wú)縫銜接
import scanpy as sc
import infercnvpy as cnv
import matplotlib.pyplot as plt
sc.settings.set_figure_params(figsize=(5, 5))
Loading the example dataset
前處理
- 應(yīng)該已經(jīng)過(guò)濾掉低質(zhì)量的細(xì)胞簿寂，并且必須對(duì)輸入數(shù)據(jù)進(jìn)行歸一化和對(duì)數(shù)轉(zhuǎn)換
- 此外宿亡，基因組位置需要存儲(chǔ)在 adata.var 中常遂。 "chromosome", "start", "end"分別保存每個(gè)基因的染色體以及在該染色體上的開始和結(jié)束位置挽荠。
- Infercnvpy provides the
infercnvpy.io.genomic_position_from_gtf()function to read these information from a GTF file and add them to.
adata = cnv.datasets.maynard2020_3k()
adata.var.loc[:, ["ensg", "chromosome", "start", "end"]].head()
sc.pl.umap(adata, color="cell_type")
Running infercnv
現(xiàn)在運(yùn)行 infercnvpy.tl.infercnv()克胳。 本質(zhì)上圈匆，該方法通過(guò)染色體和基因組位置對(duì)基因進(jìn)行分類漠另，并將基因組區(qū)域的平均基因表達(dá)與參考進(jìn)行比較。 原始的 inferCNV 方法使用 100 的窗口大小笆搓，但更大的窗口大小可能有意義，具體取決于數(shù)據(jù)集中的基因數(shù)量纬傲。
infercnv()adds a cell x genomic_region matrix to adata.obsm[“X_cnv”].
Choosing reference cells(最為關(guān)鍵的一步)
- 最常見的用例是將腫瘤與正常細(xì)胞進(jìn)行比較。 如果有關(guān)于哪些細(xì)胞正常的先驗(yàn)信息（例如叹括，來(lái)自基于轉(zhuǎn)錄組學(xué)數(shù)據(jù)的細(xì)胞類型注釋）葫录，建議將此信息提供給 infercnv()领猾。
- 可以提供的不同細(xì)胞類型越多越好米同。 一些細(xì)胞類型在生理上過(guò)度表達(dá)某些基因組區(qū)域（例如漿細(xì)胞高度表達(dá)基因組相鄰的免疫球蛋白基因）摔竿。 如果提供多種細(xì)胞類型面粮，則僅考慮與所有提供的細(xì)胞類型不同的區(qū)域受 CNV 影響。
- 如果不提供任何參考熬苍，則使用所有細(xì)胞的平均值，這可能適用于包含足夠腫瘤和正常細(xì)胞的數(shù)據(jù)集袁翁。
# We provide all immune cell types as "normal cells".
cnv.tl.infercnv(
adata,
reference_key="cell_type",
reference_cat=[
"B cell",
"Macrophage",
"Mast cell",
"Monocyte",
"NK cell",
"Plasma cell",
"T cell CD4",
"T cell CD8",
"T cell regulatory",
"mDC",
"pDC",
],
window_size=250,
)
現(xiàn)在柴底，可以按細(xì)胞類型和染色體繪制平滑的基因表達(dá)粱胜。 可以觀察到主要由腫瘤細(xì)胞組成的上皮細(xì)胞cluster似乎受到拷貝數(shù)變化的影響柄驻。
cnv.pl.chromosome_heatmap(adata, groupby="cell_type")
Clustering by CNV profiles and identifying tumor cells
為了對(duì)細(xì)胞進(jìn)行聚類和注釋焙压，infercnvpy 結(jié)合了 scanpy 工作流程鸿脓。 以下函數(shù)與它們的 scanpy 對(duì)應(yīng)函數(shù)完全一樣涯曲，除了使用 CNV 分析矩陣作為輸入野哭。 使用這些函數(shù)，可以執(zhí)行基于圖形的聚類并根據(jù) CNV 分析矩陣生成 UMAP 圖拨黔。 基于這些聚類，可以注釋腫瘤和正常細(xì)胞绰沥。
cnv.tl.pca(adata)
cnv.pp.neighbors(adata)
cnv.tl.leiden(adata)
運(yùn)行 leiden 聚類后，可以通過(guò) CNV 聚類繪制染色體熱圖揪利。 可以觀察到态兴，與底部的clusters相反疟位，頂部的clusters基本上沒有差異表達(dá)的基因組區(qū)域瞻润。 差異表達(dá)的區(qū)域可能是由于拷貝數(shù)變異甜刻，并且各自的clusters可能代表腫瘤細(xì)胞绍撞。
cnv.pl.chromosome_heatmap(adata, groupby="cnv_leiden", dendrogram=True)
UMAP plot of CNV profiles
我們可以將相同的clusters可視化為 UMAP 圖。 此外傻铣，infercnvpy.tl.cnv_score() 計(jì)算一個(gè)匯總分?jǐn)?shù)，量化每個(gè)cluster的拷貝數(shù)變異量祥绞。 它被簡(jiǎn)單地定義為每個(gè)cluster的 CNV 矩陣的絕對(duì)值的平均值鸭限。
cnv.tl.umap(adata)
cnv.tl.cnv_score(adata)
UMAP 圖由一大團(tuán)“正常”細(xì)胞和幾個(gè)具有不同 CNV 分布的較小cluster組成两踏。 除了由纖毛細(xì)胞組成的cluster“12”外，孤立的cluster都是上皮細(xì)胞梦染。 這些可能是腫瘤細(xì)胞赡麦，每個(gè)cluster代表一個(gè)單獨(dú)的亞克隆。
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(11, 11))
ax4.axis("off")
cnv.pl.umap(
adata,
color="cnv_leiden",
legend_loc="on data",
legend_fontoutline=2,
ax=ax1,
show=False,
)
cnv.pl.umap(adata, color="cnv_score", ax=ax2, show=False)
cnv.pl.umap(adata, color="cell_type", ax=ax3)
還可以在基于轉(zhuǎn)錄組學(xué)的 UMAP 圖上可視化 CNV 分?jǐn)?shù)和cluster泛粹。 同樣，可以看到存在屬于不同 CNV cluster的上皮細(xì)胞subcluster肮疗，并且這些cluster往往具有最高的 CNV 分?jǐn)?shù)。
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(
2, 2, figsize=(12, 11), gridspec_kw=dict(wspace=0.5)
)
ax4.axis("off")
sc.pl.umap(adata, color="cnv_leiden", ax=ax1, show=False)
sc.pl.umap(adata, color="cnv_score", ax=ax2, show=False)
sc.pl.umap(adata, color="cell_type", ax=ax3)
Classifying tumor cells
基于這些觀察族吻，現(xiàn)在可以將細(xì)胞分配給“腫瘤”或“正常”超歌。 為此砍艾，我們?cè)?adata.obs 中添加一個(gè)新列 cnv_status巍举。
adata.obs["cnv_status"] = "normal"
adata.obs.loc[
adata.obs["cnv_leiden"].isin(["10", "13", "15", "5", "11", "16", "12"]), "cnv_status"
] = "tumor"
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5), gridspec_kw=dict(wspace=0.5))
cnv.pl.umap(adata, color="cnv_status", ax=ax1, show=False)
sc.pl.umap(adata, color="cnv_status", ax=ax2)
cnv.pl.chromosome_heatmap(adata[adata.obs["cnv_status"] == "tumor", :])
The inferCNV method
本質(zhì)上脆荷，這個(gè)包是 infercnv 的 Python 重新實(shí)現(xiàn)懊悯。通過(guò)使用 numpy蜓谋、scipy 和稀疏矩陣，它的計(jì)算效率要高得多炭分。
Computation steps
1、從所有細(xì)胞中減去參考基因表達(dá)捧毛。 由于數(shù)據(jù)在對(duì)數(shù)空間中观堂，這有效地計(jì)算了對(duì)數(shù)倍數(shù)變化呀忧。 如果有多個(gè)類別的引用可用（即為 reference_cat 指定了多個(gè)值）师痕，則log fold change是“bounded”：
- 分別計(jì)算每個(gè)類別的平均基因表達(dá)而账。
- 在所有參考平均值的最小值和最大值范圍內(nèi)的值會(huì)收到 0 的對(duì)數(shù)倍數(shù)變化胰坟，因?yàn)樗鼈儾槐灰暈榕c背景不同。
- 從小于所有參考平均值的最小值的值中減去該最小值笔横。
- 從大于所有參考平均值的最大值的值中減去該最大值竞滓。
此過(guò)程避免了由于聚集基因區(qū)域（例如不同免疫細(xì)胞類型中的免疫球蛋白或 HLA 基因）的細(xì)胞類型特異性表達(dá)而調(diào)用假陽(yáng)性 CNV 區(qū)域。
2吹缔、Clip the fold changes at -lfc_cap and +lfc_cap.
3、通過(guò)基因組位置平滑基因表達(dá)涛菠。 計(jì)算長(zhǎng)度為 window_size 的運(yùn)行窗口的平均值。 僅計(jì)算每第 n 個(gè)窗口以節(jié)省時(shí)間和空間撇吞，其中 n = step。
4牍颈、通過(guò)從每個(gè)細(xì)胞中減去每個(gè)細(xì)胞的中位數(shù)迄薄，按細(xì)胞將平滑的基因表達(dá)居中煮岁。
5讥蔽、執(zhí)行噪聲過(guò)濾。 值 < dynamic_theshold * STDDEV 設(shè)置為 0冶伞，其中 STDDEV 是平滑基因表達(dá)的標(biāo)準(zhǔn)偏差
6、Smooth the final result using a median filter.
Preparing input data
- annadata.AnnData 對(duì)象應(yīng)該已經(jīng)被過(guò)濾為低質(zhì)量的細(xì)胞步氏。 adata.X 需要進(jìn)行規(guī)范化和對(duì)數(shù)轉(zhuǎn)換。 該方法應(yīng)該對(duì)不同的歸一化方法（scanpy.pp.normalize_total()荚醒、scran 等）相當(dāng)穩(wěn)健芋类。
- 該方法需要一個(gè)“參考”值界阁，與基因組區(qū)域的表達(dá)進(jìn)行比較侯繁。 如果數(shù)據(jù)集包含不同的細(xì)胞類型并且包括腫瘤細(xì)胞和正常細(xì)胞，則可以使用所有細(xì)胞的平均值作為參考泡躯。 這是默認(rèn)設(shè)置。
- 如果已經(jīng)知道哪些細(xì)胞是“正常的”精续，可以提供從 adata.obs 到 reference_key 的列坝锰，其中包含注釋重付。 reference_cat 在reference_key 中指定一個(gè)或多個(gè)引用正常細(xì)胞的值顷级。
鏈接在inferCNVpy
生活很好，有你更好
