參考：https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version)
jcvi安裝
可以從官網(wǎng)下載：https://codechina.csdn.net/mirrors/tanghaibao/jcvi?utm_source=csdn_github_accelerator
還有一些教程可以參考：https://blog.csdn.net/u012110870/article/details/105857278/
http://m.itdecent.cn/p/6f1a191a0c3a
當(dāng)然比較簡單的就是從conda下載啦
我是用pip下載的，教程很多都是用pip2。如果沒下pip2，但是用pip2的話，也會報錯：
ERROR:Command errored out with exit status 1: python setup.py egg_info Check the logs for command output.
如果網(wǎng)速差，前面下得挺好，后面就突然報錯：
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/
中國科技大學(xué) https://pypi.mirrors.ustc.edu.cn/simple/
華中科技大學(xué)：http://pypi.hustunique.com/
山東理工大學(xué)：http://pypi.sdutlinux.org/
豆瓣：http://pypi.douban.com/simple/
數(shù)據(jù)下載和獲取
需要的數(shù)據(jù)有：GFF3注釋文件、CDS序列文件
可下載數(shù)據(jù)來源：phytozome數(shù)據(jù)庫、ensembl plants數(shù)據(jù)庫、NCBI數(shù)據(jù)庫、文獻(xiàn).....
我用的是自己測序的廣州相思子（后續(xù)用AC代替）和在NCBI下的非洲相思子（后續(xù)用AP代替）
數(shù)據(jù)處理
將.gff3文件轉(zhuǎn)換成.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文件去重復(fù)
python -m jcvi.formats.bed uniq AC.bed
python -m jcvi.formats.bed uniq AP.bed
如果你要利用gff文件從基因組文件提取出cds或者蛋白質(zhì)序列的話，可以用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 #格式轉(zhuǎn)換
共線性分析
教程是：
python -m jcvi.compara.catalog ortholog grape peach --no_strip_names
（在這里可以加多一條“--cscore=.99”命令，過濾掉一些噪音）
因?yàn)槊罾餂]有input文件，我就很奇怪。然后看到了后續(xù)的輸出結(jié)果：
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 應(yīng)該是數(shù)據(jù)的前綴名，所以我就改為了
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文件里找不到。很奇怪，因?yàn)榻坛陶f的是在.gff3文件里把mRNA提出來轉(zhuǎn)換成.bed文件。我打算再做一次提cds序列的。
還是不行。
最終找到原因：.bed文件和.cds文件前綴要一樣，ortholog命令會從.cds文件中找.bed文件里的序列名，所以兩個文件的序列名字要一樣，要不然就會報錯。
最后的結(jié)果是
畫共線性線圖
需要三個文件：
- seqids: 需要展現(xiàn)哪些序列
- layout: 不同物種的在圖上的位置
- .simple: 從.anchors文件創(chuàng)建的更簡化格式
.simple文件
目的：從anchor文件中提取一些子集，用一種簡潔的形式概括。
python -m jcvi.compara.synteny screen --minspan=30 --simple AC.AP.anchors AC.AP.anchors.new
得到結(jié)果：
Before: 571 blocks, After: 286 blocks
seqids文件
目的：告訴軟件應(yīng)畫多少組染色體。在這里可以移除小scaffolds。第一、二行分別為AC、AP的染色體數(shù)
touch seqids # 創(chuàng)建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。前各列分別表示位置、旋轉(zhuǎn)角度、顏色、物種名稱、垂直對齊、bed文件。軌道0現(xiàn)在是AC，軌道1現(xiàn)在是AP。下一節(jié)指定在軌道之間畫什么邊。要求使用文件中的信息繪制軌道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)
