CellphoneDB 方法 3（差异表达）

在这个示例中，我们使用方法 3 (degs_analysis_method) 来研究在滋养层细胞 (Trophoblast) 分化和侵入母体子宫过程中，血管周细胞 (PVs) 和滋养层细胞之间的细胞-细胞相互作用如何变化。此方法将检索至少有一个参与相互作用的基因（参与相互作用的基因）差异表达的相互作用。

差异表达的基因应该由用户预先计算，使用他们首选的方法（例如，Limma、DEseq2、Seurat、TradeSeq 等），并将显著结果以预定义的表格格式传递给 CellphoneDB。

此笔记本假设您已经知道如何下载 CellphoneDB 数据库或创建自己的数据库。如果不是这种情况，请查看 T0_BuildDBfromFiles.ipynb 或 T0_BuildDBfromRelease.ipynb。在此笔记本中，我们将解释如何运行 CellphoneDB 的差异表达方法。

检查 Python 版本

import pandas as pd
import sys
import os

pd.set_option('display.max_columns', 100)
# 定义分析的基本目录
os.chdir('/home/jovyan/cpdb_tutorial')
# 检查环境是否包含 CellphoneDB 要求的 Python >= 3.8
print(sys.version)

安装 CellphoneDB

在当前 conda 环境中安装 CellphoneDB 的最新版本。
如果您想查看详细的安装过程，请移除 --quiet 标志。

pip install --quiet cellphonedb

---

输入文件

差异表达方法接受 6 个输入文件（其中 4 个是必需的）。

• cpdb_file_path：（必需）数据库的路径。
• meta_file_path：（必需）链接细胞条形码到集群标签的元文件的路径。
• counts_file_path：（必需）标准化计数文件的路径（不是 z 转换的），可以是文本格式或 h5ad（推荐）。
• degs_file_path：（必需）指示每个集群中差异表达基因的 DEG 文件的路径。仅应包括显著差异表达的基因。
• microenvs_file_path：（可选）将细胞类型/集群按微环境分组的微环境文件的路径。提供微环境文件时，CellphoneDB 将限制相互作用在微环境内的细胞之间。
• active_tf_path：（可选）活跃转录因子的路径。

degs_file_path 和 microenvs_file_path 的内容将根据研究者想要回答的生物学问题而定。

在这个例子中，我们研究了随着滋养层细胞分化和侵入母体子宫过程中，血管周细胞 (PVs) 和滋养层细胞之间的细胞-细胞相互作用如何变化。因此，degs_file_path 只包含沿着滋养层分化谱系差异表达的基因，因为我们对那些在此分化过程中变化的细胞-细胞通信过程感兴趣。microenvs_file_path 只包含在胎盘特定解剖区域（即在同一时空邻域中的 PVs 和滋养层细胞状态）存在的细胞。meta_file_path 和 counts_file_path 包含了我们感兴趣的所有 PVs 和滋养层细胞。

CellphoneDB 将检索在 PVs 和滋养层细胞之间发生的所有相互作用，其中：(i) 所有蛋白质在相应的细胞类型中表达，并且 (ii) 至少有一个基因由滋养层细胞差异表达。

cpdb_file_path = 'db/v5/cellphonedb.zip'
meta_file_path = 'data/metadata.tsv'
counts_file_path = 'data/normalised_log_counts.h5ad'
microenvs_file_path = 'data/microenvironment.tsv'
degs_file_path = 'data/DEGs_inv_trophoblast.tsv'
active_tf_path = 'data/active_TFs.tsv'
out_path = 'results/method3_withScore'

检查输入文件

1) 元数据文件由两列组成：

• barcode_sample：此列指示实验中每个细胞的条形码。
• cell_type：此列表示分配的细胞标签。

metadata = pd.read_csv(meta_file_path, sep='\t')
metadata.head(3)

2) 计数文件是来自 scanpy 的 h5ad 对象。此对象的维度和顺序必须与元数据文件的维度一致，即两个文件中的细胞数量必须相同。

import anndata

adata = anndata.read_h5ad(counts_file_path)
adata.shape

# 检查元数据和计数文件中的条形码是否一致
list(adata.obs.index).sort() == list(metadata['barcode_sample']).sort()

3) 差异表达基因文件是一个两列文件，指示在某个细胞类型中上调（或特定）的基因。第一列对应于集群名称（这些名称与元数据文件中的名称匹配），第二列是上调基因。其余列被 CellphoneDB 忽略。此文件中存在的所有基因都会被考虑，因此用户必须在此文件中仅提供那些被认为上调或对分析重要的基因。

pd.read_csv(degs_file_path, sep='\t').head(3)

4) 微环境定义了属于特定微环境的细胞类型。CellphoneDB 仅计算属于给定微环境的细胞之间的相互作用。在此文件中，我们仅定义了一个微环境。

microenv = pd.read_csv(microenvs_file_path, sep='\t')
microenv.head(3)

# 显示按微环境分组的细胞
microenv.groupby('microenvironment')['cell_type'].apply(lambda x: list(x.value_counts().index))

5) 活跃转录因子定义了在特定细胞类型中活跃的转录因子。

pd.read_csv(active_tf_path, sep='\t').head(3)

---

使用差异分析方法运行 CellphoneDB（方法 3）

此方法的输出将保存到 output_path，并分配到预定义的变量中。

from cellphonedb.src.core.methods import cpdb_degs_analysis_method

cpdb_results = cpdb_degs_analysis_method.call(
    cpdb_file_path=cpdb_file_path,                            # 必需：CellphoneDB 数据库 zip 文件。
    meta_file_path=meta_file_path,                            # 必需：定义条形码到细胞标签的 tsv 文件。
    counts_file_path=counts_file_path,                        # 必需：标准化计数矩阵 - 计数文件的路径或内存中的 AnnData 对象。
    degs_file_path=degs_file_path,                            # 必需：包含 DEG 的 tsv 文件。
    counts_data='hgnc_symbol',                                # 定义计数矩阵中的基因注释。
    active_tfs_file_path=active_tf_path,                      # 可选：定义细胞类型及其活跃的 TFs。
    microenvs_file_path=microenvs_file_path,                  # 可选（默认：无）：定义每个微环境中的细胞。
    score_interactions=True,                                  # 可选：是否对相互作用进行评分。
    threshold=0.1,                                            # 定义用于分析的基因在细胞中表达的最小百分比。
    result_precision=3,                                       # 设置 significant_means 中均值值的舍入。
    separator='|',                                            # 设置用于在结果数据框中分隔细胞的字符串 "cellA|CellB"。
    debug=False,                                              # 以 pkl 格式保存分析过程中使用的所有中间表。
    output_path=out_path,                                     # 保存结果的路径
    output_suffix=None,                                       # 替换输出文件中的时间戳为用户定义的字符串（默认：无）
    threads=25
)

结果保存为文件，并作为字典保存在 cpdb_results 变量中。

list(cpdb_results.keys())

---

输出文件描述

大多数输出文件共享常见列：

• id_cp_interaction：每个存储在数据库中的相互作用的唯一 CellphoneDB 标识符。
• interacting_pair：用“|”分隔的相互作用对的名称。
• partner A 或 B：第一个相互作用伙伴 (A) 或第二个 (B) 的标识符。可以是：UniProt（前缀 simple:）或复合体（前缀 complex:）。
• gene A 或 B：第一个相互作用伙伴 (A) 或第二个 (B) 的基因标识符。标识符将取决于用户输入的列表。
• secreted：如果其中一个伙伴是分泌的则为真。
• Receptor A 或 B：如果第一个相互作用伙伴 (A) 或第二个 (B) 在我们的数据库中被注释为受体则为真。
• annotation_strategy：如果相互作用由 CellphoneDB 开发人员注释则为 Curated。否则，为相互作用下载的数据库名称。
• is_integrin：如果其中一个伙伴是整合素则为真。
• directionality：指示相互作用的方向性和相互作用伙伴的特征。
• classification：相互作用伙伴的通路分类。

相关相互作用字段：

• cell_a|cell_b：如果检测到相互作用为显著，则为 1，否则为 0。

cpdb_results['relevant_interactions'].head(2)

显著均值字段：

• significant_mean：所有相互作用伙伴的显著均值计算。如果相互作用被认为是显著的，则该值为均值。否则，值设置为 0。

cpdb_results['significant_means'].head(2)

均值字段：

• means：所有相互作用伙伴的均值：均值指的是相应相互作用细胞类型对中各个伙伴平均表达值的总均值。如果其中一个均值为 0，则总均值设置为 0。

cpdb_results['means'].head(2)

去卷积字段：

• gene_name：参与 “means.csv” 文件中定义的相互作用的一个亚基的基因标识符。标识符将取决于用户输入的列表。
• uniprot：参与 “means.csv” 文件中定义的相互作用的一个亚基的 UniProt 标识符。
• is_complex：如果亚基是复合体的一部分，则为真。如果不是，则为单一，如果是复合体，则为复合体。
• protein_name：参与 “means.csv” 文件中定义的相互作用的一个亚基的蛋白质名称。
• complex_name：如果亚基是复合体的一部分，则为复合体名称。如果不是，则为空。
• id_cp_interaction：每个存储在数据库中的相互作用的唯一 CellphoneDB 标识符。
• mean：每个集群中相应基因的平均表达。

cpdb_results['deconvoluted'].head(2)

相互作用评分字段：

• scores：范围从 0 到 100 的评分。分数越高，相互作用越具体。

cpdb_results['interaction_scores'].head(2)

CellSign字段：

• interacting_pair：相互作用对的名称。
• partner A 或 B：第一个相互作用伙伴 (A) 或第二个 (B) 的标识符。可以是：UniProt（前缀 simple:）或复合体（前缀 complex:）。
• gene A 或 B：第一个相互作用伙伴 (A) 或第二个 (B) 的基因标识符。标识符将取决于用户输入的列表。
• Receptor A 或 B：如果第一个相互作用伙伴 (A) 或第二个 (B) 在我们的数据库中被注释为受体则为真。
• TF relationship：如果下游的 TF 活跃则为 1。

cpdb_results['CellSign_active_interactions'].head(2)

探索 CellphoneDB 结果

此模块允许通过指定：细胞类型对、基因或特定相互作用来筛选 CellphoneDB 结果。

我们可以指定两个细胞类型列表 (query_cell_type_1 和 query_cell_type_2)，该方法将对每个列表中定义的细胞类型进行成对比较，并将相互作用子集化为这些细胞对。如果我们有兴趣筛选在给定细胞与数据集中所有其他细胞之间发生的相互作用，我们可以定义 query_cell_type_1 = 'All' 和 query_cell_type_2 = ['cellA', 'cellB', ...]。参数 genes 允许用户筛选基因参与的相互作用，用户还可以根据相互作用名称 query_interactions 定义特定相互作用。

在这个例子中，我们将筛选一部分滋养层细胞与 PV 细胞相互作用的结果，此外，我们还将选择基因 TGFBR1 参与的相互作用以及名为 CSF1_CSF1R 的特定相互作用。

此方法筛选显著均值文件的行和列；NaN 值对应于未发现显著的相互作用对，non-NaN 值对应于发现显著的相互作用对。

from cellphonedb.utils import search_utils

search_results = search_utils.search_analysis_results(
    query_cell_types_1=['EVT_1', 'EVT_2', 'GC', 'eEVT', 'iEVT'],  # 细胞列表 1，将与细胞列表 2 配对（列表或 'All'）。
    query_cell_types_2=['PV MMP11', 'PV MYH11', 'PV STEAP4'],     # 细胞列表 2，将与细胞列表 1 配对（列表或 'All'）。
    query_genes=['TGFBR1'],                                       # 基于参与的基因筛选相互作用（列表）。
    query_interactions=['CSF1_CSF1R'],                            # 基于相互作用名称筛选相互作用（列表）。
    significant_means=cpdb_results['significant_means'],          # CellphoneDB 生成的 significant_means 文件。
    deconvoluted=cpdb_results['deconvoluted'],                    # CellphoneDB 生成的 deconvoluted 文件。
    separator='|',                                                # 用于拆分细胞的分隔符（默认：|）。
    long_format=True,                                             # 将输出转换为宽表格，移除非显著相互作用。
    query_classifications=['Signaling by Transforming growth factor']
)

search_results.head(3)

基本分析和绘图

在本节中，我们将对 CellphoneDB 返回的结果进行基本分析和绘图。用户熟悉数据的生物学背景对于从 CellphoneDB 结果中获得有意义的结果至关重要。

细胞-细胞相互作用数据可以通过多种方式绘制；描述细胞对之间相关相互作用数量的弦图/热图或点图，表示特定伙伴和相互作用细胞等。在本节中，我们选择点图，因为这些提供了相互作用和相互作用者的详细信息。

首先，我们将 relevant_interactions 和 means 的结果转换为长表格，因为这更方便工作，并且是大多数绘图库的输入格式。

# 获取注释列和相互作用
annotation = list(cpdb_results['relevant_interactions'].columns[:13])
interaction = list(cpdb_results['relevant_interactions'].columns[13:])

# 将 relevant_interactions 文件从宽表格转换为长表格
relevant_interactions_long = pd.melt(cpdb_results['relevant_interactions'],
                                     id_vars=annotation,
                                     var_name='Interacting_cell',
                                     value_vars=interaction,
                                     value_name='Relevance')

relevant_interactions_long[['Cell_a', 'Cell_b']] = relevant_interactions_long['Interacting_cell'] \
    .str.split('|', expand=True) \
    .rename(columns={0: 'Cell_a', 1: 'Cell_b'})

relevant_interactions_long.head(3)

# 将 means 文件从宽表格转换为长表格
means_long = pd.melt(cpdb_results['means'],
                     id_vars=annotation,
                     var_name='Interacting_cell',
                     value_vars=interaction,
                     value_name='Mean')

means_long.head(3)

按复发排序相互作用

为了方便选择感兴趣的相互作用，我们将计算每个相互作用的复发（即在所有相互作用伙伴中相互作用出现的次数）。在这里，我们假设我们对特定的相互作用感兴趣（即在有限数量的细胞之间发现的相互作用）。此表格与先前知识、其他数据集（如 visium、成像等）的验证一起，可以选择/验证最相关的相互作用。

# 创建一个字典，用于记录任何给定相互作用的复发
id_cp_dict = relevant_interactions_long.groupby('id_cp_interaction')['Relevance'].sum().to_dict()

# 添加新列以指示相互作用的复发
relevant_interactions_long['Recurrence'] = relevant_interactions_long['id_cp_interaction'].map(id_cp_dict)

# 按复发排序
relevant_interactions_long = relevant_interactions_long.sort_values(['Recurrence'], ascending=True)

relevant_inter

relevant_interactions_long.head(3)

基于我们自己的知识和与实验数据的验证，我们选择了一组相互作用，其 ID 如下。

interactions_to_plot = ['CPI-SC0F36DD437', 'CPI-SC06C3BA86F', 'CPI-SS0AB1F04A0', 'CPI-SC051FC855D', 'CPI-SC0DABCC4B5']

idx = [id in interactions_to_plot for id in relevant_interactions_long['id_cp_interaction']]
relevant_interactions_plot = relevant_interactions_long[idx].copy()

现在我们从之前转换为长格式的 means 表中添加相互作用伙伴的平均表达值。由于不同相互作用伙伴的平均值不能直接与其他相互作用伙伴的值进行比较，我们建议在每个相互作用内进行归一化。在这种情况下，我们应用 z 评分归一化。

from scipy.stats import zscore

# 添加相互作用伙伴的平均值
relevant_interactions_plot = relevant_interactions_plot.merge(means_long[['id_cp_interaction', 'Interacting_cell', 'Mean']],
                                                              on=['id_cp_interaction', 'Interacting_cell'],
                                                              how='inner')

在应用 z 评分归一化之前，我们将移除两个相互作用伙伴都是滋养层的相互作用。

import itertools

list_trophoblast = ['EVT_1', 'EVT_2', 'GC', 'VCT_CCC', 'eEVT', 'iEVT']
interact_rm = ['|'.join(i) for i in itertools.product(list_trophoblast, repeat=2)]

# 保留的相互作用
idx_keep = [i not in interact_rm for i in relevant_interactions_plot['Interacting_cell']]
relevant_interactions_plot = relevant_interactions_plot[idx_keep]

# 在相互作用内进行归一化
relevant_interactions_plot['Mean_scaled'] = relevant_interactions_plot.groupby('id_cp_interaction', group_keys=False)['Mean'].transform(lambda x: zscore(x, ddof=1))
relevant_interactions_plot = relevant_interactions_plot.sort_values('Interacting_cell')

relevant_interactions_plot.head(3)

绘制所选候选相互作用的点图

import seaborn as sns

g = sns.relplot(
    data=relevant_interactions_plot,
    x="Interacting_cell",
    y="interacting_pair",
    hue="Relevance",
    size="Mean_scaled",
    palette="vlag",
    hue_norm=(-1, 1),
    height=4,
    aspect=3,
    sizes=(0, 200)
)
g.set_xticklabels(rotation=90)

或者，我们可以使用 Kelvin 的实现（kt-plots tutorial）来显示 CellphoneDB 的输出。在这个示例中，我们将绘制由 TGFB2 和 CSF1R 介导的 PVs 和滋养层细胞之间的相互作用。

import os
import anndata as ad
import pandas as pd
import ktplotspy as kpy
import matplotlib.pyplot as plt
%matplotlib inline

kpy.plot_cpdb_heatmap(pvals=cpdb_results['relevant_interactions'],
                      degs_analysis=True,
                      figsize=(5, 5),
                      title="Sum of significant interactions")

kpy.plot_cpdb(
    adata=adata,
    cell_type1="PV MYH11|PV STEAP4|PV MMPP11",
    cell_type2="EVT_1|EVT_2|GC|iEVT|eEVT|VCT_CCC", 
    means=cpdb_results['means'],
    pvals=cpdb_results['relevant_interactions'],
    celltype_key="cell_labels",
    genes=["TGFB2", "CSF1R"],
    figsize=(11,4),
    title="Interactions between PV and trophoblast",
    max_size=5,
    highlight_size=0.75,
    degs_analysis=True,
    standard_scale=True,
    interaction_scores=cpdb_results['interaction_scores'],
    scale_alpha_by_interaction_scores=True,
)

绘制考虑了 CellSign 结果（受体下游的活跃 TFs）的相互作用。

kpy.plot_cpdb(
    adata=adata,
    cell_type1="PV MYH11",
    cell_type2="EVT_1|EVT_2|GC|iEVT|eEVT|VCT_CCC",
    means=cpdb_results['means'],
    pvals=cpdb_results['relevant_interactions'],
    celltype_key="cell_labels",
    figsize=(9, 4),
    title="Interactions between\nPV and trophoblast with\ndownstream significance",
    max_size=6,
    highlight_size=0.75,
    degs_analysis=True,
    standard_scale=True,
    cellsign=cpdb_results['CellSign_active_interactions'],
    filter_by_cellsign=True,
    scale_alpha_by_cellsign=True
)

相互作用也可以按通路分组绘制。

from plotnine import facet_wrap

p = kpy.plot_cpdb(
    adata=adata,
    cell_type1="PV MYH11",
    cell_type2="EVT_1|EVT_2|GC|iEVT|eEVT|VCT_CCC",
    means=cpdb_results['means'],
    pvals=cpdb_results['relevant_interactions'],
    celltype_key="cell_labels",
    genes=["TGFB2", "CSF1R", "COL1A1"],
    figsize=(12, 8),
    title="Interactions between PV and trophoblasts\n grouped by classification",
    max_size=6,
    highlight_size=0.75,
    degs_analysis=True,
    standard_scale=True
)

p + facet_wrap("~ classification", ncol=1)