5.3 Gene set enrichment and pathway analysis¶

Motivation¶

单细胞RNA测序提供了前所未有的洞见，能够洞察不同条件、组织类型、物种和个体之间的细胞类型变化。单细胞数据的差异基因表达分析几乎总是随后进行基因集富集分析，其目标是鉴定在实验条件相比对照或其他条件中过表达的基因程序，如生物过程、基因本体或调控通路，这一判断基于差异表达（DE）基因。

为了确定在两种条件下细胞类型特异性富集的通路，首先选择一个相关的基因集签名集合，其中每个基因集定义一个生物过程（例如上皮细胞到间充质细胞的转变、代谢等）或通路（例如MAPK信号传导）。对于集合中的每个基因集，使用该基因集中的DE基因获得一个测试统计量，然后用它来评估基因集的富集情况。根据选择的富集测试类型，基因表达测量可能会或可能不会用于计算测试统计量。

在本章中，我们首先提供不同类型的基因集富集测试的概览，介绍一些常用的基因签名集合，并讨论通路富集和功能富集分析的最佳实践。我们通过展示三种基因集富集分析的分析方法来结束本章。请注意，在本章中，我们将路径分析、路径富集分析、基因集富集分析和功能分析这几个术语交替使用。

Pathway and gene set collections¶

基因集是一系列经过筛选的基因名称（或基因ID），这些基因已知通过之前的研究和/或实验参与某个生物过程。分子签名数据库（MSigDB）【引用：subramanian2005gene, liberzon2011molecular】是包含9个基因集合集的最全面数据库。一些常用的集合包括C5，即基因本体（GO）集合，C2集合包括来自已发表研究的策划基因签名，这些签名通常是特定于上下文（如组织、条件）的，但也包括KEGG和REACTOME基因签名。对于癌症研究，通常使用Hallmark集合，而对于免疫学研究，则常选择C7集合。请注意，这些签名主要来源于Bulk-seq测量，并测量连续表型。最近，随着scRNA-seq数据集的广泛可用，数据库已经演变为提供来自已发表的单细胞研究的策划标记列表，这些标记定义了各种组织和物种中的细胞类型。这些数据库包括CellMarker【引用：zhang2019cellmarker】和PanglaoDB【引用：franzen2019panglaodb】。策划的标记列表不仅限于数据库中提供的，也可以由个人策划。

Null hypotheses in gene set enrichment analysis¶

基因集测试可以是竞争性的或自包含的，如Goeman和Buhlmann（2007）【引用：goeman2007analyzing】所定义。竞争性基因集测试检测集合中的基因在差异表达方面相对于不在该集合中的基因是否排名较高。这里的抽样单位是基因，因此可以使用单一样本（即单样本GSEA）进行测试。此测试需要不在集合中的基因（即背景基因）。在自包含基因集测试中，抽样单位是受试者，因此每组需要多个样本，但不需要不在集合中的基因。自包含基因集测试检查测试集中的基因是否有差异表达，而不考虑数据集中测量的任何其他基因。这两种零假设之间的区别使得对基因集富集结果的解释有所不同。请注意，在生物数据中存在基因间相关性，即同一通路中的基因表达是相关的。只有少数测试能够适应基因间相关性。我们稍后将讨论这些方法。关于各种基因集测试的详细解释可以在 limma 用户手册中找到。

Gene set tests and pathway analysis¶

在单细胞RNA测序(scRNA-seq)数据分析中，基因集富集通常在单个细胞簇或细胞类型上进行。使用超几何测试或Fisher精确测试（如 Enrichr 中所用）等简单方法，根据在簇或细胞类型中差异表达的基因来识别所选集合中过度表达的基因集。这类测试不需要实际的基因表达测量和读数计数来计算富集统计，因为它们依赖于测试一个基因集中的 $X$ 个基因在实验中差异表达的显著性与该集合中非DE基因的数量比较。

fgsea {cite}korotkevich2021fast 是更常用的基因集富集测试工具。fgsea 是广泛应用的基因集富集分析(GSEA)算法 {cite}subramanian2005gene 的计算效率更高的实现，它根据一些预排列的基因级别测试统计来计算富集统计。fgsea 使用基因集中的基因的一些带符号统计量（例如 t 统计量、对数折叠变化（logFC）或差异表达测试的 p 值）来计算富集得分。通过一些相同大小的随机基因集计算富集得分的经验（根据数据估计）零分布，并计算 p 值以确定富集得分的显著性。然后对 p 值进行多重假设测试的调整。GSVA {cite}hanzelmann2013gsva 是另一个预排列基因集富集方法的例子。应注意，预排列基因集测试不仅限于单细胞数据集，也适用于Bulk-seq测定。

测试细胞群（即相同类型的簇或细胞）中基因集富集的另一种方法是从单细胞创建伪大样本并使用为Bulk RNA-seq开发的基因集富集方法。limma {cite}ritchie2015limma 实现了几种自包含和竞争性基因集富集测试，即 fry 和 camera，这些测试通过线性模型和实证贝叶斯调整测试统计 {cite}smyth2005limma 与差异基因表达分析框架兼容。线性模型可以通过设计矩阵适应复杂的实验设计（例如受试者、干扰、批次、嵌套对比、交互作用等）。此外，limma中实现的camera和roast基因集测试考虑了基因间相关性。limma 中的基因集测试也可以应用于（经适当转化和标准化的）单细胞测量数据，而无需生成伪大样本。然而，目前还没有评估直接应用于单细胞时这些方法基因集测试结果准确性的基准。

测试	整体或单细胞	零假设类型	输入
超几何测试	两者	竞争性	基因计数
费舍尔精确测试	两者	竞争性	基因计数
GSEA$^*$	整体	竞争性	基因排名
GSVA$^*$	整体	竞争性	基因排名
fgsea	两者	竞争性	基因排名
fry$^*$	整体	自包含	表达矩阵
camera$^*$	整体	竞争性	表达矩阵
roast$^*$	整体	自包含	表达矩阵

表格：基因集测试、适用的检测类型及其测试的零假设

$^*$ 这些测试在实际中可适用于单细胞数据集，尽管将其应用于单细胞可能不是常见做法。

Gene set test vs. pathway activity inference¶

基因集测试检查某一通路在一种条件下是否被富集，换句话说，就是在一种条件下相比其他条件是否过度表达，例如，在单核细胞群体中健康捐赠者相比重症COVID-19患者。另一种方法是简单地评估单个细胞中某一通路或基因签名的活性，以绝对意义上而非比较不同条件下的差异活性。一些广泛用于推断单个细胞中的基因集活性（包括通路活性）的工具包括 VISION {cite}detomaso2019functional、AUCell {cite}aibar2017scenic、使用Pagoda2 {cite}fan2016characterizing, lake2018integrative进行通路过度分散分析和简单的组合z得分 {cite}lee2008inferring。

DoRothEA {cite}garcia2019benchmark 和 PROGENy {cite}schubert2018perturbation 是开发用来推断转录因子（TF）- 目标活性的功能分析工具，最初用于Bulk RNA数据。Holland等人 {cite}holland2020robustness 发现在模拟的scRNA-seq数据中，Bulk RNA-seq方法 DoRothEA 和 PROGENy 表现优越，甚至在单细胞数据中的数据丢失事件和低库大小的情况下，部分优于专为scRNA-seq分析设计的工具。Holland等人还得出结论，通路和TF活性推断对基因集的选择比统计方法更敏感。这一观察结果可能是特定于功能富集分析，并且可以通过TF-目标关系是上下文特定的事实来解释；也就是说，一个细胞类型中的TF-目标关联可能实际上与另一细胞类型或组织中的TF-目标关联不同。

与Holand等人的结论相反，Zhang等人 {cite}zhang2020benchmarking 发现，特别是Pagoda2等基于单细胞的工具，在准确性、稳定性和可扩展性的三个不同方面优于基于批量的方法。应注意的是，通路和基因集活性推断工具本质上不考虑批处理效应或除了感兴趣的生物变异之外的生物变异。因此，数据分析师需要确保差异基因表达分析步骤已正确执行。

此外，虽然这里提到的工具在单个细胞中评分每个基因集，但它们不能选择在所有评分基因集中最具生物学相关性的基因集。scDECAF (https://github.com/DavisLaboratory/scDECAF) 是一个基因集活性推断工具，它允许数据驱动的选择最有信息量的基因集，从而帮助解析有意义的细胞异质性。

Technical considerations¶

Filtering out the gene sets with low number of genes¶

常见的做法是在预处理步骤中排除那些与数据或高变异基因（HVG）重叠的少数基因的基因集。张等人 {cite}zhang2020benchmarking 发现，随着基因覆盖率（即通路/基因集中的基因数量）的减少，基于单细胞和基于批量的方法的性能都会下降。Holland等人 {cite}holland2020robustness 也发现，小尺寸的基因集对单细胞数据上的Bulk-seq DoRothEA 和 PROGENy 的性能产生了负面影响。这些报告共同支持了一个观点，即在通路分析中过滤掉基因计数较低的基因集（例如集合中少于10或15个基因）是有益的。Damian & Gorfine (2004) {cite}damian2004statistical 将这一现象归因于基因集中基因数量较少时基因方差更可能较大，而较大的基因集中的基因方差则倾向于较小。这影响了用于测试富集的测试统计的准确性。张等人还发现，通路分析对应用于基因表达测量的标准化程序很敏感。

Data normalization¶

在单细胞实验中，读数通常在预处理流程的早期就进行归一化，以确保不同库大小的细胞间的测量结果可比较。张等人 {cite}zhang2020benchmarking 发现，通过 SCTransform {cite}hafemeister2019normalization 和 scran {cite}lun2016step 进行归一化通常可以提高单细胞和基于批量的通路评分工具的性能。他们发现，AUCell（一种基于排名的方法）和 z-score（转换为零均值，单位标准差）的性能尤其受到不同归一化方法的影响。

Case study: Pathway enrichment analysis and activity level scoring in human PBMC single cells¶

Prepare and explore the data¶

首先，我们下载25K PBMC数据，并遵循scanpy的标准工作流程对读数进行归一化，并对高变异基因进行子集划分。该数据集包含未经处理和经IFN-$\beta$刺激的人类PBMC细胞 {cite}gspa:kang2018。我们通过UMAP表示4000个高变异基因来探索数据中的变异模式。

In [1]:

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

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

In [2]:

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

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

In [3]:

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

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

In [4]:

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

%load_ext rpy2.ipython

anndata2ri.activate()

# Setting up R dependencies import anndata2ri import rpy2 from rpy2.robjects import r import random %load_ext rpy2.ipython anndata2ri.activate()

In [5]:

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

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

In [6]:

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

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

Out[6]:

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'

In [7]:

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

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

In [8]:

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

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

In [9]:

adata

adata

Out[9]:

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'

While the current object comes with UMAP and PCA embeddings, these have been corrected for stimulation condition, which we don't want for this analysis. Instead we will recompute these.

In [10]:

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

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

In [11]:

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

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

我们通常建议按照“差异基因表达”章节的指导确定差异表达基因。为了简化操作，在这里我们使用 scanpy 中的 rank_genes_groups 函数，通过 t 检验来对基因进行排名，根据它们的差异表达测试统计量来评估：

In [12]:

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

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

In [13]:

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

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

要提取对IFN刺激反应在CD16单核细胞（FCGR3A+单核细胞）群体中差异表达的基因排名，我们可以使用这些排名和REACTOME的基因集，在decoupler中实现的GSEA中找出与所有其他细胞群体相比在此细胞群体中富集的基因集。

In [14]:

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

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

In [15]:

# 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

# 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

Out[15]:

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

Cluster-level gene set enrichment analysis with decoupler¶

现在我们将使用Python包decoupler {cite}badia2022decoupler来对我们的数据执行GSEA富集测试。

Retrieving gene sets¶

要下载并读取MSigDB的C2集合中注释的REACTOME通路的gmt文件，请访问MSigDB网站，找到并下载相应的gmt文件。

In [16]:

# Downloading reactome pathways
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

# Downloading reactome pathways 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

In [17]:

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

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

In [18]:

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

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

另外，我们可以直接从omnipath查询这些资源。

然而，为了本教程的稳定性，我们使用了一个固定版本的基因集合。

# 通过Python检索
msigdb = decoupler.get_resource("MSigDB")

# 获取reactome通路
reactome = msigdb.query("collection == 'reactome_pathways'")
# 过滤重复项
reactome = reactome[~reactome.duplicated(("geneset", "genesymbol"))]

In [19]:

reactome

reactome

Out[19]:

	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

Running GSEA¶

首先，我们将准备我们的基因集。默认情况下，decoupler 不会像 fgsea 那样按最大尺寸过滤基因集。相反，我们将手动过滤基因集，使其包含至少15个基因，最多500个基因。

In [20]:

# Filtering genesets to match behaviour of fgsea
geneset_size = reactome.groupby("geneset").size()
gsea_genesets = geneset_size.index[(geneset_size > 15) & (geneset_size < 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值。

In [21]:

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

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

我们将绘制一个条形图，显示在受IFN刺激的CD16单核细胞中与所有其他细胞类型相比显著富集的前20个通路。这个条形图将显示每个通路的富集得分或p值，从而直观地展示哪些生物通路在这些特定细胞中被激活或上调。这种可视化是理解细胞响应机制的有力工具。

In [22]:

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

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

Out[22]:

在上述条形图中，通路名称位于y轴上。x轴表示 $-\log_{10}$ 调整后的p值。因此，条形的高度越长，表示该通路的显著性越高。通路按显著性排序。大多数与干扰素相关的通路确实位于最显著富集的前20个通路之中。一些与IFN相关的通路包括REACTOME_INTERFERON_SIGNALING（排名第二）、REACTOME_INTERFERON_GAMMA_SIGNALING（排名第三）和REACTOME_INTERFERON_ALPHA_BETA_SIGNALING（排名第四）。总体而言，GSEA 在识别已知与干扰素信号传导相关的通路方面做得相当不错，考虑到我们事先知道IFN相关的通路应该是排名最高的术语。

要查看 decoupler.run_gsea 的原始输出，请直接检查或打印这个函数执行后返回的数据结构。通常，这将是一个包含了各种统计结果和基因集详细信息的DataFrame。

In [23]:

gsea_results.head(10)

gsea_results.head(10)

Out[23]:

	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 分别是富集得分和标准化富集得分。值得注意的是，富集得分是有符号的。因此，负得分表明通路被下调，正得分则表明通路或基因集中的基因被上调。这些得分提供了对基因集在特定细胞条件下表达状态变化的定量洞察。

Cell-level pathway activity scoring using AUCell¶

与之前我们按簇（或更准确地说是细胞类型）评估基因集富集的方法不同，我们可以使用像 AUCell 这样的活性评分工具，在每个单独的细胞中评分通路和基因集的活性水平，这是基于细胞中的绝对基因表达，而不考虑其他细胞中的基因表达。这样的方法允许我们获得对每个单个细胞在特定通路或基因集上的活性状态的直接了解。

类似于 GSEA，我们将使用 decoupler 中实现的 AUCell 来执行这一任务。

In [24]:

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

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

In [25]:

adata

adata

Out[25]:

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上为每个细胞标注这些通路的活动水平：

In [26]:

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

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

ifn_pathways = [ "REACTOME_INTERFERON_SIGNALING", "REACTOME_INTERFERON_ALPHA_BETA_SIGNALING", "REACTOME_INTERFERON_GAMMA_SIGNALING", ] adata.obs[ifn_pathways] = adata.obsm["aucell_estimate"][ifn_pathways]

Plot the scores on the umap

In [27]:

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

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

AUCell 在评分中将那些众所周知与干扰素信号传导相关的通路在IFN刺激的细胞中得分较高，而在对照条件下的细胞通常对这些通路的得分较低，表明使用 AUCell 进行基因集评分是成功的。同时请注意，通常在 GSEA 的基因集富集测试结果中排名较高的术语得分也较大。鉴于我们事先知道与IFN相关的通路应该是排名最高的术语，AUCell 的通路活动得分与 GSEA 的基因集富集测试的一致性是令人鼓舞的。此外，这个数据集中IFN刺激的效应非常大，这也有助于这里方法的表现。

Gene set enrichment for complex experimental designs using limma-fry and pseudo-bulks¶

在簇级t检验方法中，通过将某一簇与所有其他簇进行比较来找到差异表达的基因，这在本例中包括对照和刺激细胞。线性模型允许我们仅将刺激条件下的细胞与对照组的细胞进行比较，从而更准确地识别对刺激作出反应的基因。实际上，线性模型可以适应复杂的实验设计，例如，比较治疗1下的细胞类型A与治疗2下的细胞类型A中富集的基因集；即，跨干扰和跨细胞类型效应的同时，调整批次效应、个体间变异、小鼠模型中的性别和品系差异等。

在接下来的部分，我们将演示一个针对实际数据分析例程的 limma-fry 工作流程，例如用于单细胞案例对照研究。我们首先为每种细胞类型和条件创建伪大体复制品（每种条件-细胞类型组合3个复制品）。然后我们找出在刺激细胞与对照细胞之间富集的基因集。我们还评估两种受刺激细胞类型群体之间的基因集富集情况，以发现信号通路的差异。

Create pseudo-bulk samples and explore the data¶

In [28]:

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(),
    )

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(), )

In [29]:

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

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

Out[29]:

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'

In [30]:

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

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

让我们对这些数据进行归一化处理并快速查看。由于样本量显著减少，我们这里不会使用邻居嵌入。

In [31]:

pb_data.layers["counts"] = pb_data.X.copy()

pb_data.layers["counts"] = pb_data.X.copy()

In [32]:

sc.pp.normalize_total(pb_data)
sc.pp.log1p(pb_data)
sc.pp.pca(pb_data)

sc.pp.normalize_total(pb_data) sc.pp.log1p(pb_data) sc.pp.pca(pb_data)

In [33]:

sc.pl.pca(pb_data, color=["cell_type", "condition", "lib_size"], ncols=1, size=250)

sc.pl.pca(pb_data, color=["cell_type", "condition", "lib_size"], ncols=1, size=250)

PC1现在捕捉到了淋巴细胞（T细胞、NK细胞、B细胞）和髓样细胞（单核细胞、树突状细胞）群体之间的差异，而第二主成分捕捉到了由于刺激引起的变化（即对照和受刺激伪复制品之间的差异）。理想情况下，我们感兴趣的变异应该在伪大体数据的前几个主成分中可以检测到。

在这种情况下，因为我们确实对每种细胞类型的刺激效应感兴趣，所以我们继续进行基因集测试。我们重申，绘制主成分的目的是探索数据中不同的变异轴，并发现可能对测试结果产生重大影响的不需要的变异。如果用户对其数据中的变异感到满意，可以继续进行后续分析。

为 limma 和 fry 进行设置¶

在接下来的分析部分中，我们将使用 Bioconductor 包 limma 及其方法 fry。

我们首先设置设计矩阵和对比矩阵。设计矩阵是群体成员关系的数学表示（即样本所属的群体或条件），而对比矩阵是差异测试中感兴趣的比较的数学表示。

In [34]:

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

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

In [35]:

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

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

In [36]:

%%R
colnames(design)

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

In [37]:

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

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

Index the genes annotated in each pathway in our data as follows:

In [38]:

log_norm_X = pb_data.to_df().T

log_norm_X = pb_data.to_df().T

In [39]:

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

%%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():

As done in the gsea method, let's remove gene sets with less than 15 genes

In [40]:

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

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

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

fry test for Stimulated vs Control¶

In [41]:

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

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

Taking a look at the top ranked pathways we'll see some familiar names:

In [42]:

fry_results.head()

fry_results.head()

Out[42]:

	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

In [43]:

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

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

Out[43]:

fry test for the comparison between two stimulated cell types¶

In [44]:

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

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

In [45]:

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

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

Out[45]:

如上所示，limma-fry 可以适应具有复杂实验设计的数据集和研究问题的基因集富集测试。gsea 和 fry 都提供了富集方向的见解（gsea 中的正负分数和 fry 中的方向字段）。它们都可以应用于细胞簇或伪大体样本。然而，与 gsea 不同，fry 可以进行更灵活的测试。此外，fry 可以揭示通路中的基因在实验条件之间是否发生了变化，但变化方向是一致还是不一致的。基因变化方向一致的通路在 FDR < 0.05 时被识别出来。基因在条件之间存在差异表达但变化方向不同、不一致的通路可以在 FDR > 0.05，但 FDR.Mixed < 0.05 时被识别出来（假设 0.05 是期望的显著性水平）。fry 是双向的，适用于任意设计，并且在样本数量较少的情况下也能很好地工作（尽管这在单细胞中可能不是问题）。因此，fry 的结果在生物学上可能更有意义。

关于过滤低表达基因的影响¶

如前所述，理想情况下，我们感兴趣的变异应该在伪大体数据的前几个主成分中可以检测到。 让我们删除数据中低表达的基因，应用 $\log_2$CPM 转换并重复PCA图：

In [46]:

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

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

In [47]:

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

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

In [48]:

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

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

In [49]:

pb_data.uns["logCPM_FLE"] = logCPM.T  # FLE for filter low exprs

pb_data.uns["logCPM_FLE"] = logCPM.T # FLE for filter low exprs

In [50]:

pb_data.obsm["logCPM_FLE_pca"] = sc.pp.pca(logCPM.T.to_numpy(), return_info=False)

pb_data.obsm["logCPM_FLE_pca"] = sc.pp.pca(logCPM.T.to_numpy(), return_info=False)

In [51]:

sc.pl.embedding(pb_data, "logCPM_FLE_pca", color=pb_data.obs, ncols=1, size=250)

sc.pl.embedding(pb_data, "logCPM_FLE_pca", color=pb_data.obs, ncols=1, size=250)

在这里，“logCPM_FLE”表示过滤低表达基因后进行的 $\log_2$CPM 转换。当低表达基因被移除并且通过 $\log_2$CPM 转换调整了库大小的差异后，我们现在可以清楚地观察到PC1捕捉到了细胞类型效应，而PC2捕捉到了处理效应。

在这个案例研究中，我们确实对每种细胞类型的刺激效应感兴趣，并且这种变异在基因过滤之前得到了更好的保留，因此我们展示了未过滤数据的富集测试结果。

在实践中，过滤低丰度基因和通过 edgeR::calcNormFactors 计算归一化因子是批量RNA测序分析工作流程的标准部分。如果我们对IFN刺激的全局效应感兴趣，我们应该使用过滤后的数据。此外，需要注意的是，design <- model.matrix(~ 0 + lineage + group) 会考虑髓样和淋巴系之间的差异（即基线表达差异），从而通过IFN刺激改善伪大体样本的分离，可能沿着PC1分离。在这个案例研究中，我们感兴趣的是细胞类型特异的效应，因此我们使用了一个数据模型，其中PC1的变异是由细胞类型引起的。设计矩阵的选择必须仔细考虑，以与感兴趣的生物学问题保持一致。

关于基因集之间的冗余性以及预排列基因集测试和fry基因集测试性能的说明¶

通常，密切相关的基因集之间可能有很大的重叠。这种重叠会影响基因集在富集结果中的排名，并可能影响最终的解释。例如，Kang等人的细胞被IFN-$\beta$处理。因此，预期会看到术语REACTOME_INTERFERON_ALPHA_BETA_SIGNALING作为排名最高的术语。虽然在 fry 的输出中这一术语确实是排名最高的术语，但在 GSEA 的输出中，REACTOME_INTERFERON_SIGNALING 是排名最高的术语。这个术语包含更多的基因（52个）相比于REACTOME_INTERFERON_ALPHA_BETA_SIGNALING（24个基因），并且这两个术语之间的大多数基因是共享的。这说明了预排列基因集测试（如 GSEA 和 fry）之间的另一个差异，即防止较大的基因集主导富集结果。fry 的更好表现是由于更准确的基因表达方差估计，从而获得更敏感的DE基因结果。

Key Takeaways¶

• 使用标准的scRNA-seq归一化方法对数据进行归一化，并在通路分析前过滤数据中基因覆盖率低的基因集。
• 了解不同类型的基因集测试（即竞争性 vs 自包含）并使用适合您应用的测试。
• 了解基因集富集和基因集活性推断之间的差异。GSEA 是单细胞研究中广泛使用的基因集测试；发现Pagoda 2在通路活性评分工具中表现优越。如果您的数据集具有复杂的实验设计，考虑使用在 limma 中实现的基因集测试进行伪大体分析，因为它们兼容线性模型框架，还可以考虑基因间相关性。

Quiz¶

1. What is the difference between gene set enrichment tests and activity scoring?
2. Describe examples of settings where gene set tests should be used? Can you outline examples of settings where pathway activity scoring methods are applicable?
3. What are the two types of Null Hypotheses in gene set enrichment tests. Explain the difference between the two types.
4. What is the most important preprocessing step in pathway analysis? What are the consequences if it is not conducted properly?
5. Name one gene set testing and one gene set activity scoring algorithm and explain it briefly.

Session info¶

In [52]:

%%R
sessionInfo()

%%R sessionInfo()

R version 4.1.2 (2021-11-01)
Platform: x86_64-apple-darwin13.4.0 (64-bit)
Running under: macOS Big Sur 11.6.8

Matrix products: default
LAPACK: /Users/isaac/miniconda3/envs/pathway/lib/libopenblasp-r0.3.21.dylib

locale:
[1] C/UTF-8/C/C/C/C

attached base packages:
[1] stats4    tools     stats     graphics  grDevices utils     datasets 
[8] methods   base     

other attached packages:
 [1] SingleCellExperiment_1.16.0 SummarizedExperiment_1.24.0
 [3] Biobase_2.54.0              GenomicRanges_1.46.1       
 [5] GenomeInfoDb_1.30.1         IRanges_2.28.0             
 [7] S4Vectors_0.32.4            BiocGenerics_0.40.0        
 [9] MatrixGenerics_1.6.0        matrixStats_0.63.0         

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.9             locfit_1.5-9.6         edgeR_3.36.0          
 [4] lattice_0.20-45        bitops_1.0-7           grid_4.1.2            
 [7] zlibbioc_1.40.0        XVector_0.34.0         limma_3.50.1          
[10] Matrix_1.5-3           statmod_1.4.37         RCurl_1.98-1.9        
[13] DelayedArray_0.20.0    compiler_4.1.2         GenomeInfoDbData_1.2.7

In [53]:

session_info.show()

session_info.show()

Out[53]:

Click to view session information

-----
anndata             0.8.0
anndata2ri          0.0.0
decoupler           1.3.1
numpy               1.23.5
pandas              1.5.2
rpy2                3.5.1
scanpy              1.9.1
seaborn             0.12.1
session_info        1.0.0
-----

Click to view modules imported as dependencies

PIL                         9.2.0
appnope                     0.1.3
asttokens                   NA
backcall                    0.2.0
beta_ufunc                  NA
binom_ufunc                 NA
cffi                        1.15.1
colorama                    0.4.6
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.2
debugpy                     1.6.4
decorator                   5.1.1
dunamai                     1.15.0
entrypoints                 0.4
executing                   1.2.0
get_version                 3.5.4
h5py                        3.7.0
hypergeom_ufunc             NA
importlib_metadata          NA
ipykernel                   6.17.1
jedi                        0.18.2
jinja2                      3.1.2
joblib                      1.2.0
kiwisolver                  1.4.4
llvmlite                    0.39.1
markupsafe                  2.1.1
matplotlib                  3.6.2
matplotlib_inline           0.1.6
mpl_toolkits                NA
natsort                     8.2.0
nbinom_ufunc                NA
ncf_ufunc                   NA
numba                       0.56.4
packaging                   21.3
parso                       0.8.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
platformdirs                2.5.2
prompt_toolkit              3.0.33
psutil                      5.9.4
ptyprocess                  0.7.0
pure_eval                   0.2.2
pycparser                   2.21
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.9.1
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.13.0
pynndescent                 0.5.8
pyparsing                   3.0.9
pytz                        2022.6
pytz_deprecation_shim       NA
scipy                       1.9.3
setuptools                  65.5.1
six                         1.16.0
sklearn                     1.1.3
skmisc                      0.1.4
stack_data                  0.6.2
statsmodels                 0.13.5
threadpoolctl               3.1.0
tornado                     6.2
tqdm                        4.64.1
traitlets                   5.6.0
typing_extensions           NA
tzlocal                     NA
umap                        0.5.3
wcwidth                     0.2.5
zipp                        NA
zmq                         24.0.1
zoneinfo                    NA

-----
IPython             8.7.0
jupyter_client      7.4.8
jupyter_core        5.1.0
-----
Python 3.9.15 | packaged by conda-forge | (main, Nov 22 2022, 08:55:37) [Clang 14.0.6 ]
macOS-11.6.8-x86_64-i386-64bit
-----
Session information updated at 2022-12-11 20:31