1.基于python的单细胞数据预处理-降维可视化

参考：
[1] https://github.com/Starlitnightly/single_cell_tutorial
[2] https://github.com/theislab/single-cell-best-practices

降维的背景

虽然特征选择已经减少了维数，但为了可视化，我们需要更直观的降维方法。降维将高维数据嵌入到低维空间中。低维表示仍然捕获数据的基本结构，同时尽可能少地持有维度。比如下图，我们将三维对象可视化为投影到二维空间中。

过去的研究独立比较了10种不同的降维方法的稳定性，准确性和计算成本。他们建议使用t-分布随机邻居嵌入（t-SNE），因为它产生了最佳的性能。统一流形逼近和投影（UMAP）显示出最高的稳定性，并且最好地分离了原始细胞群体。在这种情况下，值得一提的另一种降维方法是主成分分析（PCA），它仍然被广泛使用。

首先我们加载来自特征选择处理的数据：

import omicverse as ov
import scanpy as sc

ov.ov_plot_set()

adata = sc.read("./data/s4d8_preprocess.h5ad")
print(adata)
print(adata.X.max()) # 10.989398

在标准流程中，按照特征选择出的特征切片得到新矩阵后，我们还需要对新矩阵的计数进行scale（为了利于降维）。在omicverse中，我们将scale后的值存放进adata.layer，而不是像scanpy一样默认取代adata.X。但是注意，没有任何的证据表明，数据经过scale后会取得更好的差异基因分析结果，若盲目地使用scale后的计数值，还可能会导致使用dotplot或者violinplot中忽略了基因自身的特征信息，差异基因分析最好是使用缩放前的X。

比如我们关注的一个基因A的表达值在17-20区间，而基因B的表达值在0-3的区间，经过scale后，由于平均值被缩放成了0，基因A和基因B都在-2-2的区间范围内，这一定程度上失去了基因A表达量高的信息。故scale对差异基因分析似乎是不利的。

scale如下：

# 缩放
ov.pp.scale(adata,max_value=10)
print(adata)

"""
AnnData object with n_obs × n_vars = 14814 × 2000
    obs: 'n_genes_by_counts', 'log1p_n_genes_by_counts', 'total_counts', 'log1p_total_counts', 'pct_counts_in_top_20_genes', 'total_counts_mt', 'log1p_total_counts_mt', 'pct_counts_mt', 'total_counts_ribo', 'log1p_total_counts_ribo', 'pct_counts_ribo', 'total_counts_hb', 'log1p_total_counts_hb', 'pct_counts_hb', 'outlier', 'mt_outlier', 'scDblFinder_score', 'scDblFinder_class'
    var: 'gene_ids', 'feature_types', 'genome', 'mt', 'ribo', 'hb', 'n_cells_by_counts', 'mean_counts', 'log1p_mean_counts', 'pct_dropout_by_counts', 'total_counts', 'log1p_total_counts', 'n_cells', 'percent_cells', 'robust', 'mean', 'var', 'residual_variances', 'highly_variable_rank', 'highly_variable_features'
    uns: 'hvg', 'layers_counts', 'log1p'
    layers: 'counts', 'soupX_counts', 'scaled'
"""

可以发现adata的layers层多出了一个scaled，这就是我们经过scale后的数据。

PCA

基于scaled的数据进行PCA：

ov.pp.pca(adata,layer='scaled',n_pcs=50)
print(adata)

可以发现，adata.obsm层里多出了一个scaled|original|X_pca，这代表了我们使用的是layers中的scaled层数据进行的pca计算，当然我们也可以使用counts进行pca计算，效果如下：

ov.pp.pca(adata,layer='counts',n_pcs=50)
print(adata)

使用embedding函数，来对比基于两种不同的layers计算所得出的pca的差异：

import matplotlib.pyplot as plt
fig,axes=plt.subplots(1,2,figsize=(8,4))
ov.utils.embedding(adata,
                basis='scaled|original|X_pca',frameon='small',
                color='MS4A1',show=False,ax=axes[0])
ov.utils.embedding(adata,
                basis='counts|original|X_pca',frameon='small',
                color='MS4A1',show=False,ax=axes[1])
plt.savefig('./result/2-6_pca.png')

我们会发现基于scaled的pca结果，第一主成分和第二主成分有着相似的数量级，而基于counts的pca结果，第一主成分和第二主成分的数量级则有所差异，这对于后续的2维投影（比如t-sne和umap）可能会有显著的影响。

t-sne

t-SNE 是一种基于图的、非线性的降维技术，它将高维数据投影到 2D 或 3D 分量上。该方法基于数据点之间的高维欧几里得距离定义了一个高斯概率分布。随后，使用 t 分布在低维空间中重建概率分布，其中嵌入通过梯度下降进行优化。

分别对scaled的pca和counts的pca进行tsne：

sc.tl.tsne(adata, use_rep="scaled|original|X_pca")
# tsne函数默认是存放在adata.obsm['X_tsne']中的，我们将其存放在adata.obsm['X_tsne_scaled']中来区分counts的结果
adata.obsm['X_tsne_scaled']=adata.obsm['X_tsne']
sc.tl.tsne(adata, use_rep="counts|original|X_pca")
adata.obsm['X_tsne_counts']=adata.obsm['X_tsne']

# 可视化
import matplotlib.pyplot as plt
fig,axes=plt.subplots(1,2,figsize=(8,4))
ov.utils.embedding(adata,
                basis='X_tsne_scaled',frameon='small',
                color='MS4A1',show=False,ax=axes[0])
ov.utils.embedding(adata,
                basis='X_tsne_counts',frameon='small',
                color='MS4A1',show=False,ax=axes[1])
plt.savefig('./result/2-6_tsne.png')

UMAP

UMAP 是一种基于图的、非线性的降维技术，原理上与 t-SNE 类似。它构建了数据集的高维图表示，并优化低维图表示，使其在结构上尽可能地与原始图相似。

我们首先基于PCA的结果，在单细胞数据上构建一个邻域图再运行umap：

sc.pp.neighbors(adata, n_neighbors=15, n_pcs=50,
               use_rep='scaled|original|X_pca')
sc.tl.umap(adata)
# umap函数默认是存放在adata.obsm['X_umap']中的，我们将其存放在adata.obsm['X_umap_scaled']中来区分counts的结果
adata.obsm['X_umap_scaled']=adata.obsm['X_umap']

sc.pp.neighbors(adata, n_neighbors=15, n_pcs=50,
               use_rep='counts|original|X_pca')
sc.tl.umap(adata)
adata.obsm['X_umap_counts']=adata.obsm['X_umap']

# 可视化
import matplotlib.pyplot as plt
fig,axes=plt.subplots(1,2,figsize=(8,4))
ov.utils.embedding(adata,
                basis='X_umap_scaled',frameon='small',
                color='MS4A1',show=False,ax=axes[0])
ov.utils.embedding(adata,
                basis='X_umap_counts',frameon='small',
                color='MS4A1',show=False,ax=axes[1])
plt.savefig('./result/2-6_umap.png')

检查质量控制中的指标

现在我们可以在scaled数据的UMAP中检查之前的质量控制指标，一般检查三个指标：total_counts，n_genes_by_counts，pct_counts_mt。

如果前面使用的是omicverse.pp.qc，那么我们将直接得到nUMIs，detected_genes，mito_perc三个变量，如果使用的是scanpy进行的质控，那么得到的将是total_counts，n_genes_by_counts 和 pct_counts_mt 三个变量。

我们使用numpy中的log2对数化，将数据可视化区间缩小，同时，我们定义最大最小值来衡量我们的数据质量，我们希望total_counts大于250，250的对数值是7.96，所以我们最小值设定为8，最大值则定义为30,000，即15，而线粒体的比例则在0-100的范围内。

import numpy as np
adata.obs['log2_nUMIs']=np.log2(adata.obs['total_counts'])
adata.obs['log2_nGenes']=np.log2(adata.obs['n_genes_by_counts'])
ov.utils.embedding(adata,
                basis='X_umap_scaled',
                color=['log2_nUMIs','log2_nGenes','pct_counts_mt'],
                title=['log2#(nUMIs)','log2#(nGenes)','Mito_Perc'],
                vmin=[0,0,0],
                vmax=[15,15,100],show=False,frameon='small',)
plt.savefig('./result/2-6_qc.png')

理想情况如下图，total_counts（全部高亮），n_genes_by_counts（全部高亮），pct_counts_mt（全部不高亮）。

然后保存数据用于后续分析：

adata.write_h5ad('./data/s4d8_dimensionality_reduction.h5ad', compression='gzip')