參考: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://www.reibang.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.
如果網(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/
中國科技大學 https://pypi.mirrors.ustc.edu.cn/simple/
華中科技大學:http://pypi.hustunique.com/
山東理工大學: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ù)庫、文獻.....
我用的是自己測序的廣州相思子(后續(xù)用AC代替)和在NCBI下的非洲相思子(后續(xù)用AP代替)
數(shù)據(jù)處理
將.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文件,我就很奇怪届垫。然后看到了后續(xù)的輸出結果:
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 應該是數(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文件里找不到。很奇怪装处,因為教程說的是在.gff3文件里把mRNA提出來轉換成.bed文件误债。我打算再做一次提cds序列的。
還是不行妄迁。
最終找到原因:.bed文件和.cds文件前綴要一樣寝蹈,ortholog命令會從.cds文件中找.bed文件里的序列名,所以兩個文件的序列名字要一樣登淘,要不然就會報錯箫老。
最后的結果是
畫共線性線圖
需要三個文件:
- 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
得到結果:
Before: 571 blocks, After: 286 blocks
seqids文件
目的:告訴軟件應畫多少組染色體耍鬓。在這里可以移除小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峭判。前各列分別表示位置开缎、旋轉角度、顏色林螃、物種名稱奕删、垂直對齊、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)