共线性分析可以很好地解释进化关系和多倍化事件。今天主要测试的是Python版McScan（jcvi工具包），这个包很强大，是从MCScanx升级而来的基因组共线性分析软件。

一、安装

conda create -y -c bioconda -c conda-forge -n jcvi jcvi   //为了怕依赖包冲突，新建了一个jcvi的环境
source activate jcvi
pip install git+git://github.com/tanghaibao/jcvi.git

或者 
pip install jcvi

然后，安装额外的依赖环境：Kent tools，BEDTOOLS，EMBOSS，LAST，LaTex。

Two external command-line tools are needed:

• LASTAL. Compile and then move lastal and lastdb on your PATH.

其中，BEDTOOLS, EMBOSS 和 LAST 可以用 conda 来安装。

conda install -y -n jcvi -c bioconda bedtools emboss last
pip install latex


Kent自己编译一下：
wget http://hgdownload.cse.ucsc.edu/admin/jksrc.zip
unzip jksrc.zip
cd kent/src/lib
make

二、实例

Let's assume that you want to perform synteny comparison between grape and peach. This workflow is also available as an online Jupyter notebook (credits: Wayne Decatur @fomightez).

Download data

For MCscan to work, we'll need sequences (FASTA format) and coordinates (BED format). For this example, let's get these data from Phytozome.

The first time you try to run the fetch command, the script will ask for your Phytozome login (jcvi will not store or send your Phytozome login anywhere).

$ python -m jcvi.apps.fetch phytozome
Phytozome Login: xxxxxxxx
Phytozome Password:

If you don't have a login, register one here. Then you'll be able to see the current list of genomes.

$ python -m jcvi.apps.fetch phytozome
...
             ZmaysPH207                      Zmays
                Zmarina                  Vvinifera
               Vcarteri                  Tpratense
                 Tcacao                   Sviridis
             Stuberosum                  Spurpurea
             Spolyrhiza            Smoellendorffii
          Slycopersicum                   Sitalica
                Sfallax                   Sbicolor
              Rcommunis                  Pvulgaris
              Pvirgatum               Ptrichocarpa
               Ppersica                    Ppatens
                Phallii                  Othomaeum
         OsativaKitaake                    Osativa
           Olucimarinus                Mtruncatula
              MspRCC299           MpusillaCCMP1545
            Mpolymorpha                  Mguttatus
             Mesculenta                 Mdomestica
...
$ python -m jcvi.apps.fetch phytozome Vvinifera,Ppersica
...
$ ls
Ppersica_298_v2.1.cds.fa.gz              Vvinifera_145_Genoscope.12X.cds.fa.gz
Ppersica_298_v2.1.gene.gff3.gz           Vvinifera_145_Genoscope.12X.gene.gff3.gz

1. gff文件转bed格式

$ python -m jcvi.formats.gff bed --type=mRNA --key=Name Vvinifera_145_Genoscope.12X.gene.gff3.gz -o grape.bed
$ python -m jcvi.formats.gff bed --type=mRNA --key=Name Ppersica_298_v2.1.gene.gff3.gz -o peach.bed

--type：gff文件中第三列
--key：第九列信息前缀

I also like to reformat the Phytozome FASTA files as well.

$ python -m jcvi.formats.fasta format Vvinifera_145_Genoscope.12X.cds.fa.gz grape.cds
$ python -m jcvi.formats.fasta format Ppersica_298_v2.1.cds.fa.gz peach.cds

注意：偶尔会有某些基因组注释包含大量的isoforms。MCscan不知道这些其实是相同的基因，而是把它们当作不同的基因，看起来像是tandem gene duplicates。通常情况下，这不会是一个问题，因为synsyny管道会移除串联基因复制体中除了一个基因外的所有基因，并且默认在20个基因的大窗口中寻找synteny，但这种假设会随着太多的isoforms而失效。因此，如果isoforms过多，我们建议在上面的BED生成命令中增加一个选项——--primary_only，每个基因只保留一个isoform。例如:

$ python -m jcvi.formats.gff bed --type=mRNA --key=Name --primary_only Vvinifera_145_Genoscope.12X.gene.gff3.gz -o grape.bed

deduplication
如果没有事先准备全长转录本的数据，也可以通过如下方式：

python -m jcvi.formats.bed uniq ath.bed
python -m jcvi.formats.bed uniq osa.bed
seqkit grep -f <(cut -f 4 ath.uniq.bed )Arabidopsis_thaliana.TAIR10.cdna.all.fa.gz | seqkit seq -i > ath.cds
seqkit grep -f <(cut -f 4 ath.uniq.bed )Arabidopsis_thaliana.TAIR10.pep.all.fa.gz | seqkit seq -i > ath.pep
seqkit grep -f <(cut -f 4 osa.uniq.bed)  Oryza_sativa.IRGSP-1.0.cdna.all.fa.gz| seqkit seq -i  > osa.cds
seqkit grep -f <(cut -f 4 osa.uniq.bed )Oryza_sativa.IRGSP-1.0.pep.all.fa.gz | seqkit seq -i  > osa.pep
mkdir -p cds && cd cds
ln -s ../ath.cds ath.cds
ln -s ../ath.uniq.bed ath.bed
ln -s ../osa.cds osa.cds
ln -s ../osa.uniq.bed osa.bed

2. Pairwise synteny search jcvi.compara.catalog [ 共线性分析 ]

$ ll *.???
grape.bed 
grape.cds 
peach.bed 
peach.cds


$ python -m jcvi.compara.catalog     ortholog         grape peach         --no_strip_names
23:34:43 [synteny] Assuming --qbed=grape.bed --sbed=peach.bed
23:34:43 [base] Load file `grape.bed`
23:34:43 [base] Load file `peach.bed`
23:34:44 [blastfilter] Load BLAST file `grape.peach.last` (total 515965 lines)
23:34:44 [base] Load file `grape.peach.last`
23:34:46 [blastfilter] running the cscore filter (cscore>=0.70) ..
23:34:46 [blastfilter] after filter (379462->50584) ..
23:34:46 [blastfilter] running the local dups filter (tandem_Nmax=10) ..
23:34:47 [blastfilter] after filter (50584->21696) ..
...
Stats: Min=4 Max=1002 N=687 Mean=67.1863173216885 SD=110.64978224380248 Median=27.0 Sum=46157
NR stats: Min=4 Max=356 N=687 Mean=26.595342066957787 SD=43.0459201862291 Median=11.0 Sum=18271

它调用 LAST 进行比较，过滤 LAST 输出以删除串联重复（tandem duplications）和弱命中（weak hits）。在LAST输出上执行single linkage clustering，将 cluster anchors into synteny blocks。在运行结束时，您将看到synteny块的汇总统计信息。

$ ll grape.peach.*
grape.peach.anchors   
grape.peach.lifted.anchors    
grape.peach.last     
grape.peach.last.filtered  
grape.peach.pdf

• .last 基于LAST的比对结果（原始的last输出）；
• .last.filtered LAST的比对结果过滤串联重复和低分比对；
• .anchors is the seed synteny blocks (high quality), 高质量的共线性区块
• .lifted.anchors recruits additional anchors to form the final synteny blocks. 增加了额外的锚点，形成最终的共线性区块
• .pdf点图
  前两列是LAST鉴定的homologs，第三列是LAST的分数，它们在.anchor之前的步骤中使用，比如基于C-score的筛选，或者在串联重复序列中的一系列匹配中对匹配进行优先排序。

有时，当.anchors文件是“liftover”的输出，它丰富了synteny signal，如lift .anchors。例如，Certain pairs的第三列以L结尾，第一列为一个 基因组基因ID，第二列为另一个基因组的基因ID文件，即是两个基因组共线性区域内基因对应关系。不同的共线性区块用“###” 隔开。

GSVIVT01012008001 ppa002064m 1180L

L只是强调了一个事实，这些低质量的anchors接近高质量的anchors。

Pairwise synteny visualization [ 可视化 ]
如果没有得到grape.peach.pdf点图，可以通过如下命令得到：

$ python -m jcvi.graphics.dotplot grape.peach.anchors

得到grape.peach.pdf点图.

image.png

要理解点图的含义，可以水平或垂直地观察，并统计blocks数量。

• 水平方向上，每个 peach 区，最多有3个 grape 共线性区。
• 在垂直方向上，每个grape 区，最多有3个peach 共线性区。

它有助于 理解葡萄和桃共享一个genome triplication event 基因组三重增殖事件(“gamma”)事件，从而产生这种3:3模式 。

如果你仔细观察，这三个区域中的一个通常更强，它对应着两个基因组之间的同源区域( orthologous regions)。如果我们只想得到1:1的同源区域呢?我们只是重复之前的操作，添加一个选项 --cscore=.99。

C-score is defined by the ratio of LAST hit to the best BLAST hits to either the query and hit.

A C-score cutoff of .99 effectively filters the LAST hit to contain reciprocal best hit (RBH).

$ rm grape.peach.last.filtered 
$ python -m jcvi.compara.catalog ortholog grape peach --cscore=.99 --no_strip_names
$ python -m jcvi.graphics.dotplot grape.peach.anchors

image.png

我们也可以通过运行下面的命令．快速测试synteny模式是否确实是1:1，:

$ python -m jcvi.compara.synteny depth --histogram grape.peach.anchors

image.png

3. Macrosynteny visualization 宏观的可视化

看下别人paper中的效果（2011年的Nature Genetics上A.lyrata）：

image

用JCVI的画图模块实现这种效果，只不过还需要一些其它文件，创建如下三个文件

• seqids: 需要展现哪些序列；
• layout: 不同物种的在图上的位置；
• .simple: 从.anchors文件创建的更简化格式；

现在，让我们使用相同的synteny输出文件grape.peach.anchors进行另一种可视化。除了我们已经拥有的BED文件和synteny文件之外，我们还需要准备两个额外的文件。

（1）seqids文件

它告诉 plotter 选择对哪些染色体进行画图。这里，我们移除了unplaced and small scaffolds。第一行包含19条葡萄染色体，第二行包含8条桃染色体。

chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19
Pp01,Pp02,Pp03,Pp04,Pp05,Pp06,Pp07,Pp08

（2）layout文件

它告诉plotter在哪里绘制什么。

• First, three columns specify the position of the track.
• Then rotation, color, label, vertical alignment (va), and then the genome BED file. Track 0 is now grape, track 1 is now peach.
• The next stanza specifies what edges to draw between the tracks. e, 0, 1 asks to draw edges between track 0 and 1, using information from the .simple file.轨道0现在是葡萄，轨道1现在是桃子。下一节指定在轨道之间画什么边。' e, 0,1 '使用.simple文件的信息绘制轨道0和1之间的边。

# y, xstart, xend, rotation, color, label, va,  bed
 .6,     .1,    .8,       0,      , Grape, top, grape.bed
 .4,     .1,    .8,       0,      , Peach, top, peach.bed
# edges
e, 0, 1, grape.peach.anchors.simple

其中：

• 前三列分别指定了物种 1 和物种 2 染色体在图上 Y 轴的位置，X 轴的起始和终止位置，即我们共线性图在画布上的位置，整个画布的x轴是0-1　 y轴是0-1。；
• 第四列颜色可以默认，不填写；
• 第五列是旋转角度，默认零度则是画两条平行的；
• 第六列是标记物种名；
• 第七列是染色体 ID 标记位置；
• 第八列是 bed文件。
• 最后一行是画共线性的数据文件，其中 0 和 1 代表前述第一个物种和第二个物种，进行共线性的绘制。后面跟上的是相关的数据。

前两行展示了两个物种的基础绘图信息，如果画三个物种，那么可以写三行，以此类推。

（3）simple文件

这个命令生成一个.simple文件，它的形式比.anchors文件更简洁。

$ python -m jcvi.compara.synteny screen --minspan=30 --simple grape.peach.anchors grape.peach.anchors.new 

基因组共线性简化结果文件中，一行为一个共线性区块，前两列表示一个基因组中两个基因组之间的区域，与后两列（3,4列）另一个基因组中的两个基因之间的区域存在共线性关系。第5列物区域跨度，最后一列为: +为正向，-为反向

image.png

准备好所有输入文件后，我们绘图。

$ python -m jcvi.graphics.karyotype seqids layout

这就产生了我们的 karyotype 图。

image.png

美化

Further customization is possible。例如，改变seqids中染色体的顺序可以使其在视觉上更吸引人。同时，调整layout文件中的位置、颜色、标签等。

image

其中还有一些参数，例如--shadestyle=line（共线性区域是用曲线还是线），font是字体，--diverge是track色系的调整，等等。总之，一个好看的图要不断的进行调整。

突出显示一个特定的block
我们应该去.simple件中，找到相关的block。请注意，.simple文件中的每一行都是一个synteny块，有葡萄基因的start 和 stop，然后桃子基因的start 和 stop，最后两列是分数和方向。

$ vi grape.peach.anchors.simple 
GSVIVT01012228001   GSVIVT01012173001   Prupe.1G281700.1    Prupe.1G288300.1    72  -
g*GSVIVT01012028001 GSVIVT01000604001   Prupe.1G299800.1    Prupe.1G340600.1    518 +
GSVIVT01000603001   GSVIVT01000557001   Prupe.1G239100.1    Prupe.1G242900.2    51  +    
...
(save the file)


$ python -m jcvi.graphics.karyotype seqids layout

颜色值为也可以为16进制，*隔开； 没有添加颜色值的行默认灰色。然后，我们高亮显示了图中特定的synteny块(颜色为“g”绿色)。

image.png

选择 shade 的样式

$ python -m jcvi.graphics.karyotype seqids layout --shadestyle=line

注意，阴影现在变成直线，而不是以前的Bezier曲线。

image.png

添加第三个物种

我们可以在 karyotype 图上进行有创意的画图!首先，我们可以添加任意多的基因组。让我们添加可可。我们只需要重复下载和格式化可可基因组　(提取cacao.cds和cacao.bed)。

$ python -m jcvi.apps.fetch phytozome Tcacao
$ python -m jcvi.formats.gff bed --type=mRNA --key=Name Tcacao_233_v1.1.gene.gff3.gz -o cacao.bed
$ python -m jcvi.formats.fasta format Tcacao_233_v1.1.cds.fa.gz cacao.cds

假设我们对桃子和可可比较感兴趣。

$ python -m jcvi.compara.catalog ortholog peach cacao --cscore=.99 --no_strip_names
$ python -m jcvi.compara.synteny screen --minspan=30 --simple peach.cacao.anchors peach.cacao.anchors.new

现在我们有cacao.bed和peach.cacao.anchors.simple。很简单，我们需要将它们添加到 layout中。请注意这两行——第4行和第7行。Section# edges表示我们应该将轨道0 (grape)与1 (peach)连接起来，轨道1 (peach)与2 (cocoa)连接起来。

# y, xstart, xend, rotation, color, label, va,  bed
 .7,     .1,    .8,      15,      , Grape, top, grape.bed
 .5,     .1,    .8,       0,      , Peach, top, peach.bed
 .3,     .1,    .8,     -15,      , Cacao, bottom, cacao.bed
# edges
e, 0, 1, grape.peach.anchors.simple
e, 1, 2, peach.cacao.anchors.simple

记得在 seqids 中添加可可染色体(第3行)。记住，基因组在seqids中的顺序需要与layout中的顺序匹配。

chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19
scaffold_1,scaffold_2,scaffold_3,scaffold_4,scaffold_5,scaffold_6,scaffold_7,scaffold_8
scaffold_1,scaffold_2,scaffold_3,scaffold_4,scaffold_5,scaffold_6,scaffold_7,scaffold_8,scaffold_9,scaffold_10r

我们还可以自定义旋转，所以我们不再局限于枯燥的染色体水平排列。

$ python -m jcvi.graphics.karyotype seqids layout

image.png

修改配置，是可以控制每个track的位置的。

image

这样就可以三角性状显示。

image

4. Microsynteny visualization

如果我们想看看 local synteny 呢?局部共线性主要集中在基因水平上，显示匹配的区域和对齐的基因模型。为此，我们需要计算基因水平上匹配的布局:

What if we want to look at local synteny? Local synteny mostly focuses on gene-level, showing matching regions along with aligned gene models. For this, we'll need to compute the layout for the gene-level matchings:

$ python -m jcvi.compara.synteny mcscan grape.bed grape.peach.lifted.anchors --iter=1 -o grape.peach.i1.blocks

请注意，选项—— --iter=1 表示我们正在提取与每个葡萄区最佳匹配区。如果你看一下结果文件 grape.peach.i1.blocks ，它葡萄作为第一列，桃子作为第二列。如果选项——iter被设置为2，那么将有2个peach区域，以此类推。具体来说，这将有助于绘制由基因组重复产生的区域。

OK. The file grape.peach.i1.blocks 包含大量的局部区域!我们可以选择任意的区域来进行可视化。

$ head -50 grape.peach.i1.blocks > blocks

最后，我们需要一个布局文件，就像macro-synteny图一样。blocks.layout:

# x,   y, rotation,   ha,     va,   color, ratio,            label
0.5, 0.6,        0, left, center,       m,     1,       grape Chr1
0.5, 0.4,        0, left, center, #fc8d62,     1, peach scaffold_1
# edges
e, 0, 1

通过下面的命令，我们可以生成一个local synteny 图:

$ cat grape.bed peach.bed > grape_peach.bed
$ python -m jcvi.graphics.synteny blocks grape_peach.bed blocks.layout

image.png

同样，我们也可以选择不同的阴影风格:

$ python -m jcvi.graphics.synteny blocks grape_peach.bed blocks.layout --shadestyle=line

image.png

可选的“arrow”字形样式:

$ python -m jcvi.graphics.synteny blocks grape_peach.bed blocks.layout --glyphstyle=arrow

image.png

有时我们想要根据 blocks 文件中显示的同源性来匹配基因的颜色。match the colors of the genes with respect to their homology as indicated in the blocks file.

$ python -m jcvi.graphics.synteny blocks grape_peach.bed blocks.layout --glyphcolor=orthogroup

image.png

最后，出于调试目的，我们还可以显示基因名称。对于密集区来说不是很漂亮，但有时也很方便。注意——genelabelsize 6设置标签字体大小为6。

$ python -m jcvi.graphics.synteny blocks grape_peach.bed blocks.layout --genelabelsize 6

image.png

Microsynteny getting fancy

也可以做两个以上的匹配区域! 类似于macro-synteny图，首先使用python -m jcvi.compara.synteny mcscan构造 multi-synteny块，然后修改blocks.layout文件，以指示更多区域以及区域之间的边缘。

让我们继续以葡萄、桃子和可可为例。现在让我们假设我们想要构建一个包含三个基因组的blocks文件。如果有一个单一的“ 'reference'就好了，例如，让我们挑选葡萄，并将桃子和可可分别与葡萄对齐。我们已经做了桃子，现在可可也做了类似的。

$ python -m jcvi.compara.catalog ortholog grape cacao --cscore=.99
$ python -m jcvi.compara.synteny mcscan grape.bed grape.cacao.lifted.anchors --iter=1 -o grape.cacao.i1.blocks

组合 blocks 文件。

$ python -m jcvi.formats.base join grape.peach.i1.blocks grape.cacao.i1.blocks --noheader | cut -f1,2,4,6 > grape.blocks
$ head -50 grape.blocks > blocks2

最终的blocks 文件

GSVIVT01012261001   .   .
GSVIVT01012259001   .   .
GSVIVT01012258001   .   .
GSVIVT01012257001   .   .
GSVIVT01012255001   Prupe.1G290900.1    Thecc1EG011472t1
GSVIVT01012253001   Prupe.1G290800.2    Thecc1EG011473t1
GSVIVT01012252001   Prupe.1G290700.1    Thecc1EG011474t1
GSVIVT01012250001   Prupe.1G290600.1    Thecc1EG011475t1
GSVIVT01012249001   Prupe.1G290500.1    Thecc1EG011478t1
GSVIVT01012248001   Prupe.1G290400.1    Thecc1EG011482t1

现在我们差不多准备好了。我们的新布局文件 blocks2.layout:

# x,   y, rotation,     ha,     va, color, ratio,            label
0.5, 0.6,        0, center,    top,      ,     1,       grape Chr1
0.3, 0.4,        0, center, bottom,      ,    .5, peach scaffold_1
0.7, 0.4,        0, center, bottom,      ,    .5, cacao scaffold_2
# edges
e, 0, 1
e, 0, 2

注意，我们已经将桃子和可可的ratio改为0.5，让它们适合画布。对于高质量图片，您很可能需要调整这个布局文件。就像macro synteny部分中的布局文件一样，edges stanza将grape(列0)与peach(列1)以及grape(列0)与cocoa(列2)连接起来。现在我们运行:

$ cat grape.bed peach.bed cacao.bed > grape_peach_cacao.bed
$ python -m jcvi.graphics.synteny blocks2 grape_peach_cacao.bed blocks2.layout

image.png

https://github.com/tanghaibao/jcvi/
https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version)
https://www.jianshu.com/p/7b6699153ecf