两个样本的优势碱基不同,为什么 inStrain 仍然报告群体层面的一致性为一?理解这个问题,需要同时查看每个位点保留了哪些等位基因、各自的支持数量,以及两个样本究竟比较了多少位置。
本文安装 inStrain 1.12.0,生成四组配对比对数据,实际完成四次 profile 和两次 compare。实验把参考长度、读段数量与位置固定下来,只改变一个位点的碱基组成或碱基质量。最终输出可以逐位复算,帮助读者理解不同指标描述的含义。

两类一致性指标在回答什么问题
宏基因组读段往往来自同一类生物的多个细胞,因此在一个参考位置上可能同时观察到不同碱基。以数量最多的碱基构建共识,可以得到一个简化表示,但这个表示会略去仍然存在的次要等位基因。
inStrain 同时提供共识层面与群体层面的比较。共识层面关注两个样本的优势碱基是否相同,群体层面进一步考虑能够通过条件的等位基因是否共享。阅读结果时,将两种指标与实际碱基计数放在一起,通常更容易解释差异。官方核心概念说明
工具输入通常包含参考序列以及读段比对到该参考后得到的 BAM。参考名称、参考内容与比对头必须对应。多个样本要比较同一批位置,也需要保持参考体系一致。官方输入与运行说明
四组输入只在一个位置上有所改变
本次教学参考长两千碱基,使用固定随机种子生成。每组数据包含五百对读段,一共一千条记录,每条长一百碱基。配对成员相邻排列,在整个参考上形成五十倍均匀覆盖。
参考的第五十五个位置是碱基 T。我们在甲样本中安排四十条 T 和十条 A,在乙样本中安排十条 T 和四十条 A,在丙样本中安排五十条 C。三个样本的其他位置全部与参考一致。
第四组保留甲的序列与比对位置,只把该位点十条 A 的碱基质量从四十降到二十。其他碱基质量保持四十,读段比对质量统一为六十。这样能够独立观察位点质量过滤的影响。
| 输入样本 | 目标位置 T | 目标位置 A | 目标位置 C | A 的碱基质量 |
|---|---|---|---|---|
| 甲 | 40 | 10 | 0 | Q40 |
| 乙 | 10 | 40 | 0 | Q40 |
| 丙 | 0 | 0 | 50 | 不适用 |
| 甲的低质量对照 | 40 | 10 | 0 | Q20 |
这些是直接构造的教学比对记录,没有经过真实读段搜索。它们用来验证已知输入下的过滤、变异统计和比较行为。真实项目还需要考虑错配位置、重复区域和错误比对等额外来源。
安装与运行时需要固定哪些条件
本次在独立的 Python 3.10 环境中通过官方提供的安装方式安装固定版本,配套 pysam 为 0.24.1。读段记录、排序和索引由 pysam 实际生成,随后运行正常的命令行流程。项目安装说明
python -m pip install inStrain==1.12.0
inStrain profile --version
inStrain profile A.bam reference.fna \
-o A.IS -p 2 -s genomes.stb \
--skip_plot_generation \
--skip_mm_profiling
本次关闭绘图,并关闭按读段错配数量逐级展开的分析,以集中检查最终保留读段的结果。位点覆盖门槛、变异频率门槛、碱基质量门槛和其他过滤参数保持本次版本默认值。
安装后读取实际帮助输出也很有价值。本版默认最低碱基质量为三十,默认最小覆盖为五,默认最小变异频率为百分之五。是否检出变异还会受到统计模型条件影响,因此仅达到频率门槛并不自动保证某个低计数变异被保留。
版本同样属于分析条件。作者的更新记录说明,第一点十二版调整了变异检出的阈值行为,复现更早结果时需要关注相应兼容参数。本次采用这一版本的正常默认行为。官方版本变更记录
先核对读段是否真正进入了统计
四组输入的配对过滤结果都保留五百对读段,输入记录数量均为一千。甲、乙、丙的平均覆盖均为五十,覆盖广度为一,达到位点覆盖门槛的比例也为一。
这些数值为后续比较提供了明确基础。五十倍覆盖乘以两千碱基,得到十万个比对碱基,等于一千条一百碱基读段的总长度。读段数量、参考长度与覆盖统计在本次设计中能够直接闭合。
查看自己的结果时,也可以先做类似核对。若大量配对记录被过滤,应检查配对关系、插入距离、读段一致性和比对质量,再解释后面的变异数量。过滤后的支持量发生改变,会影响可比较位置以及等位基因的检出。
本次低质量对照的读段对也全部保留,但某个位置上的部分碱基没有进入位点计数。这说明读段层面通过过滤,与每个碱基都能参与变异计算,是两个需要分别检查的环节。
甲乙优势碱基不同,为什么 popANI 仍为一
甲的目标位置以 T 为主,占百分之八十,A 占百分之二十。乙的目标位置以 A 为主,占百分之八十,T 占百分之二十。实际位点输出完整保存了这两组计数。
接着运行正常比较命令,并要求保存差异位置,以便从汇总结果退回具体位点。
inStrain compare \
-i A.IS B.IS C.IS \
-o compare_all -p 2 \
-s genomes.stb \
--skip_plot_generation \
--store_mismatch_locations
甲乙实际比较了全部两千个位点。目标位置的共识碱基分别为 T 和 A,因此共识差异计数为一,对应共识一致性为一减去两千分之一,也就是零点九九九五。
群体比较中,甲和乙都保留了 T 与 A,存在共享的等位基因。因此该位置没有被计为群体差异,群体差异计数为零,群体一致性为一。
compared_bases_count = 2000
consensus_SNPs = 1
population_SNPs = 0
conANI = 0.9995
popANI = 1.0
这个结果与两组读段组成不同同时成立。群体一致性为一,说明按照本次可比较位置和检出条件,没有观察到不共享等位基因的差异。它没有把两个样本的等位基因比例规定为相同。
如果研究问题涉及菌株比例变化,可以进一步比较位点频率、多个连锁位点以及其他证据。单个一致性数值保留的信息有限,需要回到读段支持来回答比例是否改变。
加入完全不同的等位基因后,指标如何响应
丙样本在目标位置只保留 C,与甲乙的 T 和 A 都不共享。实际位点输出显示五十条 C,参考碱基 T 的计数为零。该位点被记录为相对参考的固定替换。
甲与丙、乙与丙的比较都覆盖全部两千个位点。两组比较的共识差异和群体差异均为一,因此两个一致性指标都等于零点九九九五。
| 比较 | 实际比较位置数 | 共识差异 | 群体差异 | conANI | popANI |
|---|---|---|---|---|---|
| 甲与乙 | 2000 | 1 | 0 | 0.9995 | 1 |
| 甲与丙 | 2000 | 1 | 1 | 0.9995 | 0.9995 |
| 乙与丙 | 2000 | 1 | 1 | 0.9995 | 0.9995 |
三组比较的共识一致性相同,群体一致性却能区分甲乙之间共享等位基因的情况。这正是保存两类指标的价值。只保留其中一个字段,会失去另一层可以解释群体组成的信息。
本次程序的聚类输出把甲乙放在同一组,丙放在另一组。这个结果可以作为软件在指定条件下的分类记录,但教学参考只有两千碱基,不能据此给真实项目设定普遍适用的菌株判定界限。
碱基质量为什么会让一个变异消失
低质量对照与甲拥有相同的读段序列,只把十条 A 在目标位置的质量降低到二十。默认质量门槛为三十,这些碱基不再计入该位置的变异支持,四十条高质量 T 仍然保留。
实际输出中,低质量对照的变异表为空,变异位点计数为零。平均覆盖由五十降到四十九点九九五,恰好相当于全长累计覆盖减少十个碱基,再除以两千。
原有累计覆盖 = 50 × 2000 = 100000
减少的有效碱基 = 10
新的平均覆盖 = 99990 / 2000 = 49.995
尽管全局平均覆盖只下降了零点零零五,目标位置的次要等位基因已经完全失去有效支持。因此,对于关心的变异,应该检查位点质量和碱基计数,单看全长平均覆盖可能看不出影响。
甲与低质量对照的两个一致性指标都为一,因为它们的共识仍是 T,群体比较也共享 T。但甲保留一个多态位点,低质量对照没有检出该位点的多态性。再次说明,不同输出字段关注的是不同层面的问题。
调整质量门槛时,还应同时考虑对应错误模型和测序平台的质量含义。把门槛调低能够纳入更多碱基,也改变了如何区分真实变异与错误的条件,需要结合实际数据进行判断。
多样性数值也能从本例复算
甲乙的目标位置都存在百分之八十与百分之二十两种碱基。用一减去各碱基频率平方和计算该位点的多样性,得到零点三二。其他一千九百九十九个位点只有一种碱基,贡献为零。
将这个位点的贡献分摊到全部两千个覆盖位置,得到零点零零零一六。程序为甲乙报告的普通核苷酸多样性都接近这一数值,浮点表示导致末位略有差别。
丙只有单一 C,虽然相对参考发生固定替换,群体内部在这个位置仍然只有一种碱基,因此本例的核苷酸多样性为零。固定替换数量与群体内部多态程度可以发生不同变化,解释时应分别查看。
结果还包含经过抽样处理的多样性字段。这些字段涉及额外随机过程,复现时可能出现小幅变化。比较本例的确定性结论,应首先核对普通多样性、实际位点计数和可比较位置,再单独解释抽样指标。
在队列分析中,可以为每一对样本建立一行核对记录,同时保存参考版本、有效比较位置数、覆盖比例、两类差异计数和两个一致性指标。这样在某一对样本突然表现得特别接近时,可以先确认它们是否只是共同覆盖区域较少,再检查位点层面的证据。
比较条件也应随这张记录保存。两份结果若使用不同的质量门槛或变异频率要求,共享等位基因集合就可能不同。可以先在统一条件下重建需要比较的结果,再观察样本差异是否仍然存在。对于临界支持的位点,保留过滤前后计数,比单独保留最终判定更有利于后续复查。
从实际输出中优先检查哪些文件
第一类是读段过滤表,回答多少输入通过了各项条件。第二类是片段或基因组汇总表,回答平均覆盖、覆盖广度和变异数量。第三类是位点表,回答具体哪个位置出现哪些碱基、各有多少支持。
跨样本比较时,再查看比较位置数量和成对差异表。本次所有主要比较都有两千个有效位置,数字很容易复算。真实数据中,如果只比较到参考的一小部分,即使一致性很高,也应明确这只是已覆盖区域内的结果。
位置编号也要核对。本次位点表中的位置写作五十四,对应从一开始计数的第五十五个碱基。把它转换到其他浏览工具或注释表时,应先确认对方采用的坐标规则,避免因差一位而读取错误碱基。
inStrain 能把覆盖、群体多态性和精细比较连接起来,适合需要理解读段层面差异的分析。相应地,参考选择、过滤条件和有效覆盖都会影响结果。先检查这些证据,再把一致性与研究问题连接起来,会更容易得出有依据的解释。
下载包包含四组比对的生成脚本、正常分析命令及本次实际输出记录。读者可以改变等位基因比例或质量,观察哪些字段发生变化,并用本文的计算方式进行核对。下载实跑脚本与本次结果
当研究需要对多批样本持续进行覆盖、微多样性和菌株群体比较时,可了解云生信的生物信息学分析计算服务。
