单细胞RNA-seq数据缺失元素补全算法

崔 璐, 刘桂锋

(1. 吉林大学 中日联谊医院医疗保险管理部, 长春 130033; 2. 吉林大学 中日联谊医院放射线科, 长春 130033)

单细胞DNA和RNA测序技术能帮助人们理解细胞功能和细胞异质性, 是生物学和医学领域的前沿技术[1]. 随着单细胞测序数据的迅速增长, 数据分析与计算方法面临新的挑战[2]. 由于生物实验捕获率较低, 导致单细胞(RNA-seq)数据存在较高比例的缺失值问题. 文献[3]研究表明, 在小鼠胚胎单细胞测序数据中缺失元素比例达30%. 如果直接删除高比例的缺失值, 则将导致重要信息的缺失.

为了解决上述问题并更好地实现细胞分型和功能分析, 文献[2]提出了一种两阶段机器学习识别估计方法, 该方法能有效识别出缺失值的位置并估计缺失值大小, 但该方法计算复杂度较高. 文献[4]将缺失值定义为假阴性误差, 并提出了一种基于Bayes方法的单细胞差异表达数据分析方法, 该方法需先假设缺失值满足的分布, 但找到缺失值的准确分布是一个难题[5]. 因此, 本文提出一种基于非负矩阵分解的单细胞RNA-seq数据缺失值补全算法. 该方法通过迭代求解法寻找最佳重构基因表达矩阵.

1 单细胞RNA-seq数据缺失模型

给定一个单细胞RNA-seq基因表达矩阵E∈n×m, 其中包含n个基因和m个细胞样本. 设该表达矩阵中缺失值集合为L, 0≤i

(i,j)=F(li,j，E),ei,j=0,

(1)

E′(i,j)=Q(L(i,j),E),

(2)

其中: (i,j)表示E中的任一缺失值;F(·)表示缺失值识别函数;Q(·) 表示缺失值补全函数;E′为补全后的单细胞RNA-seq基因表达矩阵. 上述模型与传统数据补全模型的不同: 单细胞RNA-seq基因表达矩阵中的缺失值位置和数量均未知, 因此, 对该矩阵进行补全, 传统缺失值补全方法已不再适用. 单细胞RNA-seq数据缺失模型需同时解决3个问题, 即矩阵中缺失值元素的位置、 数量和补全值.

2 迭代重构基因表达矩阵算法

为解决单细胞RNA-seq数据中的缺失值问题, 文献[2]利用式(1)和式(2)先后调用不同的机器学习模型求解, 较好地解决了缺失值问题, 但上述方法求解过程繁琐且计算复杂度较高. 本文将非负矩阵分解算法引入到单细胞RNA-seq数据补全问题中求解. 对给定的单细胞RNA-seq基因表达矩阵E, 求解非负矩阵W∈n×t和H∈t×m为min{En×m-Wn×t×Ht×m}. 因此, 该求解过程是一个迭代求解W和H的过程. 其中W和H在非负矩阵分解求解过程中的更新规则[6]为

(3)

求出W和H后, 二者的矩阵点积为非负矩阵分解算法得到的结果矩阵. 结果矩阵中会估计出原始基因表达矩阵中缺失元素的表达值.

经过非负矩阵分解算法补全后的重构基因表达矩阵计算方法如下：

(4)

其中, 缺失值由非负矩阵分解结果矩阵补全, 而非缺失值位置的基因表达值E(u,v)保留原始表达值.

经典非负矩阵分解求解前需结合人工经验输入rank值, 该值大小会直接影响结果矩阵的误差. 通过前期实验表明, 即使rank值发生很小的变化, 都会对结果产生较大影响. 针对该问题, 本文提出一种迭代重构基因表达矩阵算法, 算法步骤如下.

1) 算法输入: 输入单细胞RNA-seq基因表达矩阵E和初始误差阈值ε0;

2) 初始化参数: 根据数据设定rank值搜索范围(2～θ)(50≤θ≤m)及当前误差ε=ε0;

3) 迭代求解非负矩阵: 根据式(3)和当前rankt值, 迭代求解W和H;

4) 计算误差εt+1, 若εt+1<εt则转步骤3), 否则执行步骤5);

5) 输出结果矩阵: 根据式(4)计算E′和缺失值元素集合L;

6) 用t-SNE(t-distributed stochastic neighbor embedding)方法绘制细胞分型图谱.

3 实验和分析

为检验算法的有效性, 本文收集并整理了慢性粒细胞白血病单细胞测序数据. 该数据来源于美国国家生物信息中心(NCBI)的基因表达数据库, 数据检索号为GSE76312[7]. 该原始数据共包含2 287个干细胞和23 384个基因, 本文选择其中1 102个不含络氨酸激酶抑制剂的细胞, 并利用伪发现率和差异倍数选取前234个差异表达基因. 将该数据分为5类： 急变期慢性粒细胞白血病细胞(BC-CML)、 慢性期慢性髓性白血病细胞(CP-CML)、 人红白血病细胞系(k562)、 正常造血干细胞(normal)和前急变期慢性粒细胞白血病细胞(preBC).

图1为慢性粒细胞白血病单细胞测序原始数据分型与迭代重构基因表达矩阵数据细胞分析比较结果. 用t-SNE方法进行分型数据的可视化, 该方法是一种非线性降维和可视化方法[8], 目前已被广泛应用于肿瘤亚型分析[9]、 雷达目标识别[10]及人物动作识别[11]等领域, 并取得了较好的可视化效果. 由图1(A)可见, 该数据经过细胞和差异表达基因的过滤筛选, 呈现了数据的可分性, 但5种类型的细胞簇被划分的较分散, 且正常造血干细胞簇被分割开, 图1(A)中左下侧蓝色椭圆簇类中包含了4种不同类别的细胞. 该分类结果表明, 慢性期慢性髓性白血病细胞、 前急变期慢性粒细胞白血病细胞和正常造血干细胞被混淆在一起. 由图1(B)可见, 大部分正常造血干细胞聚集在右上角的蓝色椭圆簇类中, 且该簇类中的前急变期慢性粒细胞白血病细胞和慢性期慢性髓性白血病细胞都明显少于图1(A). 因此, 迭代重构基因表达矩阵数据细胞分型的结果更准确.

图1 原始数据与迭代重构基因表达矩阵数据细胞分型结果比较Fig.1 Comparison of cell typing results between original data and iterative reconstruction matrix of gene expression data

综上所述, 针对单细胞RNA-seq数据中存在缺失值的问题, 本文提出了一种单细胞RNA-seq数据缺失元素补全算法. 首先定义了单细胞RNA-seq基因表达矩阵数据缺失模型, 然后提出了迭代重构基因表达矩阵算法, 并将该算法应用到慢性粒细胞白血病单细胞测序原始数据中, 取得了更准确的细胞分型结果.