Ikarus代码实操：运用机器学习鉴定肿瘤细胞

前言

在之前的一篇推文：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出的图，无论是配色还是构图都看起来高达上很多，小伙伴门赶紧实操起来吧。

好啦，本期分享到这里就结束啦，我们下期再会~~