预处理与分群的步骤大体和Ｒ包 Ｓｅｕｒａｔ（https://satijalab.org/seurat/v3.2/pbmc3k_tutorial.html）类似哦，
数据的下载和目录的创建可以参考一下代码
[1]
# !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
# !mkdir write
引入包
[2]
import pandas as pd
import scanpy as sc
import scanpy as sc
设置和查看当前环境
[3]
sc.settings.verbosity = 3 # verbosity: errors (0), warnings (1), info (2), hints (3)
sc.logging.print_header()
sc.settings.set_figure_params(dpi=80, facecolor='white')
指定结果路径
[4]
results_file = 'write/pbmc3k.h5ad' # the file that will store the analysis results
读入单细胞测序文件为 AnnData 对象，它包括许多注释和代表鼠的slots,它有自己的hdf5格式
读入数据
[5]
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
[6]
adata.var_names_make_unique() # this is unnecessary if using `var_names='gene_ids'` in `sc.read_10x_mtx`
[7]
adata
预处理
显示所有细胞中每个单个细胞中计数最高的基因。
[8]
sc.pl.highest_expr_genes(adata, n_top=20, )
normalizing counts per cell finished (0:00:00)
基本的数据过滤
[9]
sc.pp.filter_cells(adata, min_genes=200)
sc.pp.filter_genes(adata, min_cells=3)
filtered out 19024 genes that are detected in less than 3 cells
我们还需要收集一些有关线粒体的信息，因为这些信息对于质量控制很重要
高比例线粒体基因表明细胞质量较差（https://master.bioconductor.org/packages/release/workflows/html/simpleSingleCell.html#examining-gene-level-metrics），这可能是因为细胞穿孔后细胞质RNA丢失。线粒体大于单个转录物分子，不太可能通过细胞膜逸出。
使用 pp.calculate_qc_metrics，我们可以非常有效地计算许多指标
[10]
adata.var['mt'] = adata.var_names.str.startswith('MT-') # annotate the group of mitochondrial genes as 'mt'
sc.pp.calculate_qc_metrics(adata, qc_vars=['mt'], percent_top=None, log1p=False, inplace=True)
一些计算质量测度的小提琴图：
count 数中表达的基因数量
每个细胞的总count数
线粒体基因count数的百分比
[11]
sc.pl.violin(adata, ['n_genes_by_counts'])
sc.pl.violin(adata, ['total_counts'])
sc.pl.violin(adata, ['pct_counts_mt'])
(按官网给的代码画出来是横着的，为什么呢？有一样的小伙伴留言，不知道是不是环境包版本不同的原因)
需要删除表达太多线粒体基因或总数过多的细胞。
[12]
sc.pl.scatter(adata, x='total_counts', y='pct_counts_mt')
sc.pl.scatter(adata, x='total_counts', y='n_genes_by_counts')
实际上是通过 AnnData对象 的slots来进行过滤的。
[13]
adata = adata[adata.obs.n_genes_by_counts < 2500, :]
adata = adata[adata.obs.pct_counts_mt < 5, :]
总count数归一化（文库大小正确）数据矩阵 X 为每个细胞最多10,000reads，因此细胞的count数就可比。
[14]
sc.pp.normalize_total(adata, target_sum=1e4)
normalizing counts per cell
finished (0:00:00)
对数据进行对数转换
[15]
sc.pp.log1p(adata)
鉴定高度表达变化的基因
[16]
sc.pp.highly_variable_genes(adata, min_mean=0.0125, max_mean=3, min_disp=0.5)
[17]
sc.pl.highly_variable_genes(adata)
将.AnnData对象的.raw 属性设置为归一化和对数化的原始基因表达，以供之后在差异测试和基因表达的可视化中使用。相当于是冻结了AnnData对象的状态。
可以通过调用.raw.to_adata() 来获取AnnData对象.raw中的对象
[18]
adata.raw = adata
如果在后面没有继续使用进行校正数据做 sc.pp.regress_out和sc.pp.scale 缩放它，也可以完全不用使用.raw。
先前的高可变基因检测结果将作为注释存储在.var.highly_variabl中，并会被PCA ，sc.pp.neighbors 以及随后的流形/图形工具自动检测。在这种情况下，不需要执行下面的过滤步骤。
[19]
adata = adata[:, adata.var.highly_variable]
淘汰每个细胞的总计数和表达的线粒体基因百分比的影响。将数据缩放为单位方差。
[20]
sc.pp.regress_out(adata, ['total_counts', 'pct_counts_mt'])
将每个基因缩放至单位方差，剪辑超过标准偏差10的值
[21]
sc.pp.scale(adata, max_value=10)
主成分分析
通过运行主成分分析（PCA）来降低数据的维数，该分析可以揭示变化的主轴并对数据进行去噪。
[22]
sc.tl.pca(adata, svd_solver='arpack')
我们可以在PCA坐标中绘制散点图
[23]
sc.pl.pca(adata, color='CST3')
查看单个PC对数据总方差的贡献。这为我们提供了有关为计算细胞的群关系时（例如，在聚类函数sc.tl.louvain()或tSNE sc.tl.tsne() 中使用）应考虑的PC数量。根据经验，通常可以粗略估计PC的数量。
[24]
sc.pl.pca_variance_ratio(adata, log=True)
保存结果
[25]
adata.write(results_file)
[26]
adata
计算邻域图
使用数据矩阵的PCA表示来数据矩阵的聚类图。可以在此处简单地使用默认值。为了再现此数据Seurat的结果，我们采用以下值。
[27]
sc.pp.neighbors(adata, n_neighbors=175, n_pcs=40)
(版本不一样，使用同样参数结果也会不一样呢，n_neighbors指的是每个点的邻近点的数量，neighbors的个数越多，聚类数会越少。需要自己调整 参数 对比seurat 的分类结果参考是否正确，我也是调了很多次参数才细胞分成类似的哦)
pc的数量依赖于上面所做图，一般是选在拐点处的的主成分，只需要一个粗略值，不同的pc数量所产生的聚类形状也不同
嵌入邻域图
建议使用UMAP将图形嵌入二维中（[McInnes et al，2018] https://arxiv.org/abs/1802.03426。它比tSNE更忠实于流形的全局连通性，即它可以更好地保存轨迹。在某些情况下，可能会发现群集断开连接和类似的连接冲突。通常可以通过运行以下命令来补救它们：
tl.paga(adata)
pl.paga(adata, plot=False) # remove `plot=False` if you want to see the coarse-grained graph
tl.umap(adata, init_pos='paga')
[28]
sc.tl.umap(adata)
[29]
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'])
当设置adata的.raw属性时，先前的图显示了“raw”（标准化，对数化但未校正）的基因表达。您还可以通过明确声明不想使用来绘制缩放和校正后的基因表达。
[30]
sc.pl.umap(adata, color=['CST3', 'NKG7', 'PPBP'], use_raw=False)
聚类邻域图
与Seurat和其他一样，scanpy 推荐Traag 等人（2018）的Leiden图聚类方法（基于优化模块化的社区检测）。请注意，Leiden聚类直接将聚类邻域图细胞，这已在上一节中进行了计算。
[31]
sc.tl.leiden(adata)
绘制群集，这与Seurat的结果非常吻合
[32]
sc.pl.umap(adata, color=['leiden', 'CST3', 'NKG7'])
保存结果
[33]
adata.write(results_file)
寻找基因Marker
让我们计算每个簇中高度差异化基因的排名。为此，默认情况下AnnData属性.raw 被使用如果之前已被初始化。最简单，最快的方法是t检验。
[34]
sc.tl.rank_genes_groups(adata, 'leiden', method='t-test')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
[35]
sc.settings.verbosity = 2 # reduce the verbosity
Wilcoxon rank-sum (Mann-Whitney-U)
的结果的测试是非常相似的,建议在出版物中使用后者，例如，参见Sonison & Robinson (2018).
。您可能还会考虑功能更强大的差异测试程序包，例如MAST，limma，DESeq2，对于python，则是最新的diffxpy。
[36]
sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
保存结果
[37]
adata.write(results_file)
或者，还可以使用逻辑回归对基因进行排名。例如，Natranos et al. (2018).
已经提出了这一点。本质上的区别在于，在这里，使用多变量方法，而传统的差异检验是单变量的， Clark et al. (2014)有更多细节。
[38]
sc.tl.rank_genes_groups(adata, 'leiden', method='logreg')
sc.pl.rank_genes_groups(adata, n_genes=25, sharey=False)
除 IL7R，其仅由t检验发现，FCER1A，仅由其他两个方法发现，其它的所有的标记基因在所有的方法中发现。
| Louvain Group | Markers | Cell Type |
|---|---|---|
| 0 | IL7R | CD4 T cells |
| 1 | CD14, LYZ | CD14+ Monocytes |
| 2 | MS4A1 | B cells |
| 3 | CD8A | CD8 T cells |
| 4 | FCGR3A, MS4A7 | FCGR3A+ Monocytes |
| 5 | GNLY, NKG7 | NK cells |
| 6 | FCER1A, CST3 | Dendritic Cells |
| 7 | PPBP | Megakaryocytes |
定义标记基因的列表，以供参考
[39]
marker_genes = ['IL7R', 'CD79A', 'MS4A1', 'CD8A', 'CD8B', 'LYZ', 'CD14',
'LGALS3', 'S100A8', 'GNLY', 'NKG7', 'KLRB1',
'FCGR3A', 'MS4A7', 'FCER1A', 'CST3', 'PPBP']
重新加载已保存的Wilcoxon Rank-Sum 测试结果的对象。
在数据框中显示每个聚类0、1，...，7的前10个排名最高的基因
[40]
pd.DataFrame(adata.uns['rank_genes_groups']['names']).head(5)
[41]
获取带有评分和组的表。
[42]
result = adata.uns['rank_genes_groups']
groups = result['names'].dtype.names
pd.DataFrame(
{group + '_' + key[:1]: result[key][group]
for group in groups for key in ['names', 'pvals']}).head(5)
与单个群集进行比较
[43]
sc.tl.rank_genes_groups(adata, 'leiden', groups=['0'], reference='1', method='wilcoxon')
sc.pl.rank_genes_groups(adata, groups=['0'], n_genes=20)
如果想要某个组的更详细的视图，使用 sc.pl.rank_genes_groups_violin。
[44]
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
与其余组进行比较，重新加载对象计算差异表达式
[45]
adata = sc.read(results_file)
[46]
sc.pl.rank_genes_groups_violin(adata, groups='0', n_genes=8)
如果要跨组比较某个基因，可以使用
[47]
sc.pl.violin(adata, ['CST3', 'NKG7', 'PPBP'], groupby='leiden')
实际标记细胞类型
[48]
new_cluster_names = [
'CD4 T', 'CD14 Monocytes',
'B', 'CD8 T',
'FCGR3A Monocytes','NK',
'Dendritic', 'Megakaryocytes']
adata.rename_categories('leiden', new_cluster_names)
[49]
sc.pl.umap(adata, color='leiden',size=50, legend_loc='on data',legend_fontsize=6, title='', frameon=False, save='.pdf')
对比官网，基本一样了
注释了细胞类型可视化标记基因。
[50]
sc.pl.dotplot(adata, marker_genes, groupby='leiden');
还有一个非常紧凑的小提琴图
[51]
sc.pl.stacked_violin(adata, marker_genes, groupby='leiden', rotation=90);
在此分析过程中，AnnData累积了以下注释
[52]
adata
[53]
adata.write(results_file, compression='gzip') # `compression='gzip'` saves disk space, but slows down writing and subsequent reading
使用下来 可以大致了解该文件h5ls，该文件有很多选项-有关更多详细信息，请参见此处。文件格式将来可能仍会受到进一步优化。但是，所有读取功能都将保持向后兼容。
如果要与只想使用它进行可视化的人共享该文件，则减小文件大小的一种简单方法是删除密集缩放和校正的数据矩阵。该文件仍包含的可视化文件使用的原始数据为adata.raw。
[54]
adata.raw.to_adata().write('./write/pbmc3k_withoutX.h5ad')
[55]
如果要导出到“ csv”，则可以使用以下选项：
# Export single fields of the annotation of observations
# adata.obs[['n_counts', 'louvain_groups']].to_csv(
# './write/pbmc3k_corrected_louvain_groups.csv')
# Export single columns of the multidimensional annotation
# adata.obsm.to_df()[['X_pca1', 'X_pca2']].to_csv(
# './write/pbmc3k_corrected_X_pca.csv')
# Or export everything except the data using `.write_csvs`.
# Set `skip_data=False` if you also want to export the data.
# adata.write_csvs(results_file[:-5], )
