20. 基因集富集和通路分析#
关键要点
在做通路分析之前，使用标准的 scRNA-seq 归一化方法对数据做归一化，并过滤掉在你的数据中基因覆盖率低的基因集。
注意区分基因集富集（gene set enrichment）与基因集活性推断（gene set activity inference）。GSEA 是单细胞研究中广泛使用的基因集检验；Pagoda 2 被发现优于其他通路活性评分工具。如果你的数据集具有复杂的实验设计，可以考虑做 pseudo-bulk 分析，并使用 limma中实现的基因集检验——因为它们与线性模型框架兼容，还能额外考虑基因间的相关性。
环境设置
安装 conda:
在创建环境之前，请确保 conda 已安装在您的系统中。
保存 yml 内容:
将 yml 选项卡中的内容复制到名为
environment.yml的文件中。
创建环境:
打开终端或命令提示符。
运行以下命令：
conda env create -f environment.yml
激活环境:
创建好环境后，使用以下命令激活它：
conda activate <environment_name>
替换
<environment_name>，名称就是在environment.yml文件中指定的那个。在 yml 文件里，它看起来像这样：name: <environment_name>
验证安装:
通过运行以下命令，检查环境是否创建成功：
conda env list
name: pathway
channels:
- conda-forge
- bioconda
dependencies:
- conda-forge::python=3.12.11
- conda-forge::r-base=4.4.3
- conda-forge::pip
- conda-forge::rpy2=3.6
- conda-forge::scanpy=1.11.5
- conda-forge::scikit-misc
- bioconda::bioconductor-singlecellexperiment=1.28.0
- bioconda::bioconductor-edger=4.4.0
- bioconda::bioconductor-limma=3.62.1
- bioconda::anndata2ri=1.1
- conda-forge::r-statmod
- conda-forge::cffi
- pip:
- decoupler
20.1. 动机#
单细胞 RNA-seq 为我们提供了前所未有的洞察，去了解不同条件、组织类型、物种和个体之间细胞类型的变化。对单细胞数据做差异基因表达分析之后，几乎总会接着进行 基因集富集分析，其目的是根据差异表达（DE）基因，找出那些在某个实验条件下（相比对照或其他条件）过度代表（over-represented）的基因程序，例如生物学过程、基因本体（gene ontology）或调控通路。
为了确定在两种条件之间、以细胞类型特异方式富集的通路，首先要选定一组相关的基因集签名（gene set signature），其中每个基因集定义一个生物学过程（例如上皮-间质转化、代谢等）或一条通路（例如 MAPK 信号通路）。对于集合中的每个基因集，会用基因集中出现的 DE 基因计算一个检验统计量，再用它来评估该基因集的富集程度。视所选富集检验的类型而定，计算检验统计量时可能会、也可能不会用到基因表达的测量值。
在本章中，我们首先概述不同类型的基因集富集检验，介绍一些常用的基因签名集合，并讨论通路富集和功能富集分析的总体最佳实践。最后，我们通过演示三种基因集富集分析的分析方法来结束本章。请注意，本章中我们把通路分析、通路富集分析、基因集富集分析和功能分析这几个术语互换使用。
20.2. 通路与基因集集合#
基因集是一份经整理的基因名称（或基因 ID）列表，这些基因通过以往的研究和/或实验已知参与某个生物学过程。分子签名数据库（MSigDB） [Liberzon et al., 2011, Subramanian et al., 2005] 是最全面的数据库，由 9 个基因集集合组成。一些常用的集合是：C5，即基因本体（GO）集合；C2，来自已发表研究的、经整理的基因签名集合，通常是上下文（如组织、条件）特异的，但也包含 KEGG 和 REACTOME 基因签名。对于癌症研究，常用 Hallmark 集合；对于免疫学研究，C7 集合是常见选择。请注意，这些签名主要来自 Bulk-seq 测量，刻画的是连续表型。最近，随着 scRNA-seq 数据集的广泛普及，又出现了一些数据库，它们提供从已发表单细胞研究中整理出的标记列表，用以定义各种组织和物种中的细胞类型。这些包括 CellMarker [Zhang et al., 2019] 和 PanglaoDB [Franzén et al., 2019]。经整理的标记列表并不限于数据库中提供的那些，也可以自己整理。
20.3. 基因集富集分析中的零假设#
基因集检验可以是 competitive 或 self-contained （根据 Goeman 和 Bühlmann（2007）的定义） [Goeman and Bühlmann, 2007]。竞争性（competitive）基因集检验，检验的是：相对于不在集合中的基因，集合中的基因在差异表达方面是否排名靠前。这里的抽样单位是基因，因此可以用单个样本来做检验（即单样本 GSEA）。该检验需要不在集合中的基因（即背景基因）。而在自包含（self-contained）基因集检验中，抽样单位是受试对象，因此每组需要多个样本，但不需要有不在集合中的基因。自包含基因集检验检验的是：检验集合中的基因是否差异表达，而不考虑数据集中测量的任何其他基因。这两种零假设之间的区别，会影响对基因集富集结果的解读。注意，在生物学数据中存在基因间相关性，即同一通路中基因的表达是相关的。只有少数检验能够容纳基因间相关性，我们稍后会讨论这些方法。关于各种基因集检验的详细解释见 limma 用户手册。
20.4. 基因集检验与通路分析#
在 scRNA-seq 数据分析中，基因集富集通常是在细胞聚类或细胞类型上、一次一个地进行的。在某个聚类或细胞类型中差异表达的基因，被用来从选定的集合中识别过度代表的基因集，使用的是简单的超几何检验或 Fisher 精确检验（如 Enrichr [Chen et al., 2013]）等。这类检验不需要实际的基因表达测量值和读数来计算富集统计量，因为它们依赖于检验这样一件事有多显著： \(X\) 即与集合中非 DE 基因的数量相比，集合中差异表达的基因数量（在实验中）有多少。
fgsea [Korotkevich et al., 2021] 是更常用的基因集富集检验工具。fgsea 是对成熟的 基因集富集分析 算法在计算上更快的一种实现 [Subramanian et al., 2005]，它根据某种预排序（preranked）的基因级检验统计量来计算富集统计量。fgsea 利用基因集中基因的某种带符号统计量（例如差异表达检验得到的 t 统计量、对数倍数变化（logFC）或 p 值）来计算一个富集分数。再用若干大小相同的随机基因集，为该富集分数计算一个经验（从数据估计得到的）零分布，并计算 p 值以确定富集分数的显著性。随后对这些 p 值做多重假设检验校正。GSVA [Hänzelmann et al., 2013] 也是预排序（preranked）基因集富集方法的另一个例子。我们应注意，预排序基因集检验并不专门针对单细胞数据集，它同样适用于 Bulk-seq 测定。
在一组细胞（即聚类或同类型的细胞）中检验基因集富集的另一种方法，是从单细胞构造 pseudo-bulk 样本，并使用为 Bulk RNA-seq 开发的基因集富集方法。有几种自包含和竞争性的基因集富集检验，即 fry 和 camera 已实现于 limma [Ritchie et al., 2015]，它们通过线性模型和对检验统计量的经验贝叶斯（Empirical Bayes）压缩，与差异基因表达分析框架兼容 [Smyth, 2005]。线性模型可以通过设计矩阵容纳复杂的实验设计（如受试对象、扰动、批次、嵌套对比、交互作用等）。此外， camera 和 roast 这两个在 limma 中实现的基因集检验考虑了基因间相关性。而 limma 中的基因集检验也可以应用于（经过适当变换和归一化的）单细胞测量，而无需生成 pseudo-bulk。不过，目前还没有基准测试评估过：当这些方法直接应用于单细胞时，基因集检验结果的准确性如何。
Test |
批量或 SC |
零假设类型 |
输入 |
|---|---|---|---|
超几何 |
两者 |
竞争性 |
基因计数 |
Fisher 精确检验 |
两者 |
竞争性 |
基因计数 |
GSEA\(^*\) |
bulk |
竞争性 |
基因排名 |
GSVA\(^*\) |
bulk |
竞争性 |
基因排名 |
fgsea |
两者 |
竞争性 |
基因排名 |
fry\(^*\) |
bulk |
自包含 |
表达矩阵 |
camera\(^*\) |
bulk |
竞争性 |
表达矩阵 |
roast\(^*\) |
bulk |
自包含 |
表达矩阵 |
表：各基因集检验、适用的测量类型，以及它们所检验的零假设
\(^*\) 这些测试实际上适用于单细胞数据集，尽管其应用于单细胞可能不是常见的做法。
20.4.1. 基因集检验与通路活性推断#
基因集检验，检验的是某条通路在某种条件下相比其他条件是否富集（换句话说，是否过度代表）——例如，在单核细胞群体中，健康供体相比重症 COVID-19 患者。另一种做法是：不检验不同条件之间活性的差异，而是直接在单细胞中、以绝对意义对一条通路或一个基因签名的活性进行评分。一些广泛用于推断单细胞中基因集活性（包括通路活性）的工具包括 VISION [DeTomaso et al., 2019], AUCell [Aibar et al., 2017]、使用 Pagoda2 [Fan et al., 2016, Lake et al., 2018] 进行的通路过离散分析（pathway overdispersion analysis），以及简单的组合 Z 分数 [Lee et al., 2008]。
DoRothEA [Garcia-Alonso et al., 2019] 和 PROGENy [Schubert et al., 2018] 属于一类功能分析工具，最初是为在 Bulk RNA 数据中推断转录因子（TF）-靶标活性而开发的。Holland 等人 [Holland et al., 2020] 发现 Bulk RNA-seq 方法 DoRothEA 和 PROGENy 在模拟的 scRNA-seq 数据中具有最佳性能，甚至在单细胞数据存在漏检（drop-out）事件和文库较小的情况下，部分超过了专门为 scRNA-seq 分析设计的工具。Holland 等人还得出结论：通路和 TF 活性推断对基因集的选择比对统计方法更敏感。不过这一观察可能是功能富集分析所特有的，并且可以这样解释：TF-靶标关系是上下文特异的，即某种细胞类型中的 TF-靶标关联，实际上可能不同于另一种细胞类型或组织。
与 Holland 等人相反，Zhang 等人 [Zhang et al., 2020] 发现基于单细胞的工具（具体来说是 Pagoda2）在准确性、稳定性和可扩展性三个方面都优于基于批量的方法。应当指出，通路和基因集活性推断工具本身并不考虑批次效应、或除所关注的生物学变异之外的其他生物学变异。因此，要由数据分析者来确保差异基因表达分析这一步已经正确完成。
此外，虽然这里提到的工具会对单细胞中的每个基因集评分，但它们无法在所有被评分的基因集中，挑选出最具生物学相关性的那些。scDECAF（DavisLaboratory/scDECAF）是一个基因集活性推断工具，它允许以数据驱动的方式选出信息量最大的基因集，从而有助于解析有意义的细胞异质性。
20.5. 技术考虑因素#
20.5.1. 过滤掉基因数较少的基因集#
一个常见的做法是：在预处理步骤中，排除任何只有少数基因与数据（或高变基因 HVG）重叠的基因集。Zhang 等人 [Zhang et al., 2020] 发现，随着基因覆盖率（即通路/基因集中基因的数量）下降，基于单细胞和基于批量的方法性能都会下降。Holland 等人 [Holland et al., 2020] 也发现，规模较小的基因集会对 Bulk-seq 的性能产生不利影响 DoRothEA 和 PROGENy （应用于单细胞数据）。这些报告共同支持：在通路分析中，过滤掉基因数较少的基因集（比如集合中少于 10 或 15 个基因）是有益的。Damian 与 Gorfine（2004） [Damian and Gorfine, 2004] 把这一点归因于：基因数较少的基因集中，基因方差更可能偏大，而较大的基因集中基因方差往往较小。这会影响为检验富集而计算的检验统计量的准确性。Zhang 等人还发现，通路分析对应用于基因表达测量的归一化方法很敏感。
20.5.2. 数据归一化#
单细胞实验中的读数计数，通常在预处理流程的早期就被归一化，以确保不同文库大小的细胞之间测量具有可比性。Zhang 等人 [Zhang et al., 2020] 发现，使用 SCTransform [Hafemeister and Satija, 2019] 和 scran [Lun et al., 2016] 通常能改善基于单细胞和基于批量的通路评分工具的性能。他们发现 AUCell （一种基于排名的方法）和 z-score （变换为零均值、单位标准差）尤其会受到不同归一化方法的影响。
20.6. 案例研究：人类 PBMC 单细胞中的通路富集分析与活性水平评分#
20.6.1. 准备并探索数据#
我们首先下载了 25K PBMC 数据并遵循了标准 scanpy 工作流程，对读数计数做归一化，并在高变基因上取子集。该数据集包含未处理的和 IFN-\(\beta\) 刺激的人类 PBMC 细胞 [Kang et al., 2018]。我们用 4000 个高变基因的 UMAP 表示来探索数据中的变异模式。
from __future__ import annotations
import anndata as ad
import decoupler
import numpy as np
import pandas as pd
import scanpy as sc
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
%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,
)
我们通常建议按 Differential gene expression 一章所述来确定差异表达基因。为简单起见，这里我们运行一次 t 检验，使用 scanpy 中的 rank_genes_groups，以根据基因的差异表达检验统计量对基因排序：
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，找出与所有其他群体相比、在这个细胞群体中富集的基因集，正如 decoupler 中所实现的那样。
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
20.6.2. 用 decoupler 进行聚类级基因集富集分析#
现在我们用 Python 包 decoupler [Badia-i-Mompel et al., 2022] 对我们的数据进行 GSEA 富集检验。
20.6.2.1. 检索基因集#
下载并读取 gmt 文件，该文件对应 MSigDB 的 C2 集合中注释的 REACTOME 通路。
# 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
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
20.6.2.2. 运行 GSEA#
我们首先准备基因集。默认情况下， decoupler 不会按最大规模过滤基因集，而像 fgsea 这样的软件包则会这么做。相反，我们只是手动过滤基因集，使其至少包含 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 统计量，对 CD16 单核细胞在 IFN 刺激后的表型相关基因进行排序，并为每条通路计算 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 轴表示 \(-\log_{10}\)校正后的 p 值。因此，条形越长，该通路就越显著。通路按显著性排序。大多数与干扰素相关的通路确实排在富集程度最高的前 20 条通路之中。一些与 IFN 相关的通路包括：REACTOME_INTERFERON_SIGNALING（排名第 2）、REACTOME_INTERFERON_GAMMA_SIGNALING（排名第 3）和 REACTOME_INTERFERON_ALPHA_BETA_SIGNALING（排名第 4）。总体而言， GSEA 在识别已知与干扰素信号相关的通路方面做得相当不错，尤其考虑到我们事先就知道与 IFN 相关的通路应当是排名最靠前的条目。
我们来看看 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 分别是富集分数（enrichment score）和归一化富集分数。注意，富集分数是带符号的。因此，负分表明该通路被下调，正分则表明该通路或基因集中的基因被上调。
20.6.3. 使用 AUCell 进行细胞级通路活性评分#
与前面那种逐聚类（或者说逐细胞类型）评估基因集富集的方法不同，我们可以对每个单独细胞中通路和基因集的活性水平进行评分，它基于该细胞中基因的绝对表达，而与其他细胞中的基因表达无关。我们可以借助活性评分工具来实现这一点，例如 AUCell.
与 GSEA类似，我们将使用 decoupler 实现的 AUCell.
%%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 通路的分数加到 obs 字段（位于 AnnData 对象中），并在 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,
)
AUCell 在 IFN 刺激的细胞中，对那些已知与干扰素信号相关的通路给出了较高的分数，而对照条件下的细胞对这些通路的分数通常较低，这表明用 AUCell 进行基因集评分是成功的。还要注意，对于在 GSEA 的基因集富集检验结果中排名较高的条目，其分数通常也更大。由 AUCell得到的通路活性分数，与由 GSEA 得到的基因集富集检验结果之间的一致性令人鼓舞，尤其考虑到我们事先就知道与 IFN 相关的通路应当是排名最靠前的条目。此外，在这个数据集中 IFN 刺激的效应非常大，这也提升了这里这些方法的表现。
20.6.4. 使用 limma-fry 和 pseudo-bulk 对复杂实验设计做基因集富集#
在聚类级 t 检验方法中，差异表达基因是通过把一个聚类与所有其他聚类比较来发现的——在本例中，这既包括对照细胞也包括受刺激细胞。线性模型让我们能够只把受刺激条件下的细胞与对照组细胞比较，从而更准确地识别出对刺激作出响应的基因。事实上，线性模型可以容纳复杂的实验设计，例如，找出在 Cell type A in treatment 1 中相比 Cell type A in treatment 2中更为富集的基因集；也就是说，跨扰动、跨细胞类型的效应，同时校正批次效应、个体间差异、小鼠模型中的性别和品系差异等。
在下一节中，我们演示一个可推广到现实数据分析流程（例如单细胞病例-对照研究）的 limma-fry 工作流程。我们首先为每种细胞类型和条件创建 pseudo-bulk 重复（每个“条件-细胞类型”组合 3 个重复）。然后，我们找出在某种细胞类型中、受刺激相比对照富集的基因集。我们还会评估两个受刺激的细胞类型群体之间的基因集富集，以发现信号通路上的差异。
20.6.4.1. 创建伪块样本并探索数据#
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)
我们来把这份数据归一化，并快速看一下。由于样本量大幅减少，这里我们不会使用近邻嵌入。
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 捕捉的是因施加刺激而产生的变异（即对照与受刺激伪重复之间的差异）。理想情况下，所关注的变异应当能在 pseudo-bulk 数据的前几个 PC 中被检测到。
在这种情况下，由于我们确实关心每种细胞类型的刺激效应，我们就接着做基因集检验。我们再次强调，绘制 PC 的目的是探索数据中的各个变异轴，并发现可能严重影响检验结果的、不需要的变异。如果用户对自己数据中的变异情况满意，就可以继续进行后续分析。
20.6.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() 来检验我们上面设定的每个对比中富集的通路：
20.6.4.3. 用 fry 检验比较受刺激 vs 对照#
%%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())
)
20.6.4.4. 用 fry 检验比较两种受刺激的细胞类型#
%%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())
)
如上所述，limma-fry 可以为具有复杂实验设计的数据集和研究问题提供基因集富集检验。两者都 gsea 和 fry 都能提供关于富集方向的信息（正分或负分见于 gsea ，而 Direction 字段见于 fry）。它们都可以应用于细胞聚类或 pseudo-bulk 样本。然而，与 gsea不同，可以用 fry进行更灵活的检验。此外， fry 能揭示某条通路中的基因是否在实验条件之间发生变化、且变化方向是否一致。基因朝一致方向变化的通路，可以通过 FDR < 0.05 来识别；而基因在条件之间差异表达、但朝不同的、不一致方向变化的通路，则可以在 FDR > 0.05、但 FDR.Mixed < 0.05 处识别（假设 0.05 为期望的显著性水平）。fry 是双向的，适用于任意设计，并且在样本数较少时也表现良好（尽管这在单细胞中可能不是问题）。因此，由以下方法得到的结果 fry 也许在生物学上更有意义。
20.6.4.4.1. 关于过滤低表达基因的效果#
如前所述，理想情况下，所关注的变异应当能在 pseudo-bulk 数据的前几个 PC 中被检测到。我们来去掉数据中表达量低的基因，应用 \(\log_2\)CPM 变换，然后重复绘制 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”表示先对低表达基因进行过滤，再做 \(\log_2\)CPM 变换。现在我们可以清楚地看到：当去掉低表达基因、并对文库大小差异加以校正后，PC1 捕捉的是细胞类型效应，PC2 捕捉的是处理效应 \(\log_2\)CPM 变换。
由于在本案例研究中我们确实关心每种细胞类型的刺激效应，而这种变异在基因过滤之前保留得更好，因此我们展示的是未过滤数据上的富集检验结果。在实践中，过滤低丰度基因、并通过以下方法计算归一化因子 edgeR::calcNormFactors 是批量 RNA-seq 分析流程的标准环节。如果我们关心的是 IFN 刺激的全局效应，那就应该使用过滤后的数据。此外，可以注意到 design <- model.matrix(~ 0 + lineage + group) 会考虑髓系（myeloid）与淋巴系（lymphoid）谱系之间的差异（即基线表达差异），从而改善 pseudo-bulk 样本按 IFN 刺激的分离，这可能体现在 PC1 上。在本案例研究中，我们关心的是细胞类型特异的效应，因此我们坚持采用这样一个数据模型：PC1 方向上的变异是按细胞类型区分的。设计矩阵的选择必须仔细斟酌，以便与所关注的生物学问题保持一致。
20.6.4.4.2. 关于基因集之间的冗余、以及 preranked 与 fry 基因集检验性能的说明#
一般来说，关系密切的基因集之间可能有很大的重叠。这种重叠会影响基因集在富集结果中的排名，并可能损害最终的解读。例如，Kang 等人研究中的细胞用 IFN-\(\beta\)处理。因此，人们会预期看到 REACTOME_INTERFERON_ALPHA_BETA_SIGNALING 这个条目排名最高。虽然在 fry的输出中这个条目确实排名第一，但在 GSEA 的输出中，REACTOME_INTERFERON_SIGNALING 才是排名第一的条目。这个条目包含的基因数（52）比 REACTOME_INTERFERON_ALPHA_BETA_SIGNALING（24 个基因）更多，而且这些基因大多在两个条目之间共享。这说明了 preranked 基因集检验（例如 GSEA 和 fry）的另一个差异：在防止较大的基因集主导富集结果方面。以下方法之所以表现更好 fry，是因为它对基因表达方差的估计更准确，因而 DE 基因结果更灵敏。
20.7. Quiz#
20.8. 会话信息#
%%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
session_info.show()
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
20.9. 参考文献#
Sara Aibar, Carmen Bravo González-Blas, Thomas Moerman, Vân Anh Huynh-Thu, Hana Imrichova, Gert Hulselmans, Florian Rambow, Jean-Christophe Marine, Pierre Geurts, Jan Aerts, and others. Scenic: single-cell regulatory network inference and clustering. Nature methods, 14(11):1083–1086, 2017.
Pau Badia-i-Mompel, Jesús Vélez Santiago, Jana Braunger, Celina Geiss, Daniel Dimitrov, Sophia Müller-Dott, Petr Taus, Aurelien Dugourd, Christian H Holland, Ricardo O Ramirez Flores, and others. Decoupler: ensemble of computational methods to infer biological activities from omics data. Bioinformatics Advances, 2(1):vbac016, 2022.
Edward Y Chen, Christopher M Tan, Yan Kou, Qiaonan Duan, Zichen Wang, Gabriela Vaz Meirelles, Neil R Clark, and Avi Ma’ayan. Enrichr: interactive and collaborative html5 gene list enrichment analysis tool. BMC bioinformatics, 14(1):1–14, 2013.
Doris Damian and Malka Gorfine. Statistical concerns about the gsea procedure. Nature genetics, 36(7):663–663, 2004.
David DeTomaso, Matthew G Jones, Meena Subramaniam, Tal Ashuach, Chun J Ye, and Nir Yosef. Functional interpretation of single cell similarity maps. Nature communications, 10(1):1–11, 2019.
Jean Fan, Neeraj Salathia, Rui Liu, Gwendolyn E Kaeser, Yun C Yung, Joseph L Herman, Fiona Kaper, Jian-Bing Fan, Kun Zhang, Jerold Chun, and others. Characterizing transcriptional heterogeneity through pathway and gene set overdispersion analysis. Nature methods, 13(3):241–244, 2016.
Oscar Franzén, Li-Ming Gan, and Johan LM Björkegren. Panglaodb: a web server for exploration of mouse and human single-cell rna sequencing data. Database, 2019.
Luz Garcia-Alonso, Christian H Holland, Mahmoud M Ibrahim, Denes Turei, and Julio Saez-Rodriguez. Benchmark and integration of resources for the estimation of human transcription factor activities. Genome research, 29(8):1363–1375, 2019.
Jelle J Goeman and Peter Bühlmann. Analyzing gene expression data in terms of gene sets: methodological issues. Bioinformatics, 23(8):980–987, 2007.
Christoph Hafemeister and Rahul Satija. Normalization and variance stabilization of single-cell rna-seq data using regularized negative binomial regression. Genome biology, 20(1):1–15, 2019.
Christian H Holland, Jovan Tanevski, Javier Perales-Patón, Jan Gleixner, Manu P Kumar, Elisabetta Mereu, Brian A Joughin, Oliver Stegle, Douglas A Lauffenburger, Holger Heyn, and others. Robustness and applicability of transcription factor and pathway analysis tools on single-cell rna-seq data. Genome biology, 21(1):1–19, 2020.
Sonja Hänzelmann, Robert Castelo, and Justin Guinney. Gsva: gene set variation analysis for microarray and rna-seq data. BMC bioinformatics, 14(1):1–15, 2013.
Hyun Min Kang, Meena Subramaniam, Sasha Targ, Michelle Nguyen, Lenka Maliskova, Elizabeth McCarthy, Eunice Wan, Simon Wong, Lauren Byrnes, Cristina M Lanata, and others. Multiplexed droplet single-cell rna-sequencing using natural genetic variation. Nature biotechnology, 36(1):89–94, 2018.
Gennady Korotkevich, Vladimir Sukhov, Nikolay Budin, Boris Shpak, Maxim N Artyomov, and Alexey Sergushichev. Fast gene set enrichment analysis. BioRxiv, pages 060012, 2021.
Blue B Lake, Song Chen, Brandon C Sos, Jean Fan, Gwendolyn E Kaeser, Yun C Yung, Thu E Duong, Derek Gao, Jerold Chun, Peter V Kharchenko, and others. Integrative single-cell analysis of transcriptional and epigenetic states in the human adult brain. Nature biotechnology, 36(1):70–80, 2018.
Eunjung Lee, Han-Yu Chuang, Jong-Won Kim, Trey Ideker, and Doheon Lee. Inferring pathway activity toward precise disease classification. PLoS computational biology, 4(11):e1000217, 2008.
Arthur Liberzon, Aravind Subramanian, Reid Pinchback, Helga Thorvaldsdóttir, Pablo Tamayo, and Jill P Mesirov. Molecular signatures database (msigdb) 3.0. Bioinformatics, 27(12):1739–1740, 2011.
Aaron TL Lun, Davis J McCarthy, and John C Marioni. A step-by-step workflow for low-level analysis of single-cell rna-seq data with bioconductor. F1000Research, 2016.
Matthew E Ritchie, Belinda Phipson, DI Wu, Yifang Hu, Charity W Law, Wei Shi, and Gordon K Smyth. Limma powers differential expression analyses for rna-sequencing and microarray studies. Nucleic acids research, 43(7):e47–e47, 2015.
Michael Schubert, Bertram Klinger, Martina Klünemann, Anja Sieber, Florian Uhlitz, Sascha Sauer, Mathew J Garnett, Nils Blüthgen, and Julio Saez-Rodriguez. Perturbation-response genes reveal signaling footprints in cancer gene expression. Nature communications, 9(1):1–11, 2018.
Gordon K Smyth. Limma: linear models for microarray data. In Bioinformatics and computational biology solutions using R and Bioconductor, pages 397–420. Springer, 2005.
Aravind Subramanian, Pablo Tamayo, Vamsi K Mootha, Sayan Mukherjee, Benjamin L Ebert, Michael A Gillette, Amanda Paulovich, Scott L Pomeroy, Todd R Golub, Eric S Lander, and others. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proceedings of the National Academy of Sciences, 102(43):15545–15550, 2005.
Xinxin Zhang, Yujia Lan, Jinyuan Xu, Fei Quan, Erjie Zhao, Chunyu Deng, Tao Luo, Liwen Xu, Gaoming Liao, Min Yan, and others. Cellmarker: a manually curated resource of cell markers in human and mouse. Nucleic acids research, 47(D1):D721–D728, 2019.
20.10. 贡献者#
我们衷心感谢以下人员的贡献：