1.基于python的单细胞数据预处理-特征选择

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

特征选择背景

现在已经获得了经过归一化的测序数据，其保留了细胞异质性，同时削弱了测量误差。统计发现，一个细胞表达的基因大约是3000个左右。这意味着测序数据中的一大部分基因是0计数。对于细胞亚型的研究，大部分0计数基因都在这些细胞亚型中，因此，预处理还包含特征选择，可以排除这些不具备分析意义的基因。

基因特征选择一般有三种方法：基于基因离散度，基于基因归一化方差，基于基因的皮尔森残差。

基于基因离散度

在传统的分析流程中，我们会采用基于基因离散度的方式去计算高变基因，一般来说，我们首先确定了单细胞数据集中变异最大的一组基因。我们计算了所有单细胞中每个基因的平均值和离散度（方差/平均值），并根据基因的平均值将基因分为 20 个箱（bins）。然后，在每个箱内，我们对箱内所有基因的离散度进行z归一化，以识别表达值高度可变的基因。

我们使用移位对数归一化后的数据：

import omicverse as ov
import scanpy as sc

ov.utils.ov_plot_set()

adata = sc.read("./data/s4d8_quality_control.h5ad")

#存储原始数据以便后续还原
ov.utils.store_layers(adata,layers='counts')
adata.layers['counts'] = adata.X.copy()

sc.pp.normalize_total(adata)
sc.pp.log1p(adata)
print(adata)

调用scanpy包里的pp.highly_variable_genes函数来计算高可变基因，由于我们使用的是基于基因离散度的方法，故设置flavor='seurat'，该方法也是默认方法。基于基因离散度的方法寻找高变基因有两个方式：

• 指定HVG数量，应用广泛，简单直接。
• 指定离散度，数据敏感，应用其实很少，还是推荐指定HVG数量。

对于指定HVG数量：

adata_dis_num=sc.pp.highly_variable_genes(
    adata,
    flavor="seurat",
    n_top_genes=2000,
    subset=False,
    inplace=False,
)
print(adata)
print(adata_dis_num)
print(adata_dis_num['highly_variable'].value_counts())

设置inplace=False，将不会改变adata的var（打印adata的视图时，var中没有出现highly_variable）。输出为：

我们发现，一共选择了2000个高可变基因，这与我们最开始的分析目标一致。

基于基因归一化方差

在seurat v3中，提出了基于基因归一化方差做特征选择，我们不再使用归一化后的数据来计算高变基因。我们首先计算每一个基因的平均值${\overline{x}}_i$与方差${\sigma }_i$，然后分别对平均值与方差进行log对数变换，然后用2次多项式，将方差作为均值的函数，进行多项式回归：$$\sigma (x)=a{x}^2+bx+c$$通过这个公式，可以获得每一个基因的预测方差，然后进行z变换：$$z_{ij}=\frac{x_{ij}-{\overline{x}}_i}{\sigma (x_i)}$$其中，$z_{ij}$是细胞$j$中基因$i$的归一化值，$x_{ij}$是细胞$j$中基因$i$的原始值，${\overline{x}}_i$是所有细胞基因$i$的平均原始值，$\sigma (x_i)$是预测的方差。对于特征选择，根据预测的方差进行排序即可。

在scanpy中，需要flavor='seurat_v3'，并指定计数矩阵是没有归一化的layer='counts'：

adata_var_num=sc.pp.highly_variable_genes(
    adata,
    flavor="seurat_v3",
    layer='counts',
    n_top_genes=2000,
    subset=False,
    inplace=False,
)
print(adata_var_num['highly_variable'].value_counts())

基于基因皮尔森近似残差

基于皮尔森近似的方法也是使用原始计数：

adata_pearson_num=sc.experimental.pp.highly_variable_genes(
    adata, 
    flavor="pearson_residuals",
    layer='counts',
    n_top_genes=2000,
    subset=False,
    inplace=False,
)
print(adata_pearson_num['highly_variable'].value_counts())

特征选择总结

对比三种不同的方法：

import matplotlib.pyplot as plt
from matplotlib_venn import venn3

adata_dis_num.index=adata.var_names.copy()
adata_var_num.index=adata.var_names.copy()
adata_pearson_num.index=adata.var_names.copy()

# 三个列表的元素
list1 = set(adata_dis_num.loc[adata_dis_num['highly_variable']==True].index.tolist())
list2 = set(adata_var_num.loc[adata_var_num['highly_variable']==True].index.tolist())
list3 = set(adata_pearson_num.loc[adata_pearson_num['highly_variable']==True].index.tolist())

# 绘制 Venn 图
venn = venn3([list1, list2, list3], set_labels=('Dis', 'Var', 'Pearson'))

# 显示图形
plt.title("Venn Diagram of Three HVGs")
plt.savefig("./result/2-5.png")

发现三种不同方法所找到的高可变基因（HVGs）仅有656个是相同的，这意味着不同的方法所寻找到的高可变基因会影响下游分析的结果一致性。如果对时间要求不严格，推荐使用皮尔森残差法来获得高可变基因。如果需要快速，推荐基于基因离散度的方法。

在omicverse中，归一化和特征选择预处理被包装好了，mode参数为normalize|HVGs，前者是归一化，后者是特征选择：

adata = sc.read("./data/s4d8_quality_control.h5ad")
#存储原始数据以便后续还原
ov.utils.store_layers(adata,layers='counts')
adata.layers['counts']=adata.X.copy()

adata=ov.pp.preprocess(adata,
                       mode='shiftlog|pearson',
                       n_HVGs=2000)
print(adata)

# 存储预处理后的数据
adata.write_h5ad('./data/s4d8_preprocess.h5ad')

在结果上，注意：与scanpy不同，omicverse计算高可变基因后，将保存为var['highly_variable_features']，而在scanpy中，HVG将保存为var['highly_variable']，都是包含bool值的Series。