全基因组复制事件
1 环境准备
conda create -n comp_genome python=32 WGDI软件安装
xxxxxxxxxxgit clone https://github.com/SunPengChuan/wgdi.gitcd wgdipython setup.py install3 WGDI共线性分析
3.1 准备输入文件*.len、*.gff
新建dotplot_input_preprocessing.sh，修改输入的基因组文件名和注释文件名，以及输出的Tfa.len、Tfa.gff等。
x -0076 - '@F=map{s/[>\r\n]//gr}@F;$id=shift @F;print $id,qq{\t},length (join q{},@F),qq{\n} if $id' . >.chr.length - "/Contig/d" .chr.length - 'next if /^#/;$count{$F[0]}++ if $F[2] eq "mRNA";END{print join qq{\t},$_,$count{$_} for sort keys %count}' . > .. - 'if($flag){print join qq{\t},$_,$count{$F[0]}}else {$count{$F[0]}=$F[1]}$flag=1 if eof(ARGV)' .. .chr.length|sort -,1 >. - "/Contig/d" . - 'next unless $F[2] eq "mRNA";/ID=([^;]+)/;push @geneInfo,[$F[0],$1,$F[3],$F[4],$F[6]];END{$preChr="";for(sort {$a->[0] cmp $b->[0] || $a->[2] <=> $b->[2]} @geneInfo){if($preChr ne $_->[0]){$c=0;$preChr=$_->[0]};print join qq{\t},@{$_},++$c}}' . > .
Tfa.gff文件：为便于后期画图，可以将染色体序号简化为1,2,3等（第一列）
xxxxxxxxxx9 Ts_LG09_G01064.t1 18893935 18895118 + 10649 Ts_LG09_G01065.t1 18897494 18902194 - 10659 Ts_LG09_G01066.t1 18905024 18910662 - 10669 Ts_LG09_G01067.t1 18911694 18915629 + 10679 Ts_LG09_G01068.t1 18916950 18936981 - 106810 Ts_LG10_G00001.t1 45027 48258 - 110 Ts_LG10_G00002.t1 266222 270097 + 2
Tsi.gff文件：
xxxxxxxxxx1 Maker00019331 40209 48750 + 11 Maker00019333 446854 447609 - 21 Maker00019336 455466 456678 - 31 Maker00019338 797494 801130 - 41 Maker00019332 820015 824148 - 51 Maker00019340 910965 913873 - 61 Maker00019339 951822 976238 + 71 Maker00019334 1050357 1051216 + 8
Tfa.len文件：为便于后期画图，可以将染色体序号简化为1,2,3等（第一列）
xxxxxxxxxx1 29788854 23922 26751989 22113 23787828 16234 21062203 13185 20971970 15866 19582046 11797 19432475 14808 19374345 13219 18956111 106810 18944402 137711 18887728 131612 18608799 1359
Tsi.len文件：
xxxxxxxxxx1 12513000 5872 19105972 6663 16665695 8864 13746830 8675 19252003 9126 17559502 9527 20082163 11408 21610000 11359 18254849 111310 18285906 113011 18492580 97312 19558692 997
3.2 准备输入文件Tsi_blastp_Tfa.txt
准备基因组蛋白质文件：Toona_ciliata_LG_CDS.pep、Toona_sinensis.Protein.fasta（需要比对的两个基因组，自身比较或两者比较）
xxxxxxxxxxmakeblastdb -in Toona_ciliata_LG_CDS.pep -dbtype prot-in Toona_ciliata_LG_CDS.pep：输入文件名，即包含蛋白质序列的FASTA文件。
-dbtype prot：指定数据库类型为蛋白质。
xxxxxxxxxxblastp -num_threads 40 -db Toona_ciliata_LG_CDS.pep -query Toona_sinensis.Protein.fasta -outfmt 6 -evalue 1e-5 -num_alignments 20 -out Tsi_blastp_Tfa.txt &blastp：BLASTP程序，用于蛋白质序列比对。
-num_threads 40：使用40个线程并行处理，以加快比对速度。
-db Toona_ciliata_LG_CDS.pep：指定数据库文件Toona_ciliata_LG_CDS.pep。
-query Toona_sinensis.Protein.fasta：指定查询文件Toona_sinensis.Protein.fasta。
-outfmt 6：输出格式为表格（格式6），包含每个比对的详细信息。
-evalue 1e-5：设置期望值（E-value）阈值为1e-5，用于过滤不显著的比对结果。
-num_alignments 20：每个查询序列最多返回20个比对结果。
-out Tsi_blastp_Tfa.txt：将比对结果输出到Tsi_blastp_Tfa.txt文件中。
&：在后台运行命令。
Tfa.blastp.txt文件：
xxxxxxxxxxTs_LG01_G00894.t1 Ts_LG01_G00894.t1 100.000 226 0 0 1 226 1 226 1.17e-168 463Ts_LG01_G00894.t1 Ts_LG02_G01385.t1 93.805 226 14 0 1 226 123 348 5.48e-157 438Ts_LG01_G00378.t1 Ts_LG01_G00378.t1 100.000 352 0 0 1 352 1 352 0.0 722Ts_LG01_G00378.t1 Ts_LG03_G00318.t1 59.838 371 101 8 7 346 10 363 1.15e-151 431
3.3 基因组共线性点阵图Dotplot
3.3.1 新建配置文件
xxxxxxxxxxwgdi -d ? >>Tfa_Vs_Tsi.conf3.3.2 利用vi命令修改配置文件
xxxxxxxxxxvi Tfa_Vs_Tsi.confTfa_Vs_Tsi.conf文件如下：
xxxxxxxxxx[dotplot]blast = Tsi_blastp_Tfa.txt # BLAST比对结果文件gff1 = Tfa.gff # 第一个基因组的GFF注释文件gff2 = Tsi.gff # 第二个基因组的GFF注释文件lens1 = Tfa.len # 第一个基因组的长度文件lens2 = Tsi.len # 第二个基因组的长度文件genome1_name = T.fargesii # 第一个基因组的名称genome2_name = T.sinensis # 第二个基因组的名称multiple = 1 # 倍数因子score = 100 # 比对的最低得分阈值evalue = 1e-5 # BLAST比对的e值阈值repeat_number = 20 # 每个基因允许的最大重复数position = order # 基因位置排序blast_reverse = false # 是否反向BLASTancestor_left = none # 左边的祖先信息ancestor_top = none # 顶部的祖先信息markersize = 0.5 # 标记点的大小figsize = 14,14 # 图的大小savefig = T.fargesii_Vs_T.sinensis.dotplot.pdf # 保存图的文件名3.3.3 运行并绘制共线性点阵图
xxxxxxxxxxwgdi -d Tfa_Vs_Tsi.conf3.4 获取共线性（collinearity）区块的信息
3.4.1 新建配置文件
xxxxxxxxxxwgdi -icl ? >>Tfa_Vs_Tsi.collinearity.conf3.4.2 利用vi命令修改配置文件
xxxxxxxxxxvi Tfa_Vs_Tsi.collinearity.confTfa_Vs_Tsi.collinearity.conf文件如下：
xxxxxxxxxx[collinearity]gff1 = Tfa.gff # 第一个基因组的GFF注释文件gff2 = Tsi.gff # 第二个基因组的GFF注释文件lens1 = Tfa.len # 第一个基因组的长度文件lens2 = Tsi.len # 第二个基因组的长度文件blast = Tsi_blastp_Tfa.txt # BLAST比对结果文件blast_reverse = false # 是否反向BLASTmultiple = 1 # 倍数因子process = 40 # 处理的进程数量evalue = 1e-5 # BLAST比对的e值阈值score = 100 # 比对的最低得分阈值grading = 50,25,10 # 分级标准mg = 25,25 # 最小和最大间距pvalue = 0.2 # p值阈值repeat_number = 10 # 每个基因允许的最大重复数positon = order # 基因位置排序savefile = Tfa_Vs_Tsi.collinearity.txt # 保存结果的文件名3.4.3 运行并获取共线性区块的信息
xxxxxxxxxxwgdi -icl Tfa_Vs_Tsi.collinearity.conf4 Ka、Ks计算
4.1 两物种基因组CDS和protein文件合并
T.fa和T.si基因组CDS文件合并
xxxxxxxxxxcat Toona_ciliata_LG.fasta Toona_sinensis.CDS.fasta > Tfa_Tsi.CDS.merged.fasta T.fa和T.si基因组protein文件合并
xxxxxxxxxxcat Toona_ciliata_LG_CDS.pep Toona_sinensis.Protein.fasta > Tfa_Tsi.protein.merged.pep4.2 新建配置文件
xxxxxxxxxxwgdi -ks ? >>Tfa_Vs_Tsi.ks.conf4.3 利用vi命令修改配置文件
xxxxxxxxxxvi Tfa_Vs_Tsi.ks.confTfa_Vs_Tsi.ks.conf文件如下：
xxxxxxxxxx[ks]cds_file = Tfa_Tsi.CDS.merged.fasta#cat all cds files togetherpep_file = Tfa_Tsi.protein.merged.pep#cat all pep files togetheralign_software = musclepairs_file = Tfa_Vs_Tsi.collinearity.txtks_file = Tfa_Vs_Tsi.ks4.4 运行并获取Ka、Ks的信息
xxxxxxxxxxwgdi -ks Tfa_Vs_Tsi.ks.confKs结果文件Tfa_Vs_Tsi.ks：
xxxxxxxxxxid1 id2 ka_NG86 ks_NG86 ka_YN00 ks_YN00Ts_LG01_G00656.t1 Maker00016752 0.3139 -1.0 0.3205 3.8855Ts_LG01_G00657.t1 Maker00016757 0.228 1.4029 0.2346 1.7648Ts_LG01_G00659.t1 Maker00016716 0.2928 1.006 0.2878 1.7914
5 整合共线性区块和ks信息
5.1 新建配置文件
xxxxxxxxxxwgdi -bi ? >>Tfa_Vs_Tsi.blockinfo.conf5.2 利用vi命令修改配置文件
xxxxxxxxxxvi Tfa_Vs_Tsi.blockinfo.confTfa_Vs_Tsi.blockinfo.conf文件如下：
xxxxxxxxxx[blockinfo]blast = Tsi_blastp_Tfa.txt # BLAST比对结果文件gff1 = Tfa.gff # 第一个基因组的GFF注释文件gff2 = Tsi.gff # 第二个基因组的GFF注释文件lens1 = Tfa.len # 第一个基因组的长度文件lens2 = Tsi.len # 第二个基因组的长度文件collinearity = Tfa_Vs_Tsi.collinearity.txt # 同线性分析结果文件score = 100 # 比对的最低得分阈值evalue = 1e-5 # BLAST比对的e值阈值repeat_number = 20 # 每个基因允许的最大重复数position = order # 基因位置排序ks = Tfa_Vs_Tsi.ks # KS值文件ks_col = ks_NG86 # KS值列名称savefile = Tfa_Vs_Tsi.blockinfo.csv # 保存结果的文件名5.3 运行并获取共线性区块上所有基因对的ks值（均值，中位数）
xxxxxxxxxxwgdi -bi Tfa_Vs_Tsi.blockinfo.conf6 Ks分布拟合“单次”WGD事件
在Lynch和Conery在2000年发表在Science的论文中，他们证明了小规模基因复制的Ks分布是L型，而在L型分布背景上叠加的峰则是来自于演化历史中某个突然的大规模复制事件。 L型分布（呈指数分布, exponential distribution)的峰可能是近期的串联复制引起，随着时间推移基因丢失，形成一个向下的坡。正态分布(normal distribution)的峰则是由全基因组复制引起。
6.1 过滤并绘制ks点阵图
bk模块绘制的ks点阵图的点只包含确认了共线性的基因对，用ks值作为点的颜色信息，可以根据ks点阵图的共线性区域的颜色来区分不同时期的多倍化事件。
6.1.1 新建配置文件
xxxxxxxxxxwgdi -bk ? >>Tfa_Vs_Tsi.blockks.conf6.1.2 利用vi命令修改配置文件
xxxxxxxxxxvi Tfa_Vs_Tsi.blockks.confTfa_Vs_Tsi.blockks.conf文件如下：
xxxxxxxxxx[blockks]lens1 = Tfa.len # 第一个基因组的长度文件lens2 = Tsi.len # 第二个基因组的长度文件genome1_name = T.fargesii # 第一个基因组的名称genome2_name = T.sinensis # 第二个基因组的名称blockinfo = Tfa_Vs_Tsi.blockinfo.csv # 同线性区块信息文件pvalue = 0.05 # p值阈值tandem = true # 是否考虑串联重复tandem_length = 200 # 串联重复的长度阈值markersize = 1 # 标记点的大小area = 0,3 # KS值范围block_length = 5 # 区块最小长度figsize = 14,14 # 图的大小savefig = Tfa_Vs_Tsi.blockks.pdf # 保存图的文件名6.1.3 运行并绘制ks点阵图
xxxxxxxxxxwgdi -bk Tfa_Vs_Tsi.blockks.conf6.2 过滤并绘制ks频率分布图
通过计算共线性区块的基因对ks值，可以获得基因对复制发生的时间，如果有全基因组复制（WGD）发生，那么现有物种的基因组会留下许多ks相近的基因对，通过ks频率分布图可以看到峰，由此判断WGD的发生次数和发生时间。
6.2.1 新建配置文件
x
wgdi -kp ? >>Tfa_Vs_Tsi.kspeaks.conf6.2.2 利用vi命令修改配置文件
xxxxxxxxxxvi Tfa_Vs_Tsi.kspeaks.confTfa_Vs_Tsi.kspeaks.conf文件如下：
x
[kspeaks]blockinfo = Tfa_Vs_Tsi.blockinfo.csvpvalue = 0.05tandem = trueblock_length = 5ks_area = 0,10multiple = 1homo = -1,1fontsize = 9area = 0,3figsize = 10,6.18savefig = Tfa_Vs_Tsi.kspeaks.distri.pdfsavefile = Tfa_Vs_Tsi.ks_medain.distri.csv6.2.3 运行并绘制ks频率分布图
xxxxxxxxxxwgdi -kp Tfa_Vs_Tsi.kspeaks.conf6.3 分段（ks_area）过滤并绘制ks频率分布图
6.3.1 新建配置文件1
xxxxxxxxxxcp ./Tfa_Vs_Tsi.kspeaks.conf ./Tfa_Vs_Tsi.kspeaks1.conf6.3.2 利用vi命令修改配置文件1
xxxxxxxxxxvi Tfa_Vs_Tsi.kspeaks1.confTfa_Vs_Tsi.kspeaks1.conf文件如下：
x
[kspeaks]blockinfo = Tfa_Vs_Tsi.blockinfo.csvpvalue = 0.05tandem = trueblock_length = 5ks_area = 0,0.75multiple = 1homo = -1,1fontsize = 9area = 0,3figsize = 10,6.18savefig = Tfa_Vs_Tsi.kspeaks1.distri.pdfsavefile = Tfa_Vs_Tsi.ks_medain1.distri.csv6.3.3 运行并绘制ks频率分布图1
xxxxxxxxxxwgdi -kp Tfa_Vs_Tsi.kspeaks1.conf6.3.4 新建配置文件2
xxxxxxxxxxcp ./Tfa_Vs_Tsi.kspeaks.conf ./Tfa_Vs_Tsi.kspeaks2.conf6.3.5 利用vi命令修改配置文件2
xxxxxxxxxxvi Tfa_Vs_Tsi.kspeaks2.confTfa_Vs_Tsi.kspeaks2.conf文件如下：
x
[kspeaks]blockinfo = Tfa_Vs_Tsi.blockinfo.csvpvalue = 0.05tandem = trueblock_length = 5ks_area = 0.75,2.5multiple = 1homo = -1,1fontsize = 9area = 0,3figsize = 10,6.18savefig = Tfa_Vs_Tsi.kspeaks2.distri.pdfsavefile = Tfa_Vs_Tsi.ks_medain2.distri.csv6.3.6 运行并绘制ks频率分布图2
xxxxxxxxxxwgdi -kp Tfa_Vs_Tsi.kspeaks2.conf6.4 高斯拟合ks频率分布图的峰——pf模块
6.4.1 新建配置文件1
x
wgdi -pf ? >>Tfa_Vs_Tsi.peaksfit1.conf6.4.2 利用vi命令修改配置文件1
xxxxxxxxxxvi Tfa_Vs_Tsi.peaksfit1.confTfa_Vs_Tsi.peaksfit1.conf文件如下：
xxxxxxxxxx[peaksfit]blockinfo = Tfa_Vs_Tsi.ks_medain1.distri.csvmode = medianbins_number = 200ks_area = 0,0.75fontsize = 9area = 0,3figsize = 10,6.18shadow = truesavefig = Tfa_Vs_Tsi.peaksfit1.pdf
6.4.3 运行并绘制ks频率分布图的峰1
x
wgdi -pf Tfa_Vs_Tsi.peaksfit1.conf6.4.4 新建配置文件2
xxxxxxxxxxwgdi -pf ? >>Tfa_Vs_Tsi.peaksfit2.conf6.4.5 利用vi命令修改配置文件2
xxxxxxxxxxvi Tfa_Vs_Tsi.peaksfit2.confTfa_Vs_Tsi.peaksfit2.conf文件如下：
xxxxxxxxxx[peaksfit]blockinfo = Tfa_Vs_Tsi.ks_medain2.distri.csvmode = medianbins_number = 200ks_area = 1,2fontsize = 9area = 0,3figsize = 10,6.18shadow = truesavefig = Tfa_Vs_Tsi.peaksfit2.pdf
6.4.6 运行并绘制ks频率分布图的峰2
xxxxxxxxxxwgdi -pf Tfa_Vs_Tsi.peaksfit2.conf6.5 拟合结果作图——kf模块
6.5.1 新建拟合参数文件Tfa_Vs_Tsi.all_ks.csv：（,,,,,,代表最多有两个峰）
x
,color,linewidth,linestyle,,,,,,Tfa_Vs_Tsi,pink,1,-,1.5890660930749103,1.4696559232132298,0.25568723949817884,1.5890660930749103,1.4696559232132298,0.255687239498178846.5.2 新建配置文件
x
wgdi -kf ? >>Tfa_Vs_Tsi.ksfigure.conf6.5.3 利用vi命令修改配置文件
xxxxxxxxxxvi Tfa_Vs_Tsi.ksfigure.confTfa_Vs_Tsi.peaksfit2.conf文件如下：
xxxxxxxxxx[ksfigure]ksfit = Tfa_Vs_Tsi.all_ks.csvlabelfontsize = 15legendfontsize = 15xlabel = Synonymous nucleotide subsititution(Ks)ylabel = No.of syntenic blocks kernel densitytitle = nonearea = 0,3figsize = 10,6.18shadow = true (true/false)savefig = Tfa_Vs_Tsi.ksfigure.pdf
6.5.4 运行并绘制拟合结果
x
wgdi -kf Tfa_Vs_Tsi.ksfigure.conf