单细胞注释记不住marker怎么办--让AI帮你解释差异基因

众所周知，单细胞数据分析最磨人的环节就是确定不同的单细胞亚群的生物学名字，而且每个领域都有自己的特殊性，所以我们耗费了大量的时间做了一些亚群及其对应的marker的整理，如下所示；

• 间充质细胞的细分亚群情况
• 软骨相关单细胞亚群的标记基因
• 单细胞亚群的关键基因背诵不下来肿么办
• 内分泌器官胰岛的细分单细胞亚群
• NK细胞的单细胞层面细分亚群
• 呼吸道上皮各个单细胞亚群的层级结构
• 内皮细胞细分亚群
• 髓系免疫细胞细分亚群
• B细胞细分亚群
• 乳腺上皮细胞单细胞亚群
• 肝上皮细胞单细胞亚群
• 肺上皮细胞单细胞亚群
• 结直肠上皮细胞单细胞亚群
• 胃上皮细胞单细胞亚群
• 肾上皮细胞单细胞亚群

当前AI大语言模型逐渐改变我们的科研工作方式，甚至今年两项诺贝尔奖都给AI，我们从面向搜索引擎（谷歌，百度）工作，到面向AI（chatGPT,kimi）工作。从日常查信息翻译，到科研编程等，都会去问下AI。

生信技能树专栏 人工智能大模型的100个生信案例 提供了许多例子来演示AI在日常科研工作中的使用。

而且发在Nature methods上的GPTCelltype向我们阐述人工智能大语言模型在细胞类型注释的潜能 一行代码就能发生信顶刊的GPTCelltype做单细胞亚群注释

但显然，细胞类型注释，我们还可以通过大模型得到更多信息，例如:为什么注释到这个细胞类型？哪些gene markers 支持你注释这个细胞类型？因此这里介绍一个我们自己开发的基于Python编程语言的工具GPTBioInsightor辅助我们完成这件事。

GPTBioInsightor

GPTBioInsightor是基于Python开发的，所以我们这样安装软件：

pip install gptbioinsightor

这里也是使用经典的10X Genomics PBMC进行演示，Linux上这里下载数据

mkdir data
wget http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz -O data/pbmc3k_filtered_gene_bc_matrices.tar.gz
cd data; tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz

然后在Python环境进行单细胞数据处理, 进行预处理，单细胞降维聚类，然后计算差异基因marker：

# 更详细的scanpy数据处理参考 https://scanpy.readthedocs.io/en/stable/tutorials/basics/clustering-2017.html

import scanpy as sc

## 读取数据
adata = sc.read_10x_mtx(
    "data/filtered_gene_bc_matrices/hg19/",  # the directory with the `.mtx` file
    var_names="gene_symbols",  # use gene symbols for the variable names (variables-axis index)
    cache=True,  # write a cache file for faster subsequent reading
)

## 预处理
adata.var_names_make_unique()  

sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)

# annotate the group of mitochondrial genes as "mt"
adata.var["mt"] = adata.var_names.str.startswith("MT-")
sc.pp.calculate_qc_metrics(
    adata, qc_vars=["mt"], percent_top=None, log1p=False, inplace=True
)

adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :].copy()

sc.pp.normalize_total(adata, target_sum=1e4)

sc.pp.log1p(adata)

sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
adata.raw = adata
adata = adata[:, adata.var.highly_variable]

sc.pp.regress_out(adata, ["total_counts", "pct_counts_mt"])

sc.pp.scale(adata, max_value=10)
sc.tl.pca(adata, svd_solver="arpack")
sc.pp.neighbors(adata, n_neighbors=10, n_pcs=40)

## 聚类
sc.tl.leiden(
    adata,
    resolution=0.9,
    random_state=0,
    flavor="igraph",
    n_iterations=2,
    directed=False,
)
sc.tl.umap(adata)

## 计算差异基因marker
sc.tl.rank_genes_groups(adata, "leiden", key_added="logreg_deg", method="logreg")

计算完差异基因后，就可以使用GPTBioinsightor进行细胞类型注解，这里GPTBioinsightor可以选择不同的大模型，但需要你去申请API KEY，下面的大模型都可以试试看 :

• openai
• anthropic
• groq
• aliyun
• zhipuai
• siliconflow
• deepseek

其中咱们国内的友友们应该是用不了openai和anthropic，如果想要体验一下，考虑一下第三方中转服务，例如：https://flybirdsci.com/

在这个教程，我们使用阿里云的通义千文，API_KEY教程申请：https://help.aliyun.com/zh/model-studio/developer-reference/get-api-key

### 设置大语言模型的 API KEY
import os
# https://gptbioinsightor.readthedocs.io/zh-cn/latest/models.html
## 特定aliyun API KEY
os.environ['ALIYUN_API_KEY'] = "sk-***"
# os.environ['API_KEY'] = "sk-***"# 每个人自己的api的key都不一样哦 

import gptbioinsightor as gbi
# 设置数据的背景信息
background = "Cells are PBMCs from a Healthy Donor" 

# 使用阿里的通义千问qwen-max-latest
# 也可以设置成其他模型
res = gbi.get_celltype(adata, background=background, 
                       out="gbi.qwen.celltype.md", key="logreg_deg", 
                       topnumber=15,provider="aliyun", 
                       n_jobs=4,model="qwen-max-latest")
res
# {'0': 'CD4+ T Helper Cells',
#  '1': 'B Cells',
#  '2': 'Monocytes/Macrophages',
#  '3': 'Natural Killer (NK) cells',
#  '4': 'Cytotoxic T Cells (CD8+)',
#  '5': 'Monocytes/Macrophages',
#  '6': 'Dendritic Cells',
#  '7': 'Platelets'}

输出的markdown内容如下(这里只展示了一群的细胞类型)：主要包括：

• 关键marker: Cell-specific, Context-specific:
• 注释原因
• 备用细胞类型考虑
• 其他一些输出信息辅助判断

# CellType Analysis
## cluster geneset 0

### Gene List
'''
LDHB, LTB, RGCC, IL32, NOSIP, CD3D, CD3E, TMEM123, VIM, TMEM66, FYB, JUNB, CCR7, CD27, MYL12A
'''

### Celltype Prediction
#### Optimal Celltype: T Cells (likely CD4+ T cells)
**Key Markers**:
- Cell-specific: CD3D, CD3E, CCR7, CD27, FYB
- Context-specific: LTB, JUNB, VIM

**Evidence and Reasoning**
- **PRIMARY EVIDENCE**: The presence of CD3D and CD3E, which are essential components of the T cell receptor complex, strongly indicates a T cell population. These markers are highly specific to T cells.
- **SECONDARY EVIDENCE**: CCR7 is a chemokine receptor that is typically expressed on naïve and central memory T cells, suggesting that these cells may be in a non-activated or memory state.
- **ADDITIONAL EVIDENCE**: CD27 and FYB are also known to be expressed in T cells, particularly in activated T cells. LTB (lymphotoxin beta) and JUNB (a transcription factor) are involved in T cell activation and function.

**Validation**: Other gold standard markers for T cells (not in Geneset 0) include CD4, CD8, and TCRα/β. For CD4+ T cells, additional markers like CD45RA and CD45RO can be used to distinguish between naïve and memory T cells.

#### Alternative Considerations
- **Alternative celltype1: NK cells**
    - **WHY Alternative? Key MARKERS, Evidence and Reasoning**: NK cells can express some of the markers found in this geneset, such as CCR7 and CD27, but the presence of CD3D and CD3E, which are not expressed in NK cells, makes this less likely.
    - **OTHER Gold Standard MARKERS(NOT IN Geneset 0) TO VALIDATE THE Alternative celltype1**: NK cells would typically express NKp46, KIRs, and NKG2D, which are not present in this geneset.

- **Alternative celltype2: B cells**
    - **WHY Alternative? Key MARKERS, Evidence and Reasoning**: B cells do not typically express CD3D and CD3E, which are strong T cell markers. However, some B cell markers like CD27 and CCR7 are present, which might lead to confusion. The absence of B cell-specific markers like CD19 and CD20 makes this alternative less likely.
    - **OTHER Gold Standard MARKERS(NOT IN Geneset 0) TO VALIDATE THE Alternative celltype2**: B cells would typically express CD19, CD20, and surface IgM, which are not present in this geneset.

### Novel Insights
- **NOTEWORTHY PATTERNS**: The co-expression of CCR7 and CD27 suggests a population of T cells that are either naïve or central memory T cells.
- **CELL STATE**: The expression of JUNB and LTB, along with other activation-related genes, suggests that these T cells may be in an activated or recently activated state.
- **POTENTIAL NEW FINDINGS**: The presence of VIM (vimentin), a marker often associated with mesenchymal cells, in T cells is intriguing and could indicate a unique subset of T cells or a state of T cells that have undergone some form of stress or activation leading to the upregulation of vimentin.

和基于已知gene marker的手动注释的结果进行比较

cell_type_name = {
    "0": "CD4 T",
    "1": "B",
    "2": "FCGR3A+ Monocytes",
    "3": "NK",
    "4": "CD8 T",
    "5": "CD14+ Monocytes",
    "6": "Dendritic",
    "7": "Platelet",
}

adata.obs["celltype_manual"] = adata.obs["leiden"].map(
    cell_type_name
)
adata.obs["celltypes_gbi"] = adata.obs["leiden"].map(
    res
)
sc.pl.umap(adata, color=["leiden", "celltype_manual", "celltypes_gbi"], legend_loc="on data", frameon=False)

Seurat 用法

当前GPTBioinsightor 并没有R版本，如果你是Seurat用户，仍然可以较方便的使用GPTBioinsightor 下载数据：

mkdir data
wget http://cf.10xgenomics.com/samples/cell-exp/1.1.0/pbmc3k/pbmc3k_filtered_gene_bc_matrices.tar.gz -O data/pbmc3k_filtered_gene_bc_matrices.tar.gz
cd data; tar -xzf pbmc3k_filtered_gene_bc_matrices.tar.gz

然后在R里采用Seurat进行数据分析：

# 更详细的Seurat数据处理参考 https://satijalab.org/seurat/articles/pbmc3k_tutorial.html
library(dplyr)
library(Seurat)
library(patchwork)

# Load the PBMC dataset
pbmc.data <- Read10X(data.dir = "data/filtered_gene_bc_matrices/hg19/")
# Initialize the Seurat object with the raw (non-normalized data).
pbmc <- CreateSeuratObject(counts = pbmc.data, project = "pbmc3k", min.cells = 3, min.features = 200)
pbmc

pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & nFeature_RNA < 2500 & percent.mt < 5)
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, selection.method = "vst", nfeatures = 2000)
all.genes <- rownames(pbmc)
pbmc <- ScaleData(pbmc, features = all.genes)
pbmc <- RunPCA(pbmc, features = VariableFeatures(object = pbmc))

pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

pbmc <- RunUMAP(pbmc, dims = 1:10)

pbmc.markers <- FindAllMarkers(pbmc, only.pos = TRUE)
pbmc.markers.fil <- pbmc.markers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1 & p_val_adj < 0.01)

## 保存差异基因到一个CSV文件里
write.csv(pbmc.markers.fil, "data/pbmc.markers.fil.csv")

上面导出来了每个亚群的基因列表csv文件，接下来就可以去Python环境中使用GPTBioinsightor来读上面的基因列表csv文件进行细胞类型注解

### 设置大语言模型的 API KEY
import os
os.environ['API_KEY'] = "sk-***"

import gptbioinsightor as gbi

result_dict = gbi.get_marker_from_seurat("data/pbmc.markers.fil.csv")

# 设置数据的背景信息
background = "Cells are PBMCs from a Healthy Donor" 

# 使用阿里的通义千问qwen2-72b-instruct
# 也可以设置成openai 的模型
res = gbi.get_celltype(result_dict, background=background,
                       out="Seurat.qwen.celltype.md", n_jobs=4,
                       topnumber=15, provider="aliyun", model="qwen-max-latest")
res

#{0: ' T Cells',
# 1: ' Monocytes',
# 2: ' T Helper Cells',
# 3: ' B cells',
# 4: ' Cytotoxic T Cells',
# 5: ' Monocytes',
# 6: ' NK Cells',
# 7: ' Plasmacytoid Dendritic Cells (pDCs)',
# 8: ' Platelets'}

欢迎大家使用提出意见，AI只是辅助，需要小心求证。更多软件文档参见：https://gptbioinsightor.readthedocs.io/zh-cn/latest/index.html

load(file =  'check-by-0.5/qc-_marker_cosg.Rdata')
head(marker_cosg)
## Top10 genes
library(dplyr)  
cat(paste0('cluster',0:18,':',
           unlist(apply(marker_cosg$names,2,function(x){
             paste(head(x),collapse=',')
           })),'\n'))

cluster0:CD14,LGALS2,TMEM176B,CPVL,BLVRB,TMEM176A
 cluster1:VCAN,S100A8,LYZ,MNDA,S100A12,S100A9
 cluster2:S100A8,NEAT1,RGS2,CSF3R,VCAN,HBB
 cluster3:LEF1,MAL,CCR7,TCF7,IL7R,LTB
 cluster4:GNLY,KLRF1,SPON2,FGFBP2,KLRD1,HOPX
 cluster5:CDKN1C,HES4,LYPD2,C1QA,CKB,TCF7L2
 cluster6:GZMK,CD8A,CXCR6,CD8B,LAG3,KLRG1
 cluster7:GP9,TUBB1,CMTM5,C6orf25,SDPR,GNG11
 cluster8:PPBP,NRGN,CLU,PF4,PRKAR2B,F13A1
 cluster9:IFITM1,IL32,NKG7,PTPRCAP,CALM1,GIMAP7
 cluster10:FCER1A,ENHO,CD1C,CLEC10A,PKIB,CLIC2
 cluster11:ANK3,PARP15,SYNE2,ITK,RORA,SLFN12L
 cluster12:MS4A1,LINC00926,CD79A,VPREB3,FCER2,RP11-693J15.5
 cluster13:TNFRSF17,TXNDC5,IGLL5,IGF1,POU2AF1,IGJ
 cluster14:MMP8,MMP9,LTF,LCN2,PGLYRP1,ANXA3
 cluster15:RETN,SLC39A8,CENPF,MKI67,IL1R2,C19orf59
 cluster16:LILRA4,CLEC4C,LRRC26,SERPINF1,SCT,PTPRS
 cluster17:CD34,C19orf77,EMID1,AP001171.1,PRSS57,NPR3
 cluster18:LINC00278,KCNQ1OT1,RP11-362F19.1,CCRN4L,HCAR3,SPRED1

celltype[celltype$ClusterID %in% c( 0,1,2,14,15 ),2]='cd14-mono'  
celltype[celltype$ClusterID %in% c( 5,18),2]='cd16-mono'  
celltype[celltype$ClusterID %in% c( 10 ),2]='cDC' 
celltype[celltype$ClusterID %in% c( 16 ),2]='pDC' 
celltype[celltype$ClusterID %in% c( 3 ),2]='CD4' 
celltype[celltype$ClusterID %in% c( 4 ),2]='NK' 
celltype[celltype$ClusterID %in% c( 6,9,11 ),2]='CD8' 
celltype[celltype$ClusterID %in% c( 12 ),2]='Bcells' 
celltype[celltype$ClusterID %in% c( 13 ),2]='plasma'
celltype[celltype$ClusterID %in% c( 7,8 ),2]='megakaryocyte'