好吧小槐， 我們來(lái)繼續(xù)一些轉(zhuǎn)錄組 + VDJ的聯(lián)合分析，今天分享的軟件是scirpy，很多人應(yīng)該用過(guò)须鼎，文章在Scirpy: A Scanpy extension for analyzing single-cell T-cell receptor sequencing data,發(fā)表于Bioinformatics诞外，影響因子7分澜沟，其實(shí)這個(gè)軟件也是做了一些基礎(chǔ)的聯(lián)合分析，并不能分析的很深入峡谊，但還是很值得借鑒茫虽，方法還是很不錯(cuò)的。

Analysis of 3k T cells from cancer既们，我們以官方為例濒析，看看主要的分析內(nèi)容。 原始數(shù)據(jù)集由來(lái)自 14 名未經(jīng)治療的患者的 14 萬(wàn)多個(gè) T 細(xì)胞組成啥纸，這些患者來(lái)自四種不同類型的癌癥号杏。但若干，scripy可以直接讀取10X的分析結(jié)果數(shù)據(jù)斯棒。其他的數(shù)據(jù)格式之前都分享過(guò)盾致，大家可以回看。

圖片.png

# Load the TCR data
adata_tcr = ir.io.read_10x_vdj(
    "example_data/liao-2019-covid19/GSM4385993_C144_filtered_contig_annotations.csv.gz"
)

# Load the associated transcriptomics data
adata = sc.read_10x_h5(
    "example_data/liao-2019-covid19/GSM4339772_C144_filtered_feature_bc_matrix.h5"
)

加載

import numpy as np
import pandas as pd
import scanpy as sc
import scirpy as ir
from matplotlib import pyplot as plt, cm as mpl_cm
from cycler import cycler

sc.set_figure_params(figsize=(4, 4))
sc.settings.verbosity = 2  # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
adata = ir.datasets.wu2020_3k()  ##數(shù)據(jù)

Preprocess Transcriptomics data名船，基礎(chǔ)的scanpy分析流程

sc.pp.filter_genes(adata, min_cells=10)
sc.pp.filter_cells(adata, min_genes=100)
sc.pp.normalize_per_cell(adata, counts_per_cell_after=1000)
sc.pp.log1p(adata)
sc.pp.highly_variable_genes(adata, flavor="cell_ranger", n_top_genes=5000)
sc.tl.pca(adata)
sc.pp.neighbors(adata)

接下來(lái)直接使用作者提供的UMAP坐標(biāo)

adata.obsm["X_umap"] = adata.obsm["X_umap_orig"]
mapping = {
    "nan": "other",
    "3.1-MT": "other",
    "4.1-Trm": "CD4_Trm",
    "4.2-RPL32": "CD4_RPL32",
    "4.3-TCF7": "CD4_TCF7",
    "4.4-FOS": "CD4_FOSS",
    "4.5-IL6ST": "CD4_IL6ST",
    "4.6a-Treg": "CD4_Treg",
    "4.6b-Treg": "CD4_Treg",
    "8.1-Teff": "CD8_Teff",
    "8.2-Tem": "CD8_Tem",
    "8.3a-Trm": "CD8_Trm",
    "8.3b-Trm": "CD8_Trm",
    "8.3c-Trm": "CD8_Trm",
    "8.4-Chrom": "other",
    "8.5-Mitosis": "other",
    "8.6-KLRB1": "other",
}
adata.obs["cluster"] = [mapping[x] for x in adata.obs["cluster_orig"]]

初步的繪圖

sc.pl.umap(
    adata,
    color=["sample", "patient", "cluster", "CD8A", "CD4", "FOXP3"],
    ncols=2,
    wspace=0.7,
)

圖片.png

TCR Quality Control(TCR部分的分析),雖然大多數(shù) T 細(xì)胞受體恰好具有一對(duì) α 和 β 鏈绰上，但多達(dá)三分之一的 T 細(xì)胞可以具有雙重 TCR，即來(lái)自不同等位基因的兩對(duì)受體 .

使用 script.to.chain_qc() 函數(shù)渠驼，我們可以向 adata.obs 添加有關(guān)免疫細(xì)胞受體組成的信息蜈块。 我們可以使用 script.pl.group 豐度（）對(duì)其進(jìn)行可視化。

鏈配對(duì)
Orphan chain是指具有單個(gè)α 或β 受體鏈的細(xì)胞。
額外鏈?zhǔn)侵妇哂型暾?α/β 受體對(duì)和附加鏈的細(xì)胞百揭。
多鏈?zhǔn)侵笝z測(cè)到兩個(gè)以上受體對(duì)的細(xì)胞爽哎。 這些細(xì)胞很可能是雙聯(lián)體。
受體類型和受體亞型
受體類型是指粗略分類為 BCR 和 TCR器一。 具有 BCR 和 TCR 鏈的細(xì)胞被標(biāo)記為不明確的捍掺。
受體亞型是指更具體的分類為 α/β、?/δ位迂、IG-λ 和 IG-κ 鏈構(gòu)型操软。
有關(guān)更多詳細(xì)信息，請(qǐng)參閱 scirpy.tl.chain_qc()请毛。

ir.tl.chain_qc(adata)  ##As expected, the dataset contains only α/β T-cell receptors
ax = ir.pl.group_abundance(adata, groupby="receptor_subtype", target_col="source")

圖片.png

ax = ir.pl.group_abundance(adata, groupby="chain_pairing", target_col="source")

圖片.png

事實(shí)上志鞍，在這個(gè)數(shù)據(jù)集中，約 6% 的細(xì)胞擁有超過(guò)一對(duì)的productive T-cell receptors,接下來(lái)方仿，我們?cè)?UMAP 圖上可視化多鏈細(xì)胞并將它們從下游分析中排除：

sc.pl.umap(adata, color="chain_pairing", groups="multichain")

圖片.png

adata = adata[adata.obs["chain_pairing"] != "multichain", :].copy()
##Similarly, we can use the chain_pairing information to exclude all cells that don’t have at least one full pair of receptor sequences:
adata = adata[~adata.obs["chain_pairing"].isin(["orphan VDJ", "orphan VJ"]), :].copy()
##Finally, we re-create the chain-pairing plot to ensure that the filtering worked as expected:
ax = ir.pl.group_abundance(adata, groupby="chain_pairing", target_col="source")

圖片.png

接下來(lái)重點(diǎn)固棚，Define clonotypes and clonotype clusters

Scirpy 實(shí)現(xiàn)了一種基于網(wǎng)絡(luò)的克隆型定義方法。 創(chuàng)建和可視化克隆型網(wǎng)絡(luò)的步驟類似于使用 scanpy 從轉(zhuǎn)錄組學(xué)數(shù)據(jù)構(gòu)建鄰域圖仙蚜。

圖片.png

計(jì)算CDR3鄰域圖并定義克隆型

scirpy.pp.ir_dist() 根據(jù)序列同一性或相似性計(jì)算 CDR3 核苷酸 (nt) 或氨基酸 (aa) 序列之間的距離此洲。它創(chuàng)建兩個(gè)距離矩陣：一個(gè)用于所有獨(dú)特的 VJ 序列，另一個(gè)用于所有獨(dú)特的 VDJ 序列委粉。距離矩陣被添加到 adata.uns呜师。

函數(shù) scirpy.tl.define_clonotypes() 根據(jù)細(xì)胞的 VJ 和 VDJ CDR3 序列的距離以及函數(shù)參數(shù) dual_ir 和 receptor_arms 的值來(lái)匹配細(xì)胞。最后艳丛，它檢測(cè)圖中的連接模塊并將它們注釋為克隆型匣掸。這將向 adata.obs 添加一個(gè) clone_id 和 clone_id_size 列。

dual_ir 參數(shù)定義了 scipy 如何處理具有一對(duì)以上受體的細(xì)胞氮双。默認(rèn)值是 any 碰酝，這意味著具有任何初級(jí)或次級(jí)受體鏈匹配的細(xì)胞將被認(rèn)為是相同的克隆型。

在這里戴差，根據(jù) nt 序列身份定義克隆型送爸。在后面的步驟中，我們將根據(jù)氨基酸相似性定義克隆型簇暖释。

# using default parameters, `ir_dist` will compute nucleotide sequence identity
ir.pp.ir_dist(adata)
ir.tl.define_clonotypes(adata, receptor_arms="all", dual_ir="primary_only")

為了可視化網(wǎng)絡(luò)袭厂，首先調(diào)用 scirpy.tl.clonotype_network() 來(lái)計(jì)算布局。 然后可以使用 scirpy.pl.clonotype_network() 對(duì)其進(jìn)行可視化球匕。 建議將 min_cells 參數(shù)設(shè)置為 >=2纹磺，以防止單例克隆型混亂網(wǎng)絡(luò)。

ir.tl.clonotype_network(adata, min_cells=2)

結(jié)果圖是一個(gè)網(wǎng)絡(luò)亮曹，其中每個(gè)點(diǎn)代表具有相同受體配置的細(xì)胞橄杨。 當(dāng)將克隆型定義為具有相同 CDR3 序列的細(xì)胞時(shí)秘症，每個(gè)點(diǎn)也是一個(gè)克隆型。 對(duì)于每個(gè)克隆型式矫，數(shù)字克隆型 ID 顯示在圖中乡摹。 每個(gè)點(diǎn)的大小是指具有相同受體配置的細(xì)胞數(shù)量。 分類變量可以可視化為餅圖采转。

ir.pl.clonotype_network(
    adata, color="source", base_size=20, label_fontsize=9, panel_size=(7, 7)
)

圖片.png

重新計(jì)算 CDR3 鄰域圖并定義克隆型cluster

現(xiàn)在可以根據(jù)氨基酸序列相似性重新計(jì)算克隆型網(wǎng)絡(luò)并定義克隆型簇聪廉。

為此，需要更改 set metric="alignment" 并指定一個(gè)截止參數(shù)故慈。 距離基于 BLOSUM62 矩陣板熊。 例如，距離為 10 相當(dāng)于 2 個(gè) R 變異為 N惯悠。 這個(gè)方法最初是由 Dash 等人提出的 TCRdist邻邮。

其CDR3序列之間的距離小于cutoff的所有單元都將在網(wǎng)絡(luò)中連接竣况。

ir.pp.ir_dist(
    adata,
    metric="alignment",
    sequence="aa",
    cutoff=15,
)

ir.tl.define_clonotype_clusters(
    adata, sequence="aa", metric="alignment", receptor_arms="all", dual_ir="any"
)
ir.tl.clonotype_network(adata, min_cells=3, sequence="aa", metric="alignment")

與之前的圖相比克婶，觀察到了幾個(gè)相連的點(diǎn)。 每個(gè)完全連接的子網(wǎng)絡(luò)代表一個(gè)“克隆型簇”丹泉，每個(gè)點(diǎn)仍然代表具有相同受體配置的細(xì)胞情萤。

點(diǎn)由患者著色。 觀察到摹恨，例如筋岛，克隆型 101 和 68（左上和左下）是私有的，即它們僅包含來(lái)自單個(gè)患者的細(xì)胞晒哄。 另一方面睁宰，克隆型 159（左中）是公開(kāi)的，即它在患者 Lung1 和 Lung3 之間共享寝凌。

ir.pl.clonotype_network(
    adata, color="patient", label_fontsize=9, panel_size=(7, 7), base_size=20
)

圖片.png

現(xiàn)在可以通過(guò)子集 AnnData 從特定克隆型簇中提取信息（例如 CDR3 序列）柒傻。 在提取克隆型簇 159 的 CDR3 序列時(shí)，我們檢索了具有不同細(xì)胞數(shù)量的五種不同受體配置较木，對(duì)應(yīng)于圖中的五個(gè)點(diǎn)红符。

adata.obs.loc[adata.obs["cc_aa_alignment"] == "159", :].groupby(
    [
        "IR_VJ_1_junction_aa",
        "IR_VJ_2_junction_aa",
        "IR_VDJ_1_junction_aa",
        "IR_VDJ_2_junction_aa",
        "receptor_subtype",
    ],
    observed=True,
).size().reset_index(name="n_cells")

圖片.png

在克隆型定義中包括 V 基因

使用define_clonotypes() 中的參數(shù)use_v_gene，可以強(qiáng)制克隆型（或克隆型簇）具有相同的V 基因伐债，因此具有相同的CDR1 和2 區(qū)域预侯。 尋找具有不同 V 基因的克隆型簇：

ir.tl.define_clonotype_clusters(
    adata,
    sequence="aa",
    metric="alignment",
    receptor_arms="all",
    dual_ir="any",
    same_v_gene=True,
    key_added="cc_aa_alignment_same_v",
)
# find clonotypes with more than one `clonotype_same_v`
ct_different_v = adata.obs.groupby("cc_aa_alignment").apply(
    lambda x: x["cc_aa_alignment_same_v"].nunique() > 1
)
ct_different_v = ct_different_v[ct_different_v].index.values.tolist()
ct_different_v

在這里，看到當(dāng)設(shè)置 same_v_gene 標(biāo)志時(shí)峰锁，克隆型簇 280 和 765 分別分為 (280, 788) 和 (765, 1071)萎馅。

adata.obs.loc[
    adata.obs["cc_aa_alignment"].isin(ct_different_v),
    [
        "cc_aa_alignment",
        "cc_aa_alignment_same_v",
        "IR_VJ_1_v_call",
        "IR_VDJ_1_v_call",
    ],
].sort_values("cc_aa_alignment").drop_duplicates().reset_index(drop=True)

圖片.png

Clonotype analysis

Clonal expansion

按細(xì)胞類型可視化擴(kuò)展克隆型（即由多個(gè)細(xì)胞組成的克隆型）的數(shù)量。 第一個(gè)選項(xiàng)是將帶有 scirpy.tl.clonal_expansion() 的列添加到 adata.obs 并將其覆蓋在 UMAP 圖上虹蒋。

ir.tl.clonal_expansion(adata)###clonal_expansion refers to expansion categories, i.e singleton clonotypes, clonotypes with 2 cells and more than 2 cells. The clonotype_size refers to the absolute number of cells in a clonotype.
sc.pl.umap(adata, color=["clonal_expansion", "clone_id_size"])

圖片.png

第二個(gè)選項(xiàng)是使用 scirpy.pl.clonal_expansion() 繪圖函數(shù)在堆積條形圖中顯示屬于每個(gè)類別的擴(kuò)展克隆型的細(xì)胞數(shù)量糜芳。

ir.pl.clonal_expansion(adata, groupby="cluster", clip_at=4, normalize=False)

圖片.png

相同的圖拣技，歸一化為集群大小。 克隆擴(kuò)增是對(duì)某個(gè)反應(yīng)性 T 細(xì)胞克隆進(jìn)行陽(yáng)性選擇的標(biāo)志耍目。 因此膏斤，CD8+ 效應(yīng) T 細(xì)胞具有最大比例的擴(kuò)增克隆型是有道理的。

ir.pl.clonal_expansion(adata, "cluster")

圖片.png

預(yù)計(jì) CD8+ 效應(yīng) T 細(xì)胞具有最大比例的擴(kuò)增克隆型邪驮。

與這一觀察結(jié)果一致莫辨，它們具有最低的 scirpy.pl.alpha_diversity() 克隆型。

ax = ir.pl.alpha_diversity(
    adata, metric="normalized_shannon_entropy", groupby="cluster"
)

圖片.png

Clonotype abundance

函數(shù) scirpy.pl.group_abundance() 允許為來(lái)自 obs 的任意類別創(chuàng)建條形圖毅访。 在這里沮榜，我們用它來(lái)顯示十個(gè)最大的克隆型在細(xì)胞類型簇中的分布。

ir.pl.group_abundance(adata, groupby="clone_id", target_col="cluster", max_cols=10)

圖片.png

將計(jì)數(shù)歸一化為每個(gè)樣本的細(xì)胞數(shù)以減輕由于不同樣本大小引起的偏差可能是有益的：

ir.pl.group_abundance(
    adata, groupby="clone_id", target_col="cluster", max_cols=10, normalize="sample"
)

圖片.png

按患者對(duì)條形著色可提供有關(guān)公共和私人克隆型的信息：一些克隆型是私人的喻粹，即特定于某個(gè)組織蟆融，而其他的則是公共的，即它們?cè)诓煌M織之間共享守呜。

ax = ir.pl.group_abundance(
    adata, groupby="clone_id", target_col="source", max_cols=15, figsize=(5, 3)
)

圖片.png

However, clonotypes that are shared between patients are rare:

ax = ir.pl.group_abundance(
    adata, groupby="clone_id", target_col="patient", max_cols=15, figsize=(5, 3)
)

圖片.png

Gene usage

scirpy.tl.group_abundance() 也可以給我們一些關(guān)于 VDJ 使用的信息型酥。 我們可以選擇任何 {TRA,TRB}{1,2}{v,d,j,c}_gene 列來(lái)制作堆積條形圖。 我們使用 max_col 將圖限制為 10 個(gè)最豐富的 V 基因查乒。

ax = ir.pl.group_abundance(
    adata, groupby="IR_VJ_1_v_call", target_col="cluster", normalize=True, max_cols=10
)

圖片.png

ax = ir.pl.group_abundance(
    adata[
        adata.obs["IR_VDJ_1_v_call"].isin(
            ["TRBV20-1", "TRBV7-2", "TRBV28", "TRBV5-1", "TRBV7-9"]
        ),
        :,
    ],
    groupby="cluster",
    target_col="IR_VDJ_1_v_call",
    normalize=True,
)

圖片.png

The exact combinations of VDJ genes can be visualized as a Sankey-plot using scirpy.pl.vdj_usage().

ax = ir.pl.vdj_usage(adata, full_combination=False, max_segments=None, max_ribbons=30)

圖片.png

We can also use this plot to investigate the exact VDJ composition of one (or several) clonotypes:

ir.pl.vdj_usage(
    adata[adata.obs["clone_id"].isin(["68", "101", "127", "161"]), :],
    max_ribbons=None,
    max_segments=100,
)

圖片.png

Spectratype plots

ir.pl.spectratype(adata, color="cluster", viztype="bar", fig_kws={"dpi": 120})

圖片.png

ir.pl.spectratype(
    adata,
    color="cluster",
    viztype="curve",
    curve_layout="shifted",
    fig_kws={"dpi": 120},
    kde_kws={"kde_norm": False},
)

圖片.png

ir.pl.spectratype(
    adata[
        adata.obs["IR_VDJ_1_v_call"].isin(
            ["TRBV20-1", "TRBV7-2", "TRBV28", "TRBV5-1", "TRBV7-9"]
        ),
        :,
    ],
    cdr3_col="IR_VDJ_1_junction_aa",
    color="IR_VDJ_1_v_call",
    normalize="sample",
    fig_kws={"dpi": 120},
)

圖片.png

Comparing repertoires

Repertoire simlarity and overlaps

樣本或樣本組的適應(yīng)性免疫受體庫(kù)的重疊能夠確定重要的克隆型組弥喉，并提供樣本之間相似性的度量。 運(yùn)行 Scirpy 的 repertoire_overlap() 工具創(chuàng)建一個(gè)矩陣玛迄，其中包含每個(gè)樣本中克隆型的豐度由境。 此外，它還計(jì)算樣本的（Jaccard）距離矩陣以及層次聚類的鏈接蓖议。

df, dst, lk = ir.tl.repertoire_overlap(adata, "sample", inplace=False)
df.head()

圖片.png

The distance matrix can be shown as a heatmap, while samples are reordered based on hierarchical clustering.

ir.pl.repertoire_overlap(adata, "sample", heatmap_cats=["patient", "source"])

圖片.png

A specific pair of samples can be compared on a scatterplot, where dot size corresponds to the number of clonotypes at a given coordinate.

ir.pl.repertoire_overlap(
    adata, "sample", pair_to_plot=["LN2", "LT2"], fig_kws={"dpi": 120}
)

圖片.png

Integrating gene expression data(聯(lián)合分析)

Clonotype modularity

使用克隆型模塊性虏杰，我們可以識(shí)別由在轉(zhuǎn)錄上比隨機(jī)預(yù)期更相似的細(xì)胞組成的克隆型。

克隆型模塊性分?jǐn)?shù)表示與隨機(jī)背景模型相比勒虾，細(xì)胞-細(xì)胞鄰域圖中邊數(shù)的 log2 倍變化纺阔。 具有高模塊化評(píng)分的克隆型（或克隆型簇）由具有相似分子表型的細(xì)胞組成。

ir.tl.clonotype_modularity(adata, target_col="cc_aa_alignment")
sc.pl.umap(adata, color="clonotype_modularity")

圖片.png

ir.pl.clonotype_network(
    adata,
    color="clonotype_modularity",
    label_fontsize=9,
    panel_size=(6, 6),
    base_size=20,
)

圖片.png

ir.pl.clonotype_modularity(adata, base_size=20)

圖片.png

Let’s further inspect the two top scoring candidates. We can extract that information from adata.obs["clonotype_modularity"].

clonotypes_top_modularity = list(
    adata.obs.set_index("cc_aa_alignment")["clonotype_modularity"]
    .sort_values(ascending=False)
    .index.unique().values[:2]
)

sc.pl.umap(
    adata,
    color="cc_aa_alignment",
    groups=clonotypes_top_modularity,
    palette=cycler(color=mpl_cm.Dark2_r.colors),
)

圖片.png

觀察到它們（大部分）僅限于單個(gè)集群从撼。 通過(guò)利用 scanpy 的差異表達(dá)模塊州弟，我們可以將兩種克隆型中細(xì)胞的基因表達(dá)與其余細(xì)胞進(jìn)行比較。

sc.tl.rank_genes_groups(
    adata, "clone_id", groups=clonotypes_top_modularity, reference="rest", method="wilcoxon"
)
fig, axs = plt.subplots(1, 2, figsize=(8, 4))
for ct, ax in zip(clonotypes_top_modularity, axs):
    sc.pl.rank_genes_groups_violin(adata, groups=[ct], n_genes=15, ax=ax, show=False, strip=False)

圖片.png

Clonotype imbalance among cell clusters

例如低零，使用從基因表達(dá)簇推斷的細(xì)胞類型注釋婆翔，可以比較屬于 CD8+ 效應(yīng) T 細(xì)胞和 CD8+ 組織駐留記憶 T 細(xì)胞的克隆型。

freq, stat = ir.tl.clonotype_imbalance(
    adata,
    replicate_col="sample",
    groupby="cluster",
    case_label="CD8_Teff",
    control_label="CD8_Trm",
    inplace=False,
)
top_differential_clonotypes = stat["clone_id"].tolist()[:3]

Showing top clonotypes on a UMAP clearly shows that clonotype 101 is featured by CD8+ tissue-resident memory T cells, while clonotype 68 by CD8+ effector and effector memory cells.

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 4), gridspec_kw={"wspace": 0.6})
sc.pl.umap(adata, color="cluster", ax=ax1, show=False)
sc.pl.umap(
    adata,
    color="clone_id",
    groups=top_differential_clonotypes,
    ax=ax2,
    # increase size of highlighted dots
    size=[
        80 if c in top_differential_clonotypes else 30 for c in adata.obs["clone_id"]
    ],
    palette=cycler(color=mpl_cm.Dark2_r.colors),
)

圖片.png

Repertoire overlap of cell types

就像比較樣本之間的曲目重疊一樣掏婶，Scirpy 還提供基因表達(dá)簇或細(xì)胞亞群之間的比較啃奴。 例如，顯示了上面比較的兩種細(xì)胞類型的曲目重疊雄妥。

ir.tl.repertoire_overlap(adata, "cluster")
ir.pl.repertoire_overlap(
    adata, "cluster", pair_to_plot=["CD8_Teff", "CD8_Trm"], fig_kws={"dpi": 120}
)

圖片.png

Marker genes in top clonotypes

sc.tl.rank_genes_groups(
    adata, "clone_id", groups=["101"], reference="68", method="wilcoxon"
)
sc.pl.rank_genes_groups_violin(adata, groups="101", n_genes=15)

圖片.png

其實(shí)還是集中在分析克隆的基因選擇和多樣性最蕾，不如CoNGA依溯。

生活很好，有你更好