jcvi画两物种的共线性图

参考：https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version)

jcvi安装

可以从官网下载：https://codechina.csdn.net/mirrors/tanghaibao/jcvi?utm_source=csdn_github_accelerator
还有一些教程可以参考：https://blog.csdn.net/u012110870/article/details/105857278/
https://www.jianshu.com/p/6f1a191a0c3a
当然比较简单的就是从conda下载啦
我是用pip下载的，教程很多都是用pip2。如果没下pip2，但是用pip2的话，也会报错：

ERROR:Command errored out with exit status 1: python setup.py egg_info Check the logs for command output.

报错内容

如果网速差，前面下得挺好，后面就突然报错：

pip._vendor.urllib3.exceptions.ReadTimeoutError: HTTPSConnectionPool(host='files.pythonhosted.org', port=443): Read timed out.

解决方法：用临时镜像会特别快

pip --default-timeout=100 install ./jcvi -i http://pypi.douban.com/simple/ --trusted-host pypi.douban.com #豆瓣的镜像

其他镜像：
清华：https://pypi.tuna.tsinghua.edu.cn/simple
阿里云：http://mirrors.aliyun.com/pypi/simple/
中国科技大学 https://pypi.mirrors.ustc.edu.cn/simple/
华中科技大学：http://pypi.hustunique.com/
山东理工大学：http://pypi.sdutlinux.org/
豆瓣：http://pypi.douban.com/simple/

数据下载和获取

需要的数据有：GFF3注释文件、CDS序列文件
可下载数据来源：phytozome数据库、ensembl plants数据库、NCBI数据库、文献.....
我用的是自己测序的广州相思子（后续用AC代替）和在NCBI下的非洲相思子（后续用AP代替）

数据处理

将.gff3文件转换成.bed文件

python -m jcvi.formats.gff bed --type=mRNA --key=ID AC.gff3 -o AC.bed
python -m jcvi.formats.gff bed --type=mRNA --keyID AP.gff3 -o AP.bed

然后将.bed文件去重复

python -m jcvi.formats.bed uniq AC.bed
python -m jcvi.formats.bed uniq AP.bed

如果你要利用gff文件从基因组文件提取出cds或者蛋白质序列的话，可以用gffread命令：

conda install gffread -y #安装
gffread target.gff3 -g target.fa -y target.pep.fa #翻译后蛋白序列
gffread target.gff3 -g target.fa -x target.cds.fa #获得cds序列
gffread target.gff3 -g target.fa -w target.exons.fa  #获得exon序列
gffread target.gff3 -T -o target.gtf #格式转换

共线性分析

教程是：

python -m jcvi.compara.catalog ortholog grape peach --no_strip_names

（在这里可以加多一条“--cscore=.99”命令，过滤掉一些噪音）

因为命令里没有input文件，我就很奇怪。然后看到了后续的输出结果：

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

我就猜想 species_a 和 species_b 应该是数据的前缀名，所以我就改为了

python -m jcvi.compara.catalog ortholog AC AP --no_strip_names

但是报错了

10:09:41 [blastfilter] lcl|NW_020874292.1_cds_XP_027356858.1_2215 not in AP.bed
......
ValueError: A total of 0 anchor was found. Aborted.

报错说的是CDS序列在.bed文件里找不到。很奇怪，因为教程说的是在.gff3文件里把mRNA提出来转换成.bed文件。我打算再做一次提cds序列的。

还是不行。

最终找到原因：.bed文件和.cds文件前缀要一样，ortholog命令会从.cds文件中找.bed文件里的序列名，所以两个文件的序列名字要一样，要不然就会报错。

最后的结果是

共线性散点图

画共线性线图

需要三个文件：

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

.simple文件

目的：从anchor文件中提取一些子集，用一种简洁的形式概括。

python -m jcvi.compara.synteny screen --minspan=30 --simple AC.AP.anchors AC.AP.anchors.new

得到结果：

Before: 571 blocks, After: 286 blocks

seqids文件

目的：告诉软件应画多少组染色体。在这里可以移除小scaffolds。第一、二行分别为AC、AP的染色体数

touch seqids # 创建seqids文件
vim seqids
chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11
scaffold_1, scaffold_2,scaffold_3,scaffold_4, scaffold_5, scaffold_6 ,scaffold_7, scaffold_8, scaffold_9, scaffold_10, scaffold_11
# i 键输入；Esc 键退出；输入“：wq” 保存并退出

layout文件

目的：告诉软件图在哪里、画什么。整个画布的x轴是0-1，y轴是0-1。前各列分别表示位置、旋转角度、颜色、物种名称、垂直对齐、bed文件。轨道0现在是AC，轨道1现在是AP。下一节指定在轨道之间画什么边。要求使用文件中的信息绘制轨道0和1之间的边。

# y, xstart, xend, rotation, color, label, va,  bed
 .6,     .1,    .11,       0,      , AC, top, AC.bed
 .4,     .1,    .11,       0,      , AP, top, AP.bed
# edges
e, 0, 1, AC.AP.anchors.simple

最后画图

python -m jcvi.graphics.karyotype seqids layout