简介myTAI 是一款面向进化转录组学研究的 R 语言分析工具包适合从事 evo-devo、基因表达保守性与生物学过程进化限制探索的科研人员和生物信息学学习者。它提供标准化、自动化的分析框架通过量化转录组保守模式及其潜在基因组基础帮助用户在计算机上初步筛选感兴趣生物学过程中可能存在的进化约束信号可用于基因集、表达谱与转录组数据的探索性分析。资源包共 330 个文件约 9.32MB包含 95 个 R 源码文件、74 张 png 图表、71 个 html 文档、53 个 rd 帮助文件以及 Rmd 教程、cpp 源码、rda 数据与 yml 配置等覆盖核心函数、示例数据、可视化结果与文档说明。已有 265 人学习关注。使用者可借助其完整源码与文档快速理解 TAI 计算流程、复现分析示例并参考 Drost 等 2018 年 Bioinformatics 论文的引用信息规范开展研究适合作为进化转录组学入门与进阶的实用工具包。1. myTAI 到底解决什么问题从基因集到进化转录组学的落地路径如果你手头有一批不同物种、不同发育阶段的转录组数据想回答“哪些生物学过程在进化上更保守”这类问题myTAI 就是那个绕不开的 R 包。它专门做进化转录组学evolutionary transcriptomics分析核心思路是把表达谱映射到基因集gene-sets上再计算进化保守性指标比如 TAI转录组年龄指数和表达保守性分数。很多做 evo-devo 的研究者第一次接触它是因为想比较不同物种同一发育阶段的表达模式或者想看某个生物学过程在系统发育树上的表达漂移程度。myTAI 把这类分析封装成了可复现的 R 工作流从数据预处理、基因集富集到可视化都有现成函数。适合谁适合已经能用 R 做基础转录组分析、想往进化方向延伸的从业者也适合做比较转录组但不想自己从头写统计脚本的人。它的门槛不在 R 语法而在对输入数据格式和基因集定义的理解——这两处翻车最多。2. 环境搭建与数据准备把表达矩阵喂进 myTAI 之前2.1 安装与依赖检查myTAI 在 CRAN 和 Bioconductor 都有依赖常见做法是先装 Bioconductor 管理器再装 myTAI。我一般会固定一个 R 版本避免因为 Bioconductor 版本错位导致依赖装不上。# 安装 Bioconductor 管理器 if (!requireNamespace(BiocManager, quietly TRUE)) install.packages(BiocManager) # 安装 myTAI 及其核心依赖 BiocManager::install(myTAI) # 检查是否装好并看版本 library(myTAI) packageVersion(myTAI)逻辑说明BiocManager::install会自动解析 myTAI 在 Bioconductor 上的依赖树包括edgeR、dplyr等。参数上不需要额外指定 repos除非你用了内部镜像。装完后packageVersion能返回版本号如果报错说找不到包多半是 Bioconductor 版本和 R 版本不匹配这时候要么升 R要么用BiocManager::install(version 对应版本)降级。2.2 表达矩阵的格式要求myTAI 对输入表达矩阵有明确要求行是基因列是样本且列名要能对应到发育阶段或物种。常见做法是准备一个data.frame或matrix数值为 TPM 或 FPKM不要用 raw counts 直接进除非你先做归一化。我一般会先做 log2 转换让表达分布更接近正态后续计算 TAI 时更稳。# 假设 expr_mat 是已经归一化并 log2 转换的表达矩阵 # 行名是基因 ID列名是样本名 expr_mat - as.matrix(read.csv(expression_matrix.csv, row.names 1)) # 检查维度基因数 x 样本数 dim(expr_mat) # 检查是否有 NAmyTAI 对 NA 很敏感 sum(is.na(expr_mat)) # 如果有 NA用行均值填充或直接删基因 expr_mat - expr_mat[complete.cases(expr_mat), ]逻辑说明read.csv读进来默认是 data.frame用as.matrix转成矩阵因为 myTAI 内部很多函数要求矩阵输入。dim返回的行列数要和你预期一致比如 20000 个基因、12 个样本。sum(is.na())如果大于 0后面计算 TAI 会直接报错所以要么删基因要么用apply(expr_mat, 1, function(x) ifelse(is.na(x), mean(x, na.rmTRUE), x))填充。参数上row.names 1表示第一列是基因 ID这个别搞错否则行名会变成数字。2.3 基因集文件的准备myTAI 做富集分析时需要一个基因集列表通常是list结构每个元素是一个基因 ID 向量。你可以从 GO 或 KEGG 下载也可以自己定义。常见做法是用clusterProfiler的read.gmt读入 GMT 文件再转成 list。# 读入 GMT 格式的基因集文件 gene_sets - clusterProfiler::read.gmt(c2.all.v2023.1.Hs.symbols.gmt) # 转成 myTAI 需要的 list 结构 gene_set_list - split(gene_sets$gene, gene_sets$term) # 检查每个基因集的基因数 sapply(gene_set_list, length)[1:5]逻辑说明read.gmt返回的是长表split按 term 拆成 list。sapply看前五个基因集的基因数如果某个集合只有一两个基因后续富集结果会不稳定建议过滤掉少于 10 个基因的集合。参数上GMT 文件的物种要和你的表达矩阵一致人鼠混用会直接导致基因 ID 匹配不上。3. 核心分析流程TAI 计算与表达保守性3.1 转录组年龄指数 TAI 的计算TAI 是 myTAI 的招牌指标它把每个基因按进化年龄分层再算每个样本的加权平均。年龄分层通常用系统发育树上的节点来定义比如“后口动物特有”“脊椎动物特有”等。myTAI 提供了TAI()函数但前提是你得先有一个phylo对象和基因年龄映射表。# 假设 phylo_tree 是已经读入的系统发育树 # gene_age 是一个 data.frame两列gene_id 和 age_class # age_class 是因子按进化年龄从老到新排序 # 计算 TAI tai_result - TAI(expr_mat, phylo_tree, gene_age) # 查看前几个样本的 TAI 值 head(tai_result$TAI)逻辑说明TAI()内部会先检查基因 ID 是否匹配然后按年龄类别计算每个样本的加权平均。phylo_tree必须是ape包能识别的phylo对象gene_age的age_class顺序决定了加权方向。如果报错说“gene IDs not found”多半是表达矩阵行名和 gene_age 的 gene_id 格式不一致比如一个用 Ensembl ID一个用 Symbol。参数上TAI()还有一个log参数默认TRUE表示输入已经 log 转换如果你没转要设成FALSE。3.2 表达保守性分数与可视化算完 TAI 后通常要看不同发育阶段的 TAI 变化趋势以及哪些基因集在进化上更保守。myTAI 提供了PlotSignature()和PlotEnrichment()等函数。我一般会先画 TAI 随发育阶段的箱线图再看富集结果。# 画 TAI 随发育阶段的箱线图 PlotSignature(tai_result, xlab Developmental Stage, ylab TAI, main TAI across development) # 富集分析看哪些基因集在特定阶段显著 enrich_result - PlotEnrichment(expr_mat, gene_set_list, stages c(stage1, stage2, stage3)) # 查看富集结果表格 head(enrich_result$enrichment)逻辑说明PlotSignature的xlab和ylab只是标签不影响计算。PlotEnrichment的stages参数指定要比较的发育阶段必须和表达矩阵列名一致。enrich_result$enrichment返回的是每个基因集在每个阶段的富集 p 值和效应量。如果 p 值全为 1检查基因集里的基因 ID 是否和表达矩阵行名匹配这是最常见的翻车点。3.3 参数调优与结果解读TAI 计算里有一个age_class的加权方式默认是线性加权但你可以改成按节点深度加权。常见做法是先用默认跑一遍看趋势是否合理再调整。另外PlotEnrichment的p.adjust.method默认是 BH如果你样本量小可以改成none看原始 p 值但发文章时还是得用校正后的。# 调整 TAI 的加权方式 tai_result_weighted - TAI(expr_mat, phylo_tree, gene_age, weight depth) # 调整富集分析的 p 值校正方法 enrich_result_raw - PlotEnrichment(expr_mat, gene_set_list, stages c(stage1, stage2), p.adjust.method none)逻辑说明weight depth会让进化深度更大的节点权重更高适合想强调深层同源基因的场景。p.adjust.method none只建议在探索阶段用正式结果还是 BH。参数上TAI()的weight只接受linear和depth两个值写别的会报错。4. 避坑与常见问题那些让我重跑整晚的细节4.1 基因 ID 不匹配导致 TAI 全为 NA现象TAI()返回的矩阵全是 NA或者报错说“no overlapping genes”。原因表达矩阵行名是 Ensembl ID而 gene_age 里是 Symbol或者反过来。解决统一转成同一种 ID用clusterProfiler::bitr或org.Hs.eg.db做映射映射完再检查交集。4.2 发育阶段列名顺序错乱现象PlotSignature画出来的 x 轴阶段顺序是乱的比如“stage3”跑到“stage1”前面。原因表达矩阵列名是字符R 默认按字母排序。解决在画图前把列名设成因子并手动指定 levels比如factor(colnames(expr_mat), levels c(stage1, stage2, stage3))。4.3 基因集富集结果全不显著现象PlotEnrichment返回的 p 值全大于 0.05哪怕你预期某些通路应该显著。原因基因集里的基因和表达矩阵行名匹配率太低或者基因集太大导致校正后 p 值膨胀。解决先算匹配率sum(gene_set_list[[1]] %in% rownames(expr_mat)) / length(gene_set_list[[1]])低于 0.5 就换基因集版本另外过滤掉大于 500 个基因的集合。4.4 内存溢出与运行时间过长现象跑TAI()时 R 直接崩掉或者跑了半小时没动静。原因表达矩阵太大或者基因集数量太多。解决先对表达矩阵做基因过滤只保留在至少一半样本中表达大于 1 的基因基因集只保留和你的研究相关的比如只留 GO BP 里的免疫相关通路。4.5 系统发育树与基因年龄映射不一致现象TAI()报错说“tree tips not found in gene_age”。原因树的叶节点名称和 gene_age 里的物种名不一致比如树里是“Homo_sapiens”gene_age 里是“human”。解决用gsub统一命名或者直接改树的tip.label。5. 进阶技巧把 myTAI 嵌进可复现的 R Markdown 工作流5.1 用 R Markdown 固化整个分析我现在的习惯是任何 myTAI 分析都写成一个 R Markdown 文档参数放在 YAML 里这样换数据只需改路径。下面是一个最小模板。--- title: myTAI 进化转录组学分析 output: html_document params: expr_file: expression_matrix.csv gene_set_file: c2.all.v2023.1.Hs.symbols.gmt tree_file: species_tree.nwk --- {r setup, includeFALSE} library(myTAI) library(clusterProfiler) library(ape)expr_mat - as.matrix(read.csv(params$expr_file, row.names 1)) gene_sets - read.gmt(params$gene_set_file) gene_set_list - split(gene_sets$gene, gene_sets$term) phylo_tree - read.tree(params$tree_file)tai_result - TAI(expr_mat, phylo_tree, gene_age) PlotSignature(tai_result)逻辑说明YAML 里的 params 让路径可配置read.tree 来自 ape 包读 Newick 格式的树。gene_age 需要你提前准备好可以放在一个单独的 CSV 里用 read.csv 读入。这样整个分析从数据到图一键渲染换物种只需改三个文件路径。 ### 5.2 用 purrr 批量跑多个基因集 如果你有几十个基因集要分别算富集用 purrr::map 比写循环干净。 r library(purrr) # 批量富集 enrich_list - map(gene_set_list, ~ PlotEnrichment(expr_mat, list(.x), stages c(stage1, stage2))) # 提取每个基因集的显著通路数 sig_counts - map_int(enrich_list, ~ sum(.x$enrichment$p.adjust 0.05))逻辑说明map对每个基因集单独跑PlotEnrichment注意这里把单个基因集包成list(.x)因为函数要求 list 输入。map_int返回整数向量统计每个基因集里校正后 p 值小于 0.05 的通路数。参数上stages还是要和列名一致。5.3 验证方法用模拟数据检查流程在跑真实数据前我一般会先用模拟数据走一遍确认 TAI 和富集函数能正常返回。模拟数据可以用rnorm生成基因年龄随机分配。# 模拟 1000 个基因、6 个样本的表达矩阵 set.seed(123) sim_expr - matrix(rnorm(6000, mean 5, sd 2), nrow 1000) rownames(sim_expr) - paste0(gene, 1:1000) colnames(sim_expr) - paste0(stage, 1:6) # 模拟基因年龄 sim_age - data.frame( gene_id rownames(sim_expr), age_class sample(c(old, medium, young), 1000, replace TRUE) ) # 模拟树 sim_tree - rtree(10) # 跑 TAI sim_tai - TAI(sim_expr, sim_tree, sim_age) head(sim_tai$TAI)逻辑说明rnorm生成正态分布表达值sample随机分配年龄类别rtree生成随机树。TAI能跑通说明环境没问题再换真实数据。参数上set.seed保证可复现rtree(10)生成 10 个叶节点的树实际物种数要匹配。从那以后我每次跑 myTAI 之前都强制走一遍模拟数据检查确认基因 ID 匹配、列名顺序、树和年龄表一致再上真实数据。希望帮到你。本文还有配套的精品资源点击获取