简介myTAI 是一款面向进化转录组学研究的 R 语言分析工具包适合具备一定 R 与 Bioconductor 基础的科研人员、研究生及 evo-devo 方向学者使用。它通过量化转录组保守模式及其对应的生物学过程帮助用户在计算机上初步筛选感兴趣过程中可能存在的进化限制为后续实验提供探索性依据。资源包共 330 个文件包含 95 个 R 源码、74 张 png 图表、71 个 html 文档、53 个 rd 帮助文件以及 Rmd 教程、cpp 源码、rda 数据与 yml 配置等压缩包约 9.32MB覆盖从函数实现、示例数据到文档说明的完整结构。目前已有 265 人学习。借助标准化、自动化的分析框架读者可快速上手转录组保守性检测流程理解基因集与表达谱的进化约束分析思路并参考官方文档与引文信息完成可复现的研究实践。1. myTAI:R进化转录组学到底在算什么如果你手上有一批不同物种、不同发育阶段的转录组数据想回答“哪些基因在进化上被保守地调控、哪些在某个谱系里被重新布线”myTAI 就是 R 生态里专门干这件事的包。它不负责比对、不负责定量输入是已经算好的表达矩阵加一份物种/发育阶段的映射表输出是转录组进化速率、保守性打分、阶段间表达位移这些指标。我第一次接触它是在一个跨物种胚胎发育的项目里当时用 DESeq2 做完差异分析就卡住了——差异基因列表有了但“这些差异在进化尺度上意味着什么”没人能回答。myTAI 补的就是这一段把表达矩阵映射到系统发育树上用进化模型算出每个基因的转录组进化速率TDI、TAI 这类指标再按发育阶段做富集。适合的人群很明确做比较转录组、evo-devo、发育生物学且已经能熟练跑 R 和基本 RNA-seq 流程的人。新手直接上会卡在“表达矩阵怎么对齐到发育阶段”这一步后面会细讲。2. 把表达矩阵喂进 myTAI 之前要做的三件事myTAI 对输入格式挑剔不是随便一个 counts 矩阵就能塞进去。我见过太多人卡在Error: PhyloExpressionSet报错上根因基本都是输入没对齐。这一章把数据准备讲透后面跑分析才不会反复翻车。2.1 表达矩阵的列名必须能映射到发育阶段myTAI 的核心对象叫PhyloExpressionSet它要求表达矩阵的列不是样本名而是发育阶段。也就是说如果你有 3 个物种、每个物种 6 个发育阶段、每个阶段 2 个重复你得先把重复合并成阶段均值列名变成类似stage1到stage6并且所有物种共用同一套阶段命名。常见做法是先做 TPM 或 FPKM 标准化再按阶段取均值。# 假设 expr 是基因 x 样本的表达矩阵meta 是样本元数据 # meta 至少包含 sample_id, species, stage 三列 expr - as.matrix(read.csv(expr_matrix.csv, row.names 1)) meta - read.csv(sample_meta.csv, stringsAsFactors FALSE) # 按物种和阶段聚合取均值 library(dplyr) stage_mean_list - list() for (sp in unique(meta$species)) { sub_meta - meta[meta$species sp, ] sub_expr - expr[, sub_meta$sample_id, drop FALSE] # 按阶段分组求均值 stage_mat - t(apply(sub_expr, 1, function(x) { tapply(x, sub_meta$stage, mean) })) colnames(stage_mat) - paste0(sp, _, colnames(stage_mat)) stage_mean_list[[sp]] - stage_mat } # 合并所有物种要求阶段顺序一致 # 这里假设每个物种阶段数相同且顺序已对齐 combined - do.call(cbind, stage_mean_list)这段代码的关键在最后一步do.call(cbind, ...)之前必须确认每个物种的阶段列顺序完全一致。如果物种 A 是stage1,stage2,stage3物种 B 是stage2,stage1,stage3合并后阶段就错位了myTAI 算出来的进化速率全是错的而且不会报错——这是最阴的坑。参数上tapply里的mean可以换成median但 myTAI 官方示例用的是均值我一般也跟均值除非某个阶段有极端离群样本。2.2 系统发育树要和物种名严格对应myTAI 需要一棵有分支长度的系统发育树格式是 Newick。树上的物种名必须和表达矩阵里的物种前缀完全一致大小写、下划线都不能差。我一般用ape包读树并检查。library(ape) tree - read.tree(species_tree.nwk) # 检查树上的物种名 print(tree$tip.label) # 假设表达矩阵物种前缀是 spA, spB, spC # 树上的 tip.label 必须完全一致 # 如果不一致手动改 tree$tip.label - c(spA, spB, spC) # 按实际顺序改 # 检查树是否有分支长度 if (is.null(tree$edge.length)) { stop(树没有分支长度myTAI 无法计算进化速率) }没有分支长度的树是另一个高频翻车点。myTAI 的进化速率计算依赖分支长度做加权如果树是 ultrametric 但没有edge.length得先用chronos或compute.brlen补上。我一般用compute.brlen(tree, method Grafen)Grafen 方法生成的 branch length 在 evo-devo 里比较常用但如果你有真实的分化时间直接用时间树更好。2.3 构建 PhyloExpressionSet 对象并做第一轮质控三样东西齐了阶段均值矩阵、树、阶段与物种的映射表。myTAI 的PhyloExpressionSet构造函数需要这些。library(myTAI) # 阶段映射每个物种的阶段列在 combined 里的位置 # 假设 combined 列顺序是 spA_stage1, spA_stage2, ..., spC_stage6 # 需要告诉 myTAI 哪些列属于哪个物种、哪个阶段 # 常见做法是构造一个 data.frame stage_map - data.frame( species rep(c(spA, spB, spC), each 6), stage rep(1:6, times 3), column_index 1:18 ) # 构建对象 phy_set - PhyloExpressionSet( expr_matrix combined, tree tree, stages stage_map$stage, species stage_map$species ) # 第一轮质控检查是否有 NA、负值、全零基因 sum(is.na(phy_set$expr_matrix)) sum(phy_set$expr_matrix 0, na.rm TRUE) zero_genes - which(rowSums(phy_set$expr_matrix) 0) length(zero_genes)质控这步别省。我遇到过表达矩阵里有负值的情况——上游用 log2 转换时没处理好myTAI 内部算几何均值会直接出 NaN然后整个 TDI 曲线全是缺口。全零基因也要去掉否则富集分析时会被当成显著基因结果全是噪音。参数上PhyloExpressionSet还有个groups参数如果你有多个物种分组比如脊椎动物和无脊椎动物可以传进去后面做分组比较时用得上。3. 转录组进化速率怎么算TDI、TAI 与阶段富集数据准备好了这一章讲 myTAI 最核心的三个计算TDITranscriptome Divergence Index、TAITranscriptome Age Index和阶段富集。这三个指标是 myTAI 的招牌也是审稿人最爱问的。3.1 TDI 的计算逻辑与参数选择TDI 衡量的是每个发育阶段上基因表达谱在系统发育树上的分化程度。数值越高说明这个阶段由进化上更年轻的基因主导数值越低说明保守基因占主导。myTAI 里用TDI()函数算。# 计算 TDI tdi_result - TDI(phy_set, measure mean) # 查看结果 head(tdi_result) # 画 TDI 曲线 plot(tdi_result, type l, col blue, xlab Developmental Stage, ylab TDI, main Transcriptome Divergence Index)measure参数有两个选项mean和median。我一般用mean因为 myTAI 的原始论文和大多数已发表分析都用均值。但如果你数据里某个阶段的基因表达有极端值median更稳健。另一个隐藏参数是log默认TRUE会对表达矩阵做 log 转换。如果你的数据已经是 log 尺度记得设log FALSE否则会双重 logTDI 曲线形状会完全扭曲——这个坑我踩过当时曲线看起来“太平了”查了半天才发现是 log 套 log。TDI 的统计显著性用置换检验做# 置换检验n 是置换次数 tdi_perm - TDI(phy_set, measure mean, permutations 1000) # 提取 p 值 tdi_perm$p.value置换次数 1000 是常见做法但如果你要发文章审稿人可能要求 10000。计算时间会线性增长1000 次大概几分钟10000 次就得半小时以上。我一般先跑 1000 次看趋势确定有信号后再补到 10000。3.2 TAI 与 TDI 的区别一个看年龄一个看分化TAI 和 TDI 经常被混用但逻辑不同。TAI 衡量的是每个阶段表达的基因在系统发育树上的平均进化年龄数值越低说明越保守的基因在主导。myTAI 里用TAI()算。# 计算 TAI tai_result - TAI(phy_set, measure mean) # 同时画 TDI 和 TAI 对比 par(mfrow c(1, 2)) plot(tdi_result, type l, col blue, main TDI) plot(tai_result, type l, col red, main TAI) par(mfrow c(1, 1))TAI 和 TDI 通常呈负相关保守基因主导的阶段 TAI 低、TDI 也低年轻基因主导的阶段 TAI 高、TDI 高。但如果你的数据里两者不相关甚至正相关说明进化模式不典型可能是物种选择有问题或者阶段划分太粗。我一般会把两个曲线叠在一起看如果出现“交叉”模式反而可能是生物学上有意思的信号——比如某个阶段保守基因和年轻基因同时活跃说明这个阶段经历了调控重布线。3.3 阶段富集分析哪些基因在特定阶段被“选中”算完 TDI/TAI下一步是找具体基因。myTAI 的PlotEnrichment()和TE()函数可以做阶段富集。# 阶段富集找出在每个阶段显著高表达的基因 enrich_result - PlotEnrichment(phy_set, stages 1:6) # 提取每个阶段的富集基因 stage_genes - TE(phy_set, stages 1:6) # 查看第 3 阶段的富集基因 head(stage_genes[[3]])PlotEnrichment的stages参数指定要分析的阶段可以传单个数字或向量。TE函数返回的是一个列表每个元素对应一个阶段的富集基因名。这里有个参数p.adjust.method默认是BH我一般不改但如果你阶段数少比如只有 3 个bonferroni更保守。富集结果出来后别直接拿基因名去跑 GO——先检查一下富集基因的数量如果某个阶段富集了 2000 个基因那基本是阈值太松得回去调TE的p.value参数。4. 避坑与排查myTAI 跑不通时先看这五条myTAI 的报错信息经常很模糊这一章把我遇到过的五个高频问题按“现象 → 原因 → 解决”列出来照着排查能省不少时间。4.1 现象TDI 曲线全是 NA 或缺口原因表达矩阵里有 NA、负值或全零基因。myTAI 内部算几何均值时任何非正值都会导致 NaN然后整条曲线断掉。解决构建PhyloExpressionSet之前先做三件事——sum(is.na(expr))检查 NAsum(expr 0)检查负值rowSums(expr) 0去掉全零基因。NA 可以用impute.knn补但如果是整个阶段缺失补了也没意义直接去掉那个阶段。4.2 现象报错Error in PhyloExpressionSet: tree and expression matrix species mismatch原因树上的tip.label和表达矩阵的物种前缀不一致。常见情况是树上有 5 个物种矩阵只有 3 个或者名字里一个用下划线一个用点。解决print(tree$tip.label)和unique(stage_map$species)对比手动改到完全一致。如果树上的物种比矩阵多用drop.tip()剪掉多余的。# 剪掉树上多余的物种 tree - drop.tip(tree, tip c(spD, spE))4.3 现象置换检验跑了一小时还没结束原因permutations设太大或者基因数太多比如 20000。myTAI 的置换是逐基因做的基因数乘以置换次数就是计算量。解决先过滤低表达基因保留 TPM 1 的基因通常能砍掉一半。置换次数先跑 100 次看趋势确定有信号再补到 1000。如果还是慢用parallel包开多核但 myTAI 本身不支持并行得手动分基因块跑再合并。4.4 现象TAI 和 TDI 曲线完全重合原因measure参数设错了或者表达矩阵没有做 log 转换但log参数没关。TAI 和 TDI 在数学上不同重合说明计算逻辑被覆盖了。解决检查TAI()和TDI()调用时的measure和log参数是否一致。如果数据已经 log 过两个函数都要设log FALSE。另外如果树的分支长度全是 1比如compute.brlen默认生成等长分支TAI 和 TDI 会趋同得换Grafen方法重新生成分支长度。4.5 现象富集分析结果里基因名全是gene_1、gene_2原因表达矩阵的row.names丢了。读 CSV 时如果第一列不是行名R 会自动生成gene_1这种默认名。解决读矩阵时加row.names 1或者手动赋值rownames(expr) - expr$gene_id。myTAI 的富集结果直接返回行名行名丢了后面没法做 GO 或 KEGG。# 读矩阵时保留行名 expr - read.csv(expr_matrix.csv, row.names 1, check.names FALSE)check.names FALSE也很关键否则 R 会把基因名里的特殊字符比如-、:自动替换成点后面匹配注释文件时全对不上。5. 把 myTAI 接进完整流程从原始 counts 到进化结论myTAI 不是孤立工具它得嵌在 RNA-seq 流程里。这一章讲怎么把上游定量结果接进来以及一个我常用的验证技巧。5.1 上游定量结果的接入方式如果你用 salmon 或 kallisto 做定量得到的是 TPM 矩阵可以直接用。如果是 featureCounts 的 counts得先做标准化。我一般用edgeR的 TMM 标准化再转 TPM。library(edgeR) # counts 矩阵 counts - as.matrix(read.csv(counts.csv, row.names 1)) # 构建 DGEList dge - DGEList(counts counts) # TMM 标准化 dge - calcNormFactors(dge, method TMM) # 转 TPM tpm - cpm(dge, log FALSE) # 后续按阶段聚合TMM 标准化对比较转录组很重要因为不同物种的 RNA 组成差异大不标准化直接比 TPM 会有系统偏差。calcNormFactors的method参数默认是TMM我一般不改但如果你的数据里有很多高表达基因比如 rRNA 残留可以试RLE。5.2 用 bootstrap 验证 TDI 曲线的稳健性TDI 曲线跑出来后怎么知道它稳不稳我一般做 bootstrap对每个阶段的样本做有放回重抽样重算 TDI看曲线的置信区间。# bootstrap 验证 TDI set.seed(123) n_boot - 100 tdi_boot - matrix(NA, nrow n_boot, ncol 6) for (i in 1:n_boot) { # 对每个阶段的样本索引做重抽样 boot_idx - sample(1:ncol(combined), ncol(combined), replace TRUE) boot_expr - combined[, boot_idx] # 重新构建对象并算 TDI boot_set - PhyloExpressionSet( expr_matrix boot_expr, tree tree, stages stage_map$stage[boot_idx], species stage_map$species[boot_idx] ) tdi_boot[i, ] - TDI(boot_set, measure mean) } # 计算 95% 置信区间 tdi_ci - apply(tdi_boot, 2, quantile, probs c(0.025, 0.975))set.seed必须设否则每次跑结果不一样审稿人复现不了。n_boot我一般用 100够看趋势了要发文章可以加到 1000。置信区间如果某个阶段特别宽说明那个阶段的样本异质性大得回去检查是不是有离群样本。5.3 一个我常用的技巧把 TDI 和差异表达结果叠起来看单独看 TDI 曲线只能知道“哪个阶段进化速率高”但不知道“哪些基因驱动了它”。我一般把 TDI 和差异表达结果叠起来先算每个阶段的 TDI再对每个阶段做差异分析取交集。# 假设 de_result 是某个阶段的差异基因列表 # tdi_stage 是该阶段的 TDI 值 # 找出 TDI 高且差异显著的基因 high_tdi_stage - which(tdi_result$TDI max(tdi_result$TDI)) de_genes - de_result[[high_tdi_stage]]$gene_id # 取交集 driver_genes - intersect(rownames(combined), de_genes) length(driver_genes)这个交集基因集就是“驱动该阶段进化速率”的候选。我一般会再跑一次 GO 富集看这些基因集中在哪些功能上。如果富集到“转录调控”或“信号转导”说明这个阶段的进化重布线主要在调控层面如果富集到“代谢”说明是代谢通路在变。这个技巧不复杂但能把 myTAI 的宏观指标落到具体基因上写文章时故事更完整。最后说个习惯我每次跑 myTAI 之前会先把sessionInfo()存下来因为 myTAI 依赖ape、zoo这些包版本不同结果可能有细微差异。有一次审稿人复现我的 TDI 曲线形状差了一点查到最后是ape版本不同导致树的分支长度处理有差异。从那以后我所有分析都带sessionInfo()省得后悔药没处买。希望帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/16 15:34:32 · 阅读 1.2k myTAI实战：转录组进化速率计算与阶段富集分析避坑指南
简介myTAI 是一款面向进化转录组学研究的 R 语言分析工具包适合具备一定 R 与 Bioconductor 基础的科研人员、研究生及 evo-devo 方向学者使用。它通过量化转录组保守模式及其对应的生物学过程帮助用户在计算机上初步筛选感兴趣过程中可能存在的进化限制为后续实验提供探索性依据。资源包共 330 个文件包含 95 个 R 源码、74 张 png 图表、71 个 html 文档、53 个 rd 帮助文件以及 Rmd 教程、cpp 源码、rda 数据与 yml 配置等压缩包约 9.32MB覆盖从函数实现、示例数据到文档说明的完整结构。目前已有 265 人学习。借助标准化、自动化的分析框架读者可快速上手转录组保守性检测流程理解基因集与表达谱的进化约束分析思路并参考官方文档与引文信息完成可复现的研究实践。1. myTAI:R进化转录组学到底在算什么如果你手上有一批不同物种、不同发育阶段的转录组数据想回答“哪些基因在进化上被保守地调控、哪些在某个谱系里被重新布线”myTAI 就是 R 生态里专门干这件事的包。它不负责比对、不负责定量输入是已经算好的表达矩阵加一份物种/发育阶段的映射表输出是转录组进化速率、保守性打分、阶段间表达位移这些指标。我第一次接触它是在一个跨物种胚胎发育的项目里当时用 DESeq2 做完差异分析就卡住了——差异基因列表有了但“这些差异在进化尺度上意味着什么”没人能回答。myTAI 补的就是这一段把表达矩阵映射到系统发育树上用进化模型算出每个基因的转录组进化速率TDI、TAI 这类指标再按发育阶段做富集。适合的人群很明确做比较转录组、evo-devo、发育生物学且已经能熟练跑 R 和基本 RNA-seq 流程的人。新手直接上会卡在“表达矩阵怎么对齐到发育阶段”这一步后面会细讲。2. 把表达矩阵喂进 myTAI 之前要做的三件事myTAI 对输入格式挑剔不是随便一个 counts 矩阵就能塞进去。我见过太多人卡在Error: PhyloExpressionSet报错上根因基本都是输入没对齐。这一章把数据准备讲透后面跑分析才不会反复翻车。2.1 表达矩阵的列名必须能映射到发育阶段myTAI 的核心对象叫PhyloExpressionSet它要求表达矩阵的列不是样本名而是发育阶段。也就是说如果你有 3 个物种、每个物种 6 个发育阶段、每个阶段 2 个重复你得先把重复合并成阶段均值列名变成类似stage1到stage6并且所有物种共用同一套阶段命名。常见做法是先做 TPM 或 FPKM 标准化再按阶段取均值。# 假设 expr 是基因 x 样本的表达矩阵meta 是样本元数据 # meta 至少包含 sample_id, species, stage 三列 expr - as.matrix(read.csv(expr_matrix.csv, row.names 1)) meta - read.csv(sample_meta.csv, stringsAsFactors FALSE) # 按物种和阶段聚合取均值 library(dplyr) stage_mean_list - list() for (sp in unique(meta$species)) { sub_meta - meta[meta$species sp, ] sub_expr - expr[, sub_meta$sample_id, drop FALSE] # 按阶段分组求均值 stage_mat - t(apply(sub_expr, 1, function(x) { tapply(x, sub_meta$stage, mean) })) colnames(stage_mat) - paste0(sp, _, colnames(stage_mat)) stage_mean_list[[sp]] - stage_mat } # 合并所有物种要求阶段顺序一致 # 这里假设每个物种阶段数相同且顺序已对齐 combined - do.call(cbind, stage_mean_list)这段代码的关键在最后一步do.call(cbind, ...)之前必须确认每个物种的阶段列顺序完全一致。如果物种 A 是stage1,stage2,stage3物种 B 是stage2,stage1,stage3合并后阶段就错位了myTAI 算出来的进化速率全是错的而且不会报错——这是最阴的坑。参数上tapply里的mean可以换成median但 myTAI 官方示例用的是均值我一般也跟均值除非某个阶段有极端离群样本。2.2 系统发育树要和物种名严格对应myTAI 需要一棵有分支长度的系统发育树格式是 Newick。树上的物种名必须和表达矩阵里的物种前缀完全一致大小写、下划线都不能差。我一般用ape包读树并检查。library(ape) tree - read.tree(species_tree.nwk) # 检查树上的物种名 print(tree$tip.label) # 假设表达矩阵物种前缀是 spA, spB, spC # 树上的 tip.label 必须完全一致 # 如果不一致手动改 tree$tip.label - c(spA, spB, spC) # 按实际顺序改 # 检查树是否有分支长度 if (is.null(tree$edge.length)) { stop(树没有分支长度myTAI 无法计算进化速率) }没有分支长度的树是另一个高频翻车点。myTAI 的进化速率计算依赖分支长度做加权如果树是 ultrametric 但没有edge.length得先用chronos或compute.brlen补上。我一般用compute.brlen(tree, method Grafen)Grafen 方法生成的 branch length 在 evo-devo 里比较常用但如果你有真实的分化时间直接用时间树更好。2.3 构建 PhyloExpressionSet 对象并做第一轮质控三样东西齐了阶段均值矩阵、树、阶段与物种的映射表。myTAI 的PhyloExpressionSet构造函数需要这些。library(myTAI) # 阶段映射每个物种的阶段列在 combined 里的位置 # 假设 combined 列顺序是 spA_stage1, spA_stage2, ..., spC_stage6 # 需要告诉 myTAI 哪些列属于哪个物种、哪个阶段 # 常见做法是构造一个 data.frame stage_map - data.frame( species rep(c(spA, spB, spC), each 6), stage rep(1:6, times 3), column_index 1:18 ) # 构建对象 phy_set - PhyloExpressionSet( expr_matrix combined, tree tree, stages stage_map$stage, species stage_map$species ) # 第一轮质控检查是否有 NA、负值、全零基因 sum(is.na(phy_set$expr_matrix)) sum(phy_set$expr_matrix 0, na.rm TRUE) zero_genes - which(rowSums(phy_set$expr_matrix) 0) length(zero_genes)质控这步别省。我遇到过表达矩阵里有负值的情况——上游用 log2 转换时没处理好myTAI 内部算几何均值会直接出 NaN然后整个 TDI 曲线全是缺口。全零基因也要去掉否则富集分析时会被当成显著基因结果全是噪音。参数上PhyloExpressionSet还有个groups参数如果你有多个物种分组比如脊椎动物和无脊椎动物可以传进去后面做分组比较时用得上。3. 转录组进化速率怎么算TDI、TAI 与阶段富集数据准备好了这一章讲 myTAI 最核心的三个计算TDITranscriptome Divergence Index、TAITranscriptome Age Index和阶段富集。这三个指标是 myTAI 的招牌也是审稿人最爱问的。3.1 TDI 的计算逻辑与参数选择TDI 衡量的是每个发育阶段上基因表达谱在系统发育树上的分化程度。数值越高说明这个阶段由进化上更年轻的基因主导数值越低说明保守基因占主导。myTAI 里用TDI()函数算。# 计算 TDI tdi_result - TDI(phy_set, measure mean) # 查看结果 head(tdi_result) # 画 TDI 曲线 plot(tdi_result, type l, col blue, xlab Developmental Stage, ylab TDI, main Transcriptome Divergence Index)measure参数有两个选项mean和median。我一般用mean因为 myTAI 的原始论文和大多数已发表分析都用均值。但如果你数据里某个阶段的基因表达有极端值median更稳健。另一个隐藏参数是log默认TRUE会对表达矩阵做 log 转换。如果你的数据已经是 log 尺度记得设log FALSE否则会双重 logTDI 曲线形状会完全扭曲——这个坑我踩过当时曲线看起来“太平了”查了半天才发现是 log 套 log。TDI 的统计显著性用置换检验做# 置换检验n 是置换次数 tdi_perm - TDI(phy_set, measure mean, permutations 1000) # 提取 p 值 tdi_perm$p.value置换次数 1000 是常见做法但如果你要发文章审稿人可能要求 10000。计算时间会线性增长1000 次大概几分钟10000 次就得半小时以上。我一般先跑 1000 次看趋势确定有信号后再补到 10000。3.2 TAI 与 TDI 的区别一个看年龄一个看分化TAI 和 TDI 经常被混用但逻辑不同。TAI 衡量的是每个阶段表达的基因在系统发育树上的平均进化年龄数值越低说明越保守的基因在主导。myTAI 里用TAI()算。# 计算 TAI tai_result - TAI(phy_set, measure mean) # 同时画 TDI 和 TAI 对比 par(mfrow c(1, 2)) plot(tdi_result, type l, col blue, main TDI) plot(tai_result, type l, col red, main TAI) par(mfrow c(1, 1))TAI 和 TDI 通常呈负相关保守基因主导的阶段 TAI 低、TDI 也低年轻基因主导的阶段 TAI 高、TDI 高。但如果你的数据里两者不相关甚至正相关说明进化模式不典型可能是物种选择有问题或者阶段划分太粗。我一般会把两个曲线叠在一起看如果出现“交叉”模式反而可能是生物学上有意思的信号——比如某个阶段保守基因和年轻基因同时活跃说明这个阶段经历了调控重布线。3.3 阶段富集分析哪些基因在特定阶段被“选中”算完 TDI/TAI下一步是找具体基因。myTAI 的PlotEnrichment()和TE()函数可以做阶段富集。# 阶段富集找出在每个阶段显著高表达的基因 enrich_result - PlotEnrichment(phy_set, stages 1:6) # 提取每个阶段的富集基因 stage_genes - TE(phy_set, stages 1:6) # 查看第 3 阶段的富集基因 head(stage_genes[[3]])PlotEnrichment的stages参数指定要分析的阶段可以传单个数字或向量。TE函数返回的是一个列表每个元素对应一个阶段的富集基因名。这里有个参数p.adjust.method默认是BH我一般不改但如果你阶段数少比如只有 3 个bonferroni更保守。富集结果出来后别直接拿基因名去跑 GO——先检查一下富集基因的数量如果某个阶段富集了 2000 个基因那基本是阈值太松得回去调TE的p.value参数。4. 避坑与排查myTAI 跑不通时先看这五条myTAI 的报错信息经常很模糊这一章把我遇到过的五个高频问题按“现象 → 原因 → 解决”列出来照着排查能省不少时间。4.1 现象TDI 曲线全是 NA 或缺口原因表达矩阵里有 NA、负值或全零基因。myTAI 内部算几何均值时任何非正值都会导致 NaN然后整条曲线断掉。解决构建PhyloExpressionSet之前先做三件事——sum(is.na(expr))检查 NAsum(expr 0)检查负值rowSums(expr) 0去掉全零基因。NA 可以用impute.knn补但如果是整个阶段缺失补了也没意义直接去掉那个阶段。4.2 现象报错Error in PhyloExpressionSet: tree and expression matrix species mismatch原因树上的tip.label和表达矩阵的物种前缀不一致。常见情况是树上有 5 个物种矩阵只有 3 个或者名字里一个用下划线一个用点。解决print(tree$tip.label)和unique(stage_map$species)对比手动改到完全一致。如果树上的物种比矩阵多用drop.tip()剪掉多余的。# 剪掉树上多余的物种 tree - drop.tip(tree, tip c(spD, spE))4.3 现象置换检验跑了一小时还没结束原因permutations设太大或者基因数太多比如 20000。myTAI 的置换是逐基因做的基因数乘以置换次数就是计算量。解决先过滤低表达基因保留 TPM 1 的基因通常能砍掉一半。置换次数先跑 100 次看趋势确定有信号再补到 1000。如果还是慢用parallel包开多核但 myTAI 本身不支持并行得手动分基因块跑再合并。4.4 现象TAI 和 TDI 曲线完全重合原因measure参数设错了或者表达矩阵没有做 log 转换但log参数没关。TAI 和 TDI 在数学上不同重合说明计算逻辑被覆盖了。解决检查TAI()和TDI()调用时的measure和log参数是否一致。如果数据已经 log 过两个函数都要设log FALSE。另外如果树的分支长度全是 1比如compute.brlen默认生成等长分支TAI 和 TDI 会趋同得换Grafen方法重新生成分支长度。4.5 现象富集分析结果里基因名全是gene_1、gene_2原因表达矩阵的row.names丢了。读 CSV 时如果第一列不是行名R 会自动生成gene_1这种默认名。解决读矩阵时加row.names 1或者手动赋值rownames(expr) - expr$gene_id。myTAI 的富集结果直接返回行名行名丢了后面没法做 GO 或 KEGG。# 读矩阵时保留行名 expr - read.csv(expr_matrix.csv, row.names 1, check.names FALSE)check.names FALSE也很关键否则 R 会把基因名里的特殊字符比如-、:自动替换成点后面匹配注释文件时全对不上。5. 把 myTAI 接进完整流程从原始 counts 到进化结论myTAI 不是孤立工具它得嵌在 RNA-seq 流程里。这一章讲怎么把上游定量结果接进来以及一个我常用的验证技巧。5.1 上游定量结果的接入方式如果你用 salmon 或 kallisto 做定量得到的是 TPM 矩阵可以直接用。如果是 featureCounts 的 counts得先做标准化。我一般用edgeR的 TMM 标准化再转 TPM。library(edgeR) # counts 矩阵 counts - as.matrix(read.csv(counts.csv, row.names 1)) # 构建 DGEList dge - DGEList(counts counts) # TMM 标准化 dge - calcNormFactors(dge, method TMM) # 转 TPM tpm - cpm(dge, log FALSE) # 后续按阶段聚合TMM 标准化对比较转录组很重要因为不同物种的 RNA 组成差异大不标准化直接比 TPM 会有系统偏差。calcNormFactors的method参数默认是TMM我一般不改但如果你的数据里有很多高表达基因比如 rRNA 残留可以试RLE。5.2 用 bootstrap 验证 TDI 曲线的稳健性TDI 曲线跑出来后怎么知道它稳不稳我一般做 bootstrap对每个阶段的样本做有放回重抽样重算 TDI看曲线的置信区间。# bootstrap 验证 TDI set.seed(123) n_boot - 100 tdi_boot - matrix(NA, nrow n_boot, ncol 6) for (i in 1:n_boot) { # 对每个阶段的样本索引做重抽样 boot_idx - sample(1:ncol(combined), ncol(combined), replace TRUE) boot_expr - combined[, boot_idx] # 重新构建对象并算 TDI boot_set - PhyloExpressionSet( expr_matrix boot_expr, tree tree, stages stage_map$stage[boot_idx], species stage_map$species[boot_idx] ) tdi_boot[i, ] - TDI(boot_set, measure mean) } # 计算 95% 置信区间 tdi_ci - apply(tdi_boot, 2, quantile, probs c(0.025, 0.975))set.seed必须设否则每次跑结果不一样审稿人复现不了。n_boot我一般用 100够看趋势了要发文章可以加到 1000。置信区间如果某个阶段特别宽说明那个阶段的样本异质性大得回去检查是不是有离群样本。5.3 一个我常用的技巧把 TDI 和差异表达结果叠起来看单独看 TDI 曲线只能知道“哪个阶段进化速率高”但不知道“哪些基因驱动了它”。我一般把 TDI 和差异表达结果叠起来先算每个阶段的 TDI再对每个阶段做差异分析取交集。# 假设 de_result 是某个阶段的差异基因列表 # tdi_stage 是该阶段的 TDI 值 # 找出 TDI 高且差异显著的基因 high_tdi_stage - which(tdi_result$TDI max(tdi_result$TDI)) de_genes - de_result[[high_tdi_stage]]$gene_id # 取交集 driver_genes - intersect(rownames(combined), de_genes) length(driver_genes)这个交集基因集就是“驱动该阶段进化速率”的候选。我一般会再跑一次 GO 富集看这些基因集中在哪些功能上。如果富集到“转录调控”或“信号转导”说明这个阶段的进化重布线主要在调控层面如果富集到“代谢”说明是代谢通路在变。这个技巧不复杂但能把 myTAI 的宏观指标落到具体基因上写文章时故事更完整。最后说个习惯我每次跑 myTAI 之前会先把sessionInfo()存下来因为 myTAI 依赖ape、zoo这些包版本不同结果可能有细微差异。有一次审稿人复现我的 TDI 曲线形状差了一点查到最后是ape版本不同导致树的分支长度处理有差异。从那以后我所有分析都带sessionInfo()省得后悔药没处买。希望帮到你。本文还有配套的精品资源点击获取