公众号里的文章大多数需要编程基础，如果因为代码看不懂，而跟不上正文的节奏，可以来找我学习，相当于给自己一个新手保护期。我的课程都是循环开课，点进去咨询微信↓ 生信分析直播课程(每月初开一期，春节休一个月) 生信新手保护学习小组（每月两期） 单细胞陪伴学习小组（每月两期)
1.R包准备和数据读取
# pip install scanpy==1.9.3 -i https://pypi.tuna.tsinghua.edu.cn/simple
# pip install leidenalg scikit-misc -i https://pypi.tuna.tsinghua.edu.cn/simple
import scanpy as sc
import os
import pandas as pd
from scipy.io import mmread
import anndata
path1='04_data/'
os.listdir(path1)
['TD1', 'TD2', 'TD3']
os.listdir('04_data/TD1/')
['barcodes.tsv.gz', 'features.tsv.gz', 'matrix.mtx.gz']
如上所见，工作目录下的04_data文件夹下有3个子文件夹，每个文件夹里是3个标准文件
读取单个样本数据的代码先跑通，然后才能搞个循环批量读取
adata = sc.read_10x_mtx(
path1+os.listdir(path1)[0], # 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()
# this is unnecessary if using `var_names='gene_ids'` in `sc.read_10x_mtx`
adata
AnnData object with n_obs × n_vars = 15216 × 33538
var: 'gene_ids', 'feature_types'
python里面不能直接写1:3这样的向量单独使用，只有取子集时才能这样用。单独使用连续的数字要写：
list(range(3))
[0, 1, 2]
用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]))
['TD1:done', (15216, 33538)]
['TD2:done', (18064, 33538)]
['TD3:done', (12658, 33538)]
把多个数据拼接到一起，成为一个大的anndata对象
adata=sc.concat(adata,label='batch')
adata.obs_names_make_unique()
adata
AnnData object with n_obs × n_vars = 45938 × 33538
obs: 'batch'
adata.obs.batch.value_counts() #每个样本里的细胞数量
batch
TD2 18064
TD1 15216
TD3 12658
Name: count, dtype: int64
2.质控
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)
sc.pl.violin(adata,['pct_counts_mt', 'pct_counts_ribo','pct_counts_hb'],groupby='batch',xlabel='sample',jitter=False)
sc.pl.scatter(adata,x='total_counts',y='n_genes_by_counts',color='batch')
adata = adata[adata.obs.n_genes_by_counts < 6000, :]
adata = adata[adata.obs.pct_counts_mt < 18, :].copy()
3.常规的降维聚类流程
celltypist包需要的数据是log normlize之后的数据，会使用这里的adata.raw。
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,flavor='seurat_v3',n_top_genes=2000)
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)
2024-11-15 13:02:16,497 - harmonypy - INFO - Computing initial centroids with sklearn.KMeans...
Computing initial centroids with sklearn.KMeans...
2024-11-15 13:02:18,831 - harmonypy - INFO - sklearn.KMeans initialization complete.
sklearn.KMeans initialization complete.
2024-11-15 13:02:19,329 - harmonypy - INFO - Iteration 1 of 10
Iteration 1 of 10
2024-11-15 13:02:35,549 - harmonypy - INFO - Iteration 2 of 10
Iteration 2 of 10
2024-11-15 13:02:51,986 - harmonypy - INFO - Iteration 3 of 10
Iteration 3 of 10
2024-11-15 13:03:08,223 - harmonypy - INFO - Iteration 4 of 10
Iteration 4 of 10
2024-11-15 13:03:24,727 - harmonypy - INFO - Iteration 5 of 10
Iteration 5 of 10
2024-11-15 13:03:34,032 - harmonypy - INFO - Iteration 6 of 10
Iteration 6 of 10
2024-11-15 13:03:40,914 - harmonypy - INFO - Iteration 7 of 10
Iteration 7 of 10
2024-11-15 13:03:47,042 - harmonypy - INFO - Iteration 8 of 10
Iteration 8 of 10
2024-11-15 13:03:53,211 - harmonypy - INFO - Converged after 8 iterations
Converged after 8 iterations
sc.tl.leiden(adata, flavor="igraph", n_iterations=2,resolution=0.5)
sc.pl.umap(adata,color='leiden',legend_loc='on data')
4.完成注释
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)
predictions.predicted_labels
👀 Invalid expression matrix in `.X`, expect log1p normalized expression to 10000 counts per cell; will use `.raw.X` instead
🔬 Input data has 43649 cells and 23244 genes
🔗 Matching reference genes in the model
🧬 5842 features used for prediction
⚖️ Scaling input data
🖋️ Predicting labels
✅ Prediction done!
👀 Detected a neighborhood graph in the input object, will run over-clustering on the basis of it
⛓️ Over-clustering input data with resolution set to 20
🗳️ Majority voting the predictions
✅ Majority voting done!
adata2 = predictions.to_adata()
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
AAACCCAAGACCTTTG-1 T
AAACCCAAGCCTGAAG-1 Epithelial
AAACCCAAGCTGTGCC-1 T
AAACCCAAGGCCACTC-1 T
AAACCCAAGGGTTAAT-1 T
...
TTTGTTGCAGCGACAA-1 ILC
TTTGTTGCATCCGATA-1 T
TTTGTTGTCAGCGGAA-1 Endothelial
TTTGTTGTCCAACTGA-1 ILC
TTTGTTGTCGGAACTT-1 ILC
Name: majority_voting, Length: 43649, dtype: object
sc.tl.umap(adata2)
sc.pl.umap(adata2, color = 'majority_voting',legend_loc="on data")
... storing 'majority_voting' as categorical
5.画个好看点的图
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib.colors as mcolors
import numpy as np
umap_coords = adata2.obsm['X_umap']
labels = adata2.obs['majority_voting']
plt.figure(figsize=(6,5))
unique_labels = labels.unique()
colors = plt.cm.get_cmap('Accent', len(unique_labels))
original_colors = plt.cm.Set1(np.linspace(0, 1, 9))
# 使用LinearSegmentedColormap.from_list创建一个新的颜色映射，包含100种渐变色
colors = mcolors.LinearSegmentedColormap.from_list('colors', original_colors, N=len(unique_labels))
for i, label in enumerate(unique_labels):
plt.scatter(umap_coords[labels == label, 0], umap_coords[labels == label, 1],
color=colors(i),s=1,alpha = 0.4)
for label in unique_labels:
label_coords = umap_coords[labels == label]
center_x = label_coords[:, 0].mean()
center_y = label_coords[:, 1].mean()
plt.text(center_x, center_y, label, fontsize=10, ha='center', va='center')
plt.grid(False)
plt.title('majority_voting', fontsize=14)
plt.xlabel('UMAP 1', fontsize=12)
plt.ylabel('UMAP 2', fontsize=12)
plt.show()
