今天是生信星球陪你的第1010天
公众号里的文章大多数需要编程基础，如果因为代码看不懂，而跟不上正文的节奏，可以来找我学习，相当于给自己一个新手保护期。我的课程都是循环开课，点进去咨询微信↓
生信分析直播课程(每月初开一期，春节休一个月)
生信新手保护学习小组（每月两期）
单细胞陪伴学习小组（每月两期）
1.R包和数据准备
#pip install gseapy -i https://pypi.tuna.tsinghua.edu.cn/simple
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
随便一个h5文件即可。我这里使用的是pbmc3k，scanpy推文最后生成的文件就是它。
adata = ad.read_h5ad('../1.pbmc3k/write/pbmc3k.h5ad')
adata
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'
sc.pl.umap(adata,color='leiden',size=4,legend_loc="on data")
2.获取用于评分的基因集合
基本上大家使用的各种评分的基因集，都多数来自于gsea网站，gseapy包可以帮我们下载和读取网站上的数据，如果网络不佳可能会报错。以下代码参考自：https://gseapy.readthedocs.io/en/latest/gseapy_example.html
首先是指定自己所需要的数据是哪个版本，dbver参数是https://data.broadinstitute.org/gsea-msigdb/msigdb/release/这个网页上面的文件夹名字。而category则是该文件夹下的基因集合名字，比如人类就是h和c1~c8，都小写。
msig=Msigdb()
gmt = msig.get_gmt(category='h.all', dbver="2024.1.Hs")
列出都有哪些版本文件夹
msig.list_dbver()
列出该文件夹下都有哪些基因集合
msig.list_category(dbver="2024.1.Hs")
['c1.all',
'c2.all',
'c2.cgp',
'c2.cp.biocarta',
'c2.cp.kegg_legacy',
'c2.cp.kegg_medicus',
'c2.cp.pid',
'c2.cp.reactome',
'c2.cp',
'c2.cp.wikipathways',
'c3.all',
'c3.mir.mir_legacy',
'c3.mir.mirdb',
'c3.mir',
'c3.tft.gtrd',
'c3.tft.tft_legacy',
'c3.tft',
'c4.3ca',
'c4.all',
'c4.cgn',
'c4.cm',
'c5.all',
'c5.go.bp',
'c5.go.cc',
'c5.go.mf',
'c5.go',
'c5.hpo',
'c6.all',
'c7.all',
'c7.immunesigdb',
'c7.vax',
'c8.all',
'h.all',
'msigdb']
列出可以选择的具体基因集
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']
由上可见，其实我们只是获取到了一组基因的列表。如果你已经从别处下载或得到了基因列表，也是可以和上面的gene_set一样使用。
例如我的test.txt里放了一列基因。那么我们就可以读取并转换为python列表：
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']
3.打分并画图
敲简单了。
sc.tl.score_genes(adata,gene_set)
sc.pl.umap(adata,color='score')
WARNING: genes are not in var_names and ignored: Index(['ACADL', 'ADCY6', 'ADIG', 'ADIPOQ', 'ANGPT1', 'ANGPTL4', 'APOE',
'ATP5PO', 'CAVIN1', 'CAVIN2', 'CIDEA', 'CMBL', 'COL15A1', 'COL4A1',
'CYP4B1', 'ENPP2', 'FABP4', 'FZD4', 'GPAT4', 'HSPB8', 'ITGA7', 'ITIH5',
'LAMA4', 'LEP', 'LIFR', 'LPL', 'MIGA2', 'MRAP', 'MTARC2', 'OMD', 'ORM1',
'PPARG', 'PTGER3', 'SLC25A10', 'SLC66A3', 'SNCG', 'SOWAHC', 'SPARCL1',
'SQOR', 'UQCR11'],
dtype='object')
warning是说列表里面有些基因不在我们的表达矩阵里，很正常，无需理会。
