共线性
共线性片段是指一个同一个物种内部或者两个物种之间，由基因组复制、染色体复制、大片段复制以及物种分化而产生的大片段同源现象。同源片段内部基因排序保守，意味着功能上也可能保守。
我们做共线性分析的目的就是为了鉴定物种内或者物种间的共线性片段。
synteny：两个不同物种的一组基因位点，在每个物种中位于相同的染色体上（不一定是相同的顺序）。
collinearity：两个不同物种中的一组基因位点，它们分别位于各自的同一条染色体上，并且顺序也是一致的。
MCScanX 网址：https://github.com/wyp1125/MCScanX
MCScanX的安装请参考网上的教程
数据准备
首先去下载物种的基因组、蛋白、注释文件，推荐去ensembl上面下载
提取最长转录本（以ensembl上的拟南芥数据为例）
去除gff3文件中ID部分多余字符
# 基因结构文件 Arabidopsis_thaliana.TAIR10.47.gff3
# 基因组序列文件 Arabidopsis_thaliana.TAIR10.dna.toplevel.fa
# cds序列文件 Arabidopsis_thaliana.TAIR10.cds.all.fa
# 蛋白序列文件 Arabidopsis_thaliana.TAIR10.pep.all.fa
## 去除gff3文件中ID部分多余字符
cp Arabidopsis_thaliana.TAIR10.47.gff3 Ath.gff3
sed -i 's/=gene:/=/g' Ath.gff3
sed -i 's/=transcript:/=/g' Ath.gff3
sed -i 's/=CDS:/=/g' Ath.gff3
基于gff3提取最长cds序列ID，并过滤gff3文件
perl gff_longest.pl Ath.gff3 Ath.id Ath_final.gff3
#gff_longest.pl
#!/usr/bin/env perl
use strict;
my $a = scalar(@ARGV);
print "$a---";
die "perl $0 <ensembl.gff3> <out.longestID> <out.ensembl_longest.gff3>\n" unless (@ARGV == 3);
my $gff = shift;
my $otab = shift;
my $ogff = shift;
## read GFF and get ID info
my %mRNA2cds;
my %gene2mRNA;
my %cdslen;
open GFF, $gff or die $!;
while(<GFF>){
chomp;
next if /^#/;
my @t = split /\t/;
next unless ($t[8] =~ /ID=/ && $t[8] =~ /Parent=/);
$t[8] =~ /ID=([^;]+)/;
my $id1 = $1;
$t[8] =~ /Parent=([^;]+)/;
my $id2 = $1;
if($t[2] eq "mRNA"){
$gene2mRNA{$id2}{$id1} = 1;
}
if($t[2] eq "CDS"){
$mRNA2cds{$id2} = $id1;
## get cds length;
$cdslen{$id1} += $t[4] - $t[3] + 1;
}
}
close GFF;
## get longest cds and mRNA ID;
open O1, ">$otab" or die $!;
my %lmRNA;
foreach my $geneID (keys %gene2mRNA) {
my @mRNA = keys %{$gene2mRNA{$geneID}};
my $len = 0;
my $lmRNA;
foreach my $mRNA (@mRNA){
my $cds = $mRNA2cds{$mRNA};
next if (!exists $cdslen{$cds});
my $lcds = $cdslen{$cds};
if($lcds > $len){
$len = $lcds;
$lmRNA = $mRNA;
}
}
next if $len == 0;
print O1 "$geneID\t$lmRNA\t$mRNA2cds{$lmRNA}\n";
$lmRNA{$lmRNA} = $mRNA2cds{$lmRNA};
}
close O1;
#需要完整版请联系小编
提取最长cds ID列表
awk '{print $3}' Ath.id > Ath_mRNA.id
基于最长cds ID信息提取蛋白质序列文件
seqtk subseq Arabidopsis_thaliana.TAIR10.pep.all.fa Ath_mRNA.id > Ath.pep.fasta
MCScanX共线性分析(种内)
# 构建diamond数据库
diamond makedb --in Ath.pep.fasta --db Ath.pep.fasta
#使用diamond blastp进行比对
diamond blastp --query Ath.pep.fasta --db Ath.pep.fasta --out At.blast --outfmt 6 --max-target-seqs 5 --evalue 1e-10
## 过滤gff3文件
awk '$3=="mRNA"' At_final.gff3 | awk -F ";" '{print $1}'|awk '{print $1"\t"$9"\t"$4"\t"$5}'|sed 's/ID=//' > At.gff
## 种内基因复制模式分析,软件安装的位置
/pub/software/MCScanX/duplicate_gene_classifier ./At > At.dup_class.log
## 种内共线性分析
/pub/software/MCScanX/MCScanX ./At
生成文件Ath.collinearity
MCScanX共线性分析(种间)
diamond makedb --in B.pep.fasta --db B.pep.fasta
diamond makedb --in Ath.pep.fasta --db Ath.pep.fasta
## 使用diamond blastp进行比对
diamond blastp --query Ath.pep.fasta --db B.pep.fasta --out tmp.Ath_B.blast --outfmt 6 --max-target-seqs 5 --evalue 1e-10
diamond blastp --query B.pep.fasta --db Ath.pep.fasta --out tmp.B_Ath.blast --outfmt 6 --max-target-seqs 5 --evalue 1e-10
cat tmp.*.blast > B-Ath.blast
## 过滤gff3文件
awk '$3=="mRNA"' B_final.gff3 | awk -F ";" '{print $1}'|awk '{print $1"\t"$9"\t"$4"\t"$5}'|sed 's/ID=//' > B.gff
awk '$3=="mRNA"' Ath_final.gff3 | awk -F ";" '{print $1}'|awk '{print $1"\t"$9"\t"$4"\t"$5}'|sed 's/ID=//' > Ath.gff
cat B.gff Ath.gff > B-Ath.gff
## 种间共线性分析
/pub/software/MCScanX/MCScanX ./B-Ath
生成文件B-Ath.collinearity
可视化（使用MCScanX 自带功能）
dot_plotter
java dot_plotter -g B-Ath.gff -s B-Ath.collinearity -c B-Ath.dot.ctl -o B-Ath.dot.png
#B-Ath.dot.ctl
800 //dimension (in pixels) of x axis
800 //dimension (in pixels) of y axis
LG1,LG2,LG3,LG4,LG5 //chromosomes in x axis
chr1 chr2 chr3 chr4 chr5 //chromosomes in y axis
synteny_plotter
java dual_synteny_plotter -g B-Ath.gff -s B-Ath.collinearity -c B-Ath.dual_synteny.ctl -o B-Ath.dual_synteny.png
#B-Ath.dual_synteny.ctl
200 //plot width (in pixels)
800 //plot height (in pixels)
LG1,LG2,LG3,LG4,LG5 //chromosomes in the left column
chr1,chr2,chr3,chr4,chr5 //chromosomes in the right column
circle_plotter
java circle_plotter -g B-Ath.gff -s B-Ath.collinearity -c B-Ath.circle.ctl -o B-Ath.circle.png
#B-Ath.circle.ctl
800 //plot width and height (in pixels)
LG1,LG2,LG3,LG4,LG5,chr1,chr2,chr3,chr4,chr5 //chromosomes in the circle
bar_plotter
java bar_plotter -g B-Ath.gff -s B-Ath.collinearity -c B-Ath.bar.ctl -o B-Ath.bar.png
#B-Ath.bar.ctl
800 //dimension (in pixels) of x axis
800 //dimension (in pixels) of y axis
LG1,LG2,LG3,LG4,LG5 //reference chromosomes
chr1,chr2,chr3,chr4,chr5 //target chromosomes
详细信息请参考MCScanX 官方文档：https://github.com/wyp1125/MCScanX
