JCVI(MCscan Python version)使用记录

文章仅是记录自己的学习使用,有错误请指出,我立刻改正!

官方说明
https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version)
更多说明:
https://www.jianshu.com/p/f7971dbf5f85
https://www.jianshu.com/p/d25fb05ada97
https://www.jianshu.com/p/6f1a191a0c3a
https://zhuanlan.zhihu.com/p/353490154
https://www.omicsclass.com/article/481
https://github.com/tanghaibao/jcvi/issues/387

  • 关于JCVI的使用,唐老师在github上已经有很详细的介绍,网络上也有很多教程。
  • 这篇文章的主要是记录自己在安装和使用JCVI中产生的问题。

一、JCVI的安装

  • JCVI依赖的东西非常多,不管是直接编译还是用conda都需要额外安装别的软件,建议还是用conda安装。
  • 我自己是在ubuntu的虚拟机上安装成功了,因此后面的演示也都是在虚拟机上。
conda create -y -c bioconda -n JCVI jcvi
conda activate JCVI
conda install -c bioconda bedtools
conda install seqkit emboss
conda install last

(一)Latex安装

  • 我在安装Latex过程中遇到了各种问题,后来还是主要参考了https://www.jianshu.com/p/6f1a191a0c3a
  • 用yum安装Latex后,我在使用JCVI过程中发现字体缺少类型报错,如果只是运行- JCVI查看共线性结果的话,这个报错是没什么影响的。但如果是要可视化的话,还需要额外下载字体文件。用非root安装的Latex没有报错。

1、有root权限

# ubuntu
sudo apt-get install -y texlive texlive-latex-extra texlive-latex-recommended
# centos
sudo yum install -y  texlive texlive-latex texlive-xetex texlive-collection-latexrecommended

2、无root权限

详见https://www.jianshu.com/p/a5e9ca5886a4

二、JCVI的运行

  • 从网上下载了Arabidopsis_halleri和Arabidopsis_lyrata两种拟南芥的cdna序列和gff3文件
  • fig1.png

(一)GFF3转BED文件

python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id \
Arabidopsis_halleri.Ahal2.2.44.gff3 > halleri_all.bed
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id \
Arabidopsis_lyrata.v.1.0.44.gff3 > lyrata_all.bed
#根据cDNA文件的序列名以及类型调整type和key
#转换完成的bed文件和fasta文件中的序列名要保持一致

(二)去重复

python -m jcvi.formats.bed uniq halleri_all.bed
python -m jcvi.formats.bed uniq lyrata_all.bed
#仅保留最长且唯一序列

(三)提取CDS序列

seqkit grep -f <(cut -f 4 halleri_all.uniq.bed ) \
Arabidopsis_halleri.Ahal2.2.cdna.all.fa | seqkit seq -i > halleri.cds
seqkit grep -f <(cut -f 4 lyrata_all.uniq.bed ) \
Arabidopsis_lyrata.v.1.0.cdna.all.fa | seqkit seq -i > lyrata.cds
  • 这样就获得了每个物种唯一的cds和对应的bed文件
  • fig2.png

(四)提取共线信息

  • 这里首先要把之前获得的四个文件的名字改成以下形式,使JCVI能够识别
  • fig3.png
python -m jcvi.compara.catalog ortholog --no_strip_names --cscore=.99 halleri lyrata
#cscore=.99:这个选项能够有效地筛选出best hit,建议加上
#C分数截止点为0.99,有效过滤LAST命中,一直到包含互为最佳命中(RBH)。得到含有1:1基因对的直系同源区域
#1:1的基因对不是绝对的
#no_strip_names:是必加的参数
#dbtype:序列类型nucl还是prot
  • 主要产生四个结果文件以及一个PDF文件的点图

  • fig4.png

    1、.last文件是原始的 LAST 输出
    2、.last.filtered是过滤的 LAST 输出
    3、.anchors是种子同线性块(高质量)
    4、.lifted.anchors招募额外的锚来形成最终的同线性块

  • fig5.官方文档的配图.png

(五)测试同线性模式是否确实是1:1

python -m jcvi.compara.synteny depth --histogram halleri.lyrata.anchors
  • fig6.png

三、染色体水平的局部同线性图绘制

进行这一步需要准备三个文件

(一)生成一个.simple比文件更简洁的.anchors文件

python -m jcvi.compara.synteny screen --minspan=30 --simple \
halleri.lyrata.anchors \
halleri.lyrata.anchors.new
#minspan:是规定syntenic region的总长度(i.e. 基因个数)最低阈值,默认是0
#minsize:是规定syntenic region的anchors总个数(i.e.有共线性关系的基因对数目)最低阈值,默认是0

(二)seqids文件,记录需要展示的染色体

  • 这里由于我的示例数据中有一个不是染色体水平,因此我只选取了几条contigs

FJVB01000001.1,FJVB01000002.1,FJVB01000003.1,FJVB01000004.1,FJVB01000005.1,FJVB01000006.1,FJVB01000007.1,FJVB01000008.1,FJVB01000009.1,FJVB01000010.1,FJVB01000015.1
1,2,3,4,5,6,7,8

(三)layout文件,记录绘图的各种信息

  • y, xstart, xend:两个物种染色体在x轴或y轴的位置如y轴的0.6,0.4,如果是0.5,0.5两个物种的染色体就会重叠
  • rotation:旋转角度,0代表水平
y, xstart, xend, rotation, color, label, va,  bed
 .6,     .1,    .8,       0,      , halleri, top, halleri.bed
 .4,     .1,    .8,       0,      , lyrata, top, lyrata.bed
edges
e, 0, 1, halleri.lyrata.anchors.simple
  • 然后就是运行绘图模块
python -m jcvi.graphics.karyotype seqid.txt layout.txt
  • fig7.png

四、基因水平的局部同线性图

(一)计算基因级匹配的布局

python -m jcvi.compara.synteny mcscan \
halleri.bed \
halleri.lyrata.lifted.anchors \
--iter=1 -o halleri.lyrata.blocks
#iter=1表示我们正在提取与halleri每个区匹配的对应lyrata的唯一最佳block
#查看生成的.blocks文件,它的第一列为halleri的区块,第二为lyrata区块
#如果选项--iter设置为 2,则halleri每个区将匹配对应lyrata的2个block,依此类推

(二)合并bed文件

cat halleri.bed lyrata.bed > halleri_lyrata.bed

(三)选取可视化区域

head -150 halleri.lyrata.blocks > blocks

(四)准备layout2文件

# x,   y, rotation,   ha,     va,   color, ratio,  label
0.5, 0.6,        0, leftalign, center,       m,      1,    halleri
0.5, 0.4,      0, rightalign, center, #fc8d62,  1,    lyrata
# edges
e, 0, 1

(五)运行

python -m jcvi.graphics.synteny \
--glyphcolor=orientation \
--format=png --glyphstyle=box  --shadestyle=line \
--font=Helvetica --style=darkgrid --diverge=RdYlBu \
--figsize=8x7 \
blocks \
halleri_lyrata.bed \
layout2.txt
  • fig8.png
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 204,732评论 6 478
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 87,496评论 2 381
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 151,264评论 0 338
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,807评论 1 277
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,806评论 5 368
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,675评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,029评论 3 399
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,683评论 0 258
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 41,704评论 1 299
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,666评论 2 321
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,773评论 1 332
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,413评论 4 321
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,016评论 3 307
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,978评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,204评论 1 260
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 45,083评论 2 350
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,503评论 2 343