比较基因组学分析的多功能分析套件JCVI

实验室牛马

详细使用说明https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version)

文章链接DOI: https://doi.org/10.1002/imt2.211

安装

conda install jcvi last

使用

1.gff转成bed格式

python -m jcvi.formats.gff bed --type=mRNA --key=ID A.gff3 > ath.bed

2.去重复

python -m jcvi.formats.bed uniq ath.bed
seqkit grep -f <(cut -f 4 a.uniq.bed ) A.cdna.all.fa | seqkit seq -i > a.cds

3.共线性分析

python -m jcvi.compara.catalog ortholog --no_strip_names ath osa
#--no_strip_names 识别区分.1

4.宏观共线性分析

需要准备两个额外的文件seqids和layout

首先是seqids文件，它告诉绘图仪要包含哪些染色体集。第一系包含19条葡萄染色体，第二系包含8条桃染色体。

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，整个画布在x轴上为0-1，在y轴上为0-1。首先，前三列指定不同基因组的位置。然后旋转，颜色，标签，垂直对齐（va），然后是基因组bed文件。第一行现在是葡萄，第二行现在是桃。edges是指定要在基因组间绘制共线性关系。e、 0，1要求使用.simple文件中的信息在第一行和第二行两个基因组之间绘制共线性关系。

# y, xstart, xend, rotation, color, label, va,  bed
 .6,     .1,    .8,       0,      , Grape, top, grape.bed
 .4,     .1,    .8,       0,      , Peach, top, peach.bed

# edges
e, 0, 1, grape.peach.anchors.simple
e, 1, 2, grape.peach.anchors.simple
e, 0, 2, grape.peach.anchors.simple
#最后一行不要留空行，可能会报错

最后，将共线性分析得到的synteny 文件转化为.simple的格式

python -m jcvi.compara.synteny screen --simple grape.peach.anchors grape.peach.anchors.new

绘图

python -m jcvi.graphics.karyotype seqids layout

如果我们想突出显示一个特定的块呢？我们应该进入.simple文件，找到相关的块。请注意.simple文件中的每一行都是一个synneny块，带有start和stop grape gene，然后是start和stop peach gene，最后两列是score和orientation。

$ vi grape.peach.anchors.simple 
GSVIVT01012228001   GSVIVT01012173001   Prupe.1G281700.1    Prupe.1G288300.1    72  -
g*GSVIVT01012028001 GSVIVT01000604001   Prupe.1G299800.1    Prupe.1G340600.1    518 +
GSVIVT01000603001   GSVIVT01000557001   Prupe.1G239100.1    Prupe.1G242900.2    51  +  
...
(save the file)
$ python -m jcvi.graphics.karyotype seqids layout

5.微共线性分析

如果我们想看看局部共线呢？局部共线性主要集中在基因水平，显示匹配区域以及对齐的基因模型。为此，我们需要计算基因水平匹配的布局：

python -m jcvi.compara.synteny mcscan grape.bed grape.peach.lifted.anchors --iter=1 -o grape.peach.i1.blocks

注意，选项--iter=1表示我们正在提取一个匹配每个葡萄区域的最佳区域。查看结果文件grape.peach.i1.blocks，包含葡萄作为第一列，桃为第二列。如果选项--iter设置为2，那么将有2个桃子的区域，依此类推。具体而言，这将有助于绘制基因组复制所产生的区域。

sed -n '427,500p' grape.peach.i1.blocks >blocks #利用sed命令提取想要展示的局部blocks信息

准备一个blocks.layout文件

# x,   y, rotation,   ha,     va,   color, ratio,            label
0.5, 0.6,        0, left, center,       m,     1,       grape Chr1
0.5, 0.4,        0, left, center, #fc8d62,     1, peach scaffold_1
# edges
e, 0, 1

合并两个物种的bed文件并绘图

cat grape.bed peach.bed > grape_peach.bed
python -m jcvi.graphics.synteny blocks grape_peach.bed blocks.layout

也可以做两个以上的染色体的共线性关系，需要合并blocks信息

python -m jcvi.compara.catalog ortholog --no_strip_names grape cacao && python -m jcvi.compara.synteny mcscan Dv.bed Dv.1A_7A.lifted.anchors --iter=1 -o Dv.1A_7A.i1.blocks
python -m jcvi.formats.base join grape.peach.i1.blocks grape.cacao.i1.blocks --noheader | cut -f 1,2,4 > grape.blocks
#也可以同时合并3个
#python -m jcvi.formats.base join Dv.1A_7A.i1.blocks Dv.1B_7B.i1.blocks Dv.1D_7D.i1.blocks --noheader| cut -f 1,2,4,6 > Dv.blocks
sed -n  '5,7p' Dv.1A_7A.i1.blocks >blocks #将2-5行提取到b文件中
python -m jcvi.graphics.synteny blocks Dv.1A_7A.bed blocks.layout

微信扫一扫分享文章