引言
本系列讲解 使用Scanpy
分析单细胞(scRNA-seq)数据 教程,持续更新,欢迎关注,转发!
重新评估质控
正如之前提到的,现在将通过 UMAP
可视化不同的质量控制指标,重新审视的细胞过滤策略。
sc.pl.umap(
adata,
color=["leiden", "predicted_doublet", "doublet_score"],
# increase horizontal space between panels
wspace=0.5,
size=3,
)
AI 代码解读
sc.pl.umap(
adata,
color=["leiden", "log1p_total_counts", "pct_counts_mt", "log1p_n_genes_by_counts"],
wspace=0.5,
ncols=2,
)
AI 代码解读
手动细胞类型注释
细胞类型注释是一项繁琐且重复的工作,通常需要多次进行亚聚类和重新注释。
目前已经筛选出了一批质量较好的细胞,接下来可以对它们进行已知细胞类型的注释。通常,会使用特定细胞类型独有的基因(即标记基因)来进行注释,这些基因能够帮助区分数据中不同类型的细胞群体。此前的研究已经将各种标记基因整理到了一些资源库中,例如 CellMarker、TF-Marker 和 PanglaoDB。此外,Cellxgene 基因表达工具也非常实用,可以查看某个基因在多个数据集中被哪些细胞类型表达。
通常情况下,细胞类型注释是在细胞聚类完成后,通过标记基因来进行的。因此,需要生成一系列聚类结果,然后利用这些结果来完成细胞类型的注释。在这里,采用 Leiden 聚类算法,从最近邻图中提取细胞群落。
for res in [0.02, 0.5, 2.0]:
sc.tl.leiden(
adata, key_added=f"leiden_res_{res:4.2f}", resolution=res, flavor="igraph"
)
AI 代码解读
需要注意的是,设定的聚类数量并没有固定的标准,同样,用来调整聚类数量的分辨率参数也是主观选择的。因此,最终的聚类数量还是要根据能够区分出的稳定且具有生物学意义的细胞群体来确定,这通常需要相关领域的专家来进行,或者参考专家整理的标记基因等先验知识。
sc.pl.umap(
adata,
color=["leiden_res_0.02", "leiden_res_0.50", "leiden_res_2.00"],
legend_loc="on data",
)
AI 代码解读
sc.pl.umap(
adata,
color=["leiden_res_0.02", "leiden_res_0.50", "leiden_res_2.00"],
legend_loc="on data",
)
AI 代码解读
虽然 UMAP 的结果不能过度解读,但可以看到,在最高分辨率下,数据被过度细分成了太多聚类;而在最低分辨率下,可能会把本应属于不同细胞类型的细胞归为一类。
标记基因集
来定义一组标记基因,用于标记在该数据集中预期会看到的主要细胞类型。
marker_genes = {
"CD14+ Mono": ["FCN1", "CD14"],
"CD16+ Mono": ["TCF7L2", "FCGR3A", "LYN"],
# Note: DMXL2 should be negative
"cDC2": ["CST3", "COTL1", "LYZ", "DMXL2", "CLEC10A", "FCER1A"],
"Erythroblast": ["MKI67", "HBA1", "HBB"],
# Note HBM and GYPA are negative markers
"Proerythroblast": ["CDK6", "SYNGR1", "HBM", "GYPA"],
"NK": ["GNLY", "NKG7", "CD247", "FCER1G", "TYROBP", "KLRG1", "FCGR3A"],
"ILC": ["ID2", "PLCG2", "GNLY", "SYNE1"],
"Naive CD20+ B": ["MS4A1", "IL4R", "IGHD", "FCRL1", "IGHM"],
# Note IGHD and IGHM are negative markers
"B cells": [
"MS4A1",
"ITGB1",
"COL4A4",
"PRDM1",
"IRF4",
"PAX5",
"BCL11A",
"BLK",
"IGHD",
"IGHM",
],
"Plasma cells": ["MZB1", "HSP90B1", "FNDC3B", "PRDM1", "IGKC", "JCHAIN"],
# Note PAX5 is a negative marker
"Plasmablast": ["XBP1", "PRDM1", "PAX5"],
"CD4+ T": ["CD4", "IL7R", "TRBC2"],
"CD8+ T": ["CD8A", "CD8B", "GZMK", "GZMA", "CCL5", "GZMB", "GZMH", "GZMA"],
"T naive": ["LEF1", "CCR7", "TCF7"],
"pDC": ["GZMB", "IL3RA", "COBLL1", "TCF4"],
}
sc.pl.dotplot(adata, marker_genes, groupby="leiden_res_0.02", standard_scale="var")
AI 代码解读
从这些标记基因的表达模式来看,可以清晰地看到一些规律,利用这些规律可以为最粗略的聚类标注大致的细胞谱系。
adata.obs["cell_type_lvl1"] = adata.obs["leiden_res_0.02"].map(
{
"0": "Lymphocytes",
"1": "Monocytes",
"2": "Erythroid",
"3": "B Cells",
}
)
sc.pl.dotplot(adata, marker_genes, groupby="leiden_res_0.50", standard_scale="var")
AI 代码解读
这种分辨率似乎能够较好地区分数据中的大多数细胞类型。因此,可以尝试结合聚类的 UMAP 图和上面的点图,手动对这些细胞类型进行注释。如果需要,还可以进一步查看每个聚类,并对其进行亚聚类分析。
差异表达基因作为标记
此外,还可以计算每个聚类的标记基因,然后查看这些基因是否与已知的生物学现象(如细胞类型或状态)相关。这通常是通过简单的统计检验(如 Wilcoxon 检验或 t 检验)来实现的,即比较每个聚类与其他聚类的差异。
# Obtain cluster-specific differentially expressed genes
sc.tl.rank_genes_groups(adata, groupby="leiden_res_0.50", method="wilcoxon")
AI 代码解读
然后,可以在点图上展示每个聚类的前 5 个差异表达基因。
sc.pl.rank_genes_groups_dotplot(
adata, groupby="leiden_res_0.50", standard_scale="var", n_genes=5
)
# WARNING: dendrogram data not found (using key=dendrogram_leiden_res_0.50). Running `sc.tl.dendrogram` with default parameters. For fine tuning it is recommended to run `sc.tl.dendrogram` independently.
AI 代码解读
通过这些基因,可以推断出正在观察的细胞类型。例如,聚类 7 表达了 NKG7 和 GNLY,这表明它们可能是自然杀伤(NK)细胞。
如果你想要自己创建图表,或者采用更自动化的分析方法,可以使用 scanpy.get.rank_genes_groups_df()
以方便的格式提取差异表达基因。
sc.get.rank_genes_groups_df(adata, group="7").head(5)
AI 代码解读
dc_cluster_genes = sc.get.rank_genes_groups_df(adata, group="7").head(5)["names"]
sc.pl.umap(
adata,
color=[*dc_cluster_genes, "leiden_res_0.50"],
legend_loc="on data",
frameon=False,
ncols=3,
)
AI 代码解读
你可能已经注意到,这里的 p 值非常低。这是因为统计检验将每个细胞都视为一个独立的样本。如果希望采用更保守的方法,可以对数据进行“伪批量”处理(例如,sc.get.aggregate(adata, by=["sample", "cell_type"], func="sum", layer="counts")),然后使用更强大的差异表达分析工具,如 pydeseq2
。