前言
在之前的一篇推文：Genome Biology：运用机器学习鉴定肿瘤细胞中，Immugent根据Ikarus的原文简单介绍了一下它的功能框架，在本期推文中将通过代码实操的方式来一步步演示如何使用Ikarus来鉴定肿瘤细胞。
不同于之前的Infercnv和CopyKAT软件，主要都是基于R语言，Ikarus是由更擅长机器学习的python所构建的。事实上，在同等运算量的情况下，python处理数据的速度更快，这是因为它可以和其它底层语言衔接的更好。当然，大家也不要有畏惧python的心理，其实使用起来和R基本差不多，不断去适应就好啦。
废话不多说，下面开始展示。。。
代码实操
首先是安装Ikarus软件，要求python>=3.8。
pip install ikarus
#也可以通过下面的代码安装
python -m pip install git+https://github.com/BIMSBbioinfo/ikarus.git
接着是导入所需要的程序，和R中的library()类似。
import urllib.request
import anndata
import pandas as pd
from pathlib import Path
from ikarus import classifier, utils, data
接下来是载入示例数据。。。
Path("data").mkdir(exist_ok=True)
if not Path("data/laughney20_lung/adata.h5ad").is_file():
Path("data/laughney20_lung").mkdir(exist_ok=True)
urllib.request.urlretrieve(
url="https://bimsbstatic.mdc-berlin.de/akalin/Ikarus/part_1/data/laughney20_lung/adata.h5ad",
filename=Path("data/laughney20_lung/adata.h5ad")
)
if not Path("data/lee20_crc/adata.h5ad").is_file():
Path("data/lee20_crc").mkdir(exist_ok=True)
urllib.request.urlretrieve(
url="https://bimsbstatic.mdc-berlin.de/akalin/Ikarus/part_1/data/lee20_crc/adata.h5ad",
filename=Path("data/lee20_crc/adata.h5ad")
)
if not Path("data/tirosh17_headneck/adata.h5ad").is_file():
Path("data/tirosh17_headneck").mkdir(exist_ok=True)
urllib.request.urlretrieve(
url="https://bimsbstatic.mdc-berlin.de/akalin/Ikarus/part_1/data/tirosh17_headneck/adata.h5ad",
filename=Path("data/tirosh17_headneck/adata.h5ad")
)
paths = [
Path("data/laughney20_lung/"),
Path("data/lee20_crc/"),
Path("data/tirosh17_headneck/")
]
names = [
"laughney",
"lee",
"tirosh"
]
adatas = {}
for path, name in zip(paths, names):
adatas[name] = anndata.read_h5ad(path / "adata.h5ad")
# Uncomment to perform preprocessing. Here, the loaded anndata objects are already preprocessed.
# adatas[name] = data.preprocess_adata(adatas[name])
关键一步，训练模型。。。
signatures_path = Path("out/signatures.gmt")
pd.read_csv(signatures_path, sep="\t", header=None)
model = classifier.Ikarus(signatures_gmt=signatures_path, out_dir="out")
train_adata_list = [adatas["laughney"], adatas["lee"]]
train_names_list = ["laughney", "lee"]
obs_columns_list = ["tier_0_hallmark_corrected", "tier_0_hallmark_corrected"]
model.fit(train_adata_list, train_names_list, obs_columns_list, save=True)
model_path = Path("out/core_model.joblib")
model = classifier.Ikarus(signatures_gmt=signatures_path, out_dir="out")
model.load_core_model(model_path)
image.png
基因的名称(见示例中的第一列)在签名中列出，GMT文件对应于它们有意义的单元格类型，并且是以制表符分隔。
截止到现在，我们的模型已经搭建完毕，接下来是对新数据集进行预测。
_ = model.predict(adatas["tirosh"], "tirosh", save=True)
import numpy as np
import scanpy as sc
import matplotlib.pyplot as plt
from sklearn import metrics
def plot_confusion_matrix(
y_true, y_pred, classes, normalize=False, title=None, cmap=plt.cm.Blues, ax=None
):
"""
This function prints and plots the confusion matrix.
Normalization can be applied by setting `normalize=True`.
"""
plt.rcParams["figure.figsize"] = [6, 4]
# print(classes)
if not title:
if normalize:
title = "Normalized confusion matrix"
else:
title = "Confusion matrix, without normalization"
# Compute confusion matrix
cm = metrics.confusion_matrix(y_true, y_pred, labels=classes)
# Only use the labels that appear in the data
# classes = classes[unique_labels(y_true, y_pred)]
if normalize:
cm = cm.astype("float") / cm.sum(axis=1)[:, np.newaxis]
if ax is None:
(fig, ax) = plt.subplots()
im = ax.imshow(cm, interpolation="nearest", cmap=cmap)
ax.figure.colorbar(im, ax=ax)
# We want to show all ticks...
ax.set(
xticks=np.arange(cm.shape[1]),
yticks=np.arange(cm.shape[0]),
# ... and label them with the respective list entries
xticklabels=classes,
yticklabels=classes,
title=title,
ylabel="True label",
xlabel="Predicted label",
)
for item in (
[ax.title, ax.xaxis.label, ax.yaxis.label]
+ ax.get_xticklabels()
+ ax.get_yticklabels()
):
item.set_fontsize(12)
# Rotate the tick labels and set their alignment.
plt.setp(ax.get_xticklabels(), rotation=45, ha="right", rotation_mode="anchor")
# Loop over data dimensions and create text annotations.
fmt = ".2f" if normalize else "d"
thresh = cm.max() / 2.0
for i in range(cm.shape[0]):
for j in range(cm.shape[1]):
ax.text(
j,
i,
format(cm[i, j], fmt),
ha="center",
va="center",
color="white" if cm[i, j] > thresh else "black",
)
return fig, ax
_ = model.get_umap(adatas["tirosh"], "tirosh", save=True)
path = Path("out/tirosh")
results = pd.read_csv(path / "prediction.csv", index_col=0)
adata = anndata.read_h5ad(path / "adata_umap.h5ad")
y = adata.obs.loc[:, "tier_0"]
y_pred_lr = results["final_pred"]
acc = metrics.balanced_accuracy_score(y, y_pred_lr)
print(metrics.classification_report(y, y_pred_lr, labels=["Normal", "Tumor"]))
fig, ax = plot_confusion_matrix(
y,
y_pred_lr,
classes=["Normal", "Tumor"],
title=f"confusion matrix \n train on laughney+lee, test on tirosh \n balanced acc: {acc:.3f}",
)
fig.tight_layout()
image.png
adata.obs.loc[:, "tier_0_pred_correctness"] = "wrongly assigned"
adata.obs.loc[
adata.obs["tier_0"] == adata.obs["final_pred"],
"tier_0_pred_correctness"
] = "correctly assigned"
adata.obs.loc[:, "tier_0_pred_wrong"] = pd.Categorical(
adata.obs["tier_0"].copy(),
categories=np.array(["Normal", "Tumor", "wrongly assigned"]),
ordered=True
)
adata.obs.loc[
adata.obs["tier_0_pred_correctness"] == "wrongly assigned",
"tier_0_pred_wrong"
] = "wrongly assigned"
plt.rcParams["figure.figsize"] = [9, 6]
colors = [
["major"],
["tier_0_pred_wrong"]
]
titles = [
["major types"],
["prediction"]
]
palettes = [
["#7f7f7f", "#98df8a", "#aec7e8", "#9467bd", "#ff0000"],
["#aec7e8", "#ff0000", "#0b559f"],
]
for color, title, palette in zip(colors, titles, palettes):
ax = sc.pl.umap(
adata, ncols=1, size=20,
color=color,
title=title,
wspace=0.25,
vmax="p99",
legend_fontsize=12,
palette=palette,
show=False
)
for item in ([ax.title, ax.xaxis.label, ax.yaxis.label] +
ax.get_xticklabels() + ax.get_yticklabels()):
item.set_fontsize(12)
ax.legend(loc="best")
plt.tight_layout()
plt.show()
image.png
从这个图可以看出Ikarus预测的结果还是相当准确的。
说在最后
Ikarus其实不仅可以预测肿瘤细胞，也可以预测其它特殊的细胞，如各种干细胞。这个主要取决于使用者是通过什么数据来训练出的模型。Immugent有一个大胆的猜想，Ikarus甚至都可以用于对各种免疫细胞的注释。但是有一个问题是，在不同模型或者疾病中，各种免疫细胞的功能状态可能差异很大，鉴于此种考虑，我们不太适合用统一的模型去定义所有同类型免疫细胞。
此外，有一说一，使用python分析数据的感觉确实比R更丝滑。而且基于python出的图，无论是配色还是构图都看起来高达上很多，小伙伴门赶紧实操起来吧。
好啦，本期分享到这里就结束啦，我们下期再会~~
