目录
3.5.3 识别在静止，可逆，均匀（SRH）条件下进化的连续位点序列和区块
3. 系统发育分析
3.1 多序列比对
使用MAFFT中的LINS-i 算法对齐1,478 个直系同源群（ortholog groups，OG）。
无法对齐的问题在多序列比对（multiple sequence alignments，MSAs)中非常常见，因此要检查每个OG对齐是否存在异常值，并对其进行重新比对或删除。
多序列比对用于描述一组序列之间的相似性关系，以便了解一个基因家族的基本特征，寻找motif，保守区域等；用于描述一个同源基因之间的亲缘关系的远近，应用到分子进化（构建进化树）分析中；其他应用，如序列拼接，构建PSSM，打分矩阵等。
————引用自知乎用户Bioinfo Note的文章：第四章 多序列比对
- 首先，记录每个 OG 对齐和转录本翻译的最大 C 端和 N 端延伸。
- 随后，计算所有参照物种氨基酸序列 BLOSUM62（55）距离的平均值、中位数和四分位数。然后，计算每个转录本与最佳互作参照类群相应序列的 BLOSUM62 距离。
- 然后，以参考物种之间 BLOSUM62 距离平均值的上四分位距的 2.25 倍作为分界值，检查该距离是小还是大。
异常值：比参照物种的 BLOSUM62 距离的最小值大的转录本被列为异常值，与参考物种序列重叠的位点少于 30 个的序列也被列为异常值。
BLOSUM 打分阵矩阵用于对进化上不同的蛋白质序列之间的比对进行评分。直接基于多序列比对分析亲缘关系较远的蛋白质，而不是用近缘的序列。利用蛋白质序列高度保守区（BLOCKS）用于构建BLOSUM矩阵。BLOSUM r：由一致性小于r%的序列建立的矩阵。例如，BLOSUM62是使用一致性小于62%的序列建立的矩阵（一致性≥62%的序列被聚类以减少相近序列带来的偏好），矩阵后的数字越大，则表示关系越近。
——引用自知乎用户Bioinfo Note的文章：第三章 两序列比对（1）
链接：https://zhuanlan.zhihu.com/p/665405311
检索得到这些氨基酸序列后，将这些异常值序列提取出来，在MAFFT中找到“-add”将每个异常转录本单独重新对齐到参考物种集。这些重新对齐的序列长度不会超出参考物种氨基酸序列的长度，因此，参照物种的多序列比对可以作为骨干，将异常转录本重新整合回 OG 对齐序列中。（若重新对齐后仍被认定为异常值，则删除其核苷酸与氨基酸序列）
剔除氨基酸MSAs中的gap位点。
随后使用 pal2nal.pl 生成与氨基酸 MSAs 相对应的核苷酸 MSAs。
其在密码子查找表中增加了单字母 IUPAC 氨基酸歧义代码 B（天冬氨酸或天冬酰胺）和 Z（谷氨酸或谷氨酰胺），以处理密码子中的歧义。
将包含 1,478 个氨基酸的 MSA 和相应核苷酸水平的 MSA 存入 Dryad Digital Repository(http://dx.doi.org/10.5061/dryad.3c0f1).
3.2 结构域识别
传统上，在基于分区模型的树重建中，基因被视为主要数据块。然而，蛋白质编码基因通常编码一个或多个蛋白质结构域，而蛋白质结构域是蛋白质结构、功能和进化的单位。因此，在分区分析中，结构域编码区应优先于基因作为进化单位。
使用 Pfam 数据库27.0 (http://pfam.sanger.ac.uk/) 以及 pfam_scan.pl 1.3 版和 HMMER 3.0 (http://hmmer.org/) 软件注释同源基因多序列比对中的蛋白质结构域和族（clan，即相关蛋白质结构域的分组）。
文章中提到了Pfam-A和Pfam-B，现版本Pfam-B已经不存在了。
做保守结构蛋白的注释Domains，需要用到PfamScan软件，以及Pfam-A的数据库。安装后的PfamScan目录下主要有三个文件ChangeLog， pfam_scan.pl，README以及一个文件夹Bio（主要存放需要的模块）。主要用到的pfam_scan.pl是一个perl脚本。
pfam_scan.pl参数如下：
一般情况下，只需要用到三个参数：
-fasta需要检索的蛋白序列的fasta文件；-dir存放Pfam-A数据库的目录；-outfile需要输出的文件名字。
已知蛋白的PF编号（左）已知蛋白的名字（右） ——引用自简书作者群体遗传学的文章
链接：https://www.jianshu.com/p/e81af62323ab
3.3 Alignment Masking（比对遮蔽）
识别每个 OG 对齐中可能存在的对齐歧义块或随机 MSA 部分（文章中使用的是Aliscore ，还有许多其他的方法能达到相同的目的）
氨基酸转录本的插入缺失标记尾端用“X”填充；核苷酸序列中则使用“N”填充。
3.4 数据块连接与选择
生成了如下数据矩阵：
超级矩阵 A（2）
- 基于直系同源基因的超级矩阵：删除每个 OG 比对中的随机比对区域后，将 MSA 连接起来
- 基于结构域的超级矩阵：在对结构域进行注释并去除随机排列的区域后，将蛋白质结构域族、结构域和空白区域的联合族数据块连接起来（所有同族的结构域被集中到一个数据块中。对于没有族分类但数据集中有一个或多个注释基因座的每个域，所有基因座都被集中到一个数据块中。最后，将每个基因的所有空白区域集中到一个数据块中。）
超级矩阵 B
超级矩阵A的子集
对于超级矩阵 A 中的每个数据块，使用 mare 软件评估潜在的系统发生信号。使用定制的 PERL 脚本，只将具有潜在系统发生学信号的数据块连接到超级矩阵 B 中。
超级矩阵 C，SOS（选取最优子集）（4)
筛选出39个分类群（附件S6）后，从超级矩阵 B 中提取了所有数据块，每个分类群中至少包含一个物种，并将这些数据块连接成一个选取最优子集（selected optimal subset，SOS）。
生成了两个SOS，一个是氨基酸序列层面的 SOS，另一个是核苷酸序列层面的相应 SOS。
根据超级矩阵 C，生成相应的核苷酸超级矩阵，一个包括所有密码子位置，另一个只包括第 2 个密码子位置。
超级矩阵 D，masked SOS
超级矩阵 D 是超级矩阵 C 的子集，其中包含的位点似乎与全局 SRH 条件下的进化假设相一致。
我们使用配对同质性检验分析了超级矩阵 C，以评估数据是否符合全局静止、可逆和同质条件（SRH）下的进化（见第 3.5.3 章）。利用滑动窗口法和对称性配对检验，我们将 3,000 个氨基酸位点的区块分为似乎符合这一假设的区块和似乎违反这一假设的区块（见第 3.5.3 章）。
3.5 序列组成异质性和缺失数据模式分析
核苷酸或氨基酸排列中序列组成异质性意味着这些序列不可能是在SRH（静止，可逆，均匀 stationary, reversible and homogeneous conditions）条件下进化而来的。但大多数基于模型的分子系统发生学方法都假定进化是在这些条件下进行的，而当数据是在非 SRH 条件下进化时，出错的概率将会更大，因此在使用比对来推断系统发生树之前对比对进行调查是非常重要的。
序列组成异质性
目前常用的建树模型都假设DNA序列在不同生物类群的进化过程中, 其碱基的使用频率是不会发生
改变的, 即进化过程中序列组分是同质的。然而在真实的情况下, 不同生物的序列组成成分存在一定的差异, 如果分析的数据集在不同物种间存在较大的序列组成异质性, 则在构树过程中往往会将序列组成相似的物种错误地聚为一支, 与真实进化关系存在偏差。 在无脊椎动物(如节肢动物门)、植物、真菌等生物类群中, 序列组成异质性现象较为常见。一般来说, 相对于第一、二位密码子, 第三位密码子比较容易出现较强的序列组成异质性; 相对于核酸序列, 蛋白质序列的组成异质性会更低. 因此, 在系统发育基因组学分析中, 可以尝试对DNA序列进行简并化处理或在蛋白质水平上进行系统发育分析, 以降低序列组成异质性可能带来的系统误差影响。——李佳璇, 梁丹, 张鹏. 系统发育基因组学方法研究进展[J]. 中国科学:生命科学, 2019, 049(004):456-471.
3.5.1 缺失数据的分布
由于缺失数据在排列序列中的不均匀分布会影响树的重建，作者使用 AliStat（http://www/csiro.au/alistat）分析了超比对，并生成了序列成对比对的缺失数据分布热图，记录了这些序列对的缺失数据程度。说明了缺失数据在类群中的分布情况。结合 "流氓类群（rogue taxa ） "分析（见 第 3.9 章），利用这些信息来确定可能存在系统发育错位的物种。
3.5.2 置换检验
氨基酸或核苷酸序列组成的异质性、替代过程的非稳态性以及缺失数据的非随机分布在物种间经常出现。这些现象违反了常见的模型假设，会使基于似然法的树重建结果产生偏差。
为了评估这些现象的潜在影响，作者消除了经验数据中的系统发生学信号，但没有消除成分异质性和非稳态替代过程，也没有使用置换方法改变缺失数据的分布。
因此，作者针对 13 个选定的系统发育假说（见第 3.10 章），使用四重映射法（quartet mapping）比较了原始实验数据（original empirical data）和置换数据的信号强度，以了解原始实验数据中的信号强度是否可以用混杂作用(confounding effects)的存在来解释。
假设信号具有可加性，作者进行了以下排列：
1. 置换方案 I：在不改变缺失特征分布的情况下，对每个元分区和序列中的所有不明确的氨基酸特征进行置换。这种排列方法完全消除了物种间的系统发生学信号。
2. 置换方案 II： 使用 WAG 替代矩阵的氨基酸频率（假定数据显示了同质的位点组成）来生成在长度和缺失数据分布方面与实验元分区相对应的数据块。与实验数据相比，该方法生成的数据块大小相等，缺失数据分布相同，但没有系统发生学信号，也不存在氨基酸序列组成异质性和非稳态替换过程。
3. 置换方案 III：生成了与置换方案 II 类似的数据块（与元分区的大小相对应），但也使用了随机缺失数据，以获得具有同质氨基酸序列组成、静态置换过程和随机分布缺失数据的矩阵。
利用软件 mare并使用完整的级联超矩阵（concatenated supermatrices） C 和 D，将这些置换数据矩阵的信息含量与原始实验数据的信息含量进行比较。
4. 通过比较 permutation I II、III 和 FcLM 的结果，检验了混杂效应对所选系统发育假说的影响：我们检验了特定系统发育假说的鲁棒性是否已经可以用混杂效应的强度来解释。我们在进行所有的置换检验时都结合了对超矩阵 C 和 D 的 FcLM 分析。
鲁棒性：Robust的音译，也就是结实、耐用的意思。指在异常和危险情况下系统生存的能力。这个名词第一眼看到的时候确实觉得不够信达雅，不过搜索了一下发现有人做了考古把「robust」翻译为「鲁棒」最早起源于哪里？ - 知乎 (zhihu.com)
作者使用 RAxML 进行了分区 FcLM 分析，对方法 1 和 2 使用了最初选定的元分区模型假设（见补充文本 Supplementary Text 第 2.6 章）。对于方法 3 和 4 ，我们在所有 FcLM 分析中使用了 WAG 替代矩阵。执行这些置换的软件是用 PERL 编写的。
3.5.3 识别在静止，可逆，均匀（SRH）条件下进化的连续位点序列和区块
为了确定序列数据是否在全局 SRH 条件下进化，作者使用了三个配对检验方法（参考文献 63、74、75）对超级矩阵 C 进行了调查（而不是使用经典的双向独立性检验，后者不适合此类评估：(67)）。接着使用 SymTest软件（http://www.csiro.au/symtest2）为每对序列生成了三个 p 值，每个 p 值分别来自上述每个检验。假设序列在全局 SRH 条件下进化，观察到的 p 值将呈均匀分布。之后将观察到的 p 值与预期 p 值绘制成 PP（概率-可能性）图。如果观察到的 p 值的分布呈 J 型，低 p 值的比例大于预期，那么说明有几对序列不可能是在全局 SRH 条件下进化的。
核苷酸或氨基酸（与超矩阵 C 相对应）的对齐均采用这三种配对测试方法进行检验。在核苷酸层面，分别对第一、第二和第三密码子位置进行测试。
为了确定哪些序列对是可以被假定在全球 SRH 条件下进化的，而哪些序列对违反了这一假定，我们根据Bowker（74）的配对检验所得到的 p 值生成了热图（注意，p 值没有进行调整以应对多重比较的问题）。
此外，作者还使用滑动窗口法分析了超级矩阵 C，以确定可假定在全局 SRH 条件下演化的连续位点块。这种统计方法的威力取决于窗口的大小和特征状态的数量。因此，对于核苷酸序列来说，相对较窄的窗口就足以识别出是否违反了在全局 SRH 条件下进化的假设，而对于氨基酸序列来说，则需要更宽的窗口。在这两种情况下，窗口越宽，检验就越有效。为了确定一个窗口大小，一方面使其足以检测非 SRH 条件下的进化，另一方面使我们能够将超矩阵 C 分成尽可能多的区块，我们使用不同的窗口大小和不同的窗口起始位置对氨基酸数据进行了检验。为了解决多重比较的问题，我们对显著性水平进行了调整（采用 Bonferroni 方法）。在全局 SRH 条件下明显偏离进化的位点块从超矩阵 C 中屏蔽，生成超矩阵 D（因此，超矩阵 D 是超矩阵 C 的子集）。
滑动窗口是一种流量控制技术，通信领域专业名词（？）
对超级矩阵 D 重复进行上述测试，此外，还对所有置换（permutation）矩阵运行了 SymTest （见第 3.5.2 章），以评估作为置换方法中混杂信号潜在来源的全局 SRH 条件下进化假设的偏差。接着使用 Bowker（74）的配对检验（见上文），根据成对序列比较得出的 p 值生成了热图（见补充文本图 S20）。
热图是用 PERL 编写的定制软件生成的。
扫一扫
8074
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
