LeafCutter 实跑，删掉一条低支持连接，为什么其他内含子的比例也变了

剪接分析中，过滤低支持连接通常被理解为删掉几条不可靠记录。但如果这条连接恰好把两个局部结构连在一起，删除它还可能拆开剪接簇，进而改变其余内含子的比例分母。此时某个比例明显上升，并不意味着对应的原始支持增加了。

本篇实际安装经典 LeafCutter 项目的聚类脚本及 RegTools，构造四份可手算的比对输入，运行连接提取和三组聚类。我们从原始计数追踪到簇结构与分母变化，解释为什么过滤规则也是分析对象定义的一部分。

先确认安装的是哪一种实现

软件名称相同，不代表包内容、入口和统计实现完全相同。本次检查时，Bioconda 的同名包已经指向新的 Python 项目，而经典文档仍介绍聚类脚本及配套 R 分析。安装前应先核对包来源，不能只比较名称。

本实验选择经典项目的一个固定提交，下载并校验聚类脚本摘要，在 Python 三点十点十二环境执行。实际没有安装或测试新的同名实现，也没有运行经典 R 差异剪接模块。

官方安装说明将只需聚类计数的使用场景与需要 R 模块的差异分析分开。本篇完成前一种工作，再用手算解释结果，不把获得计数矩阵写成已经完成显著性检验。LeafCutter 安装说明

经典项目提交
2c9907ef66adf0bfb3092f0ceb6886ee5c046fbb

聚类脚本 SHA256
4228ef4ff41625a0bc3d004cff85bd7a6eb82e44b8f64c021e997077d587293f

RegTools 实际版本
1.0.0

输入包含一条容易被忽略的桥接连接

我们设计六种内含子区间，四份输入每份一百二十二条比对记录。其中五种通过共享起点或终点构成相连关系，另一种孤立存在，没有其他候选与它组成可变选择。

前两个内含子共享左端，后两个共享右端。中间的低支持内含子同时连接前一部分的一个端点和后一部分的另一个端点，因此把原本可能分开的两个集合串了起来。这里称它为桥接连接，只描述图上的作用，不暗示它一定是伪影。

A1、A2 的使用比例在前后两对输入中不同，B1、B2 则保持相等。桥接连接每份只有两条，但合并四份以后为八条。孤立内含子的支持很高，用来检查不进入可变簇是否一定意味着覆盖不足。

这些输入不是生物学重复，也没有真实测序过程。SAM 的 CIGAR、链标签和支持条数由我们明确写入，序列使用格式测试占位值。实际运行的是这些结构记录的提取与聚类，不包含上游比对准确率评估。

连接提取实际执行了哪些条件

每条设计记录长度一百，跨一个内含子，两侧各五十碱基。我们用 Samtools 完成排序和索引，再执行 RegTools 提取。链方向来自输入中明确设置的标签，所有设计连接都在所设长度范围内。

./regtools junctions extract \
  -a 8 -m 50 -M 500000 -s XS \
  -o S1.junc S1.bam

四份输入各提取六条连接记录，计数与设计一致。输出是包含两侧锚定区段的 BED12，不能只拿外层起止坐标当作内含子边界，需要结合区段长度。

安装时，容器镜像拉取处于排队状态。我们改用校验过摘要的 Bioconda 二进制包，在固定辅助容器中运行，并显式指定已有库路径。主程序版本检查和四次提取均成功，完整驱动已经包含这项安装方式。

这个操作只影响独立实验目录，没有改写系统环境。读者重跑时直接执行驱动即可，工具包摘要、二进制摘要和运行环境都在结果记录中，便于检查实际使用的程序是否一致。

RegTools 的锚定长度、链来源与内含子长度分别控制不同条件，不能把其中一个参数当成总体质量分数。本次固定这些条件，只比较后续聚类规则。RegTools 连接提取说明

三组聚类只改变两个筛选条件

第一组使用较低的比例门槛和支持门槛。第二组只提高比例门槛，第三组则保持低比例门槛、提高支持数量门槛。输入连接文件及其顺序始终相同。

python3 leafcutter_cluster_regtools.py \
  -j junctions.txt -o toy -r permissive \
  -l 500000 -p 0.001 -m 3

python3 leafcutter_cluster_regtools.py \
  -j junctions.txt -o toy -r ratio_filter \
  -l 500000 -p 0.05 -m 3

python3 leafcutter_cluster_regtools.py \
  -j junctions.txt -o toy -r support_filter \
  -l 500000 -p 0.001 -m 10

这些数值用于放大一个可以解释的变化，不是大规模队列推荐配置。三个输出目录分别保存汇总候选、精炼后的簇、整数计数和带分母的计数文件，方便逐步定位差异。

三组实际运行都正常退出。第一组得到一个含五个内含子的簇，后两组均得到两个簇，每个包含两个内含子。孤立内含子没有进入这些最终可变簇。

第一组，一个簇里的分母为什么是八十二

低门槛组保留桥接连接，五种相连内含子组成同一个簇。每份输入在这个簇中的计数和为八十二，来自四十条 A 部分、四十条 B 部分和两条桥接支持。

因此，第一份输入的 A1 记录为三十除以八十二，A2 为十除以八十二，B1 与 B2 各为二十除以八十二。第三份输入的 A1、A2 对调，但分母仍为八十二。

低门槛组 S1
A1        30/82
A2        10/82
桥接       2/82
B1        20/82
B2        20/82

这个分母不是样本总读段数。一百二十二条输入中的四十条孤立内含子支持没有进入该簇，因此不在分母中。它也不是整个基因的转录本数，因为实验根本没有提供基因注释。

比例解释必须带着它所属的集合。只写 A1 使用比例约为百分之三十六点六，却省略簇定义，读者就无法判断这个数值与另一套分析是否在比较同一件事。

提高比例门槛，删除的连接还改变了簇结构

合并四份输入后，桥接连接总支持为八，其他四种相连内含子各为八十，总计三百二十八。桥接占比约为百分之二点四四，超过第一组的低门槛，但低于第二组的百分之五。

桥接被过滤后，前一部分与后一部分不再通过共享端点相连，于是精炼过程把它们拆成两个簇。A1 与 A2 一组，B1 与 B2 另一组。原来的四种主要内含子计数都没有改变。

这时，第一份输入的 A1 变成三十除以四十，也就是零点七五，原先约为零点三六六。看起来增长了很多，但分子一直都是三十，变化全部来自分母从八十二变成四十。

如果把两次运行当成两个生物学条件，直接比较比例，就会得到一个人为制造的变化。参数敏感性检查应该解释分析集合如何改变，不能把参数造成的分母变化包装成表达响应。

支持数量门槛也可能在簇内部起作用

第三组保持低比例门槛，把支持数量条件从三改为十。桥接连接合并后只有八条，因此也被去掉，得到与第二组相同的两个簇。

这一结果提醒我们，参数名称或帮助中的简短描述不一定足以说明全部作用位置。在本次固定脚本里，该数量门槛既参与候选簇筛选，也被传入内部连接精炼。原始簇总支持远高于十，并不能保证簇内每条连接都被保留。

我们同时核对了实际输出和固定版本源码，确认第三组的变化发生在连接筛选之后。读者若遇到总支持很高却有连接消失，应检查门槛是否还作用于单条连接，而不是直接判断程序丢数据。本次使用的经典聚类脚本

这也说明，不能只记录参数值而不记录版本。同一个选项在不同实现中的作用范围，可能需要重新核对。本篇结论限定于所列提交和实际运行，不延伸为所有同名包的统一行为。

高支持孤立内含子没有输出，不等于没有表达

孤立内含子每份有四十条支持，合计一百六十，比桥接连接高得多。RegTools 确实提取到了它，但本次默认可变聚类输出没有保留它。

原因在于它没有形成需要比较的多分支集合。最终矩阵的目的不是列出所有观测到的连接，更不是提供全部基因的表达计数。因此，最终文件中不存在某条内含子，不能反推原始数据没有支持。

定位缺失项时，应先检查连接提取文件，再看汇总候选和精炼结果。这样可以区分没有被提取、没有进入候选簇、被门槛过滤以及没有构成可变集合等不同情况。

只保留最终矩阵而删除中间结构，会让这些原因混在一起。尤其在注释之外的候选分析中，能够追踪一条连接在哪一步消失，往往比只知道最终剩下多少行更有用。

坐标转换需要看实际定义

本实验 A1 的内含子按零基制半开区间写作一千到二千。RegTools 的外层区间还包含两侧各五十碱基的锚定部分，经典聚类脚本进一步转换后，输出标识中的右端为二千零一。

这个差一碱基不是随意偏移。对应的连接可以理解为左侧外显子末端到右侧外显子起始，输出端点定义与原始内含子半开区间不同。把它直接当成普通 BED 区间使用，就可能匹配不到正确位置。

在本次脚本中，左端由外层起点加左区段长度得到，右端由外层终点减右区段长度再加一得到。实际整理注释时，应把这个转换写进代码与说明，而不是凭观察临时给所有坐标加一或减一。

链方向也必须一起保留。相同数字位置在不同链上不一定属于同一个剪接对象，跨文件关联时不能只使用一个不含方向的起止字符串。

整数计数表和比例文件各有用途

本次整数计数文件保存每个内含子在每份输入中的原始支持。另一个文件把相同计数写成分子与簇内分母的形式。前者便于核对支持是否改变，后者便于理解当前聚类下的相对使用。

不要把比例文件直接当成任意差异统计方法所需的原始计数矩阵。不同分析方法对输入格式和模型有具体要求，提前把计数除掉分母，可能丢失关于覆盖深度的重要信息。

同样的零点五，可能来自二十除以四十，也可能来自一除以二。比例相同不代表证据量相同。正式统计流程应保留计数信息，并按所用方法处理样本设计与不确定性。

本篇没有运行差异检验、FDR 计算或剪接数量性状分析。因此，四份人工输入的比例差异只用于解释计算，不代表发现了显著差异事件。LeafCutter 分析流程说明

把这次经验应用到自己的数据

比较多个样本时，应尽量使用清楚定义的共同事件集合和聚类规则。若每组单独构图，再把比例按内含子名字直接相减，分母可能来自不同集合，比较就失去一致含义。

做参数敏感性分析时，可以先找原始计数不变而比例变化的条目，回查它们所在簇的成员是否改变。这样能够快速区分真正的计数变化与图结构变化，避免只盯着某一个看起来变化很大的数字。

LeafCutter 的优点是围绕内含子切除关系组织相对使用，不必先为每条短读段确定完整异构体。局限是结果仍然依赖比对、连接过滤和簇定义，局部关系不能自动回答整条转录本的远距离组合。

低支持连接该不该保留，需要结合原始证据和研究目的。保留可能把局部集合连接得更复杂，删除可能拆分图并改变分母。理解这种后果，比把某个阈值视为永远正确更有帮助。

复现本次结果

下载本组三篇的练习与实跑包。本篇驱动会下载校验过的经典脚本与 RegTools 包，生成四份输入，实际提取连接并运行三组聚类。使用 Linux 主机、Python 和 Docker 执行。

python3 155-leafcutter/run.py results/leafcutter

重跑后先核对每份提取六种连接，再检查最终分别是一个五成员簇和两个双成员簇。最后确认 A1 的整数计数保持三十，而比例从三十除以八十二变成三十除以四十。这个对照能把过滤、结构和数值含义完整连接起来。

如果后续需要持续维护转录组剪接分析环境及队列结果，可了解云生信生物分析服务器服务，统一保存输入、聚类规则和结果版本。