python单细胞数据的基因集打分

用户11414625

发布于 2024-12-20 16:37:00

13900

代码可运行

文章被收录于专栏：生信星球520生信星球520

运行总次数：0

代码可运行

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')