转自公众号：universebiologygirl
http://mp.weixin.qq.com/s?__biz=MzU1NDkzOTk2MQ==&mid=2247485464&idx=1&sn=2172488f07d242dc012750161919141d
关于噬菌体基因组组装，你不得不知道的一些技巧（一）
写在前面的
上期推文，主要讲了BAM文件中的雪茄值含义。CIGAR是记录了reads比对到参考基因组上的详细情况如：跳过、匹配，插入、缺失等。我们以基本动态规划算法，写了一个Perl脚本提取不同比对状态的subreads。上篇推文主要以CIGAR值为基础让大家更深入了解reads的比对过程。有利于大家利于BAM文件进行后续分析。
噬菌体研究简介
本次要介绍的是我研究方向内的一些内容了。可能国内研究噬菌体的是少数。因为噬菌体属于细菌病毒，一般只有研究细菌或者病毒的实验室有可能做噬菌体研究。而病毒或者细菌群体之多，研究病毒或细菌的实验室很大一部分是研究对人类或者畜类有害的细菌或病毒。例如牛逼哄哄的Ebola啊，目前国内很恐惧的非洲猪瘟病毒啊，还有细菌界的鼠疫、炭疽啊。虽说目前研究噬菌体的实验室相比于其他实验室不多，但是我相信未来研究这个的还是会越来越多。这是为啥呢？因为超级细菌啊！现在出现了越来越多的超级细菌，对多种抗生素都免疫。甚至有的细菌牛逼到通过对抗生素进行修饰把抗生素变成自己体内代谢的物质供自己代谢。抗生素在细菌面前是越来越菜了，而且国家卫生部和农业部都意识到了这个严重性。农业部已经发布了相关政策，要求在2022年以前全面取消兽用饲料中添加抗生素。这个政策之下，会有越来越多的实验室开展对噬菌体的研究。
而生物信息学应用在噬菌体研究中，能做些啥呢？首先我们肯定得获取噬菌体基因组序列啊，序列都没有还咋搞？但是噬菌体基因组序列获取相比于细菌、动物植物这些基因组的获取要难的多。这是为什么呢？我们先从基础知识开始讲，首先，对于细菌或者动植物，我们要获得其基因组序列较为容易，对于细菌我们可以培养单克隆，通过16S测序验证是不是我们想要的细菌即可以提取DNA进行全基因组测序。对于动物或者植物，我们可以直接获取其组织样本，进行DNA提取，可以在实验环节通过实验手段去除线粒体及叶绿体DNA。但是对于噬菌体就没有这么简单了。因为噬菌体的培养是通过细菌来培养的，在实验步骤上我们会先培养细菌，然后将细菌和噬菌体的混合液涂布到平板上，等一段时间后平板上会出现很多小噬菌斑。其余地方长满细菌。我们只能通过取噬菌斑上的DNA进行测序。但是这个过程中会有较大程度的细菌污染（虽然已经有一些成熟的实验手段减少细菌DNA污染，但是在实验实施过程中这个污染并不能完全去除，而且剩余细菌DNA的数量不少）。这是因为噬菌体进入到细菌体内后才能复制，行使生命周期。当噬菌体复制数量达到一定程度后，噬菌体会与细菌细胞壁接触，通过噬菌体合成的裂解酶裂解细菌细胞壁，才能让噬菌体释放。而细菌细胞壁被裂解后细菌虽然死亡，但是其DNA依然存在，会伴随噬菌体残留在噬菌斑上。我们提取噬菌斑上的DNA就无差别地提取到了裂解后细菌的DNA。这样送样测序时细菌DNA也被测序了。由于存在细菌的DNA，这对噬菌体的组装是非常不利的。可能读到这里你会纳闷有细菌基因组下载对应细菌的参考基因组然后过滤就行了呗，其实这样做不是每次都可以的。这是因为噬菌体中有一类群教prophage，翻译为前噬菌体。这些种噬菌体是整合在细菌基因组上，且有较大部分基因组与细菌基因组很类似。如果单纯采取mapping后过滤的方式会丢失掉大片段prohage的基因组区域。由于在噬菌体领域，很多实验室噬菌体都是自己筛出来了，并不清楚这个噬菌体基因组具体多大。只能通过预估的方式估计基因组大小。然后再进行组装，估计基因组大小方式有几种：1.通过荧光定量PCR；2.通过流式细胞仪;3.通过高通量测序的reads根据kmer估计基因组大小。这次我们来讲讲如何根据kmer估计噬菌体基因组大小。
kmer估计噬菌体基因组大小
Kmer是一段短的DNA片段，通常情况下kmer是一个奇数。如果了解De novo组装的过程，应该对kmer有过了解。De novo组装是基于数学中图论知识先将reads打散成一段一段的kmer，又将原始序列反向互补再取kmer，最后根据de Bruijn graph将逐渐连接。最后拼接成更大的contigs，甚至是scaffold。由于要保证组装的链唯一性，如果kmer取偶数，在遇到一些倒置回文序列时例如ATATATATAT之类的，序列在Kmer为偶数的情况下反向互补后与之前一样，这样就不利于在组装时判断kmer到底来自于哪条链。所以这也是kmer只能取奇数的原因。那么如何用kmer估计基因组大小呢？首先我们知道在理想状况下基因组大小大致等于测序碱基总数除以测序深度。测序碱基总数这个值容易获得，但是测序深度在实际情况下没有那么好算，所以可以选择通过计算kmer的深度推测测序深度。结合我在网上查阅的资料加上自己的理解，假设一个物种基因组大小为G，用k单位长度的kmer对整个测序reads进行划分，产生nk个kmer，设dk为kmer的期望深度，db为碱基期望深度，nb为碱基总数。由于kmer深度符合泊松分布，所以kmer深度的均值也是dk，那么就有如下公式
而估计基因组大小就如上图所示，用kmer的总数除以kmer峰对应的横坐标（kmer深度）即是，以上是根据kmer估计基因组大小的原理。那么实际操作，如何估计噬菌体基因组大小呢？使用jellyfish软件来估计噬菌体基因组大小。由于我们的噬菌体测序数据中存在细菌基因组，这部分细菌基因组的加入会增大基因估计的大小。jellyfish安装
wget https://github.com/gmarcais/Jellyfish/releases/download/v2.2.10/jellyfish-2.2.10.tar.gz
tar -zxvf jellyfish-2.2.10.tar.gz
cd jellyfish-2.2.10
./configure --prefix=$HOME
make -j8
make install
echo 'PATH=$PATH:~/src/jellyfish-2.2.10/bin' >> ~/.bashrc
source ~/.bashrc
使用方法
我们噬菌体基因组测序的kmer对应深度是1000，总数为276475783，对应噬菌体基因组大小估计是276475783/(1000*1024)=269.99kb, 而实际的基因组大小为138kb，这里有2倍差异。这个2倍差异不是一个随机出现的值，是因为噬菌体和细菌基因组有类似的区域。一般对于实验者提取噬菌体基因组，通常都会污染细菌基因组，如果仅是粗提取，根据实验室的组装经验而言通常会估计基因组大小为实际基因组大小的1.5~2倍。如果是精细提取DNA,（过滤细菌DNA），那么估计值会更接近于真实值。
下期预告
如何一次组装出尽可能完整的噬菌体基因组
作者简介
熊东彦，中国科学院武汉病毒研究所。擅长转录组分析，宏基因组分析，R语言编程，Perl语言编程
往期精彩推文
【编程技巧（四）】生信编程语言的经验之谈
【长文预警】你要的circRNA知识点全在这里
《生物信息学分析实用小技巧（七）：噬菌体基因组估计》来自互联网公开内容，收录仅供学习使用，如侵权请联系删除。本文URL：https://www.ezixuan.com/1020479.html