16.1 前言

單細(xì)胞RNA-seq為了解不同條件日川、組織類型、物種和個(gè)體之間細(xì)胞類型的變化提供了前所未有的信息。單細(xì)胞數(shù)據(jù)的差異基因表達(dá)分析幾乎總是隨后進(jìn)行基因集富集分析浦旱，其目的是識(shí)別與實(shí)驗(yàn)條件相比在實(shí)驗(yàn)條件下過度表達(dá)的基因程序直奋，例如生物過程能庆、基因本體或調(diào)控通路〗畔撸基于差異表達(dá)（DE）基因的對(duì)照或其他條件搁胆。
為了確定兩種條件之間以細(xì)胞類型特異性方式富集的通路，首先選擇基因集特征的相關(guān)集合邮绿，其中每個(gè)基因集定義一個(gè)生物過程（例如上皮到間質(zhì)轉(zhuǎn)化渠旁、代謝等）或通路（例如MAPK） 信號(hào)）。對(duì)于集合中的每個(gè)基因集船逮，基因集中存在的 DE 基因用于獲得測(cè)試統(tǒng)計(jì)量一死，然后使用該統(tǒng)計(jì)量來評(píng)估基因集的富集度。根據(jù)所選富集測(cè)試的類型傻唾，基因表達(dá)測(cè)量可能會(huì)或可能不會(huì)用于測(cè)試統(tǒng)計(jì)量的計(jì)算投慈。
在本篇中承耿，我們首先概述不同類型的基因集富集測(cè)試，介紹一些常用的基因特征集合伪煤，并討論通路富集和功能富集分析的最佳實(shí)踐加袋。

16.2 通路和基因集

基因集是通過先前的研究和/或?qū)嶒?yàn)已知參與生物過程的基因名稱（或基因 ID）的精選列表。分子特征數(shù)據(jù)庫 (MSigDB) 是最全面的數(shù)據(jù)庫抱既，由9個(gè)基因集集合組成职烧。一些常用的集合是C5，它是基因本體 (GO) 集合防泵，C2 是來自已發(fā)表的研究的精選基因簽名的集合蚀之，這些研究通常是上下文（例如組織、條件）特定的捷泞，但也包括KEGG和REACTOME基因特征集合足删。對(duì)于癌癥研究，通常使用Hallmark系列锁右，而對(duì)于免疫學(xué)研究失受，C7系列是常見的選擇。請(qǐng)注意咏瑟，這些特征主要來自Bulk-seq測(cè)量并測(cè)量連續(xù)表型拂到。最近，隨著scRNA-seq數(shù)據(jù)集的廣泛使用码泞，數(shù)據(jù)庫不斷發(fā)展兄旬，提供源自已發(fā)表的單細(xì)胞研究的精選標(biāo)記列表，這些單細(xì)胞研究定義了各種組織和物種中的細(xì)胞類型余寥。其中包括CellMarker[Zhang et al., 2019]和PanglaoDB[Franzén et al., 2019]辖试。

16.3 基因集測(cè)試和通路分析

在scRNA-seq數(shù)據(jù)分析中，基因集富集通常在細(xì)胞簇或細(xì)胞類型上進(jìn)行劈狐，一次一個(gè)罐孝。例如，使用簡(jiǎn)單的超幾何檢驗(yàn)或Fisher精確檢驗(yàn)肥缔，在簇或細(xì)胞類型中差異表達(dá)的基因用于從所選集合中識(shí)別過度代表性的基因集莲兢。
fgsea[Korotkevich et al., 2021]是一種更常見的基因集富集測(cè)試工具。fgsae是成熟的基因集富集分析 (GSEA) 算法的計(jì)算速度更快的實(shí)現(xiàn)续膳，該算法根據(jù)一些預(yù)先排序的基因水平測(cè)試統(tǒng)計(jì)數(shù)據(jù)計(jì)算富集統(tǒng)計(jì)數(shù)據(jù)改艇。fgsea使用基因集中基因的一些帶符號(hào)統(tǒng)計(jì)數(shù)據(jù)來計(jì)算富集分?jǐn)?shù)，例如t統(tǒng)計(jì)量坟岔、對(duì)數(shù)倍數(shù)變化 (logFC) 或差異表達(dá)測(cè)試的p值谒兄。使用一些相同大小的隨機(jī)基因集計(jì)算富集分?jǐn)?shù)的經(jīng)驗(yàn)（根據(jù)數(shù)據(jù)估計(jì)）零分布，并計(jì)算p值以確定富集分?jǐn)?shù)的顯著性社付。然后針對(duì)多重假設(shè)檢驗(yàn)調(diào)整p值承疲。 GSVA [H?nzelmann et al., 2013]是預(yù)排序基因集富集方法的另一個(gè)例子邻耕。我們應(yīng)該注意到，預(yù)先排序的基因集測(cè)試并不特定于單細(xì)胞數(shù)據(jù)集燕鸽，也適用于批量測(cè)序分析兄世。
測(cè)試一組細(xì)胞（即相同類型的簇或細(xì)胞）中基因集富集的另一種方法是從單細(xì)胞創(chuàng)建偽批量樣本，并使用為bulk RNA-seq開發(fā)的基因集富集方法啊研。limma中實(shí)現(xiàn)了幾種獨(dú)立且競(jìng)爭(zhēng)性的基因集富集測(cè)試御滩，即Fried和Camera，這些測(cè)試通過線性模型和測(cè)試統(tǒng)計(jì)的經(jīng)驗(yàn)貝葉斯調(diào)節(jié)與差異基因表達(dá)分析框架兼容党远。線性模型可以通過設(shè)計(jì)矩陣適應(yīng)復(fù)雜的實(shí)驗(yàn)設(shè)計(jì)（例如主題削解、擾動(dòng)、批次沟娱、嵌套對(duì)比氛驮、交互等）。 limma 中的基因集測(cè)試也可以應(yīng)用于（適當(dāng)轉(zhuǎn)化和歸一化的）單細(xì)胞測(cè)量花沉，而無需偽批量生成柳爽。然而媳握，目前還沒有基準(zhǔn)可以評(píng)估當(dāng)這些方法直接應(yīng)用于單細(xì)胞時(shí)基因組測(cè)試結(jié)果的準(zhǔn)確性碱屁。

16.4 實(shí)戰(zhàn)：人PBMC單細(xì)胞的通路富集分析和評(píng)分

16.4.1 準(zhǔn)備和探索數(shù)據(jù)

我們首先下載25K PBMC數(shù)據(jù)，并遵循標(biāo)準(zhǔn)scanpy工作流程蛾找，對(duì)讀數(shù)計(jì)數(shù)進(jìn)行歸一化并對(duì)高度可變的基因進(jìn)行子集化娩脾。該數(shù)據(jù)集包含未經(jīng)處理的和IFN-β刺激人類PBMC細(xì)胞[Kang et al., 2018]。 我們利用4000個(gè)高度可變基因的UMAP表示來探索數(shù)據(jù)的變異模式打毛。

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'

雖然當(dāng)前對(duì)象帶有UMAP和PCA嵌入柿赊，但這些嵌入已針對(duì)刺激條件進(jìn)行了校正，我們將重新計(jì)算這些幻枉。

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單核細(xì)胞（FCGR3A+ 單核細(xì)胞）簇中響應(yīng)IFN刺激而差異表達(dá)的基因的排名。我們使用這些排名和REACTOME中的基因集來查找該細(xì)胞群中富集的基因集椿肩，與使用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進(jìn)行簇級(jí)基因集富集分析

現(xiàn)在我們將使用python包decoupler[Badia-i-Mompel et al., 2022]對(duì)我們的數(shù)據(jù)執(zhí)行GSEA富集測(cè)試。

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查詢這些資源。

然而厂榛，為了本教程的可重復(fù)性盖矫，我們使用基因集集合的固定版本丽惭。

# 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 運(yùn)行GSEA

首先，我們將準(zhǔn)備我們的基因集炼彪。默認(rèn)情況下吐根，decoupler不會(huì)按最大大小過濾基因集。相反辐马，我們將簡(jiǎn)單地手動(dòng)過濾基因集拷橘，使其最少有15個(gè)基因，最多有500個(gè)基因喜爷。

# 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檢驗(yàn)中的t統(tǒng)計(jì)量對(duì)IFN刺激后CD16單核細(xì)胞表型的基因進(jìn)行排序冗疮，并計(jì)算每個(gè)途徑的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")
)

我們繪制了與所有其他細(xì)胞類型相比檩帐，受刺激的CD16單核細(xì)胞顯著富集的前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調(diào)整后的p值诅挑。因此，條形的高度越長(zhǎng)泛源，路徑越重要拔妥。路徑按重要性排序。大多數(shù)干擾素相關(guān)通路確實(shí)位列前20個(gè)最顯著富集通路之列达箍。一些IFN相關(guān)途徑包括REACTOME_INTERFERON_SIGNALING（排名第 2）没龙、REACTOME_INTERFERON_GAMMA_SIGNALING（排名第 3）和 REACTOME_INTERFERON_ALPHA_BETA_SIGNALING（排名第 4）《忻担總體而言硬纤，鑒于我們先驗(yàn)地知道IFN相關(guān)通路應(yīng)該是排名最高的術(shù)語，GSEA在識(shí)別已知與干擾素信號(hào)傳導(dǎo)相關(guān)的通路方面做得不錯(cuò)赃磨。

讓我們看看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是富集測(cè)試的p值，而score和norm分別是富集分?jǐn)?shù)和歸一化富集分?jǐn)?shù)邻辉。負(fù)分?jǐn)?shù)表明該途徑被下調(diào)溪王，正分?jǐn)?shù)表明該途徑或基因集中的基因上調(diào)。

16.4.3 使用AUCell進(jìn)行細(xì)胞水平通路活性評(píng)分

與之前我們?cè)u(píng)估每個(gè)簇（或更確切地說細(xì)胞類型）的基因集富集度的方法不同恩沛，我們可以對(duì)每個(gè)單獨(dú)細(xì)胞中的通路和基因集的活性水平進(jìn)行評(píng)分在扰，這是基于細(xì)胞中的絕對(duì)基因表達(dá)。 我們可以通過AUCell等評(píng)分工具來實(shí)現(xiàn)這一點(diǎn)雷客。

與GSEA類似芒珠，我們將使用AUCell的decoupler實(shí)現(xiàn)。

%%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'

現(xiàn)在搅裙，我們將干擾素相關(guān)REACTOME通路的分?jǐn)?shù)添加到AnnData對(duì)象的obs字段皱卓，并在UMAP上的每個(gè)細(xì)胞中注釋這些通路的活性水平：

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中繪制分?jǐn)?shù)情況裹芝。

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

16.4.4 使用 limma-fry和pseudo-bulks進(jìn)行復(fù)雜實(shí)驗(yàn)設(shè)計(jì)的基因集富集

在簇水平t檢驗(yàn)方法中怜械，通過將一個(gè)簇與所有其他簇（在這種情況下包括對(duì)照細(xì)胞和刺激細(xì)胞）進(jìn)行比較來發(fā)現(xiàn)差異表達(dá)的基因。線性模型使我們能夠僅將刺激條件下的細(xì)胞與對(duì)照組中的細(xì)胞進(jìn)行比較傅事，從而更準(zhǔn)確地識(shí)別對(duì)刺激做出反應(yīng)的基因缕允。事實(shí)上，線性模型可以適應(yīng)復(fù)雜的實(shí)驗(yàn)設(shè)計(jì)蹭越，例如障本，與處理2中的細(xì)胞類型A相比，處理1中細(xì)胞類型A富集的基因組的鑒定响鹃；也就是說驾霜，跨擾動(dòng)和跨細(xì)胞類型效應(yīng)，同時(shí)調(diào)整批次效應(yīng)买置、個(gè)體差異粪糙、小鼠模型中的性別和品系差異等。
在下一節(jié)中堕义，我們將演示limma-fry工作流程猜旬，該工作流程可推廣到現(xiàn)實(shí)的數(shù)據(jù)分析例程脆栋，例如單細(xì)胞病例對(duì)照研究倦卖。我們首先為每個(gè)細(xì)胞類型和條件創(chuàng)建偽批量重復(fù)（每個(gè)條件-細(xì)胞類型組合3個(gè)重復(fù)）。然后椿争，我們發(fā)現(xiàn)在某種細(xì)胞類型中怕膛，與對(duì)照細(xì)胞相比，刺激細(xì)胞中的基因集更加豐富秦踪。我們還評(píng)估兩種受刺激的細(xì)胞類型群體之間的基因集富集褐捻，以發(fā)現(xiàn)信號(hào)通路的差異。

16.4.4.1 創(chuàng)建pseudo-bulk樣本并探索數(shù)據(jù)

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)

標(biāo)準(zhǔn)化這些數(shù)據(jù)并快速瀏覽一下椅邓。我們不會(huì)在這里使用neighbor embedding柠逞，因?yàn)闃颖玖匡@著減少。

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現(xiàn)在捕獲淋巴（T景馁、NK板壮、B）和骨髓（Mono、DC）群體之間的差異合住，而第二個(gè)PC捕獲由于執(zhí)行刺激而產(chǎn)生的變化（即對(duì)照和刺激的偽重復(fù)之間的差異）绰精。理想情況下撒璧，必須在偽批量數(shù)據(jù)的前幾個(gè)PC中檢測(cè)到感興趣的變化。
在這種情況下笨使，由于我們確實(shí)對(duì)每種細(xì)胞類型的刺激效果感興趣卿樱，因此我們繼續(xù)進(jìn)行基因組測(cè)試。我們重申硫椰，繪制PC的目的是探索數(shù)據(jù)中的各個(gè)變異軸繁调，并發(fā)現(xiàn)可能對(duì)測(cè)試結(jié)果產(chǎn)生重大影響的不需要的變異。如果用戶對(duì)數(shù)據(jù)的變化感到滿意靶草，則可以繼續(xù)進(jìn)行其余的分析涉馁。

16.4.4.2 配置limma和fry

對(duì)于下一部分的分析，我們將使用Bioconductor包limma及其方法fry爱致。
我們首先設(shè)置設(shè)計(jì)和對(duì)比度矩陣烤送。

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
)

對(duì)我們數(shù)據(jù)中每個(gè)通路中注釋的基因進(jìn)行索引，如下所示：

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個(gè)基因的基因集

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

現(xiàn)在我們已經(jīng)設(shè)置了設(shè)計(jì)和對(duì)比矩陣帮坚，并對(duì)數(shù)據(jù)中每個(gè)通路中的基因進(jìn)行了索引，我們可以調(diào)用fry()來測(cè)試我們上面設(shè)置的每個(gè)對(duì)比中的富集通路：

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])

看看排名靠前的通路互艾，我們會(huì)看到一些熟悉的名字：

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 比較兩種刺激細(xì)胞類型的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())
)

如前所述勺良，理想情況下，必須在偽批量數(shù)據(jù)的前幾個(gè)PC中檢測(cè)到感興趣的變化骄噪。讓我們刪除數(shù)據(jù)中低表達(dá)的基因尚困，應(yīng)用log2CPM轉(zhuǎn)換并重復(fù)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)

這里，“l(fā)ogCPM_FLE”表示過濾低表達(dá)基因链蕊，然后進(jìn)行l(wèi)og2CPM轉(zhuǎn)換∈绿穑現(xiàn)在，當(dāng)去除低表達(dá)基因并通過log2CPM轉(zhuǎn)換調(diào)整文庫大小之間的差異時(shí)滔韵，我們現(xiàn)在可以清楚地觀察到PC1捕獲細(xì)胞類型效應(yīng)逻谦，PC2捕獲治療效應(yīng)。

由于在本案例研究中陪蜻，我們確實(shí)對(duì)每種細(xì)胞類型的刺激效果感興趣邦马，并且這種變異在基因過濾之前得到了更好的保留，因此我們展示了未過濾數(shù)據(jù)的富集測(cè)試結(jié)果宴卖。在實(shí)踐中滋将，過濾低豐度基因并通過edgeR::calcNormFactors計(jì)算標(biāo)準(zhǔn)化因子是bulk RNA-seq分析工作流程的標(biāo)準(zhǔn)部分。如果我們對(duì)干擾素刺激的整體影響感興趣症昏，我們應(yīng)該使用過濾后的數(shù)據(jù)随闽。此外，可以注意到design <- model.matrix(~ 0 + lineage + group)將考慮骨髓譜系和淋巴譜系之間的差異（即基線表達(dá)差異）肝谭，通過IFN刺激改善偽大量樣品的分離掘宪， 可能沿著PC1。在本案例研究中攘烛，我們對(duì)細(xì)胞類型特異性效應(yīng)感興趣魏滚，因此我們保留了一個(gè)數(shù)據(jù)模型，其中PC1的變異性隨細(xì)胞類型而變化医寿。必須仔細(xì)考慮設(shè)計(jì)矩陣的選擇栏赴，以與感興趣的生物學(xué)問題保持一致蘑斧。