对于用惯了R的人，冷不丁用下python，感觉还是不太习惯，但是python的功能做分析的速度也很有独特的优势，所以本次做了个简单的scanpy笔记，以后再慢慢学习~~
1 加载需要的包，感觉python的包比较好装，缺啥补啥吧
import numpy as np
import pandas as pd
import scanpy as sc
import matplotlib.pyplot as plt
2 读入数据
data=sc.read_10x_mtx('Mouse_bonemarrow_matrix/Marrow_V2_1_matrix_10X', cache=True)
# 当然scanpy可以直接读取10Xgenomics的.h5格式数据
# data=sc.read_10x_h5("./pro.h5",genome=None,gex_only=True)
# adata.var_names_make_unique()
results_file = 'bonemarrow.h5ad' #我们把后边的数据存在这个h5ad的文件中
3 数据预处理
#显示在所有细胞中在每个单细胞中产生最高计数分数的基因
sc.pl.highest_expr_genes(data, n_top=20, )
#过滤低质量细胞样本，过滤在少于三个细胞中表达，或一个细胞中表达少于200个基因的细胞样本
sc.pp.filter_cells(data, min_genes=200)
sc.pp.filter_genes(data, min_cells=3)
#查看线粒体，基因UMI的表达情况
data.var['mt'] = data.var_names.str.startswith('mt-') # 将线粒体基因标记为 mt
sc.pp.calculate_qc_metrics(data, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True) # 查看线粒体基因和总的UMI数的关系
sc.pl.violin(data, ['n_genes_by_counts', 'total_counts', 'pct_counts_mt'],jitter=0.4, multi_panel=True) #查看UMI输一局，基因数目和线粒体基因的数目
sc.pl.scatter(data, x='total_counts', y='pct_counts_mt') #查看UMI数目和mt的散点图相关性
sc.pl.scatter(data, x='total_counts', y='n_genes_by_counts') #查看UMI数目和gene的散点图相关性
# 获取线粒体基因占比在 5% 以下的细胞样本
data = data[data.obs.pct_counts_mt < 5, :]
# 获取表达基因数在 2500 以下的细胞样本
data = data[data.obs.n_genes_by_counts < 5000, :]
image.png
4 数据标准化
# 归一化，使得不同细胞样本间可比
sc.pp.normalize_total(data, target_sum=1e4)
sc.pp.log1p(data)
data.raw = data #储存标准化后的AnnaData Object
#识别特异性基因
sc.pp.highly_variable_genes(data, min_mean=0.0125, max_mean=3, min_disp=0.5)
# 绘制特异性基因散点图
#sc.pl.highly_variable_genes(data)
#绘制只有特异性基因的数据集
data = data[:, data.var.highly_variable]
# 回归每个细胞的总计数和表达的线粒体基因的百分比的影响。
sc.pp.regress_out(data, ['total_counts', 'pct_counts_mt'])
# 将每个基因缩放到单位方差。阈值超过标准偏差 10。
sc.pp.scale(data, max_value=10)
#主成分分析
# 绘制 PCA 图
#sc.pl.pca(data, color='Actb') #在这里选择了Actb这个基因，如果用这个基因画不出来可以换个确定在自己数据中表达的基因
#sc.pl.pca_variance_ratio(data, log=True)
#data.write(results_file)
# 构建图
sc.pp.neighbors(data, n_neighbors=10, n_pcs=30)
sc.tl.umap(data)
sc.tl.leiden(data, key_added="clusters")
#help(sc.pl.umap)
sc.pl.umap(data, color=["total_counts", "clusters"],legend_loc='on data', wspace=0.2)
image.png
5 注释
画出了umap图那下一步就要注释了，由于没找到自动化注释软件，所以只能用画小提琴图的方式画marker基因，还是比较麻烦的
#查看差异基因
# 查看每个cluster中markergene
sc.tl.rank_genes_groups(data, 'clusters', method='wilcoxon')
sc.pl.rank_genes_groups(data, n_genes=25, sharey=False)
pd.DataFrame(data.uns['rank_genes_groups']['names']).head(5) #定义基因列表，输出每个cluster中的marker基因
plt.savefig("dif_gene.pdf")
res<-pd.DataFrame(data.uns['rank_genes_groups']['names'])
res.to_csv("diff.csv") #把差异基因写入表格
# 查看甘兴趣的cluster的差异基因
sc.tl.rank_genes_groups(adata, 'cluster', groups=['0'], reference='1', method='wilcoxon')
sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)
#通过小提琴图的方式看marker基因的表达情况，用这个方式注释还是比较麻烦的，可能是我用惯了简便方法的原因吧，小伙伴们有什么好用的方法也可以安利给我，感谢ing~~
sc.pl.violin(data, ['Lyz2','Vcan','Fn1','Ccr2','Ly6c2'], groupby='clusters')
#plt.savefig("marker_exp_Vio.pdf")
sc.pl.umap(data, color=['Lyz2','Vcan','Fn1','Ccr2','Ly6c2'])
#plt.savefig("marker_exp_umap.pdf")
image.png
6 确定了细胞亚型，就是要把细胞亚型映射回去
new_cluster_names = ['Neutrophils','Monocytes','Neutrophils', 'Neutrophils','Monocytes','pDCs','Tcells','GMP','Neutrophils','GMP','ProEryth','Neutrophils','Bcells','ProEryth','Neutrophils','Monocytes','Plasma cells','cDCs','Pre-B','Basophils','Bcells','GMP','HSCs','GMP','cDCs','Monocytes','HSCs','Erythroblast','pDCs','GMP']
new_cluster_names_unique=list(set(new_cluster_names))
new_cluster_names_unique.sort(key=new_cluster_names.index)
celltype = [new_cluster_names[int(i)] for i in data.obs.clusters.values.tolist()]
data.obs.clusters = celltype
data.obs.clusters = data.obs.clusters.astype('category')
data.obs.clusters.cat.set_categories(new_cluster_names_unique, inplace=True)
#data.obs.head(10) #data.obs相当于seraut的meta.data，可以看有没有注释成功
sc.pl.umap(data, color='clusters', legend_loc='on data', title='', frameon=True)
image.png
本次测试先到这里，其实这个软件的功能很强大，等我用到的时候再做更深度的测试吧~
