共线性

共线性片段是指一个同一个物种内部或者两个物种之间，由基因组复制、染色体复制、大片段复制以及物种分化而产生的大片段同源现象。同源片段内部基因排序保守，意味着功能上也可能保守。

我们做共线性分析的目的就是为了鉴定物种内或者物种间的共线性片段。

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