参考：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)