【比較基因組】McScan jcvi比較兩個(gè)基因組共線性細(xì)節(jié)記錄

軟件的安裝

Python版McScan(jcvi工具包):https://github.com/tanghaibao/jcvi

以前只有python2,現(xiàn)在已有python3版本,建議用py3。安裝可用pip:

pip install jcvi
##或開發(fā)版
pip install git+git://github.com/tanghaibao/jcvi.git

pip可能會(huì)安裝很慢。建議還是用conda,要快很多,最好新建環(huán)境。

conda install -c bioconda jcvi

這時(shí),你已經(jīng)能使用命令,表面上安裝成功了,實(shí)際上可能還缺少很多依賴。比如last,latex,dvipng等。否則在后面運(yùn)行過程,可能遇到如下錯(cuò)誤:

##未安裝last
/bin/bash: lastdb: command not found
##未安裝latex、dvipng
RuntimeError: Failed to process string with tex because latex could not be found

只有一個(gè)個(gè)解決,有的可以直接conda(如last),有些則需要編譯,若有root權(quán)限,倒也好辦。

conda install -c bioconda last
sudo yum install -y  texlive texlive-latex texlive-xetex texlive-collection-latexrecommended
sudo yum install dvipng

基因組的準(zhǔn)備

若是已知物種,直接可從公共數(shù)據(jù)庫中下載gff和cds序列,jcvi提供了下載方式:

$ python -m jcvi.apps.fetch
Usage:
    python -m jcvi.apps.fetch ACTION


Available ACTIONs:
        bisect | Determine the version of the accession by querying entrez
       ensembl | Retrieve genomes and annotations from ensembl
        entrez | Fetch records from entrez using a list of GenBank accessions
     phytozome | Retrieve genomes and annotations from phytozome
    phytozome9 | Retrieve genomes and annotations from phytozome version 9.0 (legacy)
           sra | Retrieve files from SRA via the sra-instant FTP

比如從Phytozome下載,要提前注冊(cè)好,如下命令提示輸入賬號(hào)密碼。

python -m jcvi.apps.fetch phytozome Vvinifera,Ppersica

下載后無需解壓。

自己準(zhǔn)備的基因組數(shù)據(jù)也只需gff3和cds.fa(蛋白序列也可)。

gff3只保留染色體水平的ID,如:

grep '^chr' Vvinifera_145_Genoscope.12X.gene.gff3 > apricot.filter.gff3

gff3文件轉(zhuǎn)化bed文件時(shí)注意type和key類型對(duì)應(yīng)gff中第三列和第九列信息。type一般為mRNA,但是key注意你的gff文件是取Name還是ID。如:

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

若后續(xù)作圖仍報(bào)錯(cuò),可嘗試去除fasta ID中多余的描述信息(我自己不用也可跑通)。如:

# clean headers to remove description fiedls from Phytozome FASTA files.
python -m jcvi.formats.fasta format --sep="|" Vvinifera_145_cds.fa.gz grape.cds
python -m jcvi.formats.fasta format --sep="|" Ppersica_139_cds.fa.gz peach.cds

一些細(xì)節(jié)

  • 結(jié)果文件
    last比對(duì)結(jié)果,last.filtered比對(duì)過濾串聯(lián)重復(fù)和低分比對(duì)結(jié)果,anchors: 高質(zhì)量的共線性塊,lifted.anchors增加額外錨點(diǎn)的最終共線性區(qū)塊,simple簡化的anchors文件。anchors文件中每個(gè)共線性區(qū)塊以###分隔, 第一和第二列分別是兩基因組的基因ID,第三列BLAST的bit score,越大可靠性越高。

  • 調(diào)圖細(xì)節(jié)
    兩個(gè)配置文件seqid(展示染色體),layout(序列位置)。
    seqid文件中,基因組的染色體編號(hào)與其gff3文件一致(按大小順序?qū)?,而非gff文件染色體順序,轉(zhuǎn)化bed時(shí)軟件會(huì)排序)。如:

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

layout文件繪制一些選項(xiàng),若要個(gè)性化,多多修改嘗試(尤其時(shí)三個(gè)物種比較時(shí))。如:

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

若要突出顯示某一共線性區(qū)塊,可以在anchors.simple文件對(duì)應(yīng)的區(qū)塊前添加g*(g代表綠色,也可以改成其他顏色,如紅色r)。

建議和示例

建議先用示例數(shù)據(jù)跑一遍,也很快。再換自己的數(shù)據(jù),報(bào)錯(cuò)對(duì)照著尋找原因,總能解決。

示例代碼:

# 準(zhǔn)備數(shù)據(jù)(輸入帳號(hào)密碼)
python -m jcvi.apps.fetch phytozome Vvinifera,Ppersica

#去掉chr以外的序列 
grep '^chr' Vvinifera_145_Genoscope.12X.gene.gff3 > apricot.filter.gff3  

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

#reformat fasta
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

#identify blocks
python -m jcvi.compara.catalog ortholog grape peach --no_strip_names

#plot dotplot
python -m jcvi.graphics.dotplot grape.peach.anchors

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

##prepare for seqid and layout file

#  plot synteny
python -m jcvi.graphics.karyotype seqid layout
image.png

Ref:
http://m.itdecent.cn/p/a748d3a5421d
https://www.cnblogs.com/zhanmaomao/p/12525411.html
https://sr-c.github.io/2019/01/11/jcvi-MCscan/

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡書系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容