天牛是鞘翅目多食亚目天牛科昆虫的总称，咀嚼式口器，有很长的触角，常常超过身体的长度，全世界约有超过35,000种。其中一些种类属于害虫，其幼虫生活于木材中，可能对树或建筑物造成危害。

光肩星天牛（Anoplophora glabripennis，以下简称AGLAB）原产中国，现已入侵美国，英国、澳大利亚等。初孵化幼虫先在树皮和木质部之间取食，25～30天以后开始蛀入木质部；并且向上方蛀食。虫道一般长90毫米，最长的达150毫米。幼虫蛀入木质部以后，还经常回到木质部的外边，取食边材和韧皮。详见Figure1.

本研究从头组装出一个AGLAB基因组，联合另外14个昆虫基因组（其中Agrilus planipennis和Onthophagus taurus属于i5k计划中新测的物种，鞘翅目物种5个，是系统发育树上绿色枝部分，下文简称甲虫。具体的物种拉丁名及缩写参见Figure 2a)，展开比较基因组分析，探讨了AGLAB植食性的适应性进化在基因组上的体现。

一、 基因组特征

1.     组装与注释

选取AGLAB雌性幼虫，提取DNA，测序组装，一共134X数据，最后装出来的基因组大小为710Mb，Scaffold和ContigN50分别是659kb和16.5kb。

利用MAKER流程预测出22035个蛋白编码基因，手工注释出1144个基因，合并去冗余后，得到22,253个蛋白编码基因，66个假基因。

2.     Horizontal gene transfer（HGT）鉴定

目前的研究发现，甲虫物种的基因组上有很多基因编码植物细胞壁降解酶（Plantcell wall degrading enzymes, PCWDEs）,这些酶的主要作用是降解木质植物组织中的纤维素、半纤维素或者果胶（植物细胞壁中的主要成分之一），释放的糖类，矿物质及其它的营养物种。目前的研究表明：一些甲虫的PCWDEs基因，是通过水平基因转移的方式从细菌或者真菌中获得，随后进化成多基因数目的家族。

除了糖苷水解酶（GH）家族基因，利用HGT分析流程还注释出8个HGT候选基因。其中的4个基因来自细菌，与沃尔巴克氏体（Wolbachia）的关系最近；还有2个基因与Wolbachia的序列显示出95%的相似度；最后两个基因与Wolbachia相比出现了indel，只有70-71%的相似度，推测这两个基因可能转移的比较早，后期发生了改变。

预测的这8个基因在随后的RNAseq中没有检测到表达，大概是因为RNAseq选取的是单一组织，而HGT相关的基因具有组织特异性，或者表达量很低，RNAseq检测不到。

相比而言，发生过HGT的GH基因转移时间更为古老，已经进化成功能性的多基因家族。

3.     不完整的甲基化机制

AGLAB基因组上似乎存在不完整的甲基化机制，只发现了维持性甲基化转移酶DNMT1，而缺少denovo 甲基化转移酶DNMT3。这和赤拟谷盗、果蝇比较类似。但是其他的甲虫基因组中存在着完全的甲基化机制。

二、 植物细胞壁降解酶（PCWDEs）研究

在AGLAB基因组中手工注释出86个GH家族基因。其中有57个GH1，18个GH28，11个外生/内生葡聚糖酶，包括GH9、GH45、GH48及GH5亚家族2（GH5-2）的成员。其中GH5是AGLAB基因组特有的(详见Figure3)。

接下来用RNAseq的手段研究食性依赖的GH基因家族的调控表达。做法是分别用人工食物和活的糖枫树（sugarmaple trees）作为AGLAB的食物，然后测定转录组数据(每组4个三龄幼虫个体。具体差异基因表达的趋势，参考Figure4)。

在糖枫树喂食组中，GH5和GH45基因相比人工食物组的基因表达量在2倍以上，这两个基因可能参与把纤维素转化成容易消化的纤维低聚糖。

超过30个GH1基因在糖枫树组高表达，这些基因推测是Beta-葡萄糖苷酶，转化植物细胞壁中的纤维二糖或者寡聚糖成单糖。GH1家族还有部分的Beta-糖苷酸有解毒的功能。

有12个GH28家族的基因在糖枫树组高表达，GH28在AGLAB近缘物种中的同源基因的功能是聚半乳糖醛酸酶。这样看来，植物细胞壁中的果胶组分可能是AGLAB幼虫主要的糖来源。GH35在糖枫树组也被诱导，同源比对发现和Beta-半乳糖苷酶高度相似，可能在处理植物细胞基质中的Beta-1,4半乳糖苷多聚体中起重要作用。GH30同样被诱导，主要功能可能是消化植物次级细胞壁。

以上都是计算机的预测，没有被真正验证。接下来在体外检查GH酶家族的底物特异性和功能。18个GH28家族基因中有15个在平板鉴定中表现出功能特征。大多数GH28蛋白被至少一种同聚半乳糖醛酸聚合物激活。

另外6个GH5基因，2个GH45基因，1个GH9基因也在体外得到功能验证。通过体外实验分析表明，AGLAB基因组中至少有三个纤维素和半纤维素酶（GH5-2，GH9，GH45），还有1个多聚半乳糖醛酸酶（GH28）.

通过数据挖掘，只在三个甲虫物种中发现了GH28、GH45、GH5-2，其他12个昆虫基因组中缺失。其中，GH28在三个甲虫中都检测到了。而GH45只在两个甲虫物种中检测到（AGLAB，DPOND）。GH5-2只在AGLAB中发现，是叶甲总科物种特有的。这些基因没有在白蚁基因组序列中找到，但在其肠道微生物中发现了相关的基因。目前的研究表明GH28家族基因通过HGT的方式获得，并且AGLAB和DPOND是独立获得的。

另外，GH1家族的基因编码两种蛋白，一种是有消化功能，另外一种没有。AGLAB基因组中23个GH1基因和赤拟谷盗（TCAST）中的葡萄糖硫苷酶（MYR）有44%左右的同源性。对于一些昆虫来说，MYR的作用是协同性警告或者竞争性的信息素。在某些植物中常会释放出硫代葡萄糖苷，恰好是MYR的底物，因此昆虫可以解除毒性。其中还有几个基因是生氰的Beta-糖苷酶。在某些植物中有毒性的生氰糖苷是植物的防御系统，用来对抗生物胁迫。和MYR的机制类似，生氰的Beta-糖苷酶也可以解除生氰糖苷的威胁。

AGLAB肠道中的微生物对其营养物质的生物合成和循环有重要的作用，能帮助AGLAB在营养贫瘠的环境中繁殖、发育。AGLAB编码的基因中，有8个基因有血蓝蛋白的结构域。其中3个高表达的基因在以木质植物作为食物来源的生物中存在，与肠道微生物一起协作，在多个生物反应中起氧化降解作用。

三、 对植物异种化感物的解毒机制研究

1.     CYP450家族

目前植食性动物解毒相关的基因，研究最多的就是CYP450家族基因。通过手工注释，得到106个CYP450基因和19个假基因。CYP6代谢毒性化合物，CYP4参与角质层的碳氢化合物的合成，与杀虫剂的抗性相关。糖枫树喂养的幼虫组中，25个CYP450家族基因在AGLAB的肠道中被诱导。而这些诱导的基因中，只有2个在赤拟谷盗中存在。CYP基因家族的扩张，是AGLAB适应和应对植物防御系统的结果。

2.     UGT家族

UGT（UDP-葡糖醛酸基转移酶）家族的基因参与毒性物质的解毒，外源性物种的降解及内源性物质的调控。手工注释出65个UGTs基因，包括7个假基因。UGTs在甲虫物种基因组中的扩张和其广泛的宿主植物有关，其中大部分基因聚在一起形成cluster，7个cluster就包含了50个基因，其中大部分基因通过串联复制的机制形成多样性，增加了结合次级代谢物底物的范围。

最大的UGT家族是UGT352，对AGLAB物种来说是特异的，包括21个基因，其中的14个UGTs基因同向定位在一个scaffold上（参考Figure5）。AGLAB特异扩张的UGT家族是UGT321，有7个基因，这些扩张能使AGLAB适应广泛的宿主植物的防御系统。在糖枫树喂养的幼虫组中有4个UGT基因上调，包括2个UGT321和1个UGT352基因。

虽然在RNAseq数据中只发现了一部分UGT基因，这主要是因为研究只使用了1个宿主植物喂养有关，不同的宿主，可能会激发不同的UGT基因表达。

3.     酯酶类（Esterases）

在AGLAB基因组中发现了更多的esterases酶，主要是Type-B羧酸酯酶（COesterase）的扩张。COesterase对外源性物质的代谢和植物细胞次级代谢物中木质素-半纤维素的酯键断裂很重要。一共鉴定出107个COesterase，是其它昆虫的两倍多，大部分是形成cluster，只有25%的基因是singletons。

COesterase在糖枫树幼虫组被高诱导，而COesterase基因中被高诱导的基因，大多是AGLAB特异的，在基因组上形成了串联重复。COesterase在AGLAB上潜在的功能是木质植物消化或者植物异种化感素的解毒作用。

消化性蛋白酶作用是清除来自植物细胞壁蛋白或者肠道内共生体的氮物质，用来对抗植物产生的蛋白酶抑制物。蛋白酶的直系同源基因也在AGLAB基因组上发生扩张，编码胰蛋白酶基因，可以产生分泌性的丝氨酸蛋白酶。但是呢，大部分蛋白酶对每种甲虫都是特异的，这可能意味着它们的进化发生在不同甲虫物种形成之后。

四、 AGLAB感官相关的基因研究

1.     化学感官基因

化学感官包括嗅觉和味觉。在AGLAB基因中手工注释出52个OBP（odorantbinding protein），大部分是minus-C亚家族，只有一个是plus-C亚家族，和TSACT和DPOND类似。表明OBP自从190百万年前鞘翅目扁甲系起源开始，朝plus-C 方向进化。另外鉴定出131个OR（odorant receptor）。在AGLAB中共鉴定出234个GRs（gustatoryreceptors），其中3个CO₂ 受体，10个糖受体，3个果糖受体，另有127个GRs，通过可变剪接形成218个受体基因，主要是苦味感知受体和信息素感知受体。类似于TSACT，AGLAB的OR、GR相关的基因比果蝇多。

2.     视觉基因

AGLAB有一个单独的长波长敏感的视蛋白（opsin），一个单独的紫外线敏感的视蛋白。有RH7视蛋白，但是缺C-opsin.

五、 小结

1.     光肩星天牛是一种植食性昆虫，在其基因组中进化出一套完整的酶促指令，可以降解植物宿主中的大部分多聚糖，比如纤维素、木葡聚糖、木聚糖、果胶等。解毒基因和消化蛋白酶基因的多样性组合，提供了代谢的可塑性，来克服不同宿主植物的防御机制。

2.     PCWDEs基因数目在AGLAB基因组上也发生了扩张。其中三个重要的家族GH5-2、GH9和GH45推测是纤维素酶，可以处理纤维素和半纤维素，AGLAB通过利用肠道微生物或者基因组编码的酶处理木质素。

3.     CYP450家族基因、UGTs家族基因及COesterases都在AGLAB基因组中通过串联复制发生了扩张。这些基因家族的功能主要是克服植物宿主释放出来的毒素，降解外源性物质，对抗植物释放的蛋白酶抑制剂。

4.     AGLAB通过HGT，从真菌或者细菌中获得GH基因，GH基因通过基因拷贝数扩增、功能趋异，导致其代谢指令的增加、扩张和强化，使AGLAB适应广泛的植物宿主。

无论是植食性物种还是肉食性物种，在其特定的食性环境下，受到选择，进化出一套完善的对抗指令，比如植食性物种需要对抗植物的防御机制，而肉食动物需要特定的代谢系统，来消化肉食中特定的营养物种，比如高蛋白和高胆固醇以及较低的糖类物质。所以从某种程度上讲，食性塑造了物种的基因，使食物链的上下游形成了一套行之有效的进攻与防御指令。

From：McKenna, D. D. et al. Genome of the Asian longhorned beetle (Anoplophoraglabripennis), a globally significant invasive species, reveals key functionaland evolutionary innovations at the beetle–plant interface. Genome biology 17,227, doi:10.1186/s13059-016-1088-8 (2016).

 2016/11/13

这里说的猫不是家猫，是指猫科动物，包括家猫、狮子、老虎、豹子等。

近日有研究团队，公布了远东豹基因组，并且结合已经公布基因组的猫科动物及其他哺乳动物，在比较基因组的水平，就哺乳动物的食性展开了适应性进化的大讨论。

一、 数据收集

首先是测定了一个远东豹（Amur leopard）的基因组。用Illumina 小片段文库+大片段文库，测序数据高达310x。估算基因组大小是2.45Gb，最后组装出来的版本大小为2.58Gb。Scaffold N50 高达21.7Mb, 而Contig N50 仅有21.0kb。最后用两个野生远东豹个体的TSLR (TruSeq synthetic long reads)做Gap区域的校正。然后预测出19,043个编码蛋白的基因。重复序列注释发现，转座子比例高达39.04%。

另外重测序了两个野生的远东豹（不同于先前做TruSeq的两个个体）及一个野生远东豹猫。

最后加上目前已测序的17个哺乳动物，共18个物种，根据食性分为三个大类。肉食动物组包括Leopard、Tiger、Lion、Cat、Cheetah五个猫科动物，加上Polar bear 、Killer whale、Tasmanian devil共八个物种；杂食物种包括：Human、Mouse 、Dog 、Pig、Opossum五个；草食性五个：Elephant、Rabbit、Horse、Cow、Giant panda。利用这些物种的基因组序列数据，做接下来的比较分析。

二、 进化分析

1.     基因家族扩张收缩分析

对18个物种进行基因家族鉴定及扩张收缩分析。Figure 1A 展示的是五种猫科动物共有家族及特有家族数目。Figure 1B展示的基因家族扩张收缩的数目。远东豹一枝相对最近共同祖先有188个扩张家族，313个收缩家族。猫科物种相对最近共同祖先有52个扩张基因家族，567个收缩家族。猫科物种扩张的基因家族，做GO功能富集分析，发现其主要功能集中在肌肉收缩及运动活性等方面。而收缩家族的基因主要功能集中在淀粉和蔗糖的代谢通路。在肉食动物中UDP-葡糖醛酸基转移酶（UGT）1和2家族基因明显收缩，这些基因在体内解毒功能和稳态平衡上起着重要的作用。对于食草动物而言，在食用植物或者破坏植物结构时，遭受植物毒素的攻击，UGT1和UGT2家族的解毒功效就很关键了。肉食动物对这些食物需求很少，导致这些基因家族慢慢在肉食物种中呈现收缩状态。只有猫科物种中UGT3家族呈现扩张状态，UGT3的主要功能是参与N-乙酰氨基葡糖和葡萄糖的接合。UGT8A1参与神经酰胺和胆酸与半乳糖的接合，在检查的18个哺乳动物中基因拷贝数很保守（Figure2A）。淀粉酶基因家族（AMY1和AMY2）主要作为淀粉和糖原代谢的催化剂，在肉食动物中收缩。这为家猫的唾液淀粉酶低水平含量提供了遗传上的解释。

家猫缺少合成足够量的维生素A和花生四烯酸的能力，所以这些物质必须从膳食中获取。有意思的是，参与视黄醇（维生素A）、亚油酸、花生四烯酸的分解代谢的CYP450家族基因在肉食动物枝普遍收缩。视黄酸是从视黄醇转化而来，对牙齿再矿化和骨生长是必需的。花生四烯酸的主要功能是促进身体运动后骨骼肌肉组织的修复和生长。肉食动物中CYP450基因家族收缩，则能使其体内保持较高的视黄醇和花生四烯酸，这样才能使肉食动物保持强壮的肌肉、骨头和牙齿，便于捕猎。

虽然肉食动物的能量和营养主要是由动物组织驱动的，但是它们也需要调控机制，对组织保持充足的葡萄糖供应。葡糖激酶（GCK）在生物体内重要是用来调节葡萄糖的吸收和存储，被称为膳食葡萄糖的传感器。GCKR基因编码的蛋白GKRP，可以调控GCK。GCK和GKRP近年来被作为糖尿病治疗的靶点。先前的研究表明，GCKR基因在许多哺乳动物（包括家猫）中发生了移码突变，变成了假基因。经过分析目前的基因组，发现Leopard、Tiger、Lion、Cheetah、Snow leopard、Leopard cat中也发生了移码突变导致的假基因化。有意思的是，Kill whale和Domestic ferret的GCKR基因也变成假基因了，不过是通过提前终和/或移码突变。而Polar bear和Tasmaniandevil则包含完整的GCKR基因。根据这些结果推断，肉食动物没有必要移除代谢循环中的过量葡萄糖，因为他们的食物中包含大量的蛋白，而糖类成分则很少。在非肉食动物中，Cow和Opossum的GCKR基因也变成假基因。以Cow为例，可能是因为反刍动物主要利用前胃发酵产生的挥发性脂肪酸作为主要的能量来源，没有必要移除过量的葡萄糖。因此，肉食动物中GCKR基因进化性丢失，使其参与的葡萄糖敏感的代谢通路发生适应性进化。这些研究有助于我们更好地理解糖尿病中典型的葡萄糖代谢异常的机理。

2.     正选择基因（PSG）分析

为了更好的理解物种对食性环境的适应性选择，研究者利用PAML软件的branch model和branch-sitemodel筛选PSG。在Leopard基因组中筛选出586个PSG，主要功能富集在细胞增殖和微管代谢等方面。猫科物种共有的PSG有228个，主要功能富集在多糖结合、脂肪结合和免疫应答。猫科物种作为一个超级肉食者，脂肪结合相关的基因受到选择，可能有助于其保持脂肪和胆固醇的调控稳态。

接下来对三种食性组分别做正选择分析。发现肉食动物共有的PSG主要富集功能是运动轴突引导。其中的CXCL12基因主要是引导神经元迁移和轴突生长。DMP1和PTN基因主要是在骨发育和修复中扮演重要角色。有所不同的是，杂食动物和草食动物的PSG主要富集在免疫相关的功能分类。

3.     物种特异的氨基酸改变（AACs）分析

对15个猫科动物（three leopards, three lions, a snowleopard, three tigers, two leopard cats, a cheetah, and two cats）及其他13个哺乳动物做AACs分析。分析发现，猫科动物能改变的AACs基因有1509个，主要富集在DNA损伤刺激和细胞压力应答。有意思的是，对猫科动物来说，三个参与蛋白消化和吸收的通路的基因（MEP1A、ACE2、PRCP）发生功能改变性的AACs。动物肉里面的ROS（heme-related reactive oxygen species）可能会引起DNA损伤，阻碍正常的细胞增殖。因此，肉食动物中和DNA损伤和修复相关的基因发生功能改变，有助于降低饮食带来的DNA损伤。这些研究都可以启迪人的饮食和人的健康研究。

对肉食动物和草食动物还做了趋同的AACs分析（Figure2B）。只在肉食动物中发现了一个趋同的AAC基因embigin（EMB）。在食草动物中没有发现AAC趋同的基因。这样的事实比较符合之前的观点：适应性的分子趋同能关联到表型趋同上的案例真的是很少的。目前的已有的研究表明EMB主要参与运动神经元的分支和神经肌肉连接的形成。在肉食动物特异的功能改变基因中，TMOD4和SYNC基因已知的作用与肌肉收缩有关，STAR基因与类固醇激素的合成有关。

三、 保守区域分析

高度保守区域（HCR）的分析以科为单位展开，分别是猫科(cat, tiger, lion,cheetah, leopard, snow leopard, and leopard cat, divergence time: ~15.9 millionyears ago [MYA], carnivores)、人科（human, chimpanzee,bonobo, gorilla, and orangutan, ~15.8 MYA, omnivores）和牛科（cow, goat, sheep, water buffalo, and yak, ~26 MYA, herbivores）。在猫科中找到1.13Gb的HCR，在人科中找到0.93Gb的HCR，在牛科中找到088Gb的HCR。三个科HCR共有的基因数为4,342个（Figure3），主要富集在细胞周期、癌症通路、蛋白酶体，Hedgehog信号通路等方面。猫科HCR区域的基因主要富集在神经系统相关的通路上，表明猫科物种对快速条件反射的进化适应。猫科HCR也有基因富集到糖类的生物合成，这或许和猫科的饮食模式相关，比如以肉为主的膳食通常导致低糖的摄入，需要自身再合成一些糖类营养物质。另外，牛科特有的基因主要富集在气味的感应，这表明食草物种对植物毒性建立了适应性的防御机制。

四、 猫科物种的遗传多样性和种群历史

最后对猫科物种的遗传多样性和种群历史做了调研。

遗传多样性估算，基于杂合SNV比率，猫科物种有最低的遗传多样性（平均0.00094）。遗传距离的计算是基于纯合SNV数目相对基因组的比例，猫科物种也显示了最低的遗传距离（平均0.00102）。这表明极端的膳食特性，对这些物种施加了强烈的、相似的选择压力。Leopard的遗传多样性特别低，是由于隔离和近亲交配。较小的猫科物种相对的遗传多样性较高。

像猫科这样的物种，处于食物链的金子塔尖。相对底层的物种，猫科有较小的群体数目，面临着灭绝的风险。利用PSMC model估算猫科物种的有效群体大小和种群历史（Figure4）。Leopard cat在10到2百万年前有一个群体的增加，而其它大型猫科物种的有效群体一直处于减少趋势。Leopard在200百万年前到90万年前经历了一次严重的遗传瓶颈，其他的大型猫科物种则没有。在最近的3万年，组装的Leopard基因组显示出群体的扩张，而野生的Leopard则没有。Leopard相对大的有效群体很可能反映出近期远东豹和华北豹群体之间有交配发生，并非真正的群体数量扩张。

五、 小结

这篇文章主要是组装了一个远东豹的基因组，这是目前猫科物种中组装最好的基因组。然后研究者从猫科物种及其他13个哺乳动物的基因组入手，把物种分为肉食、杂食、草食三个食性特性，展开比较分析。研究揭示了肉食动物，尤其是猫科动物对饮食的适应性进化。这些研究可以启示人的糖尿病研究、人的健康膳食研究等。

当然，这篇研究的肉食物种范围仅限于以脊椎动物作为食物的物种。而肉食动物还包括食昆虫、食非脊椎动物及吸血等食性。

特有的食性环境，使肉食者进化成身形敏捷、肌肉发达、膳食高度特化的物种。也只有这样的特征才能保证其捕猎到心仪的食物。这些生活方式也带来一些负面效应，使猫科物种，尤其是大型猫科物种，呈现出较低的遗传多样性和较小的种群数量，容易走向易危甚至灭绝。

这篇研究入选了Nature杂志的一周Highlight，标题为“Cat DNA shaped by diet”。全文如下：

Carnivores haveexperienced stronger natural selection than plant-eating animals, perhapsbecause of their limited diet.

Joo-Hong Yeo at the NationalInstitute of Biological Resources in the Republic of Korea and his colleaguescompared the genomes of 18 mammals, including carnivores such as leopards,omnivores such as humans, and herbivores such as giant pandas. They found thatcarnivore genomes have lost many genes for carbohydrate digestion, and sharechanges to genes involved in muscle strength and agility, making them goodhunters. By contrast, omnivore and herbivore genomes shared fewer adaptationswith others in their groups, suggesting that their diets have imposed weakerselection compared with that of carnivores.

Many big cats showedrecent losses in genetic diversity, suggesting that population declines may belinked to their strict diet. Genome Biol. 17, 211 (2016)

From：Kim, S. et al. Comparison of carnivore,omnivore, and herbivore mammalian genomes with a new leopard assembly. Genome biology 17, 211, doi:10.1186/s13059-016-1071-4 (2016). 

                                                                                                        2016/11/12

了解一个物种或者细胞的基因组的倍性，在进化、群体、基因组研究上有重要的意义。比如组装多倍体物种的基因组，需要了解它的倍性，评估组装的难度。但是，怎样检测或者推算基因组的倍性呢？通常在实验室用流式细胞实验就可以了。但是如果你恰好做了基因组survey，有NGS数据在手头，那你完全可以用好这些数据，做一下倍性分析。

目前已经有一些流程利用NGS数据推算基因组的倍性，比如：AbsCN-seq、CLImAT或者ConPADE。但是这些流程或者软件目前都有一些明显的短板。AbsCN-seq除了mapping数据外，还需要全外显子的数据。CLImAT基于MATLAB，不能免费使用。ConPADE是专门为高等多倍体植物开发的，对mapping的质量很敏感，容易引起估算偏差。

最近有个团队开发了一款流程ploidyNGS，专门针对二代测序数据，估算基因组的倍性。PloidyNGS的原理很简单，就是统计每个等位基因点上的reads支持比例。比如，对于单倍体来说，除了测序错误的reads外，所有的reads都支持一个allele。对于二倍体物种来说，取决于杂合度。高度纯合的话，所有的正确测序reads都会支持一个allele，杂合度高的物种，大约一半reads支持一个allele，另一半支持另外一个allele。对于三倍体的物种来说，如果是多等位基因（比如ABC）的话，每个allele的支持率为1/3，如果是双等位基因（比如AAB）的话，2/3的reads支持一个allele，剩下的1/3支持另外一个allele。以此类推，详见Table1.

知道了原理，接下来说说具体怎么做。

先是把reads mapping到参考基因组（做survey的时候，建议做一个低深度的基因组组装版本），过滤掉单allele位点及过高频率的allele（比如95%），然后按照每个allele的reads支持率排序，从最低到最高排序，分别被标记为fourth、third、second、first。然后用ggplot作直方图。

利用模拟数据和测试数据，研究者对真菌的基因组做了倍性评估。如Figure1所示，横坐标是allele的频率，纵坐标是多态性位点数目。A图是单倍体基因组统计图，有两个峰值点，分别在5%和95%左右，前者是测序错误造成的，而后者是单倍体基因组的最高频率单态峰。B图有四个峰值，分别在5%、50%、50%及95%。5%是测序错误，两个50%，是杂合多态性峰，95%则是最高频率的单态峰。这些结果和Table1的理想比例值较为一致。

最后再附上三倍体和四倍体的图像，供参考（Figure2，Figure3），看看每个峰值是不是和Table1给出的理想比例值一致呢？

PloidyNGS流程用python和R串写而成。

下载路径点这里：https://github.com/diriano/ploidyNGS

参考文献：

Renato et.al. ploidyNGS: Visually exploring ploidy with Next Generation Sequencing data.

2016/11/09

当我们做一个物种的基因组de novo组装时，手里其实已经有一些分子数据了，比如该物种曾经做过转录组de novo 研究，或者在NCBI数据库上能下载到EST序列或者已经验证过的转录本数据、蛋白数据，甚至还有大量的近缘物种的蛋白序列供使用。这个时候，我们该思考的事情是：组装时，这些数据能否派上用场？能不能辅助目前的组装工作呢？

现在给大家安利两款辅助de novo组装的软件：L_RNA_scaffolder和PEP_scaffolder，分别是借助转录组数据和蛋白数据辅助scaffold构建。两个软件是一个团队开发的，并且基本原理比较相似。

一、     利用转录本辅助基因组组装：L_RNA_scaffolder

1.  背景

我们利用短reads或者长reads序列，组装出一个初级的contig版本的基因组，这个时候，需要借助大片段文库来构建scaffold。利用转录本辅助scaffold 构建，主要是基于这样的考虑：

（1）当一段转录本不能完全的mapping到一段基因组片段时，很有可能作为引导序列辅助两段基因组序列的定位；

（2）包含大内含子片段转录本长度比例能占到基因座的40-80%以上，大内含子相关的转录本跨度大，可以形成类似大片段文库的效果，有利于辅助组装；

（3）已有证据表明，基因组是广泛被转录的；所以转录本数据充足的话，非常有利组装；

（4）目前转录本研究比较火热，可以容易地获取转录本数据。

2.  原理

Figure 1 为我们展示了该软件的组装原理。

（1）   选取引导转录本。哪些转录本是引导转录本呢？主要是那些转录本的一部分区域比对到一条基因组片段上，另一端部分区域能比对到另一条基因组片段上。这样的转录本可以作为引导序列，把两个或者多个基因组片段连接起来；

（2）   转录本和基因组片段的比对区域，通过聚类，可以划分到不同的blocks里面；

（3）   每个blocks里面选取最长的比对区域可以作为代表；

（4）   通过代表性的blocks和转录本之间的关系，可以把这些blocks的顺序给确定出来；

（5）   每个代表性的blocks又和基因组片段是关联的，根据（2）的结果，可以定出基因组片段的顺序。注意：如Figure1 中图5和图6所示。c和b之间的距离，可能是一个intron的距离。但是也有可能这个距离过大，超过一般intron长度，比如a和d之间的距离，在后续组装中，这种情况被过滤掉。因此，这一步只保留了c、b、a的结果。初步定位出C-B-A的path。当然，这只是一个转录本的支持情况，需要更多的转录本作为证据支持；

（6）   在实际操作中，一个基因组片段可能会作为多个connection的起始或者终止片段。最终选取转录本支持最多的connection。比如Figure1中图7的情况，CB和CE中，CB支持证据多，最后保留CB connection。BA和BF中，BA支持率高，选取BA connection；

（7）   最后，根据CB和BA的connection，构建出C-B-A的path。这样就把三个基因组片段连接起来。通过类似的原理，可以获得更多的path。把基因组的多个contig组装成一条较长的scaffold序列。

 3.     软件的关键参数

（1）  三个重要的参数

MIL(maximal intron length)：最大内含子长度。如果两个blocks之间的长度大于MIL设定的参数，这样的connection需要被过滤掉。MIL对组装的影响见Figure2.

MLC（minimal length coverage）：最小的长度覆盖。Length coverage = alignedlength / total transcript length. MLC值以下的比对，被认为是非完全覆盖的比对，这样的转录本有可能辅助组装，作为引导序列。MLC对组装的影响见Figure3.

MPI（minimal percent identity）：最小的一致性比率。这个最容易理解。就是能比对上的区域中的一致性比对率。显然MPI越高，能比对上转录本越少，connection和path也会相应的减小。但是MPI越小，比对的准确度也会降低。所以需要做出一个权衡来。详见Figure4.

   文章中组装斑马鱼基因组时，给出的这三个参数值的经验值。MIL, MLC和 MPI 分别被设为 100 kb, 0.95 和 0.9.

（2）  Gap估算

对于mate-pair 文库，插入片段大小通常遵循一个正态分布，选取平均的插入片段大小，可以作为gap的估算值。但是在转录本辅助组装中就不一样了。因为内含子的大小并不是一个正态分布，而是一个偏态分布，较小的intron比例相对比较高。所以通常用中位数来代表一个物种的内含子长度。具体来说，就是利用那些能完全比对到基因组片段上的转录本数据，计算出intron的中位数。仍以Figure1中图5为例。c-b 之间可能的intron长度为L（intron）≧ length（C）- end（c） + start（b）；end（c）是Block c在contig C上的终止位置，start（b）是Block b 在contig B上的起始位置。如果L(intron)小于intron中位数，插入的N的数量由L(intron)和中位值的绝对差值来确定，否则插入的N数目设置为100.

4.  面临的挑战

由于转录本经常发生可变剪接，基因融合等事件，所以这种方法也面临一些问题。

（1）  顺式可变剪接的影响

如Figure5所示，可变剪接a和b能比对到相同的contig上，这样组装出来的结果不受可变剪接影响。

如Figure6所示，主要的可变剪接a、b、c包含了一个基因座的所有外显子，所以能完整的重构引导RNA和contig顺序。并不会受到d的影响。

如Figure7所示，b和c转录本在引导上占主导地位，但是它们不包含所有的exon，这样组装出来的scaffold序列，会出现relocation的错误。

（2）  反式可变剪接及基因融合事件的影响                  

如Figure 8所示，a和b是chr1上的host 基因，而c是嵌合RNA，但是a和b起主要引导。这种情况会以host基因的a、b转录本作为引导，能得到正确的组装结果。

如Figure9所示，和Figure 8刚好相反，嵌合RNA占主导地位，最后的组装结果会出现translocation之类的错误。

（3）  除了上文提到的relocation和translocation之外，还会出现inversion之类的错误。

5.  优势及案例

（1）成功地提高scaffold N50，并有较高的准确度

文章中的案例显示，能把斑马鱼的N50从140kb增加到177kb；准确度在93%以上；把人的N50从144kb增加到432kb.

在尖吻鲈基因组文章中，scaffoldN50 从1Mb提高到1.2Mb，提高10%左右。

（2）有效提高转录本在基因组上的覆盖度

通过转录本辅助组装之后，有一个显著的好处是，能有效提高转录本在基因组上的覆盖度。以人的hg19为例，能使contig版的转录本低覆盖比例，提高到scaffold版的转录本高覆盖比例。转录本的高覆盖，有助于后续的基因集注释，及转录本定量分析。

（3）对高杂合基因组有显著的提高作用

以珍珠贝为例，基因组的2/3部分是高度多态性的。利用L_RNA_scaffolder 能使N50从14.5kb上升到62.8kb，提高了近5倍。并且使转录本比对率从76.7%提高到82.8%，对后续的基因预测有很大的帮助。

（4）组装的准确性和连续性的平衡

通过和几款主流的scaffold构建软件比较，发现其在组装准确性和连续性上达到了平衡，并且处于上游位置。

（5）没有显著的数据量饱和

L_RNA_scaffolder可以接纳更多的转录本数据(包括经典的转录本数据、EST数据、Iso-Seq全长转录组数据等)，并且没有明显的饱和效应。而且较长intron的转录本更有助于组装效果的提升。

二、     利用蛋白数据辅助scaffold定位：PEP_scaffolder

和L_RNA_scaffolder的原理比较类似，不同之处是：最开始的比对blocks确定是用蛋白序列比对到基因组序列上的。

PEP_scaffolder和目前已有的两种蛋白辅助组装软件（SWiPS 和ESPRIT）相比，具有精度高、耗时短的优势。

实践证明，越多的蛋白数据，能获得越好的组装效果。而且支持近缘物种的蛋白序列，因为相比转录本数据，蛋白更保守。更多的数据能覆盖到更多的基因组片段，有助于提高N50指标等。

另外，Swiss-Prot数据比TrEMBL数据库的蛋白数据更可靠，选择需谨慎。

考虑到蛋白数据选择性比较自由，PEP_scaffolder对组装的提高在10%-20%，甚至更高。

实践中，参数推荐为：MPI大于0.9；MIL 不大于150kb；MLC 大于0.9.

参考文献

1. Xue, W. et al. L_RNA_scaffolder: scaffoldinggenomes with transcripts. BMC genomics14, 1-14,doi:10.1186/1471-2164-14-604 (2013).

2. Zhu, B.-H. et al. PEP_scaffolder: using(homologous) proteins to scaffold genomes. Bioinformatics,doi:10.1093/bioinformatics/btw378 (2016).

3. Vij, S. et al. Chromosomal-Level Assembly ofthe Asian Seabass Genome Using Long Sequence Reads and Multi-layeredScaffolding. PLoS genetics 12, e1005954,doi:10.1371/journal.pgen.1005954 (2016).                                             

2016/9/13

先来看一段关于系统发生树的定义吧：系统发生树（英文：phylogenetic tree或evolutionary tree）是表明被认为具有共同祖先的各物种相互间演化关系的树，又被译作系统发育树、系统演化树、系统进化树、种系发生树、演化树、进化树、系统树。它用来表示系统发生研究的结果，用它描述物种之间的进化关系。

今天的议题是，如何绘出高颜值的进化树。

当我们打开论文时，常常看到别人做的进化树图是这样的：

      http://s15/mw690/001pJ93szy72JnAyPX8ae&690

Figure 3. 最简单的进化树图，犹如白纸一张

相较Figure3，Figure1和Figure2不仅为我们展示了不同物种基本的进化关系，而且还展示了分歧时间及所处的地质时代，基因家族扩张收缩的情形等重要的信息。这些精美的图片，除了用常规的工具外，还应用了高级的修图工具。如果我们不会PS，也不会SVG或者R，那么我们如何做出和Figure1,2可以媲美的进化树呢？

那么今天笔者就介绍几款有意思的进化树作图工具，也让小白的我们能做出赏心悦目的进化树图。

        【工具一】 MEGA6/7 + AI

笔者当时做青稞基因组进化图时，就是采用MEGA6³+AI组合来完成文章的Figure2a的。

最开始用流程做出来的进化树图是这样的：

      http://s12/mw690/001pJ93szy72JnJOezV8b&690Figure 4.青稞分歧时间估算图

这是流程做出来最简单的图样，几乎没有什么修改。最后用在文章的附件中，详见文章⁴图Fig.S10。

当时打算在正文放一张进化树相关的图，于是就把newick文件导入MEGA软件中，调成圆形图。然后在AI中把青稞和小麦及小麦祖先种这一枝用一个淡色的椭圆形标记出来，这样突出了本图的重点。虽然还是有点丑陋，但是比最开始的Figure4 好了很多。

http://s6/mw690/001pJ93szy72JnPD30h75&690
Figure 5.青稞基因组进化树图，原文Figure2a

MEGA是一款非常好用的系统发生分析的工具，其中关于进化树美化的部分，可以做到树形调整，标记，文字修改，图片添加等基本的修改。

笔者试着用MEGA的诸多基本功能，对Figure3的白纸做了一些调整

http://s13/mw690/001pJ93szy72JnZ4sLqfc&690
Figure 6. 用MEGA6优化Fig.2

修改的内容包括树形结构调整（For Balanced Shape功能），分组标记等。

  当然啦，MEGA美化系统发育树的功能还是很基础的，为了画出更有逼格的进化树，我们还需要请出更好用的工具。

接下来轮到iTOL登场。

【工具二】在线美化工具iTOL（Interactive Tree Of Life）

目前这款工具已经更新到version3⁵了。地址是：http://itol.embl.de/index.shtml 网站首页标明这款软件的功能是：Display, annotation and management of phylogenetic trees.

新版本拥有完全的所见即所得（what you see is what you get）的输出功能。除了支持常见的几种树形结构外，最重要的是可以上传自己的数据集，能在原始树基础上增加更多更重要的信息。V3版本中最多支持13种不同的数据集类型，同时支持多个数据集展示。

Figure 7 是笔者曾经做过的一张圆形树图。由于Taxa数目较多，只展示细节部分。

 http://s16/mw690/001pJ93szy72Jo5wlNZff&690

Figure 7. iTOL 美化进化树示例。

关于更多的细节可以查看官网的video 教程。http://itol.embl.de/video_tutorial.cgi

最后笔者介绍一款北京基因组所团队开发的在线进化树展示工具。这也是一款很好用的美化工具evolview。

【工具三】 evolview 网址： http://www.evolgenius.info/evolview/#mytrees/

今年，这款工具更新到Version2⁶了，增加了几款展示方式。

Figure 8. Evolview V2 new annotation dataset showing(1)

http://s1/mw690/001pJ93szy72Joi25oY70&690

Figure 9. Evolview V2 new annotation dataset showing(2)

http://s12/mw690/001pJ93szy72JomJBJpcb&690
Figure 10. Evolview 官网提供的show case 04

同时官网提供了16个Demo及6个Case供大家练习。

好吧。关于进化树展示优化的介绍就到这里了。

 【参考文献】

1 Yu Jiang, M. X., Wenbin Chen,Richard Talbot, Jillian F. Maddox,. The sheep genome illuminates biology of the rumen and lipid metabolism. science (2014).

2 Wang, Z. et al. The draft genomes of soft-shell turtle and green sea turtle yield insights into the development and evolution of the turtle-specific body plan. Nat Genet 45, 701-706, doi:10.1038/ng.2615 (2013).

3 Tamura, K., Stecher, G., Peterson, D., Filipski, A. & Kumar, S. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Molecular biology and evolution 30, 2725-2729, doi:10.1093/molbev/mst197 (2013).

4 Zeng, X. et al. The draft genome of Tibetan hulless barley reveals adaptive patterns to the high stressful Tibetan Plateau. Proceedings of the National Academy of Sciences of the United States of America 112, 1095-1100, doi:10.1073/pnas.1423628112 (2015).

5 Letunic, I. & Bork, P. Interactive tree of life (iTOL) v3: an online tool for the display and annotation of phylogenetic and other trees. Nucleic acids research, doi:10.1093/nar/gkw290 (2016).

6 He, Z. et al. Evolview v2: an online visualization and management tool for customized and annotated phylogenetic trees. Nucleic acids research, doi:10.1093/nar/gkw370 (2016).

                                                                                                                    2016/6/19