『比较基因组学』TBtools 30秒解决基因组共线性。

直接给我惊呆了，从点开开始到出结果一共32秒ಥ_ಥ，想当年光blastp就得跑大半天，更不要说编译MCScanX，各种奇葩问题...

先看下什么是共线性：

Syntenic* = a set of loci in two different species which is located on the same chromosome in each (not necessarily in the same order).
Collinear = a set of loci in two different species which are located on the same chromosome in each, and are conserved in the same order.
Synteny now has the same meaning as collinear, despite the different origins of the terms
Source: Genomics and comparative genomics

collinear

之前做共线性分析的时候用的是MCscanX，还有python版本的JCVI，第一步的blastp基本是要了老命了，当时为了加快blastp的速度，把query的pep.fa文件切成几十个批量跑....反正也是折腾了好久。
最近CJ发布了两个插件：OneStepMCScanX-Diamond和Quick Genome Dot Plot，从名字上就能看出来这个是怎么回事了，就是基于Diamond进行序列比对（Diamond比blast快很多倍！而且准确度也相差不是很多）其实相当于把序列比对和比对完的结果用MCScanX生成共线性文件，clt文件，包括gtf文件的合并等融合到一个步骤中了。当然现在的OneStep只是两个基因组间的比较。

具体方法

插件 | 超快！全基因组共线性分析
插件 | 超快！全基因组 Dot Plot
简洁 | 优雅地准备 比较基因组分析 文件

可能遇到的其他问题

去除scaffold

一般在共线性展示的时候，如果你的参考基因组组装到了染色体水平，scaffold或者contig等是要去掉的。如果有一点shell基础的话这个就很简单了
比如

  • 首先看下你gff文件的头文件，看看染色体编号形式和scaffold或者contig等的编号形式,这里我们看下JGI的雷蒙德氏棉的gff3
head -30 gene.Grai.JGI.gff3
##gff-version 3
##sequence-region   Chr01 6158 55678711
##sequence-region   Chr02 75 62758601
##sequence-region   Chr03 1 45758218
##sequence-region   Chr04 2059 62173553
##sequence-region   Chr05 6721 64137711
##sequence-region   Chr06 50440 51073468
##sequence-region   Chr07 1895 60974818
##sequence-region   Chr08 36764 57128522
##sequence-region   Chr09 8203 70705631
##sequence-region   Chr10 1619 62168479
##sequence-region   Chr11 3378 62659606
##sequence-region   Chr12 22670 35427486
##sequence-region   Chr13 518 58321163
##sequence-region   scaffold_102 2100 3605
##sequence-region   scaffold_107 2825 16443
##sequence-region   scaffold_108 3 13412
##sequence-region   scaffold_111 13319 15939
##sequence-region   scaffold_112 4529 5790
##sequence-region   scaffold_114 53 12661
##sequence-region   scaffold_115 3197 6044
##sequence-region   scaffold_116 268 1961
##sequence-region   scaffold_119 6 12264
##sequence-region   scaffold_122 5843 8827
##sequence-region   scaffold_127 2607 11788

再看下中棉所亚洲棉的gff3

head -30 gene.Garb.CRI.gff3 
##gff-version 3
##sequence-region   Chr01 1290 113026203
##sequence-region   Chr02 50819 99045064
##sequence-region   Chr03 942 135702384
##sequence-region   Chr04 55724 98593909
##sequence-region   Chr05 32236 97875121
##sequence-region   Chr06 5813 132229846
##sequence-region   Chr07 21338 97506195
##sequence-region   Chr08 8313 129422927
##sequence-region   Chr09 10279 85030627
##sequence-region   Chr10 64098 129481832
##sequence-region   Chr11 5334 124529361
##sequence-region   Chr12 6026 103049777
##sequence-region   Chr13 15788 123934920
##sequence-region   tig00000076 27806 38254
##sequence-region   tig00000208 71464 82652
##sequence-region   tig00000247 19374 20588
##sequence-region   tig00000248 17267 18598
##sequence-region   tig00000254 5842 784687
##sequence-region   tig00000266 95796 99132
##sequence-region   tig00000267 76876 81478
##sequence-region   tig00000312 258908 329614
##sequence-region   tig00000458 14417 14761
##sequence-region   tig00000498 12571 1498618
##sequence-region   tig00000545 29227 220599
##sequence-region   tig00000736 198220 210945
##sequence-region   tig00000744 182460 186937
##sequence-region   tig00000762 9152 12184
##sequence-region   tig00000807 27129 30536
##sequence-region   tig00000820 21055 21872

现在我们发现了再雷棉的gff3中是scaffold，而再亚洲棉的gff3中是tig，那么我们用grep -v（取反）将带有contig或者scaffold的行去掉

grep -v "tig" gene.Garb.CRI.gff3 > notigs.Garb.CRI.gff3
grep -v "scaff" gene.Grai.JGI.gff3 > notigs.Grai.CRI.gff3

检查一下：

head -30 notigs.Garb.CRI.gff3
##gff-version 3
##sequence-region   Chr01 1290 113026203
##sequence-region   Chr02 50819 99045064
##sequence-region   Chr03 942 135702384
##sequence-region   Chr04 55724 98593909
##sequence-region   Chr05 32236 97875121
##sequence-region   Chr06 5813 132229846
##sequence-region   Chr07 21338 97506195
##sequence-region   Chr08 8313 129422927
##sequence-region   Chr09 10279 85030627
##sequence-region   Chr10 64098 129481832
##sequence-region   Chr11 5334 124529361
##sequence-region   Chr12 6026 103049777
##sequence-region   Chr13 15788 123934920
Chr01   EVM gene    1290    21313   .   -   .   ID=Ga01G0001
###
Chr01   EVM mRNA    1290    21313   .   -   .   ID=Ga01G0001.1
Chr01   EVM CDS 1290    1426    .   -   2   ID=cds.Ga01G0001.1;Parent=Ga01G0001.1
Chr01   EVM CDS 21271   21313   .   -   0   ID=cds.Ga01G0001.1;Parent=Ga01G0001.1
###
Chr01   EVM gene    26354   29394   .   +   .   ID=Ga01G0002;Name=LPLAT1;description=Lysophospholipid acyltransferase 1 
###
Chr01   EVM mRNA    26354   29394   .   +   .   ID=Ga01G0002.1;Name=LPLAT1;description=Lysophospholipid acyltransferase 1 
Chr01   EVM CDS 26354   26396   .   +   0   ID=cds.Ga01G0002.1;Parent=Ga01G0002.1
Chr01   EVM CDS 26960   27184   .   +   2   ID=cds.Ga01G0002.1;Parent=Ga01G0002.1
Chr01   EVM CDS 27238   27699   .   +   2   ID=cds.Ga01G0002.1;Parent=Ga01G0002.1
Chr01   EVM CDS 28202   28287   .   +   2   ID=cds.Ga01G0002.1;Parent=Ga01G0002.1
Chr01   EVM CDS 28442   28505   .   +   0   ID=cds.Ga01G0002.1;Parent=Ga01G0002.1
Chr01   EVM CDS 28585   28703   .   +   2   ID=cds.Ga01G0002.1;Parent=Ga01G0002.1
Chr01   EVM CDS 28894   29027   .   +   0   ID=cds.Ga01G0002.1;Parent=Ga01G0002.1
###################################
head -30 notigs.Grai.CRI.gff3
##gff-version 3
##sequence-region   Chr01 6158 55678711
##sequence-region   Chr02 75 62758601
##sequence-region   Chr03 1 45758218
##sequence-region   Chr04 2059 62173553
##sequence-region   Chr05 6721 64137711
##sequence-region   Chr06 50440 51073468
##sequence-region   Chr07 1895 60974818
##sequence-region   Chr08 36764 57128522
##sequence-region   Chr09 8203 70705631
##sequence-region   Chr10 1619 62168479
##sequence-region   Chr11 3378 62659606
##sequence-region   Chr12 22670 35427486
##sequence-region   Chr13 518 58321163
Chr01   phytozomev10    gene    6158    6456    .   -   .   ID=Gorai.001G000100;Name=Trmt61a;description=tRNA (adenine(58)-N(1))-methyltransferase catalytic subunit TRMT61A 
Chr01   phytozomev10    mRNA    6158    6456    .   -   .   ID=Gorai.001G000100.1;Parent=Gorai.001G000100;Name=Gorai.001G000100.1;pacid=26824196;longest=1
Chr01   phytozomev10    CDS 6158    6415    .   -   0   ID=Gorai.001G000100.1.CDS.1;Parent=Gorai.001G000100.1;pacid=26824196
Chr01   phytozomev10    five_prime_UTR  6416    6456    .   -   .   ID=Gorai.001G000100.1.five_prime_UTR.1;Parent=Gorai.001G000100.1;pacid=26824196
###
Chr01   phytozomev10    gene    6987    7491    .   -   .   ID=Gorai.001G000200
Chr01   phytozomev10    mRNA    6987    7491    .   -   .   ID=Gorai.001G000200.1;Parent=Gorai.001G000200;Name=Gorai.001G000200.1;pacid=26820978;longest=1
Chr01   phytozomev10    three_prime_UTR 6987    7146    .   -   .   ID=Gorai.001G000200.1.three_prime_UTR.1;Parent=Gorai.001G000200.1;pacid=26820978
Chr01   phytozomev10    CDS 7147    7479    .   -   0   ID=Gorai.001G000200.1.CDS.1;Parent=Gorai.001G000200.1;pacid=26820978

现在就干净了。
如果不会shell这些操作怎么办呢，别着急，我们有另外一个工具Text Block Extract and filter(再TBtools 搜索栏直接搜索)
按照下图的操作即可

Text Block Extract and filter

发现染色体排序不符合预期

乱序的染色体排布

如果用过MCScanX 你一定会了解到ctl文件的作用。我们看下这个文件：

2500
600
Ga-Chr11,Ga-Chr10,Ga-Chr02,Ga-Chr13,Ga-Chr01,Ga-Chr12,Ga-Chr04,Ga-Chr03,Ga-Chr06,Ga-Chr05,Ga-Chr08,Ga-Chr07,Ga-Chr09
Gr-Chr11,Gr-Chr10,Gr-Chr02,Gr-Chr13,Gr-Chr01,Gr-Chr12,Gr-Chr04,Gr-Chr03,Gr-Chr06,Gr-Chr05,Gr-Chr08,Gr-Chr07,Gr-Chr09

看到了吧，罪魁祸首在这里，我们只需要手动的把这个配置文件的染色体顺序改正确就可以了

2500
600
Ga-Chr01,Ga-Chr02,Ga-Chr03,Ga-Chr04,Ga-Chr05,Ga-Chr06,Ga-Chr07,Ga-Chr08,Ga-Chr09,Ga-Chr10,Ga-Chr11,Ga-Chr12,Ga-Chr13
Gr-Chr01,Gr-Chr02,Gr-Chr03,Gr-Chr04,Gr-Chr05,Gr-Chr06,Gr-Chr07,Gr-Chr08,Gr-Chr09,Gr-Chr10,Gr-Chr11,Gr-Chr12,Gr-Chr13
改ctl之后

发散下思维 前面也可以不去掉scaffold，到时候直接这里面把没有到染色体水平的那些玩意给删了是不是也可以，我这里没有测试。

©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 203,456评论 5赞 477
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 85,370评论 2赞 381
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 150,337评论 0赞 337
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 54,583评论 1赞 273
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 63,596评论 5赞 365
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 48,572评论 1赞 281
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 37,936评论 3赞 395
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 36,595评论 0赞 258
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 40,850评论 1赞 297
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,601评论 2赞 321
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,685评论 1赞 329
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 33,371评论 4赞 318
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 38,951评论 3赞 307
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,934评论 0赞 19
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,167评论 1赞 259
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 43,636评论 2赞 349
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 42,411评论 2赞 342