#18: Py-boost预测t分数
你是否注意到，在这次竞赛中很少有真正的EDA（探索性数据分析）笔记本被发布？除了分享我的机器学习模型外，我还想提供一些有助于理解数据和Limma复杂性的观察结果。
生物学知识的整合
不要相信细胞类型！
让我们以简化的形式回顾实验过程。可以想象一个实验者站在一大锅人血细胞前。这口锅含有六种细胞类型的混合物，比例各异。CD4+ T细胞占比最大（42%），而调节性T细胞仅占2%：
实验者从大锅中取出145滴液体。每滴含有1550±240个细胞（正态分布）。如果计算每滴中每种细胞类型的数量，我们会看到多项式分布。145滴液体可能像下图所示（虚构数据，按液滴大小从小到大排序）：
接下来，实验者向145滴液体中添加145种物质，并等待24小时。24小时后分析细胞。如果再次计数细胞，我们会从竞赛的训练数据中得到如下结果（测试数据的细胞计数被隐藏）：
在这个图表中，我们首先看到一些化合物毒性很强，导致某些液滴中存活的细胞少于100个。这些液滴由条形图中最左侧的条表示。
第二个观察结果更为重要：Oprozomib和IN1451的条形图中长长的红色部分表明，这些液滴含有数百个调节性T细胞——比实验开始时多得多。其他化合物（如CGM-079）则有过多CD8+ T细胞（绿色条）。我们如何解释这一观察结果？
- IN1451会刺激调节性T细胞繁殖，使24小时后数量增加五倍？不会。
- IN1451会神奇地将NK细胞转化为调节性T细胞？不会。
- IN1451是否会影响细胞，导致它们被错误分类？有可能。
如果细胞在实验过程中改变类型，那么讨论特定细胞类型的差异基因表达就变得毫无意义。对于Kaggle竞赛而言，这意味着我们必须处理许多异常值：除了至少五种有毒化合物外，还有至少七种化合物会改变细胞类型。这些异常值的差异表达难以建模，它们使交叉验证不可靠，而且私人排行榜中的异常值甚至无法通过探测公共排行榜来预测。
细胞计数不应影响差异基因表达
细胞中的基因表达是否取决于实验中的细胞数量？理论上，不会。无论实验中有10个细胞还是10000个细胞，细胞的行为方式都相同。然而，我们期望实验结果的重要性会有所不同：拥有10000个细胞的实验应该比10个细胞的实验提供更精确的测量：随着细胞数量增加，测量的方差应减小，t分数应离零更远，p值应降低。
竞赛数据并不符合这一预期。如果我们绘制602种细胞类型-化合物组合（不包括对照化合物）的平均t分数与细胞计数的对比图，会看到线性关系：对于每种细胞类型，细胞计数较低的化合物具有正的平均t分数，而细胞计数较高的化合物具有负的平均t分数。这种细胞计数与t分数之间的相关性本不该存在，它是Limma的伪影，而非生物学效应。
你可以用中位数或方差代替平均数绘制图表——结果看起来相似。你甚至可以将细胞计数与t分数的第一个主成分进行比较，看到相同的相关性。
现在我们可以列出20种因低细胞计数而被视为异常值的化合物。注意，我们不是声明数据集的单个行为异常值，而是声明与这20种化合物相关的所有86行为异常值：
异常值列表 -------- AT13387 只有7个调节性T细胞 Alvocidib 几种细胞类型≤10个 BAY 61-3606 CD8+细胞平均t分数>2 BMS-387032 只有10个CD8+T细胞，无髓系细胞 Belinostat 细胞过多的对照化合物 CEP-18770 (Delanzomib) 几种细胞类型≤10个 CGM-097 CD8+T细胞过多 CGP 60474 几种细胞类型≤10个 Dabrafenib 细胞过多的对照化合物 Ganetespib (STA-9090) 只有4个调节性T细胞，NK细胞过多 I-BET151 CD8+T细胞过多 IN1451 几种细胞类型≤10个 LY2090314 只有6个CD8+T细胞 MLN 2238 几种细胞类型≤10个 Oprozomib (ONX 0912) 几种细胞类型≤10个 Proscillaridin A;Proscillaridin-A 几种细胞类型≤10个 Resminostat 无CD8+T细胞 Scriptaid 只有2个调节性T细胞 UNII-BXU45ZH6LI 只有6个CD8+T细胞 Vorinostat 只有1个调节性T细胞
移除异常值后，图表看起来更清晰了。细胞计数的方差仍然存在，它是干扰差异表达（和预测）正确解释的噪声源。也许如果在文库大小归一化之前均衡细胞计数，我们会得到更清晰的数据。但这相当于丢弃部分测量数据，同样不可取。
考虑到数据集规模小、噪声多以及Limma伪影（下一节将展示更多），我没有尝试将任何外部生物学数据整合到我的模型中。
问题探索
概率分布的混合
训练数据单行（用Scriptaid处理的CD8+ T细胞的18211个t分数）的直方图显示分布是多峰的。
最高模式包含269个基因，这些基因都具有-3.769的相同t分数。事实证明，这些是在CD8+ T细胞中从未表达过的269个基因，无论是阴性对照还是任何其他化合物处理下。这不是很奇怪吗？在整个实验中从未表达的基因应该具有零log倍数变化，根本不应该获得t分数（因为t分数计算涉及除以方差，而从未表达的基因的方差为零）。
对于用Foretinib处理的髓系细胞，3856个基因未表达（RNA计数为零），但其中大多数具有正的t分数。它们的最高t分数为6.228（导致p值为4e-10，log10p值为9.33）。如果RNA计数为零，对应的log倍数变化（和t分数）应该永远不会是正的。
我们可以说值的分布是两个分布的混合：
- 已表达基因（蓝色）的值具有或多或少的钟形分布。
- 未表达基因（橙色）的值具有不寻常形状的分布，更奇怪的是，当没有计数到任何RNA片段时，却报告了正的差异表达。
我们在这里看到的是Limma的伪影，它影响数据集的每一行。这表明Limma输出可能存在偏差，对于研究细胞类型翻译的差异表达并非理想选择。
在CD8+ T细胞中Scriptaid表达的基因： 13560 在CD8+ T细胞中Scriptaid未表达的基因： 4651 模式：-3.804，对应269个在CD8+ T细胞中从未表达的基因
在髓系细胞中Foretinib表达的基因： 14355 在髓系细胞中Foretinib未表达的基因： 3856 模式：6.379，对应81个在髓系细胞中从未表达的基因
理想的训练集
在竞赛概述中，组织者提问：你是否有证据表明，除了随机抽样细胞类型中的化合物外，如何开发用于细胞类型翻译的理想训练集？在保留细胞类型中测量的化合物数量与模型性能之间有什么关系？
我认为我们还没有准备好回答这些问题。我们首先需要更干净的数据（以及更多数据）：
- 必须正确分类细胞类型。这可能意味着我们将工作范围限制在不损害细胞类型分类的化合物上。
- 必须从数据集中剔除细胞太少的样本。这些样本只是在我们想找到针的地方增加了干草。
- 即使我们有许多细胞，也可能需要剔除低RNA计数的基因。否则它们会向干草堆中添加更多干草。
其次，对从未表达的基因建模奇怪的t分数是浪费时间。我们需要定义一个机器学习任务和指标，奖励生物学洞察，而不是迫使人们建模上游处理步骤产生的噪声：
- 由于t分数总是受细胞计数和方差估计影响，基于处理程度较低的数据（即log倍数变化或RNA计数而非log10p值或t分数）的指标可能会将研究引向更好的方向。
- 即使使用log倍数变化，低RNA计数的基因也比高RNA计数的基因产生更多噪声。合适的指标应考虑到这一事实。
模型设计
t分数优于log10p值
Limma执行t检验。t分数（几乎）呈正态分布，这对机器学习输入很有利。本次竞赛中，t分数被非线性转换为log10p值。这种转换将漂亮的钟形压缩成峰度更高的分布。
我的机器学习模型在预处理步骤中将log10p值转换为t分数后表现更好，预测t分数，然后再将预测结果转换回去。也许使用log倍数变化或RNA计数会更好。
模型
我开发了四种模型：
- Py-boost
- 基于岭回归的推荐系统
- 基于k近邻的推荐系统
- ExtraTrees
我首先实现了Py-boost模型，源自@alexandervc的公开笔记本。
然后我实现了ExtraTrees模型，它与@alexandervc的Py-boost模型相似。所有决策树都完全生长（即过拟合）。该模型通过故意添加到目标编码特征中的噪声获得泛化能力。
为了实现多样性，我接着实现了knn推荐系统。细胞类型和化合物被分别视为用户和项目；基因表达被视为用户对项目的评分。
ExtraTrees和k近邻共享无法外推的弱点。即使经过降维，我们的训练数据集本质上也是高维空间中的614个点，因此大多数点位于凸包上。在255个测试点中，许多点将位于训练点凸包之外，这意味着模型必须进行外推。为了使外推能力发挥作用，我实现了岭回归模型。
模型的cv分数在0.878（ExtraTrees）到0.906（Py-boost）之间。Py-boost虽然是交叉验证中最差的，但拥有最佳的公共和私人lb分数（分别为0.572和0.748）。
数据增强
一个模型（k近邻）使用了数据增强：如果我们知道两种化合物的差异表达，可以假设这两种化合物的混合物将产生单个化合物差异表达的平均值。
我还尝试了另一种数据增强：由于CD4+ T细胞的数量是髓系或B细胞的两倍多，而且我知道细胞计数会使Limma结果产生偏差，因此我减少了CD4+ T细胞的计数，进行伪批量处理，通过Limma运行，并将结果作为另一种细胞类型添加到训练数据中。这种增强提高了ExtraTrees的分数，但没有达到Py-boost的水平。也许我应该将额外的细胞类型与Py-boost结合起来...
鲁棒性
我的模型的鲁棒性通过两种方式证明：
(1) 模型完全交叉验证。交叉验证策略首先在SCP快速入门中记录，确保模型在预测cell_type-sm_name组合时进行验证，因此它只知道相同细胞类型的另外17种化合物。这种交叉验证策略比普通的分层KFold更鲁棒，后者中模型知道相同细胞类型的4/5化合物。（而且比普通训练-测试分割鲁棒得多。）
不过，我必须承认，我对cv-lb对应关系不满意。
(2) 所有模型的性能在输入t分数添加高斯噪声后进行了测试。所有模型对小幅噪声都具有鲁棒性。当噪声更强时，knn和ExtraTrees模型比Py-boost和岭回归推荐系统受损更严重。
文档和代码风格
代码在笔记本中有文档记录。
可重复性
源代码在此：
结论
让我总结一下这篇文章的四个主要信息：
- 推荐系统是开发跨细胞类型差异基因表达预测模型的一个有希望的起点。由于商业利益，推荐系统是一个被深入研究的话题，并且有大量可用信息。
- 数据增强是有用的，化合物混合物是数据增强的自然方法。
- 尽管具有数据清理、异常值移除和异常指标的Kaggle竞赛很有趣，但研究目标会从另一种设置中受益。提供干净的数据并使用易于理解的指标进行评分，将有助于参与者专注于真正的话题，而不是数据中的噪声。
- 我们已经看到，Limma在某些情况下会产生有偏见的输出。我希望专业的Limma用户意识到这些影响，并在研究中解释结果时加以考虑。