基于短序列读段的转录本异构体推断

# 基于短序列读段的转录本异构体推断 ## 1. 算法核心步骤 在算法的第17和18步，构建了一个加权集合覆盖实例。其目的是找到一个有效转录本异构体的子集，该子集的总权重最小，并且能够解释所有的短读段和起始 - 结束片段。由于涉及的实例规模较小，这个集合覆盖问题可以使用GNU软件包GLPK中实现的分支限界法来求解。 ## 2. 模拟测试结果 ### 2.1 数据准备 - 从UCSC（mm9，NCBI Build 37）下载所有小鼠基因的参考基因组序列和已知异构体。此数据集包含26,989个基因和49,409个异构体。其中，16,392个（60.7%）基因只有一个异构体，59个（0.2%）基因有超过10个异构体；5830个（21.6%）基因只有一个外显子，384个（1.4%）基因有超过40个外显子 - 内含子边界。 - 仅选取至少有两个已知异构体的基因，共得到10,595个基因。 - 提取所有的起始 - 结束片段，并随机生成每个异构体的相对表达水平。考虑三种表达水平分布： - Base10：为每个异构体生成一个遵循标准正态分布的随机数r，然后将\(10^r\)作为该异构体的相对表达水平。 - Base2：为每个异构体生成一个遵循标准正态分布的随机数r，然后将\(2^r\)作为该异构体的相对表达水平。 - Uniform：为每个异构体从[0,1]中均匀生成一个随机数r，并将其作为该异构体的相对表达水平。 - 根据异构体的相对表达水平，随机生成40M单端短读段和10M双端短读段。模拟中假设双端读段的跨度是一个服从正态分布\(N(\mu, \sigma^2)\)的随机变量。 ### 2.2 性能评估指标 使用IsoInfer根据起始 - 结束片段以及单端和双端读段来恢复所有已知异构体。模拟中，单端和双端读段的长度分别为25bps和20bps，参数\(\alpha\)设为1 RPKM，\(\beta = 7\)，\(\gamma = 1000\)。考虑以下三个性能指标： - 灵敏度（Sensitivity）：恢复的异构体数量除以所有已知异构体的数量。 - 有效灵敏度（Effective Sensitivity）：恢复的异构体数量除以其外显子 - 内含子连接点得到读段数据支持的已知异构体数量。 - 精确率（Precision）：恢复的异构体数量除以推断出的异构体数量。 ### 2.3 表达水平计算 为了评估QP公式的有效性，在已知的小鼠异构体上随机生成Base10表达水平和单端短读段，并检查是否能恢复已知异构体的正确表达水平。对于表达水平为\(x_f\)和计算得到的表达水平为\(x_f'\)的异构体f，使用相对差异\(\frac{|x_f' - x_f|}{x_f}\)来衡量计算的准确性。 将该方法与简单方法（本文简记为Uniq）以及基于最大似然估计（MLE）和重要性采样（IS）的方法进行比较，结果如下表所示： |方法|1000万读段|8000万读段| | ---- | ---- | ---- | |MLE + IS|最准确|最准确| |Uniq|最差|最差| |IsoInfer|与MLE（后接IS）性能相当|与MLE（后接IS）性能相当| MLE + IS方法最准确，还能为每个表达水平估计提供95%的置信区间，但IsoInfer计算表达水平的速度比MLE + IS快得多（在标准台式PC上处理所有小鼠基因，IsoInfer需3分钟，而MLE + IS需3小时），其高效性使寻找新的异构体成为可能。 ### 2.4 表达水平分布的影响 固定双端读段跨度的分布为正态分布\(N(300, 30^2)\)，分析表达水平分布对IsoInfer推断异构体性能的影响。Iso