单细胞RNA-seq生信分析全流程——第十四篇：差异基因表达分析
14.1 前言
本篇是细胞注释章节的更详细的延续，该章节已经介绍了差异基因表达（DGE）作为用细胞类型注释簇的工具。在这里，我们重点关注更复杂的实验设计中差异基因表达测试的更高级用例，这些实验设计涉及一种或多种条件，例如疾病、基因敲除或药物。在这种情况下，我们通常对感兴趣的条件和参考条件之间基因表达模式差异的大小和重要性感兴趣。该参考可以是一切，但通常是健康样本。这种统计测试可以应用于任意组，但在单细胞RNA-Seq的情况下通常应用于细胞类型水平。
差异基因表达分析试图推断任何比较组之间（通常在每种细胞类型的健康组和状况组之间）统计上显著过度表达或表达不足的基因。
这种分析的结果可能是影响并可能解释任何观察到的表型的基因组。然后可以更仔细地检查这些基因组，例如受影响的途径或诱导的细胞间通讯变化。
差异基因表达测试通常会返回每个比较条件下每个比较基因的log2倍数变化和调整后的p值。然后可以按p值对该列表进行排序并进行更详细的研究。
流行的学生t检验是进行此类检验的一种方法。然而，它没有考虑到一些单细胞RNA-seq的特殊性，例如来自dropout的过多零或需要复杂的实验设计。更具体地说，在不汇集跨基因信息的情况下，很少有人有足够的样本数量来准确估计方差。此外，原始计数从来都不是给定样本中特定基因表达的绝对测量值。每个基因的实际读取数取决于文库制备的效率、非编码转录本的污染量和测序深度。因此，它确实缺乏单细胞RNA-seq的敏感性和特异性，更不用说实验设计灵活性了。
因此，差异基因表达测试是一个经典的生物信息学问题，已经被许多工具解决。一般来说，目前从两个角度来解决这个问题，即样本级视图，其中表达被聚合以创建“伪批量”，然后使用最初为批量表达样本设计的方法进行分析，例如edgeR或DEseq2以及细胞水平视图，其中使用广义混合效应模型（例如 MAST[Finak et al., 2015]或glmmTMB[Brooks et al., 2017]）对细胞进行单独建模。DGE工具的跨数据集的共识和稳健性较低。如前所述，尽管单细胞数据包含技术噪声伪影，例如丢失、零膨胀和高细胞间变异性， 与专门为scRNA-seq数据设计的方法相比，为批量RNA-seq数据设计的方法表现良好。发现单细胞特异性方法特别容易将高表达基因错误地标记为差异表达。
在本篇中，我们演示如何使用两种工具进行差异表达分析：具有拟似然检验的edgeR和具有随机效应的MAST。由于edgeR和MAST都是在R中实现的，因此我们使用anndata2ri包能够同时使用Python中的AnnData对象和R中的SingeCellExperiment对象。
对于这两种方法中的每一种，我们展示了两个用例：如何对完整数据进行分析以及如何对多种细胞类型和一种特定细胞类型进行测试。
14.2 环境配置
import warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
import seaborn as sns
import scanpy as sc
import pandas as pd
import numpy as np
import random
import sc_toolbox
import pertpy
import rpy2.rinterface_lib.callbacks
import anndata2ri
import logging
from rpy2.robjects import pandas2ri
from rpy2.robjects import r
sc.settings.verbosity = 0
rpy2.rinterface_lib.callbacks.logger.setLevel(logging.ERROR)
pandas2ri.activate()
anndata2ri.activate()
%load_ext rpy2.ipython
%%R
library(edgeR)
library(MAST)
14.3 数据集准备
我们将使用Kang数据集，它是来自8名狼疮患者INF-β治疗前后6小时（总共16个样本）的10倍基于scRNA-seq的外周血单核细胞 (PBMC) [Kang et al., 2018]。 干扰素β以天然成纤维细胞或重组制剂（干扰素β-1a和干扰素β-1b）的形式使用，并发挥与干扰素α类似的抗病毒和抗增殖特性。干扰素β已被批准用于治疗复发缓解型多发性硬化症和继发进展型多发性硬化症。
首先，我们加载完整的数据集。
adata = pertpy.data.kang_2018()
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'
我们将需要.obs的label（其中包含条件标签）、replicate和cell_type列。
adata.obs[:5]
nCount_RNA nFeature_RNA tsne1 tsne2 label cluster cell_type replicate nCount_SCT nFeature_SCT integrated_snn_res.0.4 seurat_clusters
index
AAACATACATTTCC-1 3017.0 877 -27.640373 14.966629 ctrl 9 CD14+ Monocytes patient_1016 1704.0 711 1 1
AAACATACCAGAAA-1 2481.0 713 -27.493646 28.924885 ctrl 9 CD14+ Monocytes patient_1256 1614.0 662 1 1
AAACATACCATGCA-1 703.0 337 -10.468194 -5.984389 ctrl 3 CD4 T cells patient_1488 908.0 337 6 6
AAACATACCTCGCT-1 3420.0 850 -24.367997 20.429285 ctrl 9 CD14+ Monocytes patient_1256 1738.0 653 1 1
AAACATACCTGGTA-1 3158.0 1111 27.952170 24.159738 ctrl 4 Dendritic cells patient_1039 1857.0 928 12 12
我们需要使用原始计数，因此我们检查.X确实包含原始计数并将它们放入AnnData对象的counts层中。
np.max(adata.X)
3828.0
adata.layers["counts"] = adata.X.copy()
我们有8名对照患者和8名疾病患者。
print(len(adata[adata.obs["label"] == "ctrl"].obs["replicate"].cat.categories))
print(len(adata[adata.obs["label"] == "stim"].obs["replicate"].cat.categories))
8
8
我们过滤含有少于200个基因的细胞以及在少于3个细胞中发现的基因，以进行基本的质量控制。
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
adata
AnnData object with n_obs × n_vars = 24562 × 15701
obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'n_genes'
var: 'name', 'n_cells'
obsm: 'X_pca', 'X_umap'
layers: 'counts'
14.4 伪批量Pseudobulk
edgeR是在R中实现的差异基因表达测试工具，最初是为批量基因表达数据而设计的。它实现了基于负二项分布、经验贝叶斯估计、精确检验、广义线性模型 (GLM) 和拟似然检验的广泛统计方法。
在这里，我们将使用拟似然检验，因为它考虑了离散度估计的不确定性。相反，精确检验假设估计的离差是真实值，这可能会导致一些不准确。此外，拟似然GLM在实验设计方面更加灵活。
由于edgeR是作为批量数据的DE分析方法引入的，因此我们首先需要从单细胞数据集中创建伪批量样本。对于每位患者，我们通过聚集每个亚群的细胞并获取该亚群内的平均基因表达来为每种细胞类型创建1个伪批量样本。
如果您正在使用的数据没有重复，则为每个患者创建多个（例如2-3个）pseudobulks以考虑患者变异性可能会有所帮助。在这里，我们选择为每个患者创建1个伪批量，因为对于每个患者，我们有2个重复，一个对照，一个刺激。
无论我们是只想对少数细胞亚群进行分析并分别为每个细胞亚群拟合一个模型，还是为所有细胞亚群拟合一个模型，我们首先需要准备数据，定义一个函数来创建pseudobulks并运行edgeR管道。首先，我们准备数据。
由于我们需要为每个患者状况组合创建伪批量，因此我们首先需要通过连接replicate和label来创建这样的列。
adata.obs["sample"] = [
f"{rep}_{l}" for rep, l in zip(adata.obs["replicate"], adata.obs["label"])
]
我们需要清理细胞类型名称，即用下划线替换空格并删除+符号，以避免Python到R的转换问题。
adata.obs["cell_type"] = [ct.replace(" ", "_") for ct in adata.obs["cell_type"]]
adata.obs["cell_type"] = [ct.replace("+", "") for ct in adata.obs["cell_type"]]
我们需要将分类元数据设置为真正分类的才能创建伪批量。
adata.obs["replicate"] = adata.obs["replicate"].astype("category")
adata.obs["label"] = adata.obs["label"].astype("category")
adata.obs["sample"] = adata.obs["sample"].astype("category")
adata.obs["cell_type"] = adata.obs["cell_type"].astype("category")
现在，让我们定义将单个细胞聚合成伪重复所需的函数：
-
aggregate_and_filter是一个函数，它从原始单细胞AnnData对象中创建一个AnnData对象，该对象为指定子群提供一个伪复制。在这里，我们还筛选出特定群体的细胞数量少于30个的供体。 - 通过更改
replicates_per_patent参数，可以为每个样本创建几个（n）伪重复；然后细胞被分成n个大小大致相等的子集。
NUM_OF_CELL_PER_DONOR = 30
def aggregate_and_filter(
adata,
cell_identity,
donor_key="sample",
condition_key="label",
cell_identity_key="cell_type",
obs_to_keep=[], # which additional metadata to keep, e.g. gender, age, etc.
replicates_per_patient=1,
):
# subset adata to the given cell identity
adata_cell_pop = adata[adata.obs[cell_identity_key] == cell_identity].copy()
# check which donors to keep according to the number of cells specified with NUM_OF_CELL_PER_DONOR
size_by_donor = adata_cell_pop.obs.groupby([donor_key]).size()
donors_to_drop = [
donor
for donor in size_by_donor.index
if size_by_donor[donor] <= NUM_OF_CELL_PER_DONOR
]
if len(donors_to_drop) > 0:
print("Dropping the following samples:")
print(donors_to_drop)
df = pd.DataFrame(columns=[*adata_cell_pop.var_names, *obs_to_keep])
adata_cell_pop.obs[donor_key] = adata_cell_pop.obs[donor_key].astype("category")
for i, donor in enumerate(donors := adata_cell_pop.obs[donor_key].cat.categories):
print(f"\tProcessing donor {i+1} out of {len(donors)}...", end="\r")
if donor not in donors_to_drop:
adata_donor = adata_cell_pop[adata_cell_pop.obs[donor_key] == donor]
# create replicates for each donor
indices = list(adata_donor.obs_names)
random.shuffle(indices)
indices = np.array_split(np.array(indices), replicates_per_patient)
for i, rep_idx in enumerate(indices):
adata_replicate = adata_donor[rep_idx]
# specify how to aggregate: sum gene expression for each gene for each donor and also keep the condition information
agg_dict = {gene: "sum" for gene in adata_replicate.var_names}
for obs in obs_to_keep:
agg_dict[obs] = "first"
# create a df with all genes, donor and condition info
df_donor = pd.DataFrame(adata_replicate.X.A)
df_donor.index = adata_replicate.obs_names
df_donor.columns = adata_replicate.var_names
df_donor = df_donor.join(adata_replicate.obs[obs_to_keep])
# aggregate
df_donor = df_donor.groupby(donor_key).agg(agg_dict)
df_donor[donor_key] = donor
df.loc[f"donor_{donor}_{i}"] = df_donor.loc[donor]
print("\n")
# create AnnData object from the df
adata_cell_pop = sc.AnnData(
df[adata_cell_pop.var_names], obs=df.drop(columns=adata_cell_pop.var_names)
)
return adata_cell_pop
我们还需要定义一个单独的函数来拟合edgeR GLM：
-
fit_model将SingleCellExperiment对象作为输入，创建设计矩阵并输出拟合的GLM。我们还输出DGEList类的edgeR对象来进行一些探索性数据分析(EDA)。
%%R
fit_model <- function(adata_){
# create an edgeR object with counts and grouping factor
y <- DGEList(assay(adata_, "X"), group = colData(adata_)$label)
# filter out genes with low counts
print("Dimensions before subsetting:")
print(dim(y))
print("")
keep <- filterByExpr(y)
y <- y[keep, , keep.lib.sizes=FALSE]
print("Dimensions after subsetting:")
print(dim(y))
print("")
# normalize
y <- calcNormFactors(y)
# create a vector that is concatentation of condition and cell type that we will later use with contrasts
group <- paste0(colData(adata_)$label, ".", colData(adata_)$cell_type)
replicate <- colData(adata_)$replicate
# create a design matrix: here we have multiple donors so also consider that in the design matrix
design <- model.matrix(~ 0 + group + replicate)
# estimate dispersion
y <- estimateDisp(y, design = design)
# fit the model
fit <- glmQLFit(y, design)
return(list("fit"=fit, "design"=design, "y"=y))
}
现在我们定义了所需的所有函数，因此我们可以继续创建伪批量。我们可能希望稍后查看可用的元数据，因此将其保留在AnnData对象中。
obs_to_keep = ["label", "cell_type", "replicate", "sample"]
我们需要将原始计数传递给edgeR。因此，我们将.X设置为counts层，以确保为原始计数创建伪重复。
adata.X = adata.layers["counts"].copy()
接下来，我们使用pseudobulks创建AnnData对象。
# process first cell type separately...
cell_type = adata.obs["cell_type"].cat.categories[0]
print(
f'Processing {cell_type} (1 out of {len(adata.obs["cell_type"].cat.categories)})...'
)
adata_pb = aggregate_and_filter(adata, cell_type, obs_to_keep=obs_to_keep)
for i, cell_type in enumerate(adata.obs["cell_type"].cat.categories[1:]):
print(
f'Processing {cell_type} ({i+2} out of {len(adata.obs["cell_type"].cat.categories)})...'
)
adata_cell_type = aggregate_and_filter(adata, cell_type, obs_to_keep=obs_to_keep)
adata_pb = adata_pb.concatenate(adata_cell_type)
Processing B_cells (1 out of 8)...
Dropping the following samples:
['patient_1039_ctrl']
Processing donor 16 out of 16...
Processing CD14_Monocytes (2 out of 8)...
Processing donor 16 out of 16...
Processing CD4_T_cells (3 out of 8)...
Processing donor 16 out of 16...
Processing CD8_T_cells (4 out of 8)...
Dropping the following samples:
['patient_101_ctrl', 'patient_1039_ctrl', 'patient_1039_stim', 'patient_107_ctrl', 'patient_107_stim', 'patient_1244_ctrl', 'patient_1244_stim']
Processing donor 16 out of 16...
Processing Dendritic_cells (5 out of 8)...
Dropping the following samples:
['patient_1016_ctrl', 'patient_1016_stim', 'patient_101_ctrl', 'patient_1039_ctrl', 'patient_1039_stim', 'patient_107_ctrl', 'patient_107_stim']
Processing donor 16 out of 16...
Processing FCGR3A_Monocytes (6 out of 8)...
Dropping the following samples:
['patient_1039_ctrl', 'patient_1039_stim', 'patient_107_ctrl', 'patient_107_stim', 'patient_1244_stim']
Processing donor 16 out of 16...
Processing Megakaryocytes (7 out of 8)...
Dropping the following samples:
['patient_1015_ctrl', 'patient_1015_stim', 'patient_1016_ctrl', 'patient_101_ctrl', 'patient_1039_stim', 'patient_107_stim', 'patient_1244_ctrl', 'patient_1244_stim', 'patient_1256_ctrl', 'patient_1256_stim', 'patient_1488_ctrl']
Processing donor 11 out of 11...
Processing NK_cells (8 out of 8)...
Dropping the following samples:
['patient_1039_ctrl', 'patient_1039_stim']
Processing donor 16 out of 16...
差异基因表达结果的有效性很大程度上取决于统计模型中变异主轴的捕获。中间数据探索步骤，例如对伪大量样本的主成分分析（PCA）或多维标度（MDS），可以识别变异源，从而可以指导构建数据建模的相应设计和对比矩阵。
如果未能考虑包括生物重复在内的实验的生物变异性的多种来源，将会导致FDR膨胀。虽然增加每个个体的细胞数量可以提高精度，但它对检测个体差异的能力影响有限。因此，提高统计功效的最佳方法是增加独立实验样本的数量。
由于我们的数据已经生成，我们无法进一步增加独立实验样本的数量。尽管如此，我们现在将探索我们的数据以确定变化的主轴，以正确生成我们的设计矩阵。
我们对创建的伪重复执行非常基本的EDA，以检查某些患者/伪批量是否是我们需要排除的异常值，以免DE结果出现偏差。我们将原始计数保存在counts层中，然后对计数进行归一化并计算归一化伪批量计数的PCA坐标。
adata_pb.layers['counts'] = adata_pb.X.copy()
sc.pp.normalize_total(adata_pb, target_sum=1e6)
sc.pp.log1p(adata_pb)
sc.pp.pca(adata_pb)
接下来，我们通过所有可用的元数据查看在PCA图上创建的伪重复和颜色，以查看是否存在我们可能希望包含在设计矩阵中的任何混杂因素。我们还添加了lib_size和log_lib_size列来检查库大小和PC组件之间是否存在相关性。
adata_pb.obs["lib_size"] = np.sum(adata_pb.layers["counts"], axis=1)
adata_pb.obs["log_lib_size"] = np.log(adata_pb.obs["lib_size"])
sc.pl.pca(adata_pb, color=adata_pb.obs, ncols=1, size=300)
我们观察PCA图上细胞类型的分离以及刺激细胞和未刺激细胞的分离。其他协变量（批次）似乎与PCA分量没有明显相关，因此我们不将它们中的任何一个包含到我们的设计矩阵中。
如上所述，edgeR将原始计数作为输入，因此我们在继续之前将计数放回.X字段。
adata_pb.X = adata_pb.layers['counts'].copy()
14.4.1 单组分析
首先，我们展示如何准备数据、构建设计矩阵以及对一种特定细胞类型进行DE测试。
我们在CD14+单核细胞数据子集上运行管道，正如论文中所示，在该子群中未确定最多数量的DE基因。
adata_mono = adata_pb[adata_pb.obs["cell_type"] == "CD14_Monocytes"]
adata_mono
View of AnnData object with n_obs × n_vars = 16 × 15701
obs: 'label', 'cell_type', 'replicate', 'sample', 'batch', 'lib_size', 'log_lib_size'
uns: 'log1p', 'pca', 'label_colors', 'cell_type_colors', 'replicate_colors', 'sample_colors', 'batch_colors'
obsm: 'X_pca'
varm: 'PCs'
layers: 'counts'
清理样本名称以使绘图不那么拥挤。
adata_mono.obs_names = [
name.split("_")[2] + "_" + name.split("_")[3] for name in adata_mono.obs_names
]
%%time
%%R -i adata_mono
outs <-fit_model(adata_mono)
[1] "Dimensions before subsetting:"
[1] 15701 16
[1] ""
[1] "Dimensions after subsetting:"
[1] 3709 16
[1] ""
CPU times: user 2.95 s, sys: 224 ms, total: 3.18 s
Wall time: 4.27 s
%%R
fit <- outs$fit
y <- outs$y
由于我们在进入分析时并未事先假设特定基因将被上调或下调，因此我们需要可视化来理解DGE结果。MDS图可以提供高层次的概述。通常，我们期望不同条件下的样品之间存在分离，如下图所示我们的结果。
%%R
plotMDS(y, col=ifelse(y$samples$group == "stim", "red", "blue"))
生物变异系数 (BCV) 图显示了生物群体中每个基因的平均变异性，作为平均表达的函数。例如，bcv为0.3表示各组间基因表达平均存在30%的变异性。
低丰度基因通常表现出较大的BCV，因为低丰度基因的读数计数测量更不确定。另一方面，平均表达量高的基因的定量更可靠，因此它们通常具有较低的变异性，因此 BCV 也较低。使用此图可以检测异常基因或查明可能需要在设计矩阵中反映的其他实验因素。例如，出现在图右上角的一组具有高平均表达和高BCV的基因可以标记实验压力、污染等，特别是如果它们属于相似的基因家族。
在下面的 BCV 图中，我们看到一些具有高bcv的低丰度基因，但没有看到具有高bcv的高丰度基因，这表明不应在设计矩阵中建模任何进一步的实验考虑因素。
%%R
plotBCV(y)
接下来，我们进行拟似然检验以找到对照条件和刺激条件之间的DE基因。让我们检查一下我们的设计矩阵必须指定哪些列来指定要测试的正确列。
%%R
colnames(y$design)
[1] "groupctrl.CD14_Monocytes" "groupstim.CD14_Monocytes"
[3] "replicatepatient_107" "replicatepatient_1015"
[5] "replicatepatient_1016" "replicatepatient_1039"
[7] "replicatepatient_1244" "replicatepatient_1256"
[9] "replicatepatient_1488"
%%R -o tt
myContrast <- makeContrasts('groupstim.CD14_Monocytes-groupctrl.CD14_Monocytes', levels = y$design)
qlf <- glmQLFTest(fit, contrast=myContrast)
# get all of the DE genes and calculate Benjamini-Hochberg adjusted FDR
tt <- topTags(qlf, n = Inf)
tt <- tt$table
让我们检查一下这个表：对于每个没有被edgeR过滤掉的基因（在我们的例子中是 3709），该表包含了DE测试的结果。该表可以保存为.csv文件以便稍后使用并用于可视化。在对完整数据进行分析后，我们将在本节末尾展示如何使用火山图。
tt.shape
(3709, 5)
tt[:5]
logFC logCPM F PValue FDR
HESX1 8.345536 6.773420 1281.013295 1.837373e-15 2.766927e-12
CD38 7.126846 7.420668 1243.793133 2.266164e-15 2.766927e-12
NT5C3A 5.657050 8.327003 1218.102628 2.628780e-15 2.766927e-12
SOCS1 4.388247 6.943768 1191.289806 3.079524e-15 2.766927e-12
GMPR 6.943484 7.031832 1159.601183 3.730018e-15 2.766927e-12
顺便说一句，人们还可以使用glmTreat测试相对于指定倍数变化所提供的系数或对比度差异表达的基因。
%%R
tr <- glmTreat(fit, contrast=myContrast, lfc=1.5)
print(head(topTags(tr)))
Coefficient: -1*groupctrl.CD14_Monocytes 1*groupstim.CD14_Monocytes
logFC unshrunk.logFC logCPM PValue FDR
HESX1 8.345536 8.702486 6.773420 2.108737e-14 3.535749e-11
CD38 7.126846 7.219184 7.420668 2.292459e-14 3.535749e-11
NT5C3A 5.657050 5.674640 8.327003 3.149544e-14 3.535749e-11
IL1RN 6.588583 6.596787 10.358946 4.249159e-14 3.535749e-11
GMPR 6.943484 7.047504 7.031832 4.766446e-14 3.535749e-11
DEFB1 6.654201 6.696759 8.067159 7.670147e-14 4.741429e-11
最后，我们可以看到有多少基因的FDR校正值小于0.01。
涂片图（也称为MA、M值与A值图）显示基因的对数倍数变化与其平均丰度的函数关系。我们通常在低丰度范围内观察到较高的logFC，因为读数计数在低丰度范围内变化更大，导致logFC估计值较大。如果我们将loess曲线拟合到logFC和平均logCPM值，则趋势应以零为中心。任何偏离此值的情况都可能表明数据尚未正确标准化。平均表达量大且绝对值 logFC 大的基因可以标记生物学上感兴趣的基因以供研究和随访。
%%R
plotSmear(qlf, de.tags = rownames(tt)[which(tt$FDR<0.01)])
14.4.2 多组分析
接下来，我们展示如何准备数据、构建设计矩阵并使用所有可用细胞类型的完整数据集的对比来执行DE测试。
%%time
%%R -i adata_pb
outs <-fit_model(adata_pb)
[1] "Dimensions before subsetting:"
[1] 15701 90
[1] ""
[1] "Dimensions after subsetting:"
[1] 2358 90
[1] ""
CPU times: user 11.9 s, sys: 39.7 ms, total: 12 s
Wall time: 12 s
%%R
fit <- outs$fit
y <- outs$y
现在我们使用对比对每种细胞类型进行拟似然检验。由于没有直接的方法将表从R循环内的R移动到Python，因此我们手动获取每种细胞类型的结果。
%%R -i adata_pb -o de_per_cell_type
de_per_cell_type <- list()
for (cell_type in unique(colData(adata_pb)$cell_type)) {
print(cell_type)
# create contrast for this cell type
myContrast <- makeContrasts(paste0("groupstim.", cell_type, "-groupctrl.", cell_type), levels = y$design)
# perform QLF test
qlf <- glmQLFTest(fit, contrast=myContrast)
# get all of the DE genes and calculate Benjamini-Hochberg adjusted FDR
tt <- topTags(qlf, n = Inf)
# save in the list with the results for all the cell types
de_per_cell_type[[cell_type]] <- tt$table
}
[1] "B_cells"
[1] "CD14_Monocytes"
[1] "CD4_T_cells"
[1] "CD8_T_cells"
[1] "Dendritic_cells"
[1] "FCGR3A_Monocytes"
[1] "NK_cells"
现在让我们使用sc_toolbox包 (https://github.com/schillerlab/sc-toolbox) 中的 sc_toolbox.tools.de_res_to_anndata将其保存在原始adata对象的.uns中，这会保存结果，就像使用scanpy创建的rank_genes_groups()函数一样。像这样存储DE表的优点是我们现在还可以使用标准的scanpy绘图函数。我们将每种细胞类型的DEG表保存为.csv，因为我们稍后在细胞间通讯章节中将需要它们。
# get cell types that we ran the analysis for
cell_types = de_per_cell_type.keys()
# add the table to .uns for each cell type
for cell_type in cell_types:
df = de_per_cell_type[cell_type]
df["gene_symbol"] = df.index
df["cell_type"] = cell_type
sc_toolbox.tools.de_res_to_anndata(
adata,
df,
groupby="cell_type",
score_col="logCPM",
pval_col="PValue",
pval_adj_col="FDR",
lfc_col="logFC",
key_added="edgeR_" + cell_type,
)
df.to_csv(f"de_edgeR_{cell_type}.csv")
为了再次将表作为pandasDataFrame获取，我们使用sc.get.rank_genes_groups_df函数。
sc.get.rank_genes_groups_df(adata, group="CD14_Monocytes", key="edgeR_CD14_Monocytes")[
:5
]
names scores logfoldchanges pvals pvals_adj
0 FTH1 16.186098 -0.572272 0.001589 0.002712
1 MALAT1 16.025682 0.010049 0.877078 0.897245
2 B2M 15.468524 0.677266 0.0 0.0
3 TMSB4X 15.054697 -0.217131 0.004663 0.007331
4 FTL 14.984937 -0.041832 0.669951 0.710956
14.4.3 关于edgeR的tips
- 需要原始计数作为输入
- 需要来自单细胞实验的pseudobulks伪重复
- 如果单细胞实验中有多个供体，并且用户想要考虑患者的变异性，我们建议为每个患者创建2或3个伪重复，并将患者信息纳入设计矩阵中
14.5 单细胞特异算法
MAST使用两部分广义线性模型对单细胞基因表达进行建模。一部分模拟细胞中每个基因的离散表达率，而另一部分模拟条件连续表达水平。
MAST框架使用两部分广义线性模型对单细胞基因表达进行建模。MAST的一个组件模拟细胞中每个基因的离散表达率，而另一个组件则模拟条件连续表达水平（以所表达的基因为条件）。 欲了解更多详细信息，请阅读原始出版物[Finak et al., 2015]。
MAST将标准化计数作为输入，因此我们首先采用“counts”层，然后执行标准化步骤。
adata.X = adata.layers["counts"].copy()
sc.pp.normalize_total(adata, target_sum=1e6)
sc.pp.log1p(adata)
我们定义了一个小辅助函数来处理R和Python之间的一些对象类型转换问题。我们还需要过滤掉每个亚群的少量细胞（本例中为3个）中表达的基因，因为模型需要能够估计每个基因的方差。
def prep_anndata(adata_):
def fix_dtypes(adata_):
df = pd.DataFrame(adata_.X.A, index=adata_.obs_names, columns=adata_.var_names)
df = df.join(adata_.obs)
return sc.AnnData(df[adata_.var_names], obs=df.drop(columns=adata_.var_names))
adata_ = fix_dtypes(adata_)
sc.pp.filter_genes(adata_, min_cells=3)
return adata_
14.5.1 单组分析
与edgeR一样，可以使用完整数据集拟合模型，然后使用对比对每种感兴趣的细胞类型执行测试，但在本篇中，我们仅显示一种细胞类型（即CD14单核细胞）的MAST-RE管道，缩短运行时间。
adata_mono = adata[adata.obs["cell_type"] == "CD14_Monocytes"].copy()
adata_mono
AnnData object with n_obs × n_vars = 5696 × 15701
obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'n_genes', 'sample'
var: 'name', 'n_cells'
uns: 'edgeR_B_cells', 'edgeR_CD14_Monocytes', 'edgeR_CD4_T_cells', 'edgeR_CD8_T_cells', 'edgeR_Dendritic_cells', 'edgeR_FCGR3A_Monocytes', 'edgeR_NK_cells', 'log1p'
obsm: 'X_pca', 'X_umap'
layers: 'counts'
sc.pp.filter_genes(adata_mono, min_cells=3)
adata_mono
AnnData object with n_obs × n_vars = 5696 × 12268
obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'n_genes', 'sample'
var: 'name', 'n_cells'
uns: 'edgeR_B_cells', 'edgeR_CD14_Monocytes', 'edgeR_CD4_T_cells', 'edgeR_CD8_T_cells', 'edgeR_Dendritic_cells', 'edgeR_FCGR3A_Monocytes', 'edgeR_NK_cells', 'log1p'
obsm: 'X_pca', 'X_umap'
layers: 'counts'
接下来我们如上所述过滤这两个对象。
adata_mono = prep_anndata(adata_mono)
adata_mono
AnnData object with n_obs × n_vars = 5696 × 12268
obs: 'nCount_RNA', 'nFeature_RNA', 'tsne1', 'tsne2', 'label', 'cluster', 'cell_type', 'replicate', 'nCount_SCT', 'nFeature_SCT', 'integrated_snn_res.0.4', 'seurat_clusters', 'n_genes', 'sample'
var: 'n_cells'
与执行多个统计测试时的情况一样，必须使用Benjamini-Hochberg校正等方法对条件下DGE测试获得的p值进行校正。
与edgeR分析类似，我们定义一个用于分析的单独函数：
-
find_de_MAST_RE采用SingleCellExperiment对象作为输入，并使用RE管道运行 MAST。该函数的输出是表（Python中的pandas DataFrame），其中包含分析结果，例如每个基因的对数倍数变化、p值和FDR校正值。
adata_mono.obs["cell_type"] = [
ct.replace(" ", "_") for ct in adata_mono.obs["cell_type"]
]
adata_mono.obs["cell_type"] = [
ct.replace("+", "") for ct in adata_mono.obs["cell_type"]
]
%%R
find_de_MAST_RE <- function(adata_){
# create a MAST object
sca <- SceToSingleCellAssay(adata_, class = "SingleCellAssay")
print("Dimensions before subsetting:")
print(dim(sca))
print("")
# keep genes that are expressed in more than 10% of all cells
sca <- sca[freq(sca)>0.1,]
print("Dimensions after subsetting:")
print(dim(sca))
print("")
# add a column to the data which contains scaled number of genes that are expressed in each cell
cdr2 <- colSums(assay(sca)>0)
colData(sca)$ngeneson <- scale(cdr2)
# store the columns that we are interested in as factors
label <- factor(colData(sca)$label)
# set the reference level
label <- relevel(label,"ctrl")
colData(sca)$label <- label
celltype <- factor(colData(sca)$cell_type)
colData(sca)$celltype <- celltype
# same for donors (which we need to model random effects)
replicate <- factor(colData(sca)$replicate)
colData(sca)$replicate <- replicate
# create a group per condition-celltype combination
colData(sca)$group <- paste0(colData(adata_)$label, ".", colData(adata_)$cell_type)
colData(sca)$group <- factor(colData(sca)$group)
# define and fit the model
zlmCond <- zlm(formula = ~ngeneson + group + (1 | replicate),
sca=sca,
method='glmer',
ebayes=F,
strictConvergence=F,
fitArgsD=list(nAGQ = 0)) # to speed up calculations
# perform likelihood-ratio test for the condition that we are interested in
summaryCond <- summary(zlmCond, doLRT='groupstim.CD14_Monocytes')
# get the table with log-fold changes and p-values
summaryDt <- summaryCond$datatable
result <- merge(summaryDt[contrast=='groupstim.CD14_Monocytes' & component=='H',.(primerid, `Pr(>Chisq)`)], # p-values
summaryDt[contrast=='groupstim.CD14_Monocytes' & component=='logFC', .(primerid, coef)],
by='primerid') # logFC coefficients
# MAST uses natural logarithm so we convert the coefficients to log2 base to be comparable to edgeR
result[,coef:=result[,coef]/log(2)]
# do multiple testing correction
result[,FDR:=p.adjust(`Pr(>Chisq)`, 'fdr')]
result = result[result$FDR<0.01,, drop=F]
result <- stats::na.omit(as.data.frame(result))
return(result)
}
我们运行单核细胞的分析流程。
%%time
%%R -i adata_mono -o res
res <-find_de_MAST_RE(adata_mono)
[1] "Dimensions before subsetting:"
[1] 12268 5696
[1] ""
[1] "Dimensions after subsetting:"
[1] 1676 5696
[1] ""
CPU times: user 9min 35s, sys: 2.96 s, total: 9min 38s
Wall time: 9min 43s
看一下运行结果。
res[:5]
primerid Pr(>Chisq) coef FDR
1 AAED1 8.410341e-19 0.709836 1.597106e-18
2 ABI1 1.401688e-05 0.312797 1.824921e-05
3 ABRACL 3.920183e-06 0.418028 5.201004e-06
4 ACADVL 2.121110e-26 -0.751305 4.656979e-26
5 ACOT9 3.424174e-151 2.921133 2.689503e-150
我们像以前一样将结果存储在.uns中。请注意，我们不需要score列，因此我们只需将logFC传递到那里即可。
res["gene_symbol"] = res["primerid"]
res["cell_type"] = "CD14_Monocytes"
sc_toolbox.tools.de_res_to_anndata(
adata,
res,
groupby="cell_type",
score_col="coef",
pval_col="Pr(>Chisq)",
pval_adj_col="FDR",
lfc_col="coef",
key_added="MAST_CD14_Monocytes",
)
adata_copy = adata.copy()
14.6 可视化
我们可以使用热图和火山图来可视化结果。热图的每一行对应一个基因，每一列对应一个单细胞。颜色越亮，该基因在特定细胞中的表达就越高。由于我们只绘制DE基因，因此我们希望看到两种条件之间表达的明显差异。火山图通常用于可视化统计测试的结果，它们在x轴上显示表达变化（对数倍数变化），在y轴上显示统计显着性（FDR校正的p值）。我们对FDR校正p值低于0.01且对数倍变化超过1.5的基因进行颜色编码。
我们在绘制热图之前对数据进行归一化，以便更好地查看两种条件之间表达的差异。
adata.X = adata.layers["counts"].copy()
sc.pp.normalize_total(adata, target_sum=1e6)
sc.pp.log1p(adata)
接下来，我们为热图定义一个辅助绘图函数。
FDR = 0.01
LOG_FOLD_CHANGE = 1.5
def plot_heatmap(adata, group_key, group_name="cell_type", groupby="label"):
cell_type = "_".join(group_key.split("_")[1:])
res = sc.get.rank_genes_groups_df(adata, group=cell_type, key=group_key)
res.index = res["names"].values
res = res[
(res["pvals_adj"] < FDR) & (abs(res["logfoldchanges"]) > LOG_FOLD_CHANGE)
].sort_values(by=["logfoldchanges"])
print(f"Plotting {len(res)} genes...")
markers = list(res.index)
sc.pl.heatmap(
adata[adata.obs[group_name] == cell_type].copy(),
markers,
groupby=groupby,
swap_axes=True,
)
最后，我们可以使用RE绘制edgeR和MAST的CD14+单核细胞的DEG热图。
plot_heatmap(adata, "edgeR_CD14_Monocytes")
Plotting 303 genes...
plot_heatmap(adata, "MAST_CD14_Monocytes")
Plotting 436 genes...
我们观察到，MAST根据我们给定的调整p值和对数倍变化的截止值识别出436个差异基因，而edgeR识别出303个差异基因。
接下来，我们定义火山图的辅助绘图函数。
FDR = 0.01
LOG_FOLD_CHANGE = 1.5
def volcano_plot(adata, group_key, group_name="cell_type", groupby="label", title=None):
cell_type = "_".join(group_key.split("_")[1:])
result = sc.get.rank_genes_groups_df(adata, group=cell_type, key=group_key).copy()
result["-logQ"] = -np.log(result["pvals"].astype("float"))
lowqval_de = result.loc[abs(result["logfoldchanges"]) > LOG_FOLD_CHANGE]
other_de = result.loc[abs(result["logfoldchanges"]) <= LOG_FOLD_CHANGE]
fig, ax = plt.subplots()
sns.regplot(
x=other_de["logfoldchanges"],
y=other_de["-logQ"],
fit_reg=False,
scatter_kws={"s": 6},
)
sns.regplot(
x=lowqval_de["logfoldchanges"],
y=lowqval_de["-logQ"],
fit_reg=False,
scatter_kws={"s": 6},
)
ax.set_xlabel("log2 FC")
ax.set_ylabel("-log Q-value")
if title is None:
title = group_key.replace("_", " ")
plt.title(title)
plt.show()
volcano_plot(adata, "MAST_CD14_Monocytes")
volcano_plot(adata, "edgeR_CD14_Monocytes")
从热图，尤其是火山图中，我们可以看到，edgeR识别出上调基因多于下调基因（刺激与对照），而MAST识别出相似数量的上调和下调基因。
