参考: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
但是报错了
ZeroDivisionError: float division by zero
原因是:bed文件中染色体名称跟seqid中的染色体名称不一致。可以通过.gff文件查看染色体名称。
最后就画好啦~
Load file `layout`
Load file `Abrus_cantoniensis.bed`
Load file `Abrus_precatorius.bed`
Figure saved to `karyotype.pdf` (2400px x 2100px)
参考教程的图
网友评论