一：数据准备
自己去下载拟南芥的基因组，gff3，CDS序列，pep文件。
1.根据gff文件获得基因位置信息。
可以自己用awk抓，也可以写个perl,perl可以在我百度云里面找到。
$perl get_gene_position_form_gff.pl Arabidopsis_thaliana.TAIR10.41.gff3 AT.gff
查看一下结果，如下图，有的在第二列前面还有别的多余的东西，自己用sed删除一下就好了。
2.提取基因的蛋白质序列
由于一个基因有多个转录本，这里统一用第一个转录本来做blast。可以写个perl，这里讲下不用perl怎么办到，首先将AT.ff中的第二列gene_id取出来，再行尾加上.1保存为pep_id文件，再利用seqkit软件抓取匹配pepe_id的序列就是我们需要的蛋白序列。
$awk '{print $2}' AT.gff |sed 's/$/&.1/g' >pep_id
$seqkit grep -f pep_id Arabidopsis_thaliana.TAIR10.pep.all.fa >pep.fa
由于要与AT.gff文件的ID保持对应，需要将文件中的.1删掉。
sed -i 's/\.//g' pep.fa
pep.fa
3.blast比对
$makeblastdb -in pep.fa -dbtype prot -title pep.fa
$blastp -num_threads 16 -db pep.fa -query pep.fa -evalue 1e-20 -outfmt 6 -seg yes > mcsan/AT.blast
二.下载安装mcscanx
wget http://chibba.pgml.uga.edu/mcscan2/MCScanX.zip
unzip MCScanX.zip
cd MCScanX
make
1.报错
g++ struct.cc mcscan.cc read_data.cc out_utils.cc dagchainer.cc msa.cc permutation.cc -o MCScanX
msa.cc: In function ‘void msa_main(const char*)’:
msa.cc:289:22: error: ‘chdir’ was not declared in this scope
if (chdir(html_fn)<0)
^
makefile:2: recipe for target 'mcscanx' failed
make: *** [mcscanx] Error 1
这是由于MCScan不支持64位，只需要在msa.h, dissect_multiple_alignment.h, and detect_collinear_tandem_arrays.h 这三个文件 前面添加
#include <unistd.h>
再次运行make即可。
2.运行MCScanX
$cp AT.gff mcsan && cd mcsan
$/home/spider/project/yuantao/soft/MCScanX-master/MCScanX AT
结果文件
到此我们已经准备好了共线性分析了2个输入文件，我们还需要一个控制文件(.ctl)和一个基因家族的ID文件，这里我就用之前的NBS-ARC基因家族id。这里由于之前是一列，现在转成如下格式，这是mcscan官网的示例文件。当然也可以用官网的数据：
image.png
sed ':a;N;$!ba;s/\n/\t/g' nbs_gene_id.txt >gene_family_id
最后就是ctl文件，可以看官网：
需要注意的是ctl里面的第二行的染色体名字要和AT.gff保持一致
到此就准备好了绘图的四个文件：AT.gff，AT.collinearity，family.ctl，gene_family_id。就可以绘图了。
三.绘图
1.family_circle_plotter
记住，绘图是必须到MCScanX的downstream_analyses文件夹下。
$java -Djava.awt.headless=true family_circle_plotter -g /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/AT.gff -s /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/AT.collinearity -c /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/family.ctl -f /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/gene_family_id -o /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/NBS_circos.png
报错
原因是windows内核集成了gui，而linux上没有启动x server。在java运行参数上加-Djava.awt.headless=true即可解决。
$java -Djava.awt.headless=true circle_plotter -g /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/AT.gff -s /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/AT.collinearity -c /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/family.ctl -f /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/MADX_box.txt -o /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/dot.png
结果
2.circle_plotter
java circle_plotter -g /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/AT.gff -s /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/AT.collinearity -c /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/family.ctl -f /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/gene_family_id -o /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/dot.png
3.最后可以看一下NBS基因家族直接哪两个存在Segmental 复制。
$perl detect_collinearity_within_gene_families.pl -i /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/gene_family_id -d /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/AT.collinearity -o /home/spider/project/yuantao/test/Athaliana/McsacnX/mcsan/NBS_collinear.pairs
结果
以上就是NBS基因家族两两复制。
参考
https://www.sogou.com/link?url=hedJjaC291NxhTf9bAf9WgOl0lSYJUQpc4MGZlkSWzC2HE9UwMwWb6Lf2avUAKyr
