# Methods (完整版,替换原稿"Skeleton"占位)

## 1. 序列数据获取与筛选

**数据来源。** 跨物种class A GPCR序列取自UniProt(SwissProt+TrEMBL),覆盖29个物种,
从非两侧对称动物(海绵、栉水母、扁盘动物)经原口动物、无脊椎后口动物到脊椎动物
(哺乳类、鸟类、爬行类、两栖类、鱼类)。全库合并文件`all_species_combined.fasta`
含23,036条候选序列。

**参考受体集与比对方法。** 采用19个跨亚科GPCRdb参考受体(人类/牛源代表,覆盖氨基类、
肽类、趋化因子肽类、脂质类、核苷类、糖蛋白激素类、视紫红质类七大配体家族)。用Biopython
`PairwiseAligner`(BLOSUM62替换矩阵,global模式,open_gap=-11,extend_gap=-1)将每条
候选序列与19个参考受体逐一比对,取相似度最高者为`best_ref`标签,序列身份百分比(pct_identity)
同时记录。

**质量筛选阈值。** 以25%序列身份为最低可信阈值(low tier: 20-25%,moderate: 25-35%,
high: >35%,三档置信度标签保留在最终数据表`confidence_tier`列,供敏感性分析使用)。
低于20%身份的比对结果被排除。6个generic-number关键位点(1.50x50/3.50x50/5.58x58/
6.30x30/7.49x49/7.50x50)全部成功调用的条目共11,609条,是本文Results §1/§3大部分
定量分析的基准样本。若仅需3个位点(6.30/5.58 + 早期分析),样本量为12,484条(见§21节
差异说明)。

**Generic numbering(通用编号)。** 采用Ballesteros-Weinstein编号并按GPCRdb结构校正
版本转移,由参考受体的已知生成号逐位映射到候选序列的比对位置。

## 2. 系统发育树构建与预处理

**基因树来源与已知局限(本次revision新增的透明化说明)。** 本文使用的
`species_tree.nwk`是对16,620条候选序列(6位点齐全的11,609条经`keep.tip`剪枝后的
子集,或3位点集12,484条,视具体检验而定)用FastTree构建的**基因树**,而非物种树——
其默认(中点)定根产生生物学上不可信的11,525:84基部分裂,表明FastTree对这棵无根树的
定根是任意的,不代表真实的祖先分歧点。这棵基因树上,同一物种的多个不同旁系同源基因
(paralog)以独立tip的形式并存,`best_ref`功能标签在树上高度多态(MRCA纯度0.001-0.125,
即同一标签下的序列散布于几乎整棵树)——这意味着直接在此树上进行的系统发育比较检验,
把不同paralog当作独立的进化谱系处理,存在因paralog非独立性未被正确建模而导致效应量
被低估或高估的风险。

**预处理流程(与原§21.2完全一致,逐步列出)。**
1. 用`ape::keep.tip`将16,620-tip基因树剪枝到目标性状数据集的tip集合;
2. 用`ape::multi2di`对残余多分支节点做随机二值化;
3. 对长度小于1×10⁻⁵的退化分支长度做floor修复(设为1×10⁻⁶),避免下游似然计算中的
   数值奇异。

**针对paralog混合问题的方法学修正(本次revision新增)。** 为直接检验paralog混合是否
影响核心共进化结论,本次revision复用项目Direction-2阶段已建立的DISCO式基因树-物种树
调和流程(`report_section22/23`):对16,620-tip基因树逐节点判定duplication/speciation
(判据:任意两个子节点物种集合有重叠即为duplication),在duplication节点处切断,得到
13,539个直系同源群(orthogroup)——其中12,441个(91.9%)是仅含单条序列的singleton,
无法承载跨物种比较;仅1,098个(8.1%)覆盖≥2个物种,合计4,179条tip(占全部候选tip的
25.1%)。**本次revision的三组pairwise共进化检验和corHMM联合模型比较,在原始全量数据
基础上新增了限定至这1,098个真直系同源群(且三性状齐全)的重复检验**,作为paralog混合
敏感性分析(见Results §3新增段落)。

**注意事项:两套树对象不可混用。** §22另行构建了一棵28-tip的**物种级**骨架树
(`species_tree_final_rooted_support.nwk`,经DL对账最小化定根,DL cost=23,178),
用于Direction-1的性状起源时间顺序分析(表征"哪个性状更古老"这一问题),这棵树的每个
tip对应一个物种而非一条序列,**不能**、也未被用于承载本文Results §3的序列级
fitPagel/corHMM检验——请勿将两棵树的用途混淆(原稿"species tree pruned to the
sequence set"这一表述曾造成歧义,已在此明确澄清)。

## 3. 系统发育比较检验

**性状定义。**
- 6.30离子锁(`lock`):6.30x30位点为Asp或Glu(负电荷)记为1,否则记为0;
- 5.58 Tyr(`tyr558`):5.58x58位点为Tyr记为1,否则记为0;
- 钠离子口袋完整(`na_pocket`):同时满足1.50x50=Asn、7.49x49=Asn、7.50x50=Pro
  (经典NPxxY-钠离子口袋组合)记为1,否则记为0。

**Pairwise共进化检验(Pagel's method)。** 用`phytools::fitPagel`(ARD模型)对
(lock, tyr558)、(lock, na_pocket)、(tyr558, na_pocket)三组二元性状分别检验系统发育
校正后的共进化关系,零假设为两性状独立演化。

**联合8态corHMM模型比较。** 将三个二元性状复合为一个8态离散字符(2³=8种组合),
构造4个候选转移速率矩阵,用corHMM 2.8(`rate.cat=1, node.states="none"`)拟合并按
AICc比较:
- M0(完全独立):三个性状各自的转移速率仅依赖自身当前状态和方向,不依赖另外两个性状,
  共6个自由参数;
- M1(全ARD/完全耦合):任一性状的转移速率可依赖另外两个性状的完整组合状态,在"一步
  只变一位"(汉明距离1)的生物学假设下,8态立方体每节点3条边,共24个独立速率参数;
- M2(钠离子口袋为轴心):锁和Tyr的转移速率只依赖钠离子口袋当前状态(不依赖对方),
  钠离子口袋自身的转移速率依赖另外两者的完整组合,共16个参数;
- M3(Tyr为轴心):对称地以Tyr为轴心,共16个参数;
- **M4a/M4b/M4c(新增,本次revision针对候选模型集合穷尽性的扩展)**:分别检验
  "仅锁-Tyr耦合、口袋独立"、"仅锁-口袋耦合、Tyr独立"、"仅Tyr-口袋耦合、锁独立"
  三种更简约的两两耦合结构,每个10个参数,填补M0(6参数)与M2/M3(16参数)之间的
  中间地带,回应"候选模型集合未穷尽其他reduced coupling structures"这一方法学关切。

所有速率矩阵的参数索引在Python中程序化生成(而非手工绘制),避免人工转录错误;
生成代码逻辑:对每条汉明距离为1的有向转移边,按其"变化的性状位"和"转移方向"分组,
M0仅分组到(性状位,方向),M1每条边独立成组,M2/M3/M4a-c的分组规则见上文。

**参数化parametric bootstrap。** 因全量数据(11,609-16,620 tip)上单次corHMM拟合
需要2-2.5小时量级,bootstrap在缩小规模的子样本(n=300-2000,视计算资源可用性而定,
seed=1/2/3做独立子抽样)上进行:先在观测数据上拟合M0,取其收敛速率矩阵作为零假设下的
生成过程,用`phytools::sim.history`模拟null性状历史(固定根状态,默认为"000"全阴性,
另设根状态敏感性检验),对每个模拟数据集重新拟合M0和M1,记录ΔAICc(M0-M1)的null分布,
与观测ΔAICc比较,按plus-one Monte Carlo方法(Davison & Hinkley 1997)计算经验p值:
p = (1 + 观测值以上的null计数) / (1 + 总模拟次数)。

**Leave-one-clade-out jackknife。** 在同一子样本上,依次移除最大的若干个`best_ref`
参考谱系分组(本次revision从原先3个扩展到6个),重新拟合M0/M1并记录ΔAICc,检验协同
演化信号是否被任何单一过度代表的谱系驱动。

**计算环境。** 全部系统发育比较分析在远程主机(`ssh:run`)的conda环境`gpcr_phylo`
(R 4.3.3,corHMM 2.8,phytools 2.5.2,ape)中运行。

## 4. Ligand-mechanism预测力的多变量检验(本次revision新增)

为直接检验"配体家族是否预测Arg3.50保护机制"这一问题(而非仅从熵值间接推断),构建了
多变量分类模型:在n=11,059条序列(排除3.50本身非Arg的550条)上,以mechanism_top
(5类)为因变量,分别用ligand_family(7类)、gene/best_ref(18类精配对参考受体)、
clade(15个演化谱系)及其组合作为特征,用scikit-learn `LogisticRegression`(one-vs-rest,
L2正则化)拟合多分类逻辑回归,5折分层交叉验证评估log-loss。核心检验:在gene身份已知
的条件下,加入ligand_family是否带来统计显著的额外预测力——用置换检验(200次组内
[按gene分组]洗牌ligand_family标签,保留gene-mechanism关联不变)构造零分布,与观测的
log-loss降幅比较。另计算ligand_family、gene、clade与mechanism之间的归一化互信息
作为辅助描述统计量。

## 5. Arg3.50保护机制分类

对全库11,609条序列,提取3.50(3.46x46-3.54x54窗口)和6.30(6.26x26-6.34x34窗口)
周围各±4个残基的generic-number列,按以下规则逐条分类为六类之一:
1. **经典跨螺旋锁**(classic_cross_helix_lock):6.30为Asp/Glu,与3.50-Arg形成跨
   TM3-TM6盐桥;
2. **同电荷+螺旋内救援**(same_charge_intrahelix_rescue):6.30为Arg/Lys(与3.50
   同电荷),但3.50或6.30窗口内另有Asp/Glu残基就近提供电荷中和;
3. **6.30中性+局部酸性**(neutral_630_with_local_acidic):6.30本身既非酸性也非
   碱性,但窗口内仍有酸性残基;
4. **同电荷+无局部救援**(same_charge_no_local_rescue):6.30同电荷且窗口内确未
   找到任何酸性残基;
5. **6.30中性+无局部酸性**(neutral_630_no_local_acidic):6.30中性且窗口内也无
   酸性残基;
6. **3.50非Arg**(no_arg350_excluded):3.50位点本身发生了替换,单独讨论,不计入
   前5类的机制多样性统计。

## 6. 直系同源群(Orthogroup)重建

采用DISCO式基因树分解算法,对16,620-tip基因树逐节点判定duplication/speciation:
若一节点的任意两个子节点物种集合有重叠,判定为duplication节点,并在此节点处切断
树结构,得到互不重叠的直系同源群。判据通过两种独立方式验证:(a) 全树64.0% duplication
/ 36.0% speciation的比例与方向2项目独立计算的bitmask方法基准完全一致;(b) 分解后
全部13,539个直系同源群逐一验证无单拷贝违规(0 violations)。

对每个直系同源群,计算三个版本的机制多样性Shannon熵(按ligand_family分组):
(a) best_ref版(原始功能标签分组);(b) orthogroup pooled版(多序列直系同源群内
全部序列按多数ligand_family标签汇总后统一算熵);(c) orthogroup within-group版
(先对每个直系同源群单独算熵,再按多数标签取平均)。

## 7. 结构分析

**结构来源。** 实验测定结构(X-ray/cryo-EM)从RCSB PDB获取(mmCIF格式);
AlphaFold单体预测结构从EBI AlphaFold DB API获取。

**距离测量(本次revision修订:改用精确重原子距离,不再仅用侧链质心距离)。**
对每个候选intra-TM6酸性搭档案例,用gemmi解析mmCIF坐标,计算候选酸性残基
(Asp的OD1/OD2,或Glu的OE1/OE2)与Arg候选搭档(NE/NH1/NH2)之间的**最小重原子
间距离**(而非侧链质心间距离),同时报告该距离对应的具体原子对。对可靠解析的案例
额外计算salt-bridge几何角度(酸性残基CG/CD-最近O原子-Arg最近N原子夹角)作为
补充几何描述。当候选残基侧链在实验密度图中未完整解析(仅到Cβ)时,明确标注
"不可评估精确距离",不用截断侧链的代理距离冒充完整盐桥测量。

**构建体工程属性标注。** 对每个纳入Table 1的PDB结构,通过mmCIF的`_struct_ref_seq_dif`
(SEQADV)记录统计工程突变数、表达标签数、连接子(linker)数,并结合结构标题/文献
判断该结构是否为热稳定化(thermostabilized)构建体、融合蛋白嵌合体(如BRIL/T4溶菌酶
插入)、或特定构象状态(活化态复合物 vs 非活化态),明确区分"野生型静息态直接证据"
与"工程改造构建体推断证据"。

**排除标准。** 融合蛋白插入残基(编号规则通常为PDB residue number ≥1000,如T4溶菌酶、
rubredoxin、cytochrome b562插入片段)、hetero残基、backbone原子均从酸性搭档搜索中
排除。搜索半径:结构级分析12Å,全库generic-number普查15Å。

## 8. 数据与代码可用性

本文全部原始分类表、模型比较输出、bootstrap/jackknife结果、结构分辨note及图源文件
存放于Zenodo(CC-BY 4.0开放许可):https://doi.org/10.5281/zenodo.21351361 (已于
本次revision核实记录公开可解析,压缩包含14个文件,详见正文"Data availability")。
本次revision新增的orthogroup限定重跑数据、ligand预测力检验结果、精确结构距离测量、
扩展bootstrap/jackknife结果将同步更新至该Zenodo记录的新版本。