「JCVI教程」使用JCVI进行基因组共线性分析(上)
本教程借鉴https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version).
我们先从http://plants.ensembl.org/index.html选择两个物种做分析, 这里选择的就是前两个物种,也就是拟南芥和水稻(得亏没有小麦和玉米)
我们下载它的GFF文件,cdna序列和蛋白序列
#Athaliana
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/cdna/Arabidopsis_thaliana.TAIR10.cdna.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.44.gff3.gz
#Osativa
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/oryza_sativa/cdna/Oryza_sativa.IRGSP-1.0.cdna.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/oryza_sativa/pep/Oryza_sativa.IRGSP-1.0.pep.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/oryza_sativa/Oryza_sativa.IRGSP-1.0.44.gff3.gz
保证要有6个文件以便下游分析
$ ls
Arabidopsis_thaliana.TAIR10.44.gff3.gz Arabidopsis_thaliana.TAIR10.pep.all.fa.gz Oryza_sativa.IRGSP-1.0.cdna.all.fa.gz
Arabidopsis_thaliana.TAIR10.cdna.all.fa.gz Oryza_sativa.IRGSP-1.0.44.gff3.gz Oryza_sativa.IRGSP-1.0.pep.all.fa.gz
我们分析只需要用到每个基因最长的转录本就行,之前我用的是自己写的脚本,但其实我发现jcvi其实可以做到这件事情
先将gff转成bed格式,
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_thaliana.TAIR10.44.gff3.gz > ath.bed
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Oryza_sativa.IRGSP-1.0.44.gff3.gz > osa.bed
然后将bed进行去重复
python -m jcvi.formats.bed uniq ath.bed
python -m jcvi.formats.bed uniq osa.bed
最后我们得到了ath.uniq.bed
和osa.uniq.bed
, 根据bed文件第4列就可以用于提取cds序列和蛋白序列。
# Athaliana
seqkit grep -f <(cut -f 4 ath.uniq.bed ) Arabidopsis_thaliana.TAIR10.cdna.all.fa.gz | seqkit seq -i > ath.cds
seqkit grep -f <(cut -f 4 ath.uniq.bed ) Arabidopsis_thaliana.TAIR10.pep.all.fa.gz | seqkit seq -i > ath.pep
# Osativa
seqkit grep -f <(cut -f 4 osa.uniq.bed ) Oryza_sativa.IRGSP-1.0.cdna.all.fa.gz | seqkit seq -i > osa.cds
seqkit grep -f <(cut -f 4 osa.uniq.bed ) Oryza_sativa.IRGSP-1.0.pep.all.fa.gz | seqkit seq -i > osa.pep
这里用到的seqkit建议学习,非常好用
下面使用python -m jcvi.compara.catalog ortholog
进行共线性分析,这是一个非常行云流水的过程(除非你报错)
新建一个文件夹,方便在报错的时候,把全部都给删了,
mkdir -p cds && cd cds
ln -s ../ath.cds ath.cds
ln -s ../ath.uniq.bed ath.bed
ln -s ../osa.cds osa.cds
ln -s ../osa.uniq.bed osa.bed
运行代码
python -m jcvi.compara.catalog ortholog --no_strip_names ath osa
输出结果如下
$ ls ath.osa.*
ath.osa.anchors ath.osa.last ath.osa.last.filtered ath.osa.lifted.anchors ath.osa.pdf
其中我们最感兴趣都是pdf结果,不出意外没啥共线性。
我们还可以用蛋白序列做共线性分析
# 在之前输出cds,pep都文件夹操作
mkdir -p pep && cd pep
ln -s ../ath.pep ath.pep
ln -s ../ath.uniq.bed ath.bed
ln -s ../osa.pep osa.pep
ln -s ../osa.uniq.bed osa.bed
运行代码
python -m jcvi.compara.catalog ortholog --dbtype prot --no_strip_names ath osa
我之前以为他不可以基于蛋白序列分析,幸亏有人提醒。
你会发现这是一个自动化分析流程,我们只是提供了4个文件,它就完成了一些列事情。它生成的文件里除了PDF外,其他还有啥用呢?
- ath.osa.last: 基于LAST的比对结果
- ath.osa.last.filtered: LAST的比对结果过滤串联重复和低分比对
- ath.osa.anchors: 高质量的共线性区块
- ath.osa.lifted.anchors:增加了额外的锚点,形成最终的共线性区块
anchors文件特别有用,之后会写一篇介绍如何利用他进行可视化,这里介绍它的格式。
###
AT1G28395.5 Os01t0238800-02 66
AT1G28440.1 Os01t0239700-02 1360
AT1G28480.1 Os01t0241400-01 136
AT1G28510.1 Os01t0242300-01 241
###
AT1G11100.3 Os01t0779400-01 943
AT1G11125.1 Os01t0779800-01 52
AT1G11160.2 Os01t0780400-02 535
AT1G11180.1 Os01t0780500-01 483
AT1G11330.2 Os01t0784700-00 742
AT1G11360.1 Os01t0783500-01 305
AT1G11540.2 Os01t0786800-01 422
AT1G11570.3 Os01t0788200-01 162
AT1G11580.2 Os01t0788400-01 550
AT1G11630.1 Os01t0793200-01 321
每个共线性区块以###
进行分隔, 第一列是检索的基因,第二列是被检索的基因,第三列则是两个序列的BLAST的bit score,值越大可靠性越高。
用水稻和拟南芥进行了比较之后,发现后面基本上也没啥可以分析了。因此下面基于「JCVI教程」如何基于物种的CDS的blast结果绘制点图(dotplot)得到的cds和bed文件进行分析。
之前已经得到了如下四个文件
ls ???.???
aly.bed aly.cds ath.bed ath.cds
所以我们只要运行
python -m jcvi.compara.catalog ortholog --no_strip_names aly ath
就得到了一个非常好看的点图
我们可以发现,都作为Arabidopsis属的两个物种,他们之间存在很高的同源性,并且同源区比例是1:1,
这其实和2011年的Nature Genetics上Alyrata的文章的结果是相似的,只不过他不是用点图进行展示
我们也可以用JCVI的画图模块实现这种效果,只不过还需要一点额外操作,创建如下三个文件
- seqids: 需要展现哪些序列
- layout: 不同物种的在图上的位置
- .simple: 从
.anchors
文件创建的更简化格式
第一步,创建.simple
文件
python -m jcvi.compara.synteny screen --minspan=30 --simple aly.ath.anchors aly.ath.anchors.new
第二步, 创建seqid
文件,非常简单,就是需要展示的scaffold或染色体的编号
scaffold_1,scaffold_2,scaffold_3,scaffold_4,scaffold_5,scaffold_6,scaffold_7,scaffold_8
Chr1,Chr2,Chr3,Chr4,Chr5
第二步,创建layout
文件,用于设置绘制的一些选项。
# y, xstart, xend, rotation, color, label, va, bed.6, .2, .8, 0, , Alyrata, top, aly.bed.4, .2, .8, 0, , Athaliana, top, ath.bed
# edges
e, 0, 1, aly.ath.anchors.simple
注意, #edges
下的每一行开头都不能有空格
最后运行下面的命令,会得到一个karyotype.pdf
python -m jcvi.graphics.karyotype seqids layout
如何让这个图垂直呢?(导入AI里就好了)
「JCVI教程」使用JCVI进行基因组共线性分析(上)相关推荐
- 「JCVI教程」如何基于物种的CDS的blast结果绘制点图(dotplot)
这是唐海宝老师GitHub上的JCVI工具的非官方说明书. 该工具集的功能非常多,但是教程资料目前看起来并不多,因此为了能让更多人用上那么好用的工具,我就一边探索,一边写教程 这一篇文章教大家如何利用 ...
- applewatch的健身目标如何修改「苹果教程」
我们来看看如何更改Apple Watch的「体能训练」和「健身记录」目标,调整卡路里目标或运动时间. Apple Watch的「健身记录」由「站立」.「锻炼」和「移动」目标组成,每天需要达到一定数值才 ...
- “十四五”要建设的「交通强国」,会让我们都坐上自动驾驶车么?
伊瓢 发自 凹非寺 量子位 报道 | 公众号 QbitAI 自动驾驶落地,到底啥时候才能广泛开展? 很多人都把它归于一个「政策」问题. 说道政策问题,前不久正式通过的十四五规划就提到了这一点. 倒是没 ...
- 「技术人生」第2篇:学会分析事物的本质
简介: 对于研发同学而言,探究事物的本质,是最基础最核心最先需要被掌握的技能,没有之一. 作者:贺科学 技术一号位不是岗位,更多的是技术人员在公司中做事的一种心态,这个系列的文章适合所有想要对日常工作 ...
- JCVI安装以及数据下载(用于共线性分析)
JCVI安装以及数据下载(MCScan) 安装 最简单的方法是通过PyPI安装它: pip install jcvi#或者安装开发版本 pip install git+git://github.com ...
- 得出来的视差图左边有黑色补上原图_「PS教程」Photoshop使用通道快速抠图的详细教程...
抠图的方法很多,通道抠图的话就有很多刚学习PS的小伙伴不知道了,今天我就来演示一次使用通道给人物抠图的全过程,教程还算蛮详细大家可以跟着教程一起来学习一遍. 先看一下原图: 一.打开通道面版 我们首先 ...
- nextpolish安装_「三代组装」使用Pilon对基因组进行polish
软件安装 官方提供了编译好的jar包,方便使用 wget https://github.com/broadinstitute/pilon/releases/download/v1.23/pilon-1 ...
- ai如何旋转画布_「AI教程」使用AI制作3D立体文字效果
今天macdown小编要通过AI制作一种3D立体字体,主要的知识点是混合工具的使用.Shift键.Alt键的灵活运用.3D旋转命令.投影效果的创建等,学会后可以应用在需要的设计中,比如海报设计,log ...
- dreamweaver排列顺序怎么用_「化妆教程」粉底液是怎么用的?用在哪个具体步骤顺序...
今日学而仿给大家带来一节粉底液"化妆教程",详情见下文:当然你也可以购买在线教程(如下):手把手教会你(文末) 粉底液怎么使用? 教您打粉底的技巧,涂抹粉底液其实不需要任何复杂的技 ...
- host ntrip 千寻rtk_「图文教程」千寻RTK连接千寻cors账号的操作步骤
一.SR2-(HC3手簿)连接千寻cors账号教程 打开SR2接收机,信号灯全亮,随后仪器进入自检状态,信号灯正常闪烁,开始搜星. 信号灯简介图 手簿连接主机 1.打开手簿app软件,选择[设备],点 ...
最新文章
- android debug database 源码解析
- 使用idea 在springboot添加本地jar包的方法 部署的时候本地jar没有包含的解决方法
- jupyter中中文显示不正常_jupyter 中文乱码设置编码格式 避免控制台输出的解决...
- 二、stm32f103+enc28j60
- http协议下:为什么请求与响应会做到准确误的对应。不会出现请求与响应的错乱...
- Redis系统性介绍
- 视频光端机常见故障问题及处理方法大全
- 输入 3 个正数,判断能否构成一个三角形
- [转载]c#委托事件简单例子
- Ruby数据结构-数组和哈希表
- C语言动态规划——背包问题详解
- 1247:河中跳房子
- 鸡啄米:C++编程入门系列之目录和总结
- 【江枫】lvm2与powerpath的Found duplicate PV问题
- 工作后,又想读个名校的计算机硕士,该怎么做?
- Mybatis的mapper代理开发方法
- 路边停车系统无线地磁车辆传感器
- Javascript 实现城市选择控件
- The ST Intranet updater server is unknown:mcucrossselector.codex.cro.st.com
- PS(2018)更换启动界面