iLocator

此篇论文分析有部分不通透的地方可以从上一篇了解:A Framework for the Automated Analysis of Subcellular Patterns in Human Protein Atlas Images

论文材料

论文下载路径： An image-based multi-label human protein subcellular localization predictor (iLocator) reveals protein mislocalizations in cancer tissues
数据、代码、补充材料：http://www.csbio.sjtu.edu.cn/bioinf/iLocator/
作者：Ying-Ying Xu, Fan Yang, Yang Zhang, and Hong-Bin Shen
发表年份：2013
期刊： bioinformatics (IF:6.931，2022年中科院2区，TOP)

摘要及背景

  蛋白质需要在正确的时间出现在正确的亚细胞区室中发挥作用，才能使得人体机能正常运作。基于以下背景：(1) 蛋白质被错误的表达在错误的位置会引起病理性疾病，如肺癌，乳腺癌等等；(2) 一定数量的蛋白质会同时出现在多个亚细胞位置，如NSDHL同时出现在内质网和囊泡，然而用于多标签的蛋白质亚细胞位置预测的模型还需要进一步研究。因此，作者采用免疫组化图像(IHC)作为研究对象，提出了一个面向多标签多类别蛋白质亚细胞位置预测的系统，即iLocator，模型的框架如图1所示。主要工作如下：
   1. 在前人的Haralick的特征基础上，引入获取图像局部统计算子 局部二值模式(Local Binary Pattern, LBP) 描述IHC图像的局部特性；
   2. 面向多标签预测问题，作者引入链式分类器(Chain Classifier, CC) , 二元关联(Binary Relevance, BR) 处理多标签数据；
   3. 为了处理多标签数据预测输出类标决策问题，基于贝叶斯最小风险决策理论从训练集预测置信度和标签拟合出决策标签的阈值；
   4. 基于训练好的模型，应用于判定蛋白质是否能作为癌症组织中的生物标志物(biomarker)，由于病理图像没有标签，作者采用T检验的方式判决预测结果。

图1. iLocator的逻辑框图

研究方法

   如图1所示，iLocator分两个部分，首先以IHC为研究对象，采用线性光谱分离(LIN)方法解离出protein channel 和DNA channel，并用于获取Haralick, DNA和LBP特征；采用逐步判别分析(Stepwise Discriminative Analysis, SDA)用于特征筛选，并使用CC或者BR拟合数据；第二阶段将拟合好的模型用于预测癌症组织下蛋白质，并使用假设检验的方法认定生物标志物可靠。

数据收集及预处理

   作者收集了28个位于细胞质, 内质网, 高尔基体, 溶酶体, 线粒体, 细胞核, 囊泡等七个位置的3240张正常组织下染色水平为高的IHC图片，分辨率为3000*3000, 具体的每个蛋白质所属位置和数量可从补充材料Table 1中获取，其中多标签图片有827张占比25.52%，单标签累计2413张占比74.48%；表格中包括了实验中使用的癌症组织下蛋白质亚细胞的对应图片数量。为了避免数据中存在的其他染色水平下的图片，采用阈值13作为有效样本的筛选判据。然后，采用LIN的方法将IHC图片分解为Protein channel和DNA channel用于获取Haralick, DNA和LBP特征。
注：LIN的原理可以从这里了解

特征提取

   根据图1可知，蛋白质通道可以获取Haralick和LBP特征，蛋白质和DNA通道一起获取DNA分布特征。

Haralick特征

   Haralick特征是一种在灰度共生矩阵Gray Level Co-occurrence Matrix, GLCM）下基于小波变换的统计特征。首先将protein channel的灰度级别缩放到0-31(正常是0-255)；然后获取图像的GLCM，在灰度共生矩阵中可以获取13种，诸如惯量，相关性，熵等特征。图像可以获取水平，垂直，35°，135°四个方向的灰度共生矩阵，具体如图2案例所示。灰度共生矩阵是相邻像素$(i,j)$的数量统计矩阵，如水平方向上$(1,2)$相邻像素的数量即为1个，默认相邻的距离为1。因此，IHC在0-31的灰度级上的数学表示如下公式(1)所示[1]。在获取GLCM后计算相应特征，具体如[2]所述。
$$G=\left[\begin{array}{cccc}p(0,0)& p(0,1)& ...& p(0,Ng)\\ p(1,0)& p(1,1)& ...& p(1,Ng)\\ ...& ....& p(i,j)& ....\\ p(Ng,0)& p(Ng,1)& ...& p(Ng,Ng)\end{array}\right]\text{\hspace{1em}(1)}$$

这里的$Ng=31$，$p(i,j)$表示像素对$(i,j)$的数量。

图2. GLCM的一个示例

LBP特征

  局部二值模式(LBP)是一种获取中心像素与邻接像素二值化的统计信息。论文中的LBP如下图3，公式2和3所示。具体来说，在0-255灰度级的蛋白质通道中, $g_c$表示的是像素的灰度值，$(g_0,g_1,g_2,...,g_{U-1})$为$g_c$在半径为$R$下的$U$个领接像素值；如图中心像素值为$90$，对应的半径为1，邻接8个像素值为$[97,92,85,80,88,95,100,108]$； 然后进行公式2的计算，即从中心像素水平顺时针第一个像素开始进行大小比较，并返回逻辑值，结果为$[1,1,0,0,0,1,1,1]$； 接下来，将所示的逻辑结果转化为十进制数值；最后对所有的十进制数值进行直方图统计，0-255的频数即为LBP特征；因此LBP的特征维度为256。 论文中$U=8$, $R=1$。

$$s(g_k)=\{\begin{array}{rcl}1& & g_k\ge g_c\\ 0& & g_k<g_c\end{array}\text{\hspace{1em}(2)}$$
$$LB{P}_{U,R}=\sum _{k=1}^{U-1}s(g_k)2^k\text{\hspace{1em}(3)}$$

图3. LBP特征的图示

逐步判别分析(SDA)

  由特征提取，图片获得了$836_{Haralick}+4_{DNA}+256_{LBP}=1096$维度的特征。为避免特征冗杂，作者采用了逐步判别分析(SDA)进行特征筛选。SDA是一 种以特征在类内类间的差异，利用F检验(ANOVA)来判定特征是否能被选为判别性特征，如图4所示，数学关系如公式4-8所示。具体来说，原假设$H_0$：特征$m$不能够被选择，备择假设$H_1$：特征$m$可以被选择；在图4中当样本统计量的$P$值(蓝色区域的面积)落在临界值$\alpha$值(红色区域面积)右侧，则拒绝原假设$H_0$，即当前$m$特征不能被选为判别性特征；其中$\alpha$值是根据F分布的自由度查表获得，同时会对应一个临界值，论文中$\alpha =0.0001$。
  由于需要考虑的是特征在多个类别中的特异性，因而是一种多特征维度匹对多类别的假设检验问题，因而采用F检验；为了满足F检验中多数据符合中心极限定理下卡方分布的特点，采用大数据量下的类内方差样本和类间方差样本作为两种卡方分布，即公式6所示；由此即可计算得到F统计量样本下的拒绝域分布概率$P$, 意思是说只有$P$的概率大小才能发生原假设$H_0$, 因此当$P$小于0.0001时拒绝原假设$H_0$, 可以考虑选为判别性特征。实际上，根据公式4和5，利用F检验体现出的实际意义是指特征在类内和类间的关系，当$\Lambda (m)$越小，反映出的是类内差异小类间差异大的现象 [3] [4] [5]。

$$P=1-fcdf(F,df1,df2),F为F检验中的样本统计量,df1,df2为F检验的两个自由度\text{\hspace{1em}(4)}$$
$$F=(\frac{n-q-m}{q-1}(\frac{1-\Lambda (m+1)}{\Lambda (m+1)},n为特征数,q为类别数,m为当前特征数量\text{\hspace{1em}(5)}$$
$$\Lambda (m)=\frac{|W(X)|}{|T(X)|},X=[X_1,X_2,...,X_m],W是类内协方差矩阵，T是类间协方差矩阵\text{\hspace{1em}(6)}$$
$$W(i,j)=\sum _{g=1}^q\sum _{t=1}^{n^g}(X_{igt}-\overline{X_{ig}})(X_{jgt}-\overline{X_{jg}}),i,j\in 1,2,...,m\text{\hspace{1em}(7)}$$
$$T(i,j)=\sum _{g=1}^q\sum _{t=1}^{n^g}(X_{igt}-\overline{X_i})(X_{jgt}-\overline{X_j}),i,j\in 1,2,...,m\text{\hspace{1em}(8)}$$
$$X_{igt}是g类中t数据点的第i个特征,\overline{X_{ig}}是g类中第i个特征的均值,\overline{X_i}是所有特征的均值,n_g是g类的数据量。$$

图4. 假设检验的图示

多标签分类器

  对于实验中的蛋白质样本会同时出现在多个亚细胞区室中，是一种多标签多类别问题，因此普通的单标签的分类器并不适用。作者在实验中引入了链式分类器(chain classifier, CC) 和二元关联(binary relevance, BR)，并且实验结果表明多标签分类任务中CC的性能优于BR。关于两种分类器的介绍如图５所示。CC和BR的出发点是相似的，即基于one vs rest(OVR)思想，训练出与亚细胞区室类别数量相当的7个SVM分类器，训练过程中拟合其中一类为正样本其他为负样本，即每一个分类器次训练时都是一个二分类问题。区别在于，BR是单纯的遍历每一个类别，而CC从训练第二个类别开始将上一次样本的决策结果拼接到训练样本的特征中并作为下一个分类器的训练特征；测试时BR将7个分类器依次对测试样本进行预测获得7维的置信度，而CC会将上一次的预测的决策结果拼接到测试样本中进行下一次预测，最后也是7维置信度输出。 图中$x_i=[f_1^i,f_2^i,...,f_d^i]$表示训练样本$i$, 它的标签为$[1,0,1,0,0,0]$, $d$表示特征维度；$t_j=[f_1^j,f_2^j,...,f_d^j]$表示测试样本$j$, $S_1^j$表示第一个分类器对测试样本$j$的预测置信度。

图5. A, C分别是BR和CC的训练过程；B, D是BR和CC对一个样本的测试过程

决策层

  对于多标签多类别的置信度决策问题，作者基于贝叶斯最小风险判决准则拟合出一个单标签和多标签样本的判定阈值。基于所属类别置信度与非属类别置信度之间的分布差异在模型性能的驱导下会很明显，因此在大量数据和中心极限定理的假设下，可以利用高斯分布拟合出所属类别置信度与非属类别置信度的两种分布曲线，并获得满足最小风险下的决策阈值，这种风险是指图中的阴影部分，具体如图6，公式9和10所示。

d i f = m a x { s 1 , s 2 , . . . , s 7 } − s c (9) dif=max \lbrace s_1, s_2, ..., s_7 \rbrace -s_c \tag{9} dif=max{s1​,s2​,...,s7​}−sc​(9)
$$\begin{array}{rl}b& =\underset{ϵ=1,2}{\mathrm{argmax}}P(H_{ϵ}|dif)\\ & =\underset{ϵ=1,2}{\mathrm{argmax}}P(dif|H_{ϵ})\ast P(H_{ϵ})/P(dif)\\ & =\underset{ϵ=1,2}{\mathrm{argmax}}P(dif|H_{ϵ})\ast P(H_{ϵ})\end{array}\text{\hspace{1em}(10)}$$
   $s_c$是样本在$7$个分类器下的置信度；$dif$是最大置信度与每个置信度的差值；$H_1$：在$dif\ge T$时，为所属类别，$H_2$：在$dif<T$时,不是所属类别；$P(H_{ϵ}|dif)$表示在差值为dif的情况下发生$H_{ϵ}$的概率, $b$表示在两种情况下的最大值；根据贝叶斯公式可以推导公式10的结果，其中$P(dif)=1$；$P(H_{ϵ})$是$H_1$和$H_2$的后验概率，可以根据统计结果利用频数和经典概率公式计算获取，物理意义即是判定为单和多标签的置信度的数量的比例。因此，计算$P(dif|H_{ϵ})$很重要。

图6. 阈值T计算过程，右下图为LIN下db1特征训练结果的拟合过程

   计算$P(dif|H_{ϵ})$分了如图6所示的B，C，D， E四个过程。$di{f}_1$是最大值置信度与非属类别置信度最大值的差值的统计，$di{f}_2$为最大值置信度与所属类别置信度之间的差值统计；B，C图表示在两种差值下的直方图统计；D图根据B和C图直方图统计拟合的高斯曲线，由此，即可计算在$P(dif|H_{ϵ})$, 比如，$P(dif|H_1)$是指在所属类别下发生差值$dif$的概率。最后，与后验概率相乘获得E图，在此分布下，两曲线交点的值即为所求的判别是否为所属类别的阈值。论文中的拟合结果$T=0.8897$, 当然因为不同的数据T应该是不一样的，所以如果更换数据集T也是会有相应的变化。

实验结果

正常组织

  正常组织下28个蛋白质所属7个位置的3240张图片，按照训练集测试集1：1的比例设置，获得了如图7所示的结果。从图中可以获得以下信息：(1) IHC在加入LBP后提高了特征的描述能力; (2) 两种分类器中，CC对多标签数据有更好的结果; (3) 模型对单标签的数据拟合能力可以有较大的提升，单引入多标签数据结果下降了很多, 说明模型还是对单标签的数据更敏感，可能的原因在于单标签的数据量占比更大。

图7. 左图为两种分类器在不同特征下的结果，右图为db8特征下的单标签和多标签数据下的结果

癌症组织

   为了判断28个蛋白质哪些可以用于生物标记物，作者根据相同蛋白质收集了3696张IHC图片，其中多标签图片含有912张占比24.68%。由于病理组织的IHC图片不带有标签，因此不能采用标签进行对比结果，作者在这里采用了T检验的方法来判定是否接受预测结果。被认为是标记物的条件：1. 输出置信度不一样；2. 正常和癌症的对应位置的预测分数是符号相反的；除此之外还满足$P_{-values}<0.01$。同假设检验原理，由于是相同蛋白质下正常组织和癌症组织两种状态下的检验，输出置信度的均值接近的，但由于预测量的方差未知，且样本少于30，因此采用双总体t检验。其中原假设$H_0$：两个蛋白质来自同一个分布，即所谓标签可认定一致，不能认定为该位置的生物标记去；备择假设$H_1$：两个蛋白质来自不同分布，即所谓标签不一致，可以被认定为该位置的生物标记物[6]。
$$S_g^{normal/cancer}=\frac{1}{n}\sum _{i=1}^n(s_{ig}-{\bar{s}}_g{)}^2\text{\hspace{1em}(11)}$$
$$t=\frac{({\bar{s}}_g^{normal}-{\bar{s}}_g^{cancer})-(u_g^{normal}-u_g^{cancer})}{\sqrt{\frac{S_{norma{l}_g}^2}{N}+\frac{S_{cance{r}_g}^2}{M}}}\text{\hspace{1em}(12)}$$
   $s_{ig}$是第$g$类置信度下图片i的置信度分数，${\bar{s}}_g$是$g$类的置信度均值,$N$表示采样的正常组织下的图片数量，$M$表示采样的癌症组织下的图片数量。即先计算蛋白质组织的IHC样本的置信度的类别均值，总体均值和样本标准差，公式11；然后计算每个蛋白质组织的T值，公式12；并采用T分布计算获得$P_{-values}$；当$P_{-values}<0.01$时，拒绝原假设$H_0$，即认为该蛋白在正常组织和癌症组织下不处于同一个分布，不属于相同的位置；综合贝叶斯最小风险判决获取的阈值决策结果，可以获取癌症组织下蛋白质的位置；由于只有蛋白质位置发生变化才能作为生物标志物，因此认为该情况下可认为蛋白质为生物标志物。比如，CPT2蛋白质，正常组织下该蛋白质位于高尔基体，经过模型预判定癌症组织下位于细胞质($P_{value}=0.0013<0.01$)和线粒体($P_{value}=0.1741>0.01$), 因此CPT2在可以作为肺腺癌的生物标记物。

图8. T检验被用于检验蛋白质能被接受认定为生物标记物

总结与讨论

1. 引入LBP但中心像素的特征未被考虑，因此也可以考虑Local Ternary Patterns(LTP) 或者 Local Quinary Patterns(LQP)。
2. 特征选择层面SDA的方法可以高效率的获得判别性的特征，但在不同数据下可以具体讨论，其中遗传算法也可以使用, 在[5]中性能比SDA略微好一些。
3. CC比BR性能更好，更新了下一步的研究方向，可以考虑标签之间的潜在关联性，挖掘样本在多标签之间的差异，进一步探索多标签预测器及置信度的决策方法
4. 对于蛋白质亚细胞位置预测，还存在大量的序列数据，因此可以考虑图片与序列相结合的研究。

[1] https://murphylab.web.cmu.edu/publications/boland/boland_node26.html
[2] https://blog.csdn.net/weixin_49171484/article/details/126261788
[3] https://murphylab.web.cmu.edu/publications/boland/boland_node27.html
[4] https://www.bilibili.com/video/BV1dU4y147Sx?
[5] Huang K, Velliste M, Murphy R F. Feature reduction for improved recognition of subcellular location patterns in fluorescence microscope images[C]//Manipulation and Analysis of Biomolecules, Cells, and Tissues. SPIE, 2003, 4962: 307-318.
[6] https://www.bilibili.com/video/BV1Qv411W7GQ?