python单细胞学习笔记-day9(发在Science的celltypist软件注释)

生信技能树

发布于 2025-03-28 17:04:09

5500

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

前面的python学习笔记：：

今天继续学习视频：学习celltypist细胞注释，单样本和多样本~

0.环境准备

在之前的conda环境中安装新的本次上课需要的几个包：（建议先不要全都安装上，很容易后面出现包导入失败，可以运行到哪个代码缺少模块再开始安装~）

# bash终端
conda activate sc
# 安装 scikit-learn库
pip install singler -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install singlecellexperiment -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install scranpy -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install celldex -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install celltypist -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install harmonypy -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple
pip install loompy -i https://mirrors.tuna.tsinghua.edu.cn/pypi/web/simple

文件说明：

2.celltypist.ipynb：单样本的单细胞标准分析以及后面的celltypilst细胞注释
4.celltypist_multisamples.ipynb：多样本的单细胞标准分析以及后面的celltypilst细胞注释

1.单样本的celltypist注释

1.1 数据读取

使用用read_10x_mtx读取10x cellranger的三个标准文件，01_data 的文件：'barcodes.tsv', 'genes.tsv', 'matrix.mtx'

import singlecellexperiment as sce
import scanpy as sc
import os
print(os.listdir("01_data"))

# 读取数据
adata = sc.read_10x_mtx("01_data/")
print(adata.shape)
adata.X

1.2 数据质控

都是常规流程，快速看下数据分布特点决定数据过滤条件：

# 计算线粒体基因表达比例
adata.var['mt']=adata.var_names.str.startswith('MT-')
sc.pp.calculate_qc_metrics(adata,qc_vars=['mt'],log1p=False,percent_top=None,inplace=True)
# 过滤前数据分布情况
sc.pl.violin(adata,["n_genes_by_counts", "total_counts", "pct_counts_mt"],jitter=0.4, multi_panel=True)

过滤：这里这个数据我不太清楚来自哪里以及是什么组织类型，所以如果是自己的数据需要结合项目背景来看：

# 确定过滤条件并过滤
# 初始过滤
sc.pp.filter_cells(adata,min_genes=200)
sc.pp.filter_genes(adata,min_cells=3)
adata=adata[adata.obs.n_genes_by_counts>200]
adata=adata[adata.obs.n_genes_by_counts<2500]
adata=adata[adata.obs.pct_counts_mt<5].copy()
# 过滤后数据分布情况
sc.pl.violin(adata,["n_genes_by_counts", "total_counts", "pct_counts_mt"],jitter=0.4, multi_panel=True)
# 过滤后的细胞与基因数
print(adata.shape)

这好像是个过滤了的数据，小提琴过滤前后没啥变化：

1.3 数据标准化、降维聚类分群

还是scanpy的标准代码：

# 数据标准化
sc.pp.normalize_total(adata,target_sum=1e4)
sc.pp.log1p(adata)
adata.raw=adata
# 高变基因
sc.pp.highly_variable_genes(adata,n_top_genes=2000)
# 归一化
sc.pp.scale(adata)
sc.pp.pca(adata)
sc.pp.neighbors(adata,n_pcs=15)
sc.tl.leiden(adata,flavor="igraph",n_iterations=2,resolution=0.5)
sc.tl.umap(adata)
# 可视化
sc.pl.umap(adata,color='leiden')

到这里单细胞的标准分析就做完了，聚类结果如下：

1.4 celltypist自动注释

celltypist软件于2022年5月13号发表在顶刊Science杂志上，文献标题为《Cross-tissue immune cell analysis reveals tissue-specific features in humans》。

与别的自动注释算法不同，CellTypist可以自定义高精度和低精度，即直接注释出细胞的亚群。

这个软件已经训练好的参考模型：https://www.celltypist.org/models

首先下载model：

import celltypist
from celltypist import models

# 第一次使用时，会自动下载所有的参考数据
model = models.Model.load(model='Immune_All_High.pkl') #https://www.celltypist.org/models

动态日志：

🔎 No available models. Downloading...
📜 Retrieving model list from server https://celltypist.cog.sanger.ac.uk/models/models.json
📚 Total models in list: 54
📂 Storing models in /home/zhangj/.celltypist/data/models
💾 Downloading model [1/54]: Immune_All_Low.pkl
💾 Downloading model [2/54]: Immune_All_High.pkl
💾 Downloading model [3/54]: Adult_COVID19_PBMC.pkl
💾 Downloading model [4/54]: Adult_CynomolgusMacaque_Hippocampus.pkl
💾 Downloading model [5/54]: Adult_Human_MTG.pkl
💾 Downloading model [6/54]: Adult_Human_PancreaticIslet.pkl
💾 Downloading model [7/54]: Adult_Human_PrefrontalCortex.pkl
💾 Downloading model [8/54]: Adult_Human_Skin.pkl
💾 Downloading model [9/54]: Adult_Human_Vascular.pkl
💾 Downloading model [10/54]: Adult_Mouse_Gut.pkl
💾 Downloading model [11/54]: Adult_Mouse_OlfactoryBulb.pkl
💾 Downloading model [12/54]: Adult_Pig_Hippocampus.pkl
💾 Downloading model [13/54]: Adult_RhesusMacaque_Hippocampus.pkl
💾 Downloading model [14/54]: Autopsy_COVID19_Lung.pkl
💾 Downloading model [15/54]: COVID19_HumanChallenge_Blood.pkl
💾 Downloading model [16/54]: COVID19_Immune_Landscape.pkl
💾 Downloading model [17/54]: Cells_Adult_Breast.pkl
💾 Downloading model [18/54]: Cells_Fetal_Lung.pkl
💾 Downloading model [19/54]: Cells_Human_Tonsil.pkl
💾 Downloading model [20/54]: Cells_Intestinal_Tract.pkl
💾 Downloading model [21/54]: Cells_Lung_Airway.pkl
💾 Downloading model [22/54]: Developing_Human_Brain.pkl
💾 Downloading model [23/54]: Developing_Human_Gonads.pkl
💾 Downloading model [24/54]: Developing_Human_Hippocampus.pkl
💾 Downloading model [25/54]: Developing_Human_Organs.pkl
💾 Downloading model [26/54]: Developing_Human_Thymus.pkl
💾 Downloading model [27/54]: Developing_Mouse_Brain.pkl
💾 Downloading model [28/54]: Developing_Mouse_Hippocampus.pkl
💾 Downloading model [29/54]: Fetal_Human_AdrenalGlands.pkl
💾 Downloading model [30/54]: Fetal_Human_Pancreas.pkl
💾 Downloading model [31/54]: Fetal_Human_Pituitary.pkl
💾 Downloading model [32/54]: Fetal_Human_Retina.pkl
💾 Downloading model [33/54]: Fetal_Human_Skin.pkl
💾 Downloading model [34/54]: Healthy_Adult_Heart.pkl
💾 Downloading model [35/54]: Healthy_COVID19_PBMC.pkl
💾 Downloading model [36/54]: Healthy_Human_Liver.pkl
💾 Downloading model [37/54]: Healthy_Mouse_Liver.pkl
💾 Downloading model [38/54]: Human_AdultAged_Hippocampus.pkl
💾 Downloading model [39/54]: Human_Colorectal_Cancer.pkl
💾 Downloading model [40/54]: Human_Developmental_Retina.pkl
💾 Downloading model [41/54]: Human_Embryonic_YolkSac.pkl
💾 Downloading model [42/54]: Human_Endometrium_Atlas.pkl
💾 Downloading model [43/54]: Human_IPF_Lung.pkl
💾 Downloading model [44/54]: Human_Longitudinal_Hippocampus.pkl
💾 Downloading model [45/54]: Human_Lung_Atlas.pkl
💾 Downloading model [46/54]: Human_PF_Lung.pkl
💾 Downloading model [47/54]: Human_Placenta_Decidua.pkl
💾 Downloading model [48/54]: Lethal_COVID19_Lung.pkl
💾 Downloading model [49/54]: Mouse_Dentate_Gyrus.pkl
💾 Downloading model [50/54]: Mouse_Isocortex_Hippocampus.pkl
💾 Downloading model [51/54]: Mouse_Postnatal_DentateGyrus.pkl
💾 Downloading model [52/54]: Mouse_Whole_Brain.pkl
💾 Downloading model [53/54]: Nuclei_Lung_Airway.pkl
💾 Downloading model [54/54]: Pan_Fetal_Human.pkl
CellTypist model with 32 cell types and 6639 features
    date: 2022-07-16 08:53:00.959521
    details: immune populations combined from 20 tissues of 18 studies
    source: https://doi.org/10.1126/science.abl5197
    version: v2
    cell types: B cells, B-cell lineage, ..., pDC precursor
    features: A1BG, A2M, ..., ZYX

这是默认的下载目录：

下载的models：

predictions = celltypist.annotate(adata, model = model, majority_voting = True)
print(predictions.predicted_labels)

预测结果：这里使用第三列的结果作为细胞类型注释结果majority_voting

                         predicted_labels over_clustering majority_voting
AAACATACAACCAC-1          T cells              40         T cells
AAACATTGAGCTAC-1          B cells              14         B cells
AAACATTGATCAGC-1          T cells              10         T cells
AAACCGTGCTTCCG-1        Monocytes               6       Monocytes
AAACCGTGTATGCG-1              ILC              22             ILC
...                           ...             ...             ...
TTTCGAACTCTCAT-1        Monocytes               1       Monocytes
TTTCTACTGAGGCA-1     Plasma cells              14         B cells
TTTCTACTTCCTCG-1          B cells               8         B cells
TTTGCATGAGAGGC-1          B cells               8         B cells
TTTGCATGCCTCAC-1          T cells              34         T cells

整理一下预测结果：

adata2 = predictions.to_adata()
print(adata2.obs)

# 修改一下预测的细胞名字，个人习惯，也可以不改
# 去掉名字冲的cells
adata2.obs.majority_voting = adata2.obs.majority_voting.str.replace(' cells', '', regex=False)
# 把monocytes改为Mono
adata2.obs.majority_voting = adata2.obs.majority_voting.str.replace(' monocytes', ' Mono', regex=False)
adata2.obs.majority_voting

整理之后：

AAACATACAACCAC-1            T
AAACATTGAGCTAC-1            B
AAACATTGATCAGC-1            T
AAACCGTGCTTCCG-1    Monocytes
AAACCGTGTATGCG-1          ILC
                      ...    
TTTCGAACTCTCAT-1    Monocytes
TTTCTACTGAGGCA-1            B
TTTCTACTTCCTCG-1            B
TTTGCATGAGAGGC-1            B
TTTGCATGCCTCAC-1            T

umap可视化看一下：

sc.tl.umap(adata2)
sc.pl.umap(adata2, color = 'majority_voting',legend_loc="on data")

2.多样本的celltypist注释

这里与单样本的区别在于数据的读取和合并，先来读取一下：读取的数据为02_data，这里我也没有获取到数据的背景信息（数据是什么组织，如何取样，以便了解里面可能包含的细胞类型，这是做细胞类型注释必须了解的内容，后面celltypist的模型选择也会有参考），所以这里就看看代码就好了。

02_data的目录下有三个样本，'TD1', 'TD2', 'TD3'，每个样本下面是标准的三个文件：'barcodes.tsv.gz', 'features.tsv.gz', 'matrix.mtx.gz'

import scanpy as sc
import os
import pandas as pd
from scipy.io import mmread
import anndata

# 列出目录下的文件有哪些
path1='02_data/'
os.listdir(path1)
os.listdir('02_data/TD1/')

关于python读取单细胞，各种格式，可以看我们前面的帖子：python版读取不同的单细胞数据格式（单样本与多样本）

多样本读取：用for循环批量读取数据

# 用for循环批量读取数据
adata = {}
for i in range(len(os.listdir(path1))):
    data = sc.read_10x_mtx(path1 +os.listdir(path1)[i], var_names="gene_symbols", cache=True)
    data.var_names_make_unique()  
    adata[os.listdir(path1)[i]] = data
    print([f'{os.listdir(path1)[i]}:done']+list([data.shape]))
    
# 把多个数据拼接到一起，成为一个大的anndata对象
adata = sc.concat(adata,label='sample')
adata.obs_names_make_unique()
adata

# 每个样本里的细胞数量
adata.obs.sample.value_counts()

然后是质控：

sc.pp.filter_cells(adata,min_genes=300)
sc.pp.filter_genes(adata,min_cells=5)
#计算线粒体基因比例
adata.var['mt']=adata.var_names.str.startswith('MT-')
# 计算核糖体基因比例
adata.var['ribo']=adata.var_names.str.match('^RP[SL]')
# 计算红血细胞基因比例
adata.var['hb']=adata.var_names.str.match('^HB[^(P)]')
sc.pp.calculate_qc_metrics(adata,qc_vars=['mt','ribo','hb'],percent_top=None,log1p=False,inplace=True)
# 可视化细胞的上述比例情况

sc.pl.violin(adata,['n_genes_by_counts','total_counts'],groupby='batch',jitter=False)

质控前的小提琴图：

过滤细胞：

adata = adata[adata.obs.n_genes_by_counts < 6000, :]
adata = adata[adata.obs.pct_counts_mt < 18, :].copy()

又是常规的降维聚类：

sc.pp.normalize_total(adata,target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata.copy()
#celltypist包需要的数据是log normlize之后的数据，会使用这里的adata.raw。
sc.pp.highly_variable_genes(adata)
sc.pp.scale(adata)
sc.tl.pca(adata,use_highly_variable=True)
sc.external.pp.harmony_integrate(adata,'batch')
#聚类
sc.pp.neighbors(adata,use_rep='X_pca_harmony',n_pcs=15)
sc.tl.umap(adata)
sc.tl.leiden(adata, flavor="igraph", n_iterations=2,resolution=0.5)
sc.pl.umap(adata,color='leiden',legend_loc='on data')

注释：

import celltypist
from celltypist import models
model = models.Model.load(model='Immune_All_High.pkl') #https://www.celltypist.org/models
model
predictions = celltypist.annotate(adata, model = model, majority_voting = True)
print(predictions.predicted_labels)
adata2 = predictions.to_adata()
print(adata2.obs)
adata2.obs.majority_voting = adata2.obs.majority_voting.str.replace(' cells', '', regex=False)
adata2.obs.majority_voting = adata2.obs.majority_voting.str.replace(' monocytes', ' Mono', regex=False)
adata2.obs.majority_voting
sc.tl.umap(adata2)
sc.pl.umap(adata2, color = 'majority_voting',legend_loc="on data")

注释结果如下：