跨物种比较转录组数据的处理

Posted on | Post modified | In

该流程参考自Perry 等人（2012）。跨物种的全长转录本存在非同源插入/缺失（Indels）和不保守剪接，直接比对会引入物种特异性的长度偏差，掩盖真实的表达量演化信号。必须通过多序列比对强制裁剪掉所有存在序列变异或缺失的区域，仅保留跨物种100%保守的“最大直系同源区域”。

先看看作者是如何处理数据的

基于 Perry 等人（2012）和他们补充材料（Supplemental Methods）中提供的详尽细节，处理流程为：

定制化 de novo 组装与直系同源基因识别（逐物种进行）

将同一物种所有个体的 reads 合并，进行基于人类参考基因组引导的定制化组装：

  • 构建 de Bruijn 图： 使用 $k=39$ 构建图。使用 Bloom Filter 过滤掉仅出现一次的 k-mers（覆盖度为1），并去除代表测序错误的短末端 (tips)。

  • 基于种子的同源性搜索： 将人类 RefSeq 基因的所有亚型合并为单一共有序列。在图中搜索与人类基因匹配的短种子序列。种子长度和允许的错配数量根据物种与人类的系统发育距离严格设定（例如：黑猩猩种子长度27，允许2个错配；小鼠种子长度15，允许4个错配）。

  • 局部图过滤： 提取上一步识别出的同源 k-mers 的覆盖度，取其第 90 百分位数作为预期覆盖度。从图中删除覆盖度低于该预期值 10%（且绝对覆盖度不小于3）的节点，以此去除内含子序列和测序错误。

  • 路径分析与共有序列提取： 提取图中的连通分量，并去除长度差异不超过 10 bp 的短泡 (bubbles)。将保留的长度至少达到人类参考序列编码区 50% 的路径，使用 Fast Statistical Alignment (FSA) 算法与人类参考序列比对，选择匹配碱基最多的路径作为该基因的共有序列。

  • 旁系同源基因过滤： 将组装出的共有序列使用 BLAST 反向比对回人类 RefSeq 基因库，如果最佳匹配不是原始的目标基因，则将其作为旁系同源基因剔除。

Reads 映射与表达量初始计数（逐样本单独进行）

  • 严格唯一比对： 使用 BWA（默认参数）将每个样本的 reads 比对到其所属物种组装好的基因序列库上。双端 reads 被拆开独立比对，并且仅保留唯一比对 (uniquely mapped) 的 reads。

  • 跨剪接位点比对： 对于第一步未比对上的 reads，取其首尾 20 bp 独立进行 BWA 比对。如果两端比对到同一转录本且能向内部完美延伸，则计为有效的外显子连接 (exon-junction) read。

跨物种定量的核心：最大直系同源区域限制

为了防止测序偏差或部分物种特有的可变剪接外显子影响表达量估计，作者在计算最终表达量时实施了极其严格的区域截断：

  • 使用 FSA 对所有物种组装出的同源基因序列执行多序列比对。

  • 识别在所有物种中完全对齐的最大直系同源区域 (maximum orthologous region)，并从中强制排除任何非编码区 (UTRs)。

  • 仅仅统计物理位置落入这个保守同源区域内的 reads（包括上述的剪接 reads），作为该转录本的表达量衡量指标。

多维数据标准化

针对上述落在直系同源区域内的 reads 计数，必须进行连续的两步数学调整：

  • 种内调整（长度与 GC 含量）：

    • 首先除以转录本的长度。

    • 为了校正 PCR 扩增偏好，计算该转录本 reads 在个体中的比例 $y_{ij}$，并针对转录本的 GC 含量 $g_{j}$ 进行 loess 回归。利用拟合值 $f_{ij}$ 对表达量进行调整：$x_{ij} = (mean_{j}(f_{ij}) / f_{ij}) \times Z_{ij}$。

  • 跨物种调整（测序深度校正）：

    • 计算每个基因跨所有物种的平均表达值，将这些值排序形成一个合成分布，并计算其中位数 (medSyn)。

    • 对于每个样本，计算归一化因子：NormSamp = medSyn / (样本中位数)。

    • 使用该因子乘以该样本所有组装转录本的表达量，迫使每个样本的表达量中位数都对齐到 medSyn。

该标准化思路可否应用于可变剪接？

1。鉴定直系同源外显子后，对每个外显子量化跳过率； 

2。鉴定直系同源外显子后，构建他们之间的三联体，量化剪接事件的psi； 

3。仿照perry，严格提取直系同源无gap连续区域，在这个区域map并量化基因的剪接发生率； 

4。尝试组装转录本（但需要严格控制假阳性），再鉴定直系同源转录本，量化标准化表达量。

四条思路的分子生物学与进化意义阐明

思路 1：量化直系同源单外显子的跳过率 (Single Exon Skip Rate)

这正是 Perry 等人实际执行的分析逻辑。

  • 分子生物学意义： 孤立地评估一段特定编码序列（CDS）在转录本中被保留的绝对频率。它不关心这段序列的上游和下游连接了什么，只测量该外显子本身的“剪接亲和力”。

  • 进化意义（功能域的得失）： 它用于检验特定的蛋白质功能结构域是否受到了选择。例如，如果发现某个特定的外显子在食树胶的懒猴（Slow loris）中跳过率极低（几乎全部纳入），而在其他灵长类中跳过率极高，这提示该外显子编码的结构域可能在解毒或代谢树胶产物（如视黄醇代谢基因 SDR16C5）中发挥了适应性作用。

  • 局限性： 缺乏结构上下文，若物种间侧翼序列发生重组，相同的跳过率可能产生截然不同的蛋白质产物。

思路 2：构建直系同源三联体量化 PSI (Exon Triplet PSI)

  • 分子生物学意义： 将分析对象从“物理序列”升级为“剪接事件”。三联体严格控制了剪接体面临的竞争性环境（固定的上下游外显子）。PSI 测定的是剪接体在 E1​−E2​−E3​ 这一固定微环境中，做出特定连接选择的精确概率。

  • 进化意义（顺式调控元件的微调）： 它专门用于检测剪接调控网络的进化。PSI 的改变通常不源于蛋白质功能域的直接选择，而是源于cis-regulatory elements如剪接位点（SS）强度或剪接增强子/沉默子（ESE/ESS）的碱基突变。通过 EVE 模型检验三联体 PSI，可以回答：物种间顺式调控元件的突变是中性漂变，还是受到了定向选择的驱动。

思路 3：提取直系同源无 Gap 连续区域量化基因级剪接率 (Gene-level AS in Gapless Regions)

这是将 Perry 定量总表达量的逻辑 创新性地平移到 AS 领域。

  • 分子生物学意义： 忽略具体的异构体形式，直接把基因的保守核心区看作一个整体，测定该区域内产生的所有非标准剪接（包括跳跃、替代位点）的总体频率。它反映的是该基因转录过程中的“总剪接变异度”。

  • 进化意义（剪接保真度与遗传漂变）： 这种基因级别的总体变异度，通常用于评估“剪接噪声（Splicing Noise）”。Perry 等人发现克氏原狐猴（Coquerel’s sifaka）具有极高的遗传多样性，而指猴（Aye-aye）的遗传多样性极低，这与有效种群大小有关。思路3可以用来检验假说：在有效种群较小、纯化选择效率较低的物种中，基因级别的总剪接噪声是否由于遗传漂变（剪接保真度下降）而显著升高。

思路 4：组装并鉴定直系同源转录本量化表达量 (Orthologous Isoform Expression)

  • 分子生物学意义： 直接量化细胞内最终产生的功能产物（特定亚型的蛋白质产物前体）。这是分辨率最高的测量，直接关联到最终的表型输出。

  • 进化意义（完整功能模块的演化）： 它试图评估物种是否进化出了完全不同的整体调控策略。然而，正如 Perry 等人在对比16个物种后得出的结论：整个外显子结构的固定种间差异极其罕见，绝大多数变异是种内的多态性。

  • 局限性（假阳性极高）： 跨越远缘物种时，组装算法在遇到重复序列和测序深度不足时极易产生嵌合转录本。要求转录本全长保持直系同源且结构完全对应，在客观数据层面上几乎无法实现，强行比对会引入巨大的系统性技术偏差。

总结

如果研究目标是使用 EVE 模型来检验系统发育树上的适应性演化：

  • 思路 2（三联体 PSI） 提供的数据最符合系统发育方差分析的数学假设，且生物学意义（调控机制演化）最为精准排他。

  • 思路 3（核心区总剪接率） 适合用于探讨大尺度进化中的“剪接噪声”与选择压力强度的关系。

  • 思路 4（全转录本） 难以控制组装的假阳性，且演化尺度过大不一定能找到所谓的同源转录本。