基于全基因组复制（WGD）的 KEGG 功能富集及 Ka/Ks 进化分析
01 分析背景与核心目标
本分析聚焦基因复制事件中全基因组复制（WGD）类型，结合 KEGG 功能富集分析解析 WGD 基因的生物学功能特征，通过 Ka/Ks（非同义替换率 / 同义替换率）及 4DTv（四简并位点颠换率）分析揭示 WGD 基因的进化选择压力与分化特征；同时兼顾串联重复（TD）、近端重复（PD）等其他复制类型，形成多维度对比分析。
02 依赖软件与核心工具
| 软件 / 脚本 | 功能用途 |
|---|---|
| blastp | 蛋白质序列比对，为基因复制类型鉴定提供同源性依据 |
| DupGen-finder | 基因复制类型分类（核心工具，将基因分为 WGD/TD/PD/TRD/DSD/SL 六类） |
| KaKs_Calculator/ParaAT.pl | 计算基因对的 Ka、Ks 值，解析进化选择压力 |
| mafft/muscle | 序列比对（推荐 muscle，速度更快） |
| R (ggplot2/clusterProfiler) | KEGG/GO 富集分析及可视化（气泡图、小提琴图、密度图） |
| 辅助脚本 | axt2one-line.py、calculate_4DTV_correction.pl（计算 4DTv 值） |
03 基因复制类型定义
| 复制类型 | 英文全称 | 定义 |
|---|---|---|
| WGD | Whole-genome duplication | 全基因组复制（核心分析对象） |
| TD | Tandem duplication | 串联重复（相邻的两个重复基因） |
| PD | Proximal duplication | 近端重复（相隔 10 个以内基因的重复基因） |
| TRD | Transposed duplication | 转置重复（祖先和新基因座组成的重复基因） |
| DSD | Dispersed duplication | 分散重复（不相邻也不共线性的重复基因） |
| SL | Singleton | 单拷贝（无同源重复基因） |
04 输入文件要求
| 文件类型 | 文件名示例 | 格式要求 |
|---|---|---|
| 基因组注释 | Spd.gff/Ath.gff | 制表符分隔，格式：染色体ID\t基因ID\t起始位置\t终止位置（示例：Ath-Chr1 AT1G01010.1 3630 5899） |
| 蛋白质序列 | Spd.pep/Ath.pep | FASTA 格式，序列 ID 需与 gff 文件中基因 ID 完全一致（示例：>AT3G05780.1 + 氨基酸序列） |
| 扩张基因列表 | expansion.genes | 单列基因 ID，需与上述文件的基因 ID 命名规则统一 |
| 功能注释库 | GO.info/KEGG.info | 制表符分隔，GO.info：GOID\t功能描述\t基因ID\t功能分类；KEGG.info：KOID\t通路名称\t基因ID |
05 分析流程（分模块梳理，聚焦 WGD）
模块 1：基因复制类型鉴定（DupGen-finder 核心步骤）
1.1 输入文件预处理（统一格式）
# 物种1（Spd）gff文件格式化（添加物种前缀，调整列顺序） cat Spd.bed | sed 's/^/Spd-/g' | awk '{print $1"\t"$4"\t"$2"\t"$3}' > Spd.gff sed -i 's/Chr0/Chr/g' Spd.gff # 统一染色体命名格式 # 物种2（Ath）gff文件格式化（可选，跨物种对比用） cat Ath.bed | sed 's/^/Ath-Chr/g' | awk '{print $1"\t"$4"\t"$2"\t"$3}' > Ath.gff # 合并文件（跨物种分析用，单物种分析仅需Spd.gff） cat Spd.gff Ath.gff > Spd_Ath.gff1.2 蛋白质序列比对（blastp）
# 构建物种蛋白序列数据库（单物种分析仅需构建Spd数据库） makeblastdb -in Spd.pep -dbtype prot -title Spd -parse_seqids -out Spd blastp -query Spd.pep -db Spd -evalue 1e-10 -max_target_seqs 5 -outfmt 6 -out Spd.blast # 跨物种分析：构建Ath数据库并比对（单物种分析可跳过） makeblastdb -in Ath.pep -dbtype prot -title Ath -parse_seqids -out Ath blastp -query Ath.pep -db Ath -evalue 1e-10 -max_target_seqs 5 -outfmt 6 -out Ath.blast mkdir Spd_Ath && cat Spd.blast Ath.blast > Spd_Ath.blast1.3 复制类型分类（DupGen-finder，核心步骤）
# 模式1：物种内自身比较（推荐，聚焦WGD分析） DupGen_finder.pl -i $PWD -t Spd -c Spd -o ${PWD}/Spd_self/results1 # 自身对比 DupGen_finder-unique.pl -i $PWD -t Spd -c Spd -o ${PWD}/Spd_self/results2 # 严格模式（去重） # 模式2：跨物种比较（可选，Spd vs Ath） DupGen_finder.pl -i $PWD -t Spd -c Ath -o ${PWD}/Spd_Ath/results1 DupGen_finder-unique.pl -i $PWD -t Spd -c Ath -o ${PWD}/Spd_Ath/results21.4 核心输出文件（聚焦 WGD）
| 文件名 | 内容说明 |
|---|---|
| Spd.wgd.genes-unique | 严格模式下 WGD 类型的所有基因列表（KEGG 富集分析输入） |
| Spd.wgd.pairs-unique | 严格模式下 WGD 类型的基因对详细信息（Ka/Ks 分析输入） |
| Spd.collinearity | 物种内共线性文件（辅助验证 WGD 事件） |
模块 2：WGD 基因筛选（与扩张基因取交集）
扩张变多的基因大概率发生复制事件！
cd results2 # 进入严格模式结果目录 # 筛选扩张基因中属于WGD类型的基因（核心输入文件：expansion_wgd.csv） grep -f expansion.genes Spd.wgd.genes-unique | cut -f1 > expansion_wgd.csv # 其他复制类型（TD/PD/TRD/DSD）同步筛选（用于对比） grep -f expansion.genes Spd.td.genes-unique | cut -f1 > expansion_td.csv grep -f expansion.genes Spd.pd.genes-unique | cut -f1 > expansion_pd.csv grep -f expansion.genes Spd.trd.genes-unique | cut -f1 > expansion_trd.csv grep -f expansion.genes Spd.dsd.genes-unique | cut -f1 > expansion_dsd.csv模块 3：WGD 基因 KEGG 功能富集分析
3.1 核心函数
setwd("e:/Family_expansion/GO/results2") # 替换为实际工作路径 library(clusterProfiler) library(ggplot2) library(cowplot) # 定义GO/KEGG富集分析函数（输入：基因列表、复制类型前缀） get_go_kegg <- function(filename,genetype,GO_list="GO.info",KEGG_list="KEGG.info"){ # 读取基因列表 gene <- read.csv(filename,header = FALSE) # 读取GO/KEGG注释库 golist <- read.table(file=GO_list,sep = "\t",header = FALSE) kegglist <- read.table(file=KEGG_list,sep = "\t",header = FALSE) # GO富集分析（输出PDF+CSV） spo2gene <- data.frame(golist$V1,golist$V3) spo2name <- data.frame(golist$V1,golist$V2) spo_GO <- enricher(t(gene$V1),TERM2GENE = spo2gene,TERM2NAME = spo2name) p1 <- dotplot(spo_GO) ggsave(file=paste(genetype,".GO.pdf",sep=""), dpi=300) write.csv(spo_GO,file=paste(genetype,"_GO.csv",sep = "")) # KEGG富集分析（输出PDF+CSV） kegg2gene <- data.frame(kegglist$V1,kegglist$V3) kegg2name <- data.frame(kegglist$V1,kegglist$V2) spo_KEGG <- enricher(t(gene$V1),TERM2GENE = kegg2gene,TERM2NAME = kegg2name) p2 <- dotplot(spo_KEGG) ggsave(file=paste(genetype,".KEGG.pdf",sep=""), dpi=300) write.csv(spo_KEGG,file=paste(genetype,"_kegg.csv",sep="")) # 组合图输出 plot_grid(p1, p2, nrow=2, labels=c("A", "B")) ggsave(file=paste(genetype,"_go_kegg.pdf",sep=""), dpi=300) } # 运行富集分析（聚焦WGD，其他类型用于对比） get_go_kegg("expansion_wgd.csv","wgd") # WGD核心分析 get_go_kegg("expansion_td.csv","td") get_go_kegg("expansion_pd.csv","pd") get_go_kegg("expansion_trd.csv","trd") get_go_kegg("expansion_dsd.csv","dsd")3.2 KEGG 富集结果可视化（多类型对比气泡图）
# 定义数据读取函数（提取前15个显著富集的KEGG term） getdata <- function(Prefix,enrich_type="KEGG",count_num=15){ kegg_res <- read.csv(paste(Prefix,"_",enrich_type,".csv",sep = ""),header = T) pic_kegg <- head(kegg_res,count_num) pic_kegg$adjust <- -log10(pic_kegg$p.adjust) # 计算校正后P值的-log10 pic_kegg$type <- Prefix return(pic_kegg) } # 读取各复制类型的KEGG富集结果（聚焦WGD） wgd_KEGG <- getdata("wgd",enrich_type = "KEGG") td_KEGG <- getdata("td",enrich_type = "KEGG") pd_KEGG <- getdata("pd",enrich_type = "KEGG") trd_KEGG <- getdata("trd",enrich_type = "KEGG") dsd_KEGG <- getdata("dsd",enrich_type = "KEGG") # 合并数据并绘制气泡图 data_kegg_all <- rbind(wgd_KEGG,td_KEGG,pd_KEGG,trd_KEGG,dsd_KEGG) ggplot(data_kegg_all,aes(type,Description)) + geom_point(aes(fill=adjust,size=Count),alpha=0.9,pch=21) + scale_fill_gradient(low='SpringGreen',high='DeepPink') + scale_size_area() + labs(y='KEGG Pathway',x='Duplication Type',fill='-log10(P.adjust)',size='Gene Count')+ theme_bw() ggsave("geneduplication.KEGG.pdf",dpi=300)模块 4：WGD 基因 Ka/Ks 与 4DTv 进化分析（Shell+R）
4.1 Ka/Ks 计算（Shell 脚本：KaKs.sh）
#!/bin/bash # 运行方式：bash KaKs.sh Spd（Spd为物种缩写） species=$1 proc=16 # 线程数（适配大规模基因分析） # 1. 提取各复制类型的基因对（聚焦WGD） cat $species.wgd.pairs-unique | cut -f 1,3 | sed '1d' > wgd.homolog cat $species.tandem.pairs-unique | cut -f 1,3 | sed '1d' > td.homolog cat $species.proximal.pairs-unique | cut -f 1,3 | sed '1d' > pd.homolog cat $species.transposed.pairs-unique | cut -f 1,3 | sed '1d' > trd.homolog cat $species.dispersed.pairs-unique | cut -f 1,3 | sed '1d' > dsd.homolog # 2. 定义Ka/Ks计算函数（ParaAT.pl） getKaKs(){ ParaAT.pl -h $1.homolog -n $2.cds -a $2.pep -p <(echo $proc) -m mafft -f axt -g -k -o $1.result_dir # 提取Ka/Ks结果并格式化 cat $1.result_dir/*.axt.kaks | cut -f 1,3,4,5 | grep -v 'Sequence' | sort | uniq >$1.KaKs.txt cat $1.KaKs.txt | tr "\t" "," | sed '1i Sequence,Ka,Ks,Ka/Ks' >$1.KaKs.csv # 添加复制类型列（用于后续可视化） cat $1.KaKs.csv | awk -F , -v type=$1 '{print $0","type}' | sed '1d' >$1.KaKs.Item.csv } # 3. 批量计算各类型Ka/Ks（聚焦WGD） for ele in wgd td pd trd dsd; do getKaKs $ele $species done # 4. 合并所有Ka/Ks结果 cat *.KaKs.Item.csv | sed '1i Sequence,Ka,Ks,Ka/Ks,Item' > gene_dup_KaKs.csv # 5. 计算4DTv值（聚焦WGD） get4DTv(){ cd $1.result_dir # 转换axt文件为单行格式 for i in `ls *.axt`; do axt2one-line.py $i ${i}.one-line; done # 计算4DTv ls *.axt.one-line | while read id; do calculate_4DTV_correction.pl $id >${id%%one-line}4dtv; done # 合并结果 for i in `ls *.4dtv`; do awk 'NR>1{print $1"\t"$3}' $i >>$1.4DTv.txt; done sort $1.4DTv.txt | uniq > ../$1.4DTv.txt cat ../$1.4DTv.txt | tr "\t" "," | sed '1i Gene,4DTv' > ../$1.4DTv.csv # 添加复制类型列 cat ../$1.4DTv.csv | awk -F , -v type=$1 '{print $0","type}' | sed '1d' > ../$1.4DTv.Item.csv cd ../ } # 6. 批量计算4DTv并合并结果 for ele in wgd td pd trd dsd; do get4DTv $ele done cat *.4DTv.Item.csv | sed '1i Gene,4DTv,Item' > gene_dup_4DTv.csv4.2 Ka/Ks 可视化（R 脚本：kaks.R）
library(ggplot2) library(cowplot) # 1. Ka/Ks分布可视化（聚焦WGD的选择压力） data_kaks <- read.csv("gene_dup_KaKs.csv", na.strings = c("NA", "")) data_kaks <- na.omit(data_kaks) # 删除缺失值 # 小提琴图（Ka/Ks值分布，核心展示WGD的选择压力） p0 <- ggplot(data=data_kaks, aes(x=Item, y=Ka.Ks, fill=Item)) + geom_violin(alpha=0.8, width=1) + geom_boxplot(alpha=0.5, width=0.3) + coord_flip() + # 翻转坐标轴，便于查看 ylim(0, 2) + # 聚焦0-2区间（Ka/Ks>1为正选择，<1为纯化选择） labs(x='Duplication Type', y='Ka/Ks Ratio') + guides(fill=FALSE) + theme_bw() ggsave('distribute_of_KaKs.pdf', dpi=300) # 2. Ka/Ks密度图（WGD与其他类型对比） p1 <- ggplot(data=data_kaks, aes(x=Ka, group=Item, color=Item)) + geom_density(alpha=0.4) + labs(x='Ka Value', y='Density', title='Distribution of Ka') + theme_bw() p2 <- ggplot(data=data_kaks, aes(x=Ks, group=Item, color=Item)) + geom_density(alpha=0.4) + labs(x='Ks Value', y='Density', title='Distribution of Ks') + theme_bw() # 3. 4DTv密度图（WGD进化分化特征） data_4DTv <- read.csv("gene_dup_4DTv.csv", na.strings = c("NA", "")) data_4DTv <- na.omit(data_4DTv) p3 <- ggplot(data=data_4DTv, aes(x=X4DTv, group=Item, color=Item)) + geom_density(alpha=0.4) + labs(x='4DTv Value', y='Density', title='Distribution of 4DTv') + theme_bw() # 4. 组合图输出（一站式展示所有结果） p4.1 <- plot_grid(p1, p2, ncol=2, labels=c("a", "b")) p4.2 <- plot_grid(p0, p3, ncol=2, labels=c("c", "d")) p4 <- plot_grid(p4.1, p4.2, nrow=2) ggsave("Distribution_of_Ka_Ks_Kaks_4DTv.pdf", dpi=300)06 结果
WGD的分类，扩张的复制基因富集，扩张基因的kaks选择压力，全基因组事件发生时间分析
KEGG 富集结果：WGD 基因显著富集的通路反映其核心生物学功能（如次生代谢、胁迫响应、生长发育），可对比 TD/PD 等类型，揭示 WGD 在物种进化中的功能贡献；
Ka/Ks 分析：若 WGD 基因的 Ka/Ks 均值 <1，表明其受纯化选择（功能保守）；若 Ka/Ks>1，提示正选择（功能分化）；
4DTv 分析：WGD 基因的 4DTv 值分布峰值可反映全基因组复制事件的发生时间，峰值越集中，说明 WGD 事件越同步。
