python版本单细胞数据基因集打分并可视化
从我们的新课《掌握Python，解锁单细胞数据的无限可能》学习了python版本的对单细胞数据用某个基因集打分，现在就来实践一下~
单细胞数据的不同基因集打分有非常广泛的应用，比如文献《Single-cell dissection of cellular and molecular features underlying mesenchymal stem cell therapy in ischemic acute kidney injury》展示了6个基因集在不同细胞亚群中的打分小提琴图：
基因集变异分析（GSVA）进一步揭示，促纤维化TEC（肾小管上皮细胞）和受损TEC在炎症和纤维化相关通路中表现出富集，包括TNF-α信号通路、缺氧、TGF-β信号通路、IL-6-JAK-STAT3通路和纤维化，这与其促纤维化和促炎症特性一致（图2E）
example02-GSVA
0.环境准备
这一次需要用到python版本的gsea软件
代码语言：javascript代码运行次数：0运行复制# bash安装
pip install gseapy lxml scipy==1.11.4 fa2-modified
加载模块：
代码语言：javascript代码运行次数：0运行复制# conda activate sc
# pip install gseapy lxml scipy==1.11.4 fa2-modified
from gseapy import Msigdb
import pandas as pd
import numpy as np
import scanpy as sc
import matplotlib.pyplot as plt
import anndata as ad
1.获取基因集
基因集可以来自很多地方，数据库，文献等，就是一群有生物学意义的基因的集合。常见的功能数据库有GO/KEGG/msigdb。
本次就用 msigdb 数据库的，对应的python代码如下：
代码语言：javascript代码运行次数：0运行复制## 基因集：使用来自 Msigdb 数据库的
# gp为前面加载的gseapy模块的缩写 import gseapy as gp
msig = gp.Msigdb()
## `get_gmt()`函数有两个参数：
# dbver="2024.1.Hs": list msigdb version you wanna query
msig.list_dbver()
# category='h.all': list categories given dbver. 列出该数据库版本下都有哪些基因集合
msig.list_category(dbver="2024.1.Hs")
# 获取 2024.1.Hs 版本下的 h.all，即那50个非常有名的癌症相关基因集，超多文献里面都用它做一个打分来讲各种生物学故事
gmt = msig.get_gmt(category='h.all', dbver="2024.1.Hs")
msig.get_gmt 函数有两个参数：
- dbver：指定 Msigdb数据库的版本，可以使用函数
msig.list_dbver()列出来，也可以再网址这里查看得到 / - category：指定某一个版本下的基因集合类型，可以使用函数
msig.list_category(dbver="2024.1.Hs")列出来，也可以再网址这里查看得到 /
以前也介绍过：python版本的功能富集分析：GSEApy
对于获取到的gmt基因集简单探索
上面下载了 2024.1.Hs 版本中 50个 hallmark 通路，探索一下：
代码语言：javascript代码运行次数：0运行复制# 列出具体的通路名
list(gmt.keys())[0:10] #只列了前10
['HALLMARK_ADIPOGENESIS',
'HALLMARK_ALLOGRAFT_REJECTION',
'HALLMARK_ANDROGEN_RESPONSE',
'HALLMARK_ANGIOGENESIS',
'HALLMARK_APICAL_JUNCTION',
'HALLMARK_APICAL_SURFACE',
'HALLMARK_APOPTOSIS',
'HALLMARK_BILE_ACID_METABOLISM',
'HALLMARK_CHOLESTEROL_HOMEOSTASIS',
'HALLMARK_COAGULATION']
# 列出某个基因集里的基因
gene_set=gmt['HALLMARK_ADIPOGENESIS']
print(gene_set)
['ABCA1', 'ABCB8', 'ACAA2', 'ACADL', 'ACADM', 'ACADS', 'ACLY', 'ACO2', 'ACOX1', 'ADCY6', 'ADIG', 'ADIPOQ', 'ADIPOR2', 'AGPAT3', 'AIFM1', 'AK2', 'ALDH2', 'ALDOA', 'ANGPT1', 'ANGPTL4', 'APLP2', 'APOE', 'ARAF', 'ARL4A', 'ATL2', 'ATP1B3', 'ATP5PO', 'BAZ2A', 'BCKDHA', 'BCL2L13', 'BCL6', 'C3', 'CAT', 'CAVIN1', 'CAVIN2', 'CCNG2', 'CD151', 'CD302', 'CD36', 'CDKN2C', 'CHCHD10', 'CHUK', 'CIDEA', 'CMBL', 'CMPK1', 'COL15A1', 'COL4A1', 'COQ3', 'COQ5', 'COQ9', 'COX6A1', 'COX7B', 'COX8A', 'CPT2', 'CRAT', 'CS', 'CYC1', 'CYP4B1', 'DBT', 'DDT', 'DECR1', 'DGAT1', 'DHCR7', 'DHRS7', 'DHRS7B', 'DLAT', 'DLD', 'DNAJB9', 'DNAJC15', 'DRAM2', 'ECH1', 'ECHS1', 'ELMOD3', 'ELOVL6', 'ENPP2', 'EPHX2', 'ESRRA', 'ESYT1', 'ETFB', 'FABP4', 'FAH', 'FZD4', 'G3BP2', 'GADD45A', 'GBE1', 'GHITM', 'GPAM', 'GPAT4', 'GPD2', 'GPHN', 'GPX3', 'GPX4', 'GRPEL1', 'HADH', 'HIBCH', 'HSPB8', 'IDH1', 'IDH3A', 'IDH3G', 'IFNGR1', 'IMMT', 'ITGA7', 'ITIH5', 'ITSN1', 'JAGN1', 'LAMA4', 'LEP', 'LIFR', 'LIPE', 'LPCAT3', 'LPL', 'LTC4S', 'MAP4K3', 'MCCC1', 'MDH2', 'ME1', 'MGLL', 'MGST3', 'MIGA2', 'MRAP', 'MRPL15', 'MTARC2', 'MTCH2', 'MYLK', 'NABP1', 'NDUFA5', 'NDUFAB1', 'NDUFB7', 'NDUFS3', 'NKIRAS1', 'NMT1', 'OMD', 'ORM1', 'PDCD4', 'PEMT', 'PEX14', 'PFKFB3', 'PFKL', 'PGM1', 'PHLDB1', 'PHYH', 'PIM3', 'PLIN2', 'POR', 'PPARG', 'PPM1B', 'PPP1R15B', 'PRDX3', 'PREB', 'PTCD3', 'PTGER3', 'QDPR', 'RAB34', 'REEP5', 'REEP6', 'RETN', 'RETSAT', 'RIOK3', 'RMDN3', 'RNF11', 'RREB1', 'RTN3', 'SAMM50', 'SCARB1', 'SCP2', 'SDHB', 'SDHC', 'SLC19A1', 'SLC1A5', 'SLC25A1', 'SLC25A10', 'SLC27A1', 'SLC5A6', 'SLC66A3', 'SNCG', 'SOD1', 'SORBS1', 'SOWAHC', 'SPARCL1', 'SQOR', 'SSPN', 'STAT5A', 'STOM', 'SUCLG1', 'SULT1A1', 'TALDO1', 'TANK', 'TKT', 'TOB1', 'TST', 'UBC', 'UBQLN1', 'UCK1', 'UCP2', 'UQCR10', 'UQCR11', 'UQCRC1', 'UQCRQ', 'VEGFB', 'YWHAG']
也可以是自己定义的基因集：
代码语言：javascript代码运行次数：0运行复制gene_set2 = pd.read_table('test.txt',header=None)[0].tolist()
print(gene_set2)
['CD3D', 'CD3E', 'CD3G', 'CD247', 'CD4', 'CD8A', 'CD8B', 'CD8B2', 'PTPRC', 'LCK', 'FYN', 'ZAP70', 'LCP2', 'LAT', 'ITK', 'TEC', 'NCK1', 'NCK2', 'VAV3', 'VAV1', 'VAV2', 'GRAP2', 'GRB2', 'PAK1', 'PAK2', 'PAK3', 'PAK4', 'PAK5', 'PAK6', 'BUB1B-PAK6', 'RHOA', 'CDC42', 'DLG1', 'MAPK11', 'MAPK12', 'MAPK13', 'MAPK14', 'PLCG1', 'PPP3CA', 'PPP3CB', 'PPP3CC', 'PPP3R1', 'PPP3R2', 'NFATC1', 'NFATC2', 'NFATC3', 'SOS1', 'SOS2', 'RASGRP1', 'HRAS', 'KRAS', 'NRAS', 'RAF1', 'MAP2K1', 'MAP2K2', 'MAPK1', 'MAPK3', 'FOS', 'JUN', 'PRKCQ', 'CARD11', 'BCL10', 'MALT1', 'MAP3K7', 'MAP2K7', 'MAPK8', 'MAPK10', 'MAPK9', 'CHUK', 'IKBKB', 'IKBKG', 'NFKB1', 'RELA', 'NFKBIA', 'NFKBIB', 'NFKBIE', 'CD28', 'ICOS', 'CD40LG', 'PIK3R1', 'PIK3R2', 'PIK3R3', 'P3R3URF-PIK3R3', 'PIK3CA', 'PIK3CD', 'PIK3CB', 'PDPK1', 'AKT1', 'AKT2', 'AKT3', 'MAP3K8', 'MAP3K14', 'GSK3B', 'PDCD1', 'CTLA4', 'PTPN6', 'PTPN11', 'PPP2CA', 'PPP2CB', 'PPP2R1B', 'PPP2R1A', 'PPP2R2A', 'PPP2R2B', 'PPP2R2C', 'PPP2R2D', 'PPP2R3B', 'PPP2R3C', 'PPP2R3A', 'PPP2R5B', 'PPP2R5C', 'PPP2R5D', 'PPP2R5E', 'PPP2R5A', 'CBLB', 'IL2', 'IL4', 'IL5', 'IL10', 'IFNG', 'CSF2', 'TNF', 'CDK4']
2.打分计算
这里需要注意的是使用的adata对象，保证数据中有全部的基因而不是只有高变基因，不然基因集中很多基因可能就没有办法计算到打分中：
示例数据这里就用经典的pbmc3k数据集：
代码语言：javascript代码运行次数：0运行复制adata = ad.read_h5ad('pbmc3k_anno.h5ad')
adata
这里共有 13714个基因：
代码语言：javascript代码运行次数：0运行复制AnnData object with n_obs × n_vars = 2638 × 13714
obs: 'n_genes', 'n_genes_by_counts', 'total_counts', 'total_counts_mt', 'pct_counts_mt', 'leiden'
var: 'gene_ids', 'n_cells', 'mt', 'n_cells_by_counts', 'mean_counts', 'pct_dropout_by_counts', 'total_counts', 'highly_variable', 'means', 'dispersions', 'dispersions_norm', 'mean', 'std'
uns: 'dendrogram_leiden', 'hvg', 'leiden', 'leiden_colors', 'log1p', 'neighbors', 'pca', 'rank_genes_groups', 'umap'
obsm: 'X_pca', 'X_umap'
varm: 'PCs'
obsp: 'connectivities', 'distances'
数据中的 leiden 为已经做好注释的标签：
代码语言：javascript代码运行次数：0运行复制import matplotlib.pyplot as plt
fig, ax = plt.subplots(figsize=(5, 5))
sc.pl.umap(adata,color='leiden',size=12,legend_loc="on data", ax=ax)
# 显示图表
plt.show()
打分函数sc.tl.score_genes
代码语言：javascript代码运行次数：0运行复制# 列出某个基因集里的基因
gene_set1 = gmt['HALLMARK_ADIPOGENESIS']
sc.tl.score_genes(adata,gene_set1,score_name="HALLMARK_ADIPOGENESIS")
现在看一下数据，会发现多了一列：
代码语言：javascript代码运行次数：0运行复制adata
具体看一下：
代码语言：javascript代码运行次数：0运行复制adata.obs.head()
可以看基因集与数据中基因的交集个数：
代码语言：javascript代码运行次数：0运行复制import numpy as np
# 列出某个基因集里的基因
gene_set1 = gmt['HALLMARK_ADIPOGENESIS']
list1 = gene_set1
list1[1:10]
# 单细胞中的基因
list2 = adata.var.index.tolist()
len(list2)
# 取交集
inter = np.intersect1d(list1, list2).tolist()
print(inter)
# 交集个数
len(inter)
# 160个
如果你的adata之前保存只用了高变基因用于下游分析，那这里的基因交集个数会大大减少，需要注意~
多个基因集打分
上面只做了50个基因集里面的一个，如果我想对50个通路全部计算打分呢？来看看：
使用一个for循环gmt，gmt是一个字典
代码语言：javascript代码运行次数：0运行复制for index, key in enumerate(gmt):
sc.tl.score_genes(adata,gmt[key],score_name=key)
3.打分可视化
可以在ump上可视化打分：
代码语言：javascript代码运行次数：0运行复制# 单个通路
sc.pl.umap(adata,color=['HALLMARK_ADIPOGENESIS'])
小提琴可视化：
adata：AnnData对象。keys：指定要绘制的列名，这里是"score"。groupby：如果需要按某个分类变量（例如细胞类型）分组绘制小提琴图，可以指定一个列名。如果不需要分组，可以设置为None。size：小提琴图的大小。log：是否对数据取对数。cut：小提琴图的截断范围。
# 绘制小提琴图
sc.pl.violin(adata, keys="HALLMARK_APOPTOSIS", groupby="leiden", jitter=0.2,size=2, log=False, cut=0,rotation=90)
还可以可视化多个通路：
ump多个：
代码语言：javascript代码运行次数：0运行复制sc.pl.umap(adata,color=list(gmt.keys())[0:10])
小提琴多个：
代码语言：javascript代码运行次数：0运行复制# 绘制小提琴图
sc.pl.violin(adata, keys=['HALLMARK_ADIPOGENESIS','HALLMARK_APOPTOSIS'], groupby="leiden", jitter=0.1, size=2, log=False, cut=0,rotation=90)
堆叠小提琴：
代码语言：javascript代码运行次数：0运行复制sc.pl.stacked_violin(adata, list(gmt.keys())[0:10], groupby="leiden", swap_axes=False,
figsize=(7, 5)) # 调整画布大小为 10x6 英寸
气泡图：
代码语言：javascript代码运行次数：0运行复制sc.pl.dotplot(adata, list(gmt.keys())[0:10], groupby='leiden', dendrogram=False,
figsize=(7, 5) )
热图：
代码语言：javascript代码运行次数：0运行复制sc.pl.matrixplot(adata, var_names=list(gmt.keys())[0:10], groupby='leiden', cmap='RdBu_r')
热图：
代码语言：javascript代码运行次数：0运行复制sc.pl.heatmap(adata=adata, var_names=list(gmt.keys())[0:10], groupby='leiden', cmap='RdBu_r', dendrogram=True)
发布评论