WGD分类进阶--随笔021

基于全基因组复制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 基因复制类型定义复制类型英文全称定义WGDWhole-genome duplication全基因组复制核心分析对象TDTandem duplication串联重复相邻的两个重复基因PDProximal duplication近端重复相隔 10 个以内基因的重复基因TRDTransposed duplication转置重复祖先和新基因座组成的重复基因DSDDispersed duplication分散重复不相邻也不共线性的重复基因SLSingleton单拷贝无同源重复基因04 输入文件要求文件类型文件名示例格式要求基因组注释Spd.gff/Ath.gff制表符分隔格式染色体ID\t基因ID\t起始位置\t终止位置示例Ath-Chr1 AT1G01010.1 3630 5899蛋白质序列Spd.pep/Ath.pepFASTA 格式序列 ID 需与 gff 文件中基因 ID 完全一致示例AT3G05780.1 氨基酸序列扩张基因列表expansion.genes单列基因 ID需与上述文件的基因 ID 命名规则统一功能注释库GO.info/KEGG.info制表符分隔GO.infoGOID\t功能描述\t基因ID\t功能分类KEGG.infoKOID\t通路名称\t基因ID05 分析流程分模块梳理聚焦 WGD模块 1基因复制类型鉴定DupGen-finder 核心步骤1.1 输入文件预处理统一格式# 物种1Spdgff文件格式化添加物种前缀调整列顺序 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 # 统一染色体命名格式 # 物种2Athgff文件格式化可选跨物种对比用 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 事件模块 2WGD 基因筛选与扩张基因取交集扩张变多的基因大概率发生复制事件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模块 3WGD 基因 KEGG 功能富集分析3.1 核心函数setwd(e:/Family_expansion/GO/results2) # 替换为实际工作路径 library(clusterProfiler) library(ggplot2) library(cowplot) # 定义GO/KEGG富集分析函数输入基因列表、复制类型前缀 get_go_kegg - function(filename,genetype,GO_listGO.info,KEGG_listKEGG.info){ # 读取基因列表 gene - read.csv(filename,header FALSE) # 读取GO/KEGG注释库 golist - read.table(fileGO_list,sep \t,header FALSE) kegglist - read.table(fileKEGG_list,sep \t,header FALSE) # GO富集分析输出PDFCSV 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(filepaste(genetype,.GO.pdf,sep), dpi300) write.csv(spo_GO,filepaste(genetype,_GO.csv,sep )) # KEGG富集分析输出PDFCSV 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(filepaste(genetype,.KEGG.pdf,sep), dpi300) write.csv(spo_KEGG,filepaste(genetype,_kegg.csv,sep)) # 组合图输出 plot_grid(p1, p2, nrow2, labelsc(A, B)) ggsave(filepaste(genetype,_go_kegg.pdf,sep), dpi300) } # 运行富集分析聚焦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_typeKEGG,count_num15){ 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(filladjust,sizeCount),alpha0.9,pch21) scale_fill_gradient(lowSpringGreen,highDeepPink) scale_size_area() labs(yKEGG Pathway,xDuplication Type,fill-log10(P.adjust),sizeGene Count) theme_bw() ggsave(geneduplication.KEGG.pdf,dpi300)模块 4WGD 基因 Ka/Ks 与 4DTv 进化分析ShellR4.1 Ka/Ks 计算Shell 脚本KaKs.sh#!/bin/bash # 运行方式bash KaKs.sh SpdSpd为物种缩写 species$1 proc16 # 线程数适配大规模基因分析 # 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 NR1{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.Rlibrary(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(datadata_kaks, aes(xItem, yKa.Ks, fillItem)) geom_violin(alpha0.8, width1) geom_boxplot(alpha0.5, width0.3) coord_flip() # 翻转坐标轴便于查看 ylim(0, 2) # 聚焦0-2区间Ka/Ks1为正选择1为纯化选择 labs(xDuplication Type, yKa/Ks Ratio) guides(fillFALSE) theme_bw() ggsave(distribute_of_KaKs.pdf, dpi300) # 2. Ka/Ks密度图WGD与其他类型对比 p1 - ggplot(datadata_kaks, aes(xKa, groupItem, colorItem)) geom_density(alpha0.4) labs(xKa Value, yDensity, titleDistribution of Ka) theme_bw() p2 - ggplot(datadata_kaks, aes(xKs, groupItem, colorItem)) geom_density(alpha0.4) labs(xKs Value, yDensity, titleDistribution 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(datadata_4DTv, aes(xX4DTv, groupItem, colorItem)) geom_density(alpha0.4) labs(x4DTv Value, yDensity, titleDistribution of 4DTv) theme_bw() # 4. 组合图输出一站式展示所有结果 p4.1 - plot_grid(p1, p2, ncol2, labelsc(a, b)) p4.2 - plot_grid(p0, p3, ncol2, labelsc(c, d)) p4 - plot_grid(p4.1, p4.2, nrow2) ggsave(Distribution_of_Ka_Ks_Kaks_4DTv.pdf, dpi300)06 结果WGD的分类扩张的复制基因富集扩张基因的kaks选择压力全基因组事件发生时间分析KEGG 富集结果WGD 基因显著富集的通路反映其核心生物学功能如次生代谢、胁迫响应、生长发育可对比 TD/PD 等类型揭示 WGD 在物种进化中的功能贡献Ka/Ks 分析若 WGD 基因的 Ka/Ks 均值 1表明其受纯化选择功能保守若 Ka/Ks1提示正选择功能分化4DTv 分析WGD 基因的 4DTv 值分布峰值可反映全基因组复制事件的发生时间峰值越集中说明 WGD 事件越同步。