参考：https://www.jianshu.com/p/0a14d971bd0a
require software
blastp，DupGen-finder
geneDuplication analysis
基因的复制类型分类分析，使用的是DupGen-finder，可以把所有基因按照复制类型分为5类，
WGD：全基因组复制
TD：串联重复（相邻的两个重复基因）
PD：近端重复（相隔10个以内基因的重复基因）
TRD：转置重复（祖先和新基因座组成的重复基因）
DSD：分散重复（不相邻也不共线性的重复基因）
SL：单拷贝
require input file
研究物种的gff,pep
#分析模式1（自身比较）和模式2（和其他物种比较）
####分析模式1
cat Spd.bed |sed 's/^/Spd-/g'|awk '{print $1"\t"$4"\t"$2"\t"$3}' >Spd.gff
cat Ath.bed |sed 's/^/Ath-Chr/g'|awk '{print $1"\t"$4"\t"$2"\t"$3}' >Ath.gff
sed -i 's/Chr0/Chr/g' Spd.gff
cat Spd.gff Ath.gff >Spd_Ath.gff
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
# Create a reference database
makeblastdb -in Ath.pep -dbtype prot -title Ath -parse_seqids -out Ath
# Align protein query sequences against the reference database
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.blast
#-t 是实验组 -c是外源对照组
#一般模式
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/results2
输入文件格式
要求的Ath.gff文件格式:
Ath-Chr1 AT1G01010.1 3630 5899
Ath-Chr1 AT1G01020.1 6787 9130
Ath-Chr1 AT1G01030.1 11648 13714
Ath.pep文件的格式：
>AT3G05780.1
MMPKRFNTSGFDTTLRLPSYYGFLHLTQSLTLNSRVFYGARHVTPPAIRIGSNPVQSLLL
FRAPTQLTGWNRSSRDLLGRRVSFSDRSDGVDLLSSSPILSTNPNLDDSLTVIALPLPHK
PLIPGFYMPIHVKDPKVLAALQESTRQQSPYVGAFLLKDCASTDSSSRSETEDNVVEKFK
VKGKPKKKRRKELLNRIHQVGTLAQISSIQGEQVILVGRRRLIIEEMVSEDPLTVRVDHL
>AT4G32100.1
MATNACKFLCLVLLFAFVTQGYGDDSYSLESLSVIQSKTGNMVENKPEWEVKVLNSSPCY
FTHTTLSCVRFKSVTPIDSKVLSKSGDTCLLGNGDSIHDISFKYVWDTSFDLKVVDGYIA
CS
输出文件
在results1和results2中输出结果文件基本一致，results2中的是严格模式，文件后面后缀是-unique
Spd.pairs.stats-unique #是各组复制信息的pairs数量。
Spd.wgd.genes-unique #WGD类型的所有基因
Spd.wgd.pairs-unique #WGD类型的pairs的详细信息
Spd.singletons #singletons类型的基因(即没有同源基因的基因)
Spd.tandem.pairs-unique # TD类型
Spd.proximal.pairs-unique #PD类型pairs
Spd.transposed.pairs-unique #TRD类型pairs
Spd.dispersed.pairs-unique #DSD类型pairs
Spd.collinearity #Spd物种内的共线性文件
Spd_Ath.collinearity #Spd和Ath物种的共线性文件
注意：一定要分清楚gene和gene pairs的数量并不是一一对应的关系，会有一对多。
后续分析
- GO,KEGG (使用的是5个genes-unique文件)
- Ka/Ks （使用的是5个pairs-unique文件和2个collinearity文件）
获取目标基因后，可以和扩张的基因取并集，然后富集分析这部分基因的功能。
expansion.genes是扩张基因列表
cd results2
grep -f expansion.genes Spd.wgd.genes-unique|cut -f1 >expansion_wgd.csv
grep -f expansion.genes Spd.td.genes-unique|cut -f1 >expansion_td.csv
grep -f expansion.genes Spd.trd.genes-unique|cut -f1 >expansion_trd.csv
grep -f expansion.genes Spd.pd.genes-unique|cut -f1 >expansion_pd.csv
grep -f expansion.genes Spd.dsd.genes-unique|cut -f1 >expansion_dsd.csv
expansion_wgd.genes.csv即为扩张的wgd的基因。
注意：gff文件，protein文件和expansion文件，各处使用的geneid一定要一致
可以按照这5种分类
GO、KEGG富集分析
下面代码可以直接在R-studio运行，
要求先安装ggplot2和clusterProfiler
修改工作路径，工作路径下要有以下文件：
GO.info #GO的数据库
KEGG.info #KEGG的数据库
expansion_wgd.csv
expansion_td.csv
expansion_dsd.csv
expansion_pd.csv
expansion_trd.csv
5个csv文件时扩张基因和对应复制方式的共有的基因的列表
#使用新的GO和KEGG进行富集分析
setwd("e:/Family_expansion/GO/results2")
library(clusterProfiler)
library(ggplot2)
##获取spo基因的富集分析结果，需要输入基因列表和genetype(用于输出文件名的前缀)
get_go_kegg <- function(filename,genetype,GO_list="GO.info",KEGG_list="KEGG.info"){
#定义需要分析的基因列表
gene <- read.csv(filename,header = FALSE)
genetype <- genetype
#修改GO和KEGG文件名
golist <- read.table(file=GO_list,sep = "\t",header = FALSE)
kegglist <- read.table(file=KEGG_list,sep = "\t",header = FALSE)
#下面内容基本不需要修改
spo2gene <- data.frame(golist$V1,golist$V3) #GO,Geneid
spo2name <- data.frame(golist$V1,golist$V2) #GO,GO描述
gene1 <- t(gene$V1) #获取需要分析的基因列表
#GO富集分析
spo_GO <- enricher(gene1,TERM2GENE = spo2gene,TERM2NAME = spo2name)
p1 <- dotplot(spo_GO)
p1
ggsave(file=paste(genetype,".GO.pdf",sep=""))
write.csv(spo_GO,file=paste(genetype,"_GO.csv",sep = ""))
#kegg富集
kegg2gene <- data.frame(kegglist$V1,kegglist$V3) #keggid,Geneid
kegg2name <- data.frame(kegglist$V1,kegglist$V2) #keggid,kegg描述
gene1 <- t(gene$V1) #获取需要分析的基因列表
spo_KEGG <- enricher(gene1,TERM2GENE = kegg2gene,TERM2NAME = kegg2name)
p2 <- dotplot(spo_KEGG)
p2
ggsave(file=paste(genetype,".KEGG.pdf",sep=""))
write.csv(spo_KEGG,file=paste(genetype,"_kegg.csv",sep=""))
cowplot::plot_grid(p1, p2, nrow=2, labels=c("A", "B")) #ggplot2的函数，nrow指定2行，如果是ncol则是2列。
ggsave(file=paste(genetype,"_go_kegg.pdf",sep=""))
}
#运行富集分析，第1个参数是基因列表csv格式，第2个是前缀类型，第3个是GO,第4个是KEGG，参数3和4要求是tab分割符
get_go_kegg("expansion_trd.csv","trd")
get_go_kegg("expansion_td.csv","td")
get_go_kegg("expansion_dsd.csv","dsd")
get_go_kegg("expansion_wgd.csv","wgd")
get_go_kegg("expansion_pd.csv","pd")
#GO和KEGG文件的格式
head(GO.info)
# GO:0005452 inorganic anion exchanger activity Sp09G018780.1 Molecular Function
# GO:0006820 anion transport Sp09G018780.1 Biological Process
# GO:0016021 integral component of membrane Sp09G018780.1 Cellular Component
head(KEGG.info)
# ko00010 Glycolysis / Gluconeogenesis Sp12G016220.1
# ko00010 Glycolysis / Gluconeogenesis Sp10G003910.1
# ko00010 Glycolysis / Gluconeogenesis Sp16G007040.1
#合并多组GO分析，生成气泡图
#在进行气泡图分析之前，可以手动修改上一步输出的富集分析的csv文件中的GO/KEGG term的词条，把长词条变成短词条，以方便在气泡图上展示。
#从csv文件读取富集的信息到数据框（count_num可以指定获取的GO/KEGG term的数量，默认是15;enrich_type指定富集类型GO/KEGG）
getdata <- function(Prefix,enrich_type="GO",count_num=15){
trd_GO <- read.csv(paste(Prefix,"_",enrich_type,".csv",sep = ""),header = T)
pic_trd_GO <- head(trd_GO,count_num)
pic_trd_GO$adjust <- -log10(pic_trd_GO$p.adjust)
pic_trd_GO$type <- Prefix
return(pic_trd_GO)
}
#读取GO
trd_GO <- getdata("trd")
wgd_GO <- getdata("wgd")
dsd_GO <- getdata("dsd")
pd_GO <- getdata("pd")
td_GO <- getdata("td")
#合并5组GO数据
data_GO_all <- rbind(trd_GO,wgd_GO,dsd_GO,pd_GO,td_GO)
#气泡图可视化富集结果
ggplot(data_GO_all,aes(type,Description)) +
geom_point(aes(fill=adjust,size=Count),alpha=0.9,pch=21) + #fill对应点的填充色，colour对应点的边框色
scale_fill_gradient(low='SpringGreen',high='DeepPink') + #设定颜色的变化范围
scale_size_area() + #设定点的大小比例和图例上显示的间隔
labs(y='GO term',x='type',fill='-log10(P.adjust)',size='Metabolites number')+
theme_bw()
ggsave("geneduplication.GO.pdf",dpi=300)
write.csv(data_GO_all,file="data_GO_all.csv")
#读取KEGG
trd_KEGG <- getdata("trd",enrich_type = "KEGG")
wgd_KEGG <- getdata("wgd",enrich_type = "KEGG")
dsd_KEGG <- getdata("dsd",enrich_type = "KEGG")
pd_KEGG <- getdata("pd",enrich_type = "KEGG")
td_KEGG <- getdata("td",enrich_type = "KEGG")
#合并5组GO数据
data_kegg_all <- rbind(trd_KEGG,wgd_KEGG,dsd_KEGG,pd_KEGG,td_KEGG)
#气泡图可视化富集结果
ggplot(data_kegg_all,aes(type,Description)) +
geom_point(aes(fill=adjust,size=Count),alpha=0.9,pch=21) + #fill对应点的填充色，colour对应点的边框色
scale_fill_gradient(low='SpringGreen',high='DeepPink') + #设定颜色的变化范围
scale_size_area() + #设定点的大小比例和图例上显示的间隔
labs(y='KEGG term',x='type',fill='-log10(P.adjust)',size='Metabolites number')+
theme_bw()
ggsave("geneduplication.KEGG.pdf",dpi=300)
write.csv(data_kegg_all,file="data_KEGG_all.csv")
分别分析Ka/Ks 参考
需要安装的软件 KaKs,mafft,ParaAT.pl (推荐使用muscle比对，速度更快)
当基因比较多超过1万时，比较耗时，请增加cpu数量，对应的参数时proc脚本中设置的是16.
KaKs.sh脚本如下,工作目录results2
#定义物种的缩写
species=$1
cat $species.dispersed.pairs-unique|cut -f 1,3 |sed '1d' >dsd.homolog
cat $species.tandem.pairs-unique | cut -f 1,3 |sed '1d' >td.homolog
cat $species.transposed.pairs-unique | cut -f 1,3 |sed '1d' >trd.homolog
cat $species.wgd.pairs-unique | cut -f 1,3 |sed '1d' >wgd.homolog
cat $species.proximal.pairs-unique | cut -f 1,3 |sed '1d' >pd.homolog
echo "16" >proc
function getKaKs(){
ParaAT.pl -h $1.homolog -n $2.cds -a $2.pep -p proc -m mafft -f axt -g -k -o $1.result_dir
#把KaKs输出到文件KaKs.txt
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
#添加新的列，Item用于后续合并后可视化分组(注意：awk里的$1是全局的环境的$1，awk调用外部变量，使用多重引号)
cat $1.KaKs.csv|awk -F , '{print $0",""'"$1"'"}' | sed '1d' >$1.KaKs.Item.csv
}
list=(wgd td pd trd dsd)
for ele in ${list[@]}
do
getKaKs $ele $species
done
#合并所有的KaKs输出到一个csv文件
cat td.KaKs.Item.csv wgd.KaKs.Item.csv pd.KaKs.Item.csv trd.KaKs.Item.csv dsd.KaKs.Item.csv |sed '1i Sequence,Ka,Ks,Ka.Ks,Item' >gene_dup_KaKs.csv
##计算4DTv
# 将多行axt文件转换成单行
function get4DTv(){
cd $1.result_dir
#使用axt2one-line.py合并axt为1行。https://github.com/scbgfengchao/4DTv/blob/master/axt2one-line.py
for i in `ls *.axt`;do axt2one-line.py $i ${i}.one-line;done #使用calculate_4DTV_correction.pl脚本计算4dtv值。脚本地址：https://github.com/JinfengChen/Scripts/blob/master/FFgenome/03.evolution/distance_kaks_4dtv/bin/calculate_4DTV_correction.pl
ls *.axt.one-line|while read id;do calculate_4DTV_correction.pl $id >${id%%one-line}4dtv;done
#合并所有同源基因对的4dtv
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
cd ../
#添加新的列，Item用于后续合并后可视化分组
cat $1.4DTv.csv|awk -F , '{print $0",""'"$1"'"}' | sed '1d' >$1.4DTv.Item.csv
}
list=(wgd td pd trd dsd)
for ele in ${list[@]}
do
get4DTv $ele;
done
cat td.4DTv.Item.csv wgd.4DTv.Item.csv pd.4DTv.Item.csv trd.4DTv.Item.csv dsd.4DTv.Item.csv |sed '1i Gene,4DTv,Item' >gene_dup_4DTv.csv
运行bash KaKs.sh Spd即可得到最终的总KaKs和4DTv, gene_dup_KaKs.csv和gene_dup_4DTv.csv,还有各组分组的.KaKs.csv和.4DTv.csv。
可视化Ka/Ks的分布,
kaks.R脚本代码如下，直接在上述结果目录运行Rscript kaks.R 即可得到ka/ks的分布，输出4个pdf,3个是单独的，1个是组图。
#setwd("e:/****/GeneDuplication")
library("ggplot2")
data_kaks <- read.csv("gene_dup_KaKs.csv")
#删除包含NA的行
data_kaks <- na.omit(data_kaks)
p0 <- ggplot(data=data_kaks, aes(x=Item, y=Ka.Ks, fill=Item)) +
geom_violin(alpha=0.8,width=1)+ guides(fill=F)+xlab(' ')+ylab('Ka/Ks')+
geom_boxplot(alpha=0.5,width = 0.3)+
coord_flip()+ #颠倒xy轴
ylim(0,2) #设置y轴最大值
p0
ggsave('distribute_of_KaKs.pdf',dpi=300)
p1 <- ggplot(data=data_kaks, aes(x=Ka,group=Item)) +
geom_density(alpha=0.4,aes(color = Item))+ xlab('Ka value')+ylab('Density')+
labs(title = "Distribution of Ka distances", size = 1.5)+guides()
p1
ggsave('Distribute_of_Ka.pdf',dpi=300)
p2 <- ggplot(data=data_kaks, aes(x=Ks,group=Item)) +
geom_density(alpha=0.4,aes(color = Item))+xlab('Ks value')+ylab('Density')+
labs(title = "Distribution of Ks distances", size = 1.5)
p2
ggsave('Distribute_of_Ks.pdf',dpi=300)
#可视化4DTv
data_4DTv <- read.csv("gene_dup_4DTv.csv")
#删除包含NA的行
data_4DTv <- na.omit(data_4DTv)
p3 <- ggplot(data=data_4DTv, aes(x=X4DTv,group=Item)) +
geom_density(alpha=0.4,aes(color = Item))+xlab('4DTv value')+ylab('Density')
#labs(title = "Distribution of 4DTv distances", size = 1.5)
p3
ggsave("Distribution_of_4DTv.pdf",dpi=300)
#拼图，直接出组合图
p4.1 <- cowplot::plot_grid(p1, p2, ncol=2, labels=c("a", "b")) #ggplot2的函数，nrow指定2行，如果是ncol则是2列。
p4.2<-cowplot::plot_grid(p0,p3,ncol = 2,labels = c("c","d"))
p4<-cowplot::plot_grid(p4.1,p4.2,nrow = 2)
p4
ggsave("Distribution_of_Ka_Ks_Kaks_4DTv.pdf",dpi = 300)
2022/06/07更新
注意：
- 即使使用和其他物种比较，即模式2，获取的仍然是本身物种的复制信息。所以最终获取到的仍然是自身的分化时间。但是，外群不一样，分化时间会有变化。可以根据物种的进化树来选择合适的外群。距离研究物种比较近的物种作为外群，可能会鉴定到比较少的复制基因对。因为外群模式的结果是物种和外群分开后的复制信息。
例如：rice作为研究物种，A. thaliana作为外群。这时候计算出的结果就是rice和Ath分开后，rice经历的复制。分开前共同经历的复制没有被计算。
选择物种进化比较近的物种作为外群，这样可能更方便讲故事。 - 外群的选择，可以根据研究分析目的来选择不同的外群。同时还可以比较同一个物种和不同外群分开后的复制差异。https://github.com/qiao-xin/DupGen_finder/issues/5
- 作者提供了下游分析的pipline
- Calculating Ks values for WGD-derived gene pairs pipeline 计算KaKs的值的管道
- Inferring the Ks peaks corresponding to WGD events of different ages pipeline计算WGD峰值的管道
- compute Pearson's coorelation coefficient (r) between expression profiles of duplicate gene pairs derived from different modes of gene duplication 计算不同复制模式的基因对和表达量之间的皮尔逊相关性
- detect gene conversion events检测目标物种和外群之间的 WGD 衍生基因对中的基因转换事件
