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

pip install 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

# 更详细的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")

### 设置大语言模型的 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'}

# 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.

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)

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

# 更详细的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")

### 设置大语言模型的 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'}

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'