16.1 前言

单细胞RNA-seq为了解不同条件、组织类型、物种和个体之间细胞类型的变化提供了前所未有的信息。单细胞数据的差异基因表达分析几乎总是随后进行基因集富集分析，其目的是识别与实验条件相比在实验条件下过度表达的基因程序，例如生物过程、基因本体或调控通路。基于差异表达（DE）基因的对照或其他条件。
为了确定两种条件之间以细胞类型特异性方式富集的通路，首先选择基因集特征的相关集合，其中每个基因集定义一个生物过程（例如上皮到间质转化、代谢等）或通路（例如MAPK） 信号）。对于集合中的每个基因集，基因集中存在的 DE 基因用于获得测试统计量，然后使用该统计量来评估基因集的富集度。根据所选富集测试的类型，基因表达测量可能会或可能不会用于测试统计量的计算。
在本篇中，我们首先概述不同类型的基因集富集测试，介绍一些常用的基因特征集合，并讨论通路富集和功能富集分析的最佳实践。

16.2 通路和基因集

基因集是通过先前的研究和/或实验已知参与生物过程的基因名称（或基因 ID）的精选列表。分子特征数据库 (MSigDB) 是最全面的数据库，由9个基因集集合组成。一些常用的集合是C5，它是基因本体 (GO) 集合，C2 是来自已发表的研究的精选基因签名的集合，这些研究通常是上下文（例如组织、条件）特定的，但也包括KEGG和REACTOME基因特征集合。对于癌症研究，通常使用Hallmark系列，而对于免疫学研究，C7系列是常见的选择。请注意，这些特征主要来自Bulk-seq测量并测量连续表型。最近，随着scRNA-seq数据集的广泛使用，数据库不断发展，提供源自已发表的单细胞研究的精选标记列表，这些单细胞研究定义了各种组织和物种中的细胞类型。其中包括CellMarker[Zhang et al., 2019]和PanglaoDB[Franzén et al., 2019]。

16.3 基因集测试和通路分析

在scRNA-seq数据分析中，基因集富集通常在细胞簇或细胞类型上进行，一次一个。例如，使用简单的超几何检验或Fisher精确检验，在簇或细胞类型中差异表达的基因用于从所选集合中识别过度代表性的基因集。
fgsea[Korotkevich et al., 2021]是一种更常见的基因集富集测试工具。fgsae是成熟的基因集富集分析 (GSEA) 算法的计算速度更快的实现，该算法根据一些预先排序的基因水平测试统计数据计算富集统计数据。fgsea使用基因集中基因的一些带符号统计数据来计算富集分数，例如t统计量、对数倍数变化 (logFC) 或差异表达测试的p值。使用一些相同大小的随机基因集计算富集分数的经验（根据数据估计）零分布，并计算p值以确定富集分数的显着性。然后针对多重假设检验调整p值。 GSVA [Hänzelmann et al., 2013]是预排序基因集富集方法的另一个例子。我们应该注意到，预先排序的基因集测试并不特定于单细胞数据集，也适用于批量测序分析。
测试一组细胞（即相同类型的簇或细胞）中基因集富集的另一种方法是从单细胞创建伪批量样本，并使用为bulk RNA-seq开发的基因集富集方法。limma中实现了几种独立且竞争性的基因集富集测试，即Fried和Camera，这些测试通过线性模型和测试统计的经验贝叶斯调节与差异基因表达分析框架兼容。线性模型可以通过设计矩阵适应复杂的实验设计（例如主题、扰动、批次、嵌套对比、交互等）。 limma 中的基因集测试也可以应用于（适当转化和归一化的）单细胞测量，而无需伪批量生成。然而，目前还没有基准可以评估当这些方法直接应用于单细胞时基因组测试结果的准确性。

16.4 实战：人PBMC单细胞的通路富集分析和评分

16.4.1 准备和探索数据

我们首先下载25K PBMC数据，并遵循标准scanpy工作流程，对读数计数进行归一化并对高度可变的基因进行子集化。该数据集包含未经处理的和IFN-β刺激人类PBMC细胞[Kang et al., 2018]。 我们利用4000个高度可变基因的UMAP表示来探索数据的变异模式。

from __future__ import annotations

import numpy as np
import pandas as pd

import scanpy as sc
import anndata as ad
import decoupler
import seaborn.objects as so

import session_info

sc.settings.set_figure_params(dpi=200, frameon=False)
sc.set_figure_params(dpi=200)
sc.set_figure_params(figsize=(4, 4))

# Filtering warnings from current version of matplotlib
import warnings

warnings.filterwarnings(
    "ignore", message=".*Parameters 'cmap' will be ignored.*", category=UserWarning
)
warnings.filterwarnings(
    "ignore", message="Tight layout not applied.*", category=UserWarning
)

# Setting up R dependencies
import anndata2ri
import rpy2
from rpy2.robjects import r
import random

%load_ext rpy2.ipython

anndata2ri.activate()

%%R
suppressPackageStartupMessages({
    library(SingleCellExperiment)
})

adata = sc.read(
    "kang_counts_25k.h5ad", backup_url="https://figshare.com/ndownloader/files/34464122"
)
adata

AnnData object with n_obs × n_vars = 24673 × 15706
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters'
    var: 'name'
    obsm: 'X_pca', 'X_umap'

# Storing the counts for later use
adata.layers["counts"] = adata.X.copy()
# Renaming label to condition
adata.obs = adata.obs.rename({"label": "condition"}, axis=1)

# Normalizing
sc.pp.normalize_total(adata)
sc.pp.log1p(adata)

# Finding highly variable genes using count data
sc.pp.highly_variable_genes(
    adata, n_top_genes=4000, flavor="seurat_v3", subset=False, layer="counts"
)

adata

AnnData object with n_obs × n_vars = 24673 × 15706
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'condition', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters'
    var: 'name', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm'
    uns: 'log1p', 'hvg'
    obsm: 'X_pca', 'X_umap'
    layers: 'counts'

虽然当前对象带有UMAP和PCA嵌入，但这些嵌入已针对刺激条件进行了校正，我们将重新计算这些。

sc.pp.pca(adata)
sc.pp.neighbors(adata)
sc.tl.umap(adata)

/Users/isaac/miniconda3/envs/pathway/lib/python3.9/site-packages/tqdm/auto.py:22: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html
  from .autonotebook import tqdm as notebook_tqdm

sc.pl.umap(
    adata,
    color=["condition", "cell_type"],
    frameon=False,
    ncols=2,
)

adata.obs["group"] = adata.obs.condition.astype("string") + "_" + adata.obs.cell_type

# find DE genes by t-test
sc.tl.rank_genes_groups(adata, "group", method="t-test", key_added="t-test")

让我们提取CD16单核细胞（FCGR3A+ 单核细胞）簇中响应IFN刺激而差异表达的基因的排名。我们使用这些排名和REACTOME中的基因集来查找该细胞群中富集的基因集，与使用GSEA的所有其他群体相比。

celltype_condition = "stim_FCGR3A+ Monocytes"  # 'stimulated_B',  'stimulated_CD8 T', 'stimulated_CD14 Mono'

# extract scores
t_stats = (
    # Get dataframe of DE results for condition vs. rest
    sc.get.rank_genes_groups_df(adata, celltype_condition, key="t-test")
    # Subset to highly variable genes
    .set_index("names")
    .loc[adata.var["highly_variable"]]
    # Sort by absolute score
    .sort_values("scores", key=np.abs, ascending=False)
    # Format for decoupler
    [["scores"]]
    .rename_axis(["stim_FCGR3A+ Monocytes"], axis=1)
)
t_stats

stim_FCGR3A+ Monocytes  scores
names   
IFITM3  123.019180
ISG15   119.732079
TYROBP  91.894241
TNFSF10 87.408890
S100A11 85.721817
... ...
NR1D1   -0.005578
PIK3R5  0.004145
FHL2    0.002915
CLECL1  -0.000262
ADCK4   0.000002

4000 rows × 1 columns

16.4.2 使用decoupler进行簇级基因集富集分析

现在我们将使用python包decoupler[Badia-i-Mompel et al., 2022]对我们的数据执行GSEA富集测试。

16.4.2.1 检索基因集

下载并阅读MSigDB C2集合中注释的REACTOME路径的gmt文件。

from pathlib import Path

if not Path("c2.cp.reactome.v7.5.1.symbols.gmt").is_file():
    !wget -O 'c2.cp.reactome.v7.5.1.symbols.gmt' https://figshare.com/ndownloader/files/35233771

def gmt_to_decoupler(pth: Path) -> pd.DataFrame:
    """
    Parse a gmt file to a decoupler pathway dataframe.
    """
    from itertools import chain, repeat

    pathways = {}

    with Path(pth).open("r") as f:
        for line in f:
            name, _, *genes = line.strip().split("\t")
            pathways[name] = genes

    return pd.DataFrame.from_records(
        chain.from_iterable(zip(repeat(k), v) for k, v in pathways.items()),
        columns=["geneset", "genesymbol"],
    )

reactome = gmt_to_decoupler("c2.cp.reactome.v7.5.1.symbols.gmt")

或者，我们可以只从omnipath查询这些资源。

然而，为了本教程的可重复性，我们使用基因集集合的固定版本。

# Retrieving via python
msigdb = decoupler.get_resource("MSigDB")

# Get reactome pathways
reactome = msigdb.query("collection == 'reactome_pathways'")
# Filter duplicates
reactome = reactome[~reactome.duplicated(("geneset", "genesymbol"))]

reactome

    geneset genesymbol
0   REACTOME_INTERLEUKIN_6_SIGNALING    JAK2
1   REACTOME_INTERLEUKIN_6_SIGNALING    TYK2
2   REACTOME_INTERLEUKIN_6_SIGNALING    CBL
3   REACTOME_INTERLEUKIN_6_SIGNALING    STAT1
4   REACTOME_INTERLEUKIN_6_SIGNALING    IL6ST
... ... ...
89471   REACTOME_ION_CHANNEL_TRANSPORT  FXYD7
89472   REACTOME_ION_CHANNEL_TRANSPORT  UBA52
89473   REACTOME_ION_CHANNEL_TRANSPORT  ATP6V1E2
89474   REACTOME_ION_CHANNEL_TRANSPORT  ASIC5
89475   REACTOME_ION_CHANNEL_TRANSPORT  FXYD1

89476 rows × 2 columns

16.4.2.2 运行GSEA

首先，我们将准备我们的基因集。默认情况下，decoupler不会按最大大小过滤基因集。相反，我们将简单地手动过滤基因集，使其最少有15个基因，最多有500个基因。

# Filtering genesets to match behaviour of fgsea
geneset_size = reactome.groupby("geneset").size()
gsea_genesets = geneset_size.index[(geneset_size > 15) & (geneset_size < 500)]

我们将使用t检验中的t统计量对IFN刺激后CD16单核细胞表型的基因进行排序，并计算每个途径的p值。

scores, norm, pvals = decoupler.run_gsea(
    t_stats.T,
    reactome[reactome["geneset"].isin(gsea_genesets)],
    source="geneset",
    target="genesymbol",
)

gsea_results = (
    pd.concat({"score": scores.T, "norm": norm.T, "pval": pvals.T}, axis=1)
    .droplevel(level=1, axis=1)
    .sort_values("pval")
)

我们绘制了与所有其他细胞类型相比，受刺激的CD16单核细胞显着富集的前20条途径的条形图。

(
    so.Plot(
        data=(
            gsea_results.head(20).assign(
                **{"-log10(pval)": lambda x: -np.log10(x["pval"])}
            )
        ),
        x="-log10(pval)",
        y="source",
    ).add(so.Bar())
)

在上图中，y轴给出了通路名称。x轴描述-log10调整后的p值。因此，条形的高度越长，路径越重要。路径按重要性排序。大多数干扰素相关通路确实位列前20个最显着富集通路之列。一些IFN相关途径包括REACTOME_INTERFERON_SIGNALING（排名第 2）、REACTOME_INTERFERON_GAMMA_SIGNALING（排名第 3）和 REACTOME_INTERFERON_ALPHA_BETA_SIGNALING（排名第 4）。总体而言，鉴于我们先验地知道IFN相关通路应该是排名最高的术语，GSEA在识别已知与干扰素信号传导相关的通路方面做得不错。

让我们看看decoupler.run_gsea的原始输出：

gsea_results.head(10)

    score   norm    pval
source          
REACTOME_NEUTROPHIL_DEGRANULATION   0.624770    5.587953    2.297617e-08
REACTOME_INTERFERON_SIGNALING   0.844158    5.155074    2.535313e-07
REACTOME_INTERFERON_GAMMA_SIGNALING 0.831962    4.137064    3.517788e-05
REACTOME_INTERFERON_ALPHA_BETA_SIGNALING    0.893431    4.107962    3.991655e-05
REACTOME_SIGNALING_BY_INTERLEUKINS  0.376129    3.535838    4.064833e-04
REACTOME_MHC_CLASS_II_ANTIGEN_PRESENTATION  0.701555    3.378736    7.282002e-04
REACTOME_TRANSLATION    -0.628266   -3.277846   1.046026e-03
REACTOME_RRNA_PROCESSING    -0.703607   -3.205217   1.349605e-03
REACTOME_PLATELET_ACTIVATION_SIGNALING_AND_AGGREGATION  0.475259    3.162945    1.561817e-03
REACTOME_LEISHMANIA_INFECTION   0.531540    2.964611    3.030663e-03

在上面，pval是富集测试的p值，而score和norm分别是富集分数和归一化富集分数。负分数表明该途径被下调，正分数表明该途径或基因集中的基因上调。

16.4.3 使用AUCell进行细胞水平通路活性评分

与之前我们评估每个簇（或更确切地说细胞类型）的基因集富集度的方法不同，我们可以对每个单独细胞中的通路和基因集的活性水平进行评分，这是基于细胞中的绝对基因表达。 我们可以通过AUCell等评分工具来实现这一点。

与GSEA类似，我们将使用AUCell的decoupler实现。

%%time
decoupler.run_aucell(
    adata,
    reactome,
    source="geneset",
    target="genesymbol",
    use_raw=False,
)

CPU times: user 16min 42s, sys: 5.09 s, total: 16min 47s
Wall time: 1min 9s

adata

AnnData object with n_obs × n_vars = 24673 × 15706
    obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'condition', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'group'
    var: 'name', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm'
    uns: 'log1p', 'hvg', 'pca', 'neighbors', 'umap', 'condition_colors', 'cell_type_colors', 't-test'
    obsm: 'X_pca', 'X_umap', 'aucell_estimate'
    varm: 'PCs'
    layers: 'counts'
    obsp: 'distances', 'connectivities'

现在，我们将干扰素相关REACTOME通路的分数添加到AnnData对象的obs字段，并在UMAP上的每个细胞中注释这些通路的活性水平：

ifn_pathways = [
    "REACTOME_INTERFERON_SIGNALING",
    "REACTOME_INTERFERON_ALPHA_BETA_SIGNALING",
    "REACTOME_INTERFERON_GAMMA_SIGNALING",
]

adata.obs[ifn_pathways] = adata.obsm["aucell_estimate"][ifn_pathways]

在UMAP中绘制分数情况。

sc.pl.umap(
    adata,
    color=["condition", "cell_type"] + ifn_pathways,
    frameon=False,
    ncols=2,
    wspace=0.3,
)

16.4.4 使用 limma-fry和pseudo-bulks进行复杂实验设计的基因集富集

在簇水平t检验方法中，通过将一个簇与所有其他簇（在这种情况下包括对照细胞和刺激细胞）进行比较来发现差异表达的基因。线性模型使我们能够仅将刺激条件下的细胞与对照组中的细胞进行比较，从而更准确地识别对刺激做出反应的基因。事实上，线性模型可以适应复杂的实验设计，例如，与处理2中的细胞类型A相比，处理1中细胞类型A富集的基因组的鉴定；也就是说，跨扰动和跨细胞类型效应，同时调整批次效应、个体差异、小鼠模型中的性别和品系差异等。
在下一节中，我们将演示limma-fry工作流程，该工作流程可推广到现实的数据分析例程，例如单细胞病例对照研究。我们首先为每个细胞类型和条件创建伪批量重复（每个条件-细胞类型组合3个重复）。然后，我们发现在某种细胞类型中，与对照细胞相比，刺激细胞中的基因集更加丰富。我们还评估两种受刺激的细胞类型群体之间的基因集富集，以发现信号通路的差异。

16.4.4.1 创建pseudo-bulk样本并探索数据

def subsampled_summation(
    adata: ad.AnnData,
    groupby: str | list[str],
    *,
    n_samples_per_group: int,
    n_cells: int,
    random_state: None | int | np.random.RandomState = None,
    layer: str = None,
) -> ad.AnnData:
    """
    Sum sample of X per condition.

    Drops conditions which don't have enough samples.

    Parameters
    ----------
    adata
        AnnData to sum expression of
    groupby
        Keys in obs to groupby
    n_samples_per_group
        Number of samples to take per group
    n_cells
        Number of cells to take per sample
    random_state
        Random state to use when sampling cells
    layer
        Which layer of adata to use

    Returns
    -------
    AnnData with same var as original, obs with columns from groupby, and X.
    """
    from scipy import sparse
    from sklearn.utils import check_random_state

    # Checks
    if isinstance(groupby, str):
        groupby = [groupby]
    random_state = check_random_state(random_state)

    indices = []
    labels = []

    grouped = adata.obs.groupby(groupby)
    for k, inds in grouped.indices.items():
        # Check size of group
        if len(inds) < (n_cells * n_samples_per_group):
            continue

        # Sample from group
        condition_inds = random_state.choice(
            inds, n_cells * n_samples_per_group, replace=False
        )
        for i, sample_condition_inds in enumerate(np.split(condition_inds, 3)):
            if isinstance(k, tuple):
                labels.append((*k, i))
            else:  # only grouping by one variable
                labels.append((k, i))
            indices.append(sample_condition_inds)

    # obs of output AnnData
    new_obs = pd.DataFrame.from_records(
        labels,
        columns=[*groupby, "sample"],
        index=["-".join(map(str, l)) for l in labels],
    )
    n_out = len(labels)

    # Make indicator matrix
    indptr = np.arange(0, (n_out + 1) * n_cells, n_cells)
    indicator = sparse.csr_matrix(
        (
            np.ones(n_out * n_cells, dtype=bool),
            np.concatenate(indices),
            indptr,
        ),
        shape=(len(labels), adata.n_obs),
    )

    return ad.AnnData(
        X=indicator @ sc.get._get_obs_rep(adata, layer=layer),
        obs=new_obs,
        var=adata.var.copy(),
    )

pb_data = subsampled_summation(
    adata, ["cell_type", "condition"], n_cells=75, n_samples_per_group=3, layer="counts"
)
pb_data

AnnData object with n_obs × n_vars = 42 × 15706
    obs: 'cell_type', 'condition', 'sample'
    var: 'name', 'highly_variable', 'highly_variable_rank', 'means', 'variances', 'variances_norm'

# Does PC1 captures a meaningful biological or technical fact?
pb_data.obs["lib_size"] = pb_data.X.sum(1)

标准化这些数据并快速浏览一下。我们不会在这里使用neighbor embedding，因为样本量显着减少。

pb_data.layers["counts"] = pb_data.X.copy()
sc.pp.normalize_total(pb_data)
sc.pp.log1p(pb_data)
sc.pp.pca(pb_data)
sc.pl.pca(pb_data, color=["cell_type", "condition", "lib_size"], ncols=1, size=250)

PC1现在捕获淋巴（T、NK、B）和骨髓（Mono、DC）群体之间的差异，而第二个PC捕获由于执行刺激而产生的变化（即对照和刺激的伪重复之间的差异）。理想情况下，必须在伪批量数据的前几个PC中检测到感兴趣的变化。
在这种情况下，由于我们确实对每种细胞类型的刺激效果感兴趣，因此我们继续进行基因组测试。我们重申，绘制PC的目的是探索数据中的各个变异轴，并发现可能对测试结果产生重大影响的不需要的变异。如果用户对数据的变化感到满意，则可以继续进行其余的分析。

16.4.4.2 配置limma和fry

对于下一部分的分析，我们将使用Bioconductor包limma及其方法fry。
我们首先设置设计和对比度矩阵。

groups = pb_data.obs.condition.astype("string") + "_" + pb_data.obs.cell_type

%%R -i groups
group <-  as.factor(gsub(" |\\+","_", groups))
design <- model.matrix(~ 0 + group)
head(design)

     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11] [,12] [,13] [,14]
[1,]    0    0    1    0    0    0    0    0    0     0     0     0     0     0
[2,]    0    0    1    0    0    0    0    0    0     0     0     0     0     0
[3,]    0    0    1    0    0    0    0    0    0     0     0     0     0     0
[4,]    0    0    0    0    0    0    0    0    0     1     0     0     0     0
[5,]    0    0    0    0    0    0    0    0    0     1     0     0     0     0
[6,]    0    0    0    0    0    0    0    0    0     1     0     0     0     0

%%R
colnames(design)

 [1] "groupctrl_B_cells"           "groupctrl_CD14__Monocytes"  
 [3] "groupctrl_CD4_T_cells"       "groupctrl_CD8_T_cells"      
 [5] "groupctrl_Dendritic_cells"   "groupctrl_FCGR3A__Monocytes"
 [7] "groupctrl_NK_cells"          "groupstim_B_cells"          
 [9] "groupstim_CD14__Monocytes"   "groupstim_CD4_T_cells"      
[11] "groupstim_CD8_T_cells"       "groupstim_Dendritic_cells"  
[13] "groupstim_FCGR3A__Monocytes" "groupstim_NK_cells"        

%%R 
kang_pbmc_con <- limma::makeContrasts(
    
    # the effect if stimulus in CD16 Monocyte cells
    groupstim_FCGR3A__Monocytes - groupctrl_FCGR3A__Monocytes,
    
    # the effect of stimulus in CD16 Monocytes compared to CD8 T Cells
    (groupstim_FCGR3A__Monocytes - groupctrl_FCGR3A__Monocytes) - (groupstim_CD8_T_cells - groupctrl_CD8_T_cells), 
    levels = design
)

对我们数据中每个通路中注释的基因进行索引，如下所示：

log_norm_X = pb_data.to_df().T

%%R -i log_norm_X -i reactome
# Move pathway info from python to R
pathways = split(reactome$genesymbol, reactome$geneset)
# Map gene names to indices
idx = limma::ids2indices(pathways, rownames(log_norm_X))

/Users/isaac/miniconda3/envs/pathway/lib/python3.9/site-packages/rpy2/robjects/pandas2ri.py:54: FutureWarning: iteritems is deprecated and will be removed in a future version. Use .items instead.
  for name, values in obj.iteritems():
/Users/isaac/miniconda3/envs/pathway/lib/python3.9/site-packages/rpy2/robjects/pandas2ri.py:54: FutureWarning: iteritems is deprecated and will be removed in a future version. Use .items instead.
  for name, values in obj.iteritems():

正如gsea方法中所做的那样，让我们删除少于15个基因的基因集

%%R
keep_gs <- lapply(idx, FUN=function(x) length(x) >= 15)
idx <- idx[unlist(keep_gs)]

现在我们已经设置了设计和对比矩阵，并对数据中每个通路中的基因进行了索引，我们可以调用fry()来测试我们上面设置的每个对比中的富集通路：

16.4.4.3 fry test for Stimulated vs Control

%%R -o fry_results
fry_results <- limma::fry(log_norm_X, index = idx, design = design, contrast = kang_pbmc_con[,1])

看看排名靠前的通路，我们会看到一些熟悉的名字：

fry_results.head()

    NGenes  Direction   PValue  FDR PValue.Mixed    FDR.Mixed
REACTOME_INTERFERON_ALPHA_BETA_SIGNALING    57  Up  3.836198e-24    3.410380e-21    8.018820e-39    9.504975e-38
REACTOME_INTERFERON_SIGNALING   177 Up  5.651011e-18    2.511875e-15    3.888212e-51    1.047461e-49
REACTOME_INTERFERON_GAMMA_SIGNALING 84  Up  6.080234e-13    1.801776e-10    4.268886e-61    2.710743e-59
REACTOME_MRNA_SPLICING_MINOR_PATHWAY    50  Down    8.311795e-11    1.847296e-08    5.137952e-20    2.003351e-19
REACTOME_DDX58_IFIH1_MEDIATED_INDUCTION_OF_INTERFERON_ALPHA_BETA    67  Up  1.555236e-09    2.765210e-07    9.966640e-53    3.055291e-51

(
    so.Plot(
        data=(
            fry_results.head(20)
            .assign(**{"-log10(FDR)": lambda x: -np.log10(x["FDR"])})
            .rename_axis(index="Pathway")
        ),
        x="-log10(FDR)",
        y="Pathway",
    ).add(so.Bar())
)

16.4.4.4 比较两种刺激细胞类型的fry test

%%R -o fry_results_negative_ctrl
fry_results_negative_ctrl <- limma::fry(log_norm_X, index = idx, design = design, contrast = kang_pbmc_con[,2])

(
    so.Plot(
        data=(
            fry_results_negative_ctrl.head(20)
            .assign(**{"-log10(FDR)": lambda x: -np.log10(x["FDR"])})
            .rename_axis(index="Pathway")
        ),
        x="-log10(FDR)",
        y="Pathway",
    ).add(so.Bar())
)

如前所述，理想情况下，必须在伪批量数据的前几个PC中检测到感兴趣的变化。让我们删除数据中低表达的基因，应用log2CPM转换并重复PCA图：

counts_df = pb_data.to_df(layer="counts").T

%%R -i counts_df
keep <- edgeR::filterByExpr(counts_df) # in real analysis, supply the desig matrix to the function to retain as more genes as possible
counts_df <- counts_df[keep,]
logCPM <- edgeR::cpm(counts_df, log=TRUE, prior.count = 2)

/Users/isaac/miniconda3/envs/pathway/lib/python3.9/site-packages/rpy2/robjects/pandas2ri.py:54: FutureWarning: iteritems is deprecated and will be removed in a future version. Use .items instead.
  for name, values in obj.iteritems():
R[write to console]: No group or design set. Assuming all samples belong to one group.

%%R -o logCPM
logCPM = data.frame(logCPM)

pb_data.uns["logCPM_FLE"] = logCPM.T  # FLE for filter low exprs
pb_data.obsm["logCPM_FLE_pca"] = sc.pp.pca(logCPM.T.to_numpy(), return_info=False)
sc.pl.embedding(pb_data, "logCPM_FLE_pca", color=pb_data.obs, ncols=1, size=250)

这里，“logCPM_FLE”表示过滤低表达基因，然后进行log2CPM转换。现在，当去除低表达基因并通过log2CPM转换调整文库大小之间的差异时，我们现在可以清楚地观察到PC1捕获细胞类型效应，PC2捕获治疗效应。

由于在本案例研究中，我们确实对每种细胞类型的刺激效果感兴趣，并且这种变异在基因过滤之前得到了更好的保留，因此我们展示了未过滤数据的富集测试结果。在实践中，过滤低丰度基因并通过edgeR::calcNormFactors计算标准化因子是bulk RNA-seq分析工作流程的标准部分。如果我们对干扰素刺激的整体影响感兴趣，我们应该使用过滤后的数据。此外，可以注意到design <- model.matrix(~ 0 + lineage + group)将考虑骨髓谱系和淋巴谱系之间的差异（即基线表达差异），通过IFN刺激改善伪大量样品的分离， 可能沿着PC1。在本案例研究中，我们对细胞类型特异性效应感兴趣，因此我们保留了一个数据模型，其中PC1的变异性随细胞类型而变化。必须仔细考虑设计矩阵的选择，以与感兴趣的生物学问题保持一致。