lyao222lll的个人博客分享 http://blog.sciencenet.cn/u/lyao222lll

博文

(动物)线粒体基因组二代测序组装经验分享

已有 4192 次阅读 2019-9-28 13:51 |个人分类:经验分享|系统分类:科研笔记| _基因组组装

(动物)线粒体基因组二代测序组装经验分享

 

介绍完了叶绿体的二代测序组装,顺便再讲一下线粒体的吧,相较于叶绿体,线粒体简单多了。不过我这里指的是动物线粒体,平均大小十几kb,二代测序组装挺容易得到完成图的。真菌、植物等的线粒体,我也没碰过不清楚它们的具体结构。不过考虑到真菌线粒体基因组也很小,所以下文的方法对真菌线粒体可能也适用;但是听说植物线粒体很复杂了,所以您若是研究植物线粒体的,绕道吧。

话不多说了,上正文。同样地,因为本文是经验分享,所以对以下涉及的软件,以及一些常规概念,默认大家已知了。

 


(动物)线粒体二代测序组装过程(个人常用经验)


叶绿体基因组测序时,可以使用新鲜的植物叶片提取总DNA直接测序,由于叶绿体细胞器在新鲜叶片组织的细胞中很丰富,所以测序数据中叶绿体基因组的reads含量也可观,通常无需单独提纯叶绿体。动物线粒体类似,它也是细胞中多拷贝的细胞器,含量可观,通常也无需单独提纯,直接提取组织的总DNA测序即可;再加上它基因组也很小,组装起来也比较轻松。当然如果你还是选择先分离细胞器再提取DNA测序,效果会更好(相当于目标基因组被“富集”了,组装会更方便)。如果实在懒得自己搞,可以交给测序公司顺便帮你把样本处理了,这点他们会更有经验5.gif(你信得过他们的话,我见过有人想测鱼线粒体的,然后直接生鲜发货给测序公司了一条鱼,样本处理时发现尾巴还在动的,哈哈)。

好了,下文假设我们想研究某昆虫的线粒体基因组,直接提取组织的全部DNA后,测序,而非先提纯线粒体再测序的方式。然后期望通过来自全基因组的测序结果,组装得到其线粒体基因组序列。

 

0、准备工作



测序数据clean reads

测序后,拿到手了fastq格式的测序数据clean dataclean reads)。如下所示,通常情况下,为了节省空间会将fastq压缩为gz格式存储。

1.png

 

关于参考基因组

叶绿体组装策略不同,对于动物线粒体的拼接,参考基因组并不是那么重要。首先动物线粒体没叶绿体那么保守,靠近缘物种线粒体序列在测序数据中挑选reads,再用于组装,效果一般不咋的,反而容易产生更多gap;其次线粒体序列非常小,而且同样是多拷贝,大可以直接使用全部测序数据从头拼接,然后我们再从组装结果中挑选合适的contigs/scaffolds序列。此外,对于后续contigs/scaffolds序列的定位和定向,也并非一定要有参考基因组。

详见下文具体内容。

 

NT数据库Nucleotide Sequence Database

既然不用参考基因组,那么该如何确定组装结果中哪些是线粒体序列呢?我们需要用到本地版的NCBI核酸数据库(下文简称NT库),通过后续进行本地BLAST,在组装结果中挑选出比对到线粒体的contigs/scaffolds序列。详见下文具体内容。

NT库获取链接:ftp://ftp.ncbi.nih.gov/blast/db/。数据库很大,慢慢下。

 

1、基因组的二代初步组装


1步,直接拿所有的clean reads开始作组装。

基因组二代组装工具挺多的,比如SOAPdenovoABySSSPAdesA5-miseq等。大家根据实际情况选择吧。这一步多尝试几个拼接软件(至少两个吧,只用一个可能不太够用,原因后面会提到),分别得到不同的拼接结果。

例如分别使用SPAdesA5-miseq作组装。

##1、基因组组装,使用 clean data 直接进行基因组的二代组装
#使用 SPAdes
spades.py -1 test_mitochondria_cleandata_R1.fastq.gz -2 test_mitochondria_cleandata_R2.fastq.gz -o ./spades -t 4 --cov-cutoff auto --careful
 
#使用 A5-miseq
a5_pipeline.pl --threads=4 test_mitochondria_cleandata_R1.fastq.gz test_mitochondria_cleandata_R2.fastq.gz a5miseq

 

2、在组装序列中挑选出目标线粒体序列


由于clean data数据中,绝大多数是核基因组序列,因此直接使用全部clean reads拼接的话,几乎所有得到的contigs/scaffolds序列都是核基因组的,同时还不可避免的存在一些细菌的基因组序列片段,真正属于线粒体的组装序列只有几条而已。第2步就是要从中找到这些属于线粒体的contigs/scaffolds序列。

下面是常用的3种挑选方法,结合使用。

 

BLAST比对数据库,初步获取线粒体序列

NT库就在这个时候用的。通过将所有得到的contigs/scaffolds序列与NT库中收录的核酸序列做BLAST比对,定位目标序列。

##2、BLAST 比对,确定组装序列的来源
#例如我们以上述 SPAdes 组装结果 fasta 文件“scaffolds.fasta”为例
 
#blastn 核酸比对,指定 NT 库路径,默认对每条序列输出一条最佳 hits 
blastn -db /database/nt/nt -query spades/scaffolds.fasta -out spades_blast -num_threads 4 -num_descriptions 1 -num_alignments 1 -dust no
 
#对 blast 结果格式作个转化,perl 脚本获取链接:https://pan.baidu.com/s/1-HkUh_C9JgYH9q-J2R7ZDA
perl blast_trans.pl spades_blast spades_blast.txt
 
#根据注释描述,提取其中命中到“线粒体”的序列比对结果
grep 'mitocho' spades_blast.txt > spades_blast.select.txt

查看“spades_blast.select.txt”,该文件中只保留了能够比对至数据库中已知线粒体序列的结果,即可大致确定哪些contigs/scaffolds序列是来自线粒体的。

2.png

 

根据长度挑选序列

继续查看“spades_blast.select.txt”中的内容,BLAST结果中也记录了各序列的长度。通常动物线粒体长度十几、二十kb左右,我们先把里面较长的序列找到。

3.png

 

根据“覆盖度”挑选序列

此外,由于线粒体是多拷贝的,所以在测序数据中,它的reads丰度相对还高一些。某些组装工具,如SPAdes等,还能给出组装序列的一个“覆盖度”指标。对于这个“覆盖度”,有的工具是指k-mer平均深度,有的是指reads平均深度,无论哪种,我们都可以判断出,线粒体contigs/scaffolds序列的“覆盖度”是会高一些的,所以也可以作为挑选依据。

如果你的测序样本不是直接提取的组织总DNA,而是先提取纯化了线粒体后,再提取的DNA,那么这个“覆盖度”指标更是一个靠谱的选择依据。

4.png

 

3、基因组环化


综上一步的contigs/scaffolds识别结果,把符合条件的序列挑选出来。对于contigs/scaffolds序列的挑选,可使用SAMtools等工具,或者自写脚本提取,就算手动打开fasta文件选择序列也不麻烦。

之后查找这些contigs/scaffolds序列之间是否有重叠(overlap)区域,并根据这些重叠区域连接序列。

有这么几点要注意的:

 

先试一下最长的那条序列,是否自己成环了

就我们的示例而言,上图中我们看到获得的最长contigs/scaffolds序列约接近17kb,差不多就是一条动物线粒体基因组的平均长度,那么,这个是否就已经组装完整了呢?

如果是的话,真的就很省事了。二代测序组装工具普遍存在这么一个问题,它们无法判断环状序列,所以有时对于已经组装成环的序列,会多出来那么一小截。这个就可以成为我们判断序列成环的依据。

5.png

事实上,对于动物线粒体基因组的组装,一半以上的情况就是这样的,最长的那条contigs/scaffolds序列就自己成环了。这种情况下,参考基因组真的就可直接忽略了。

可以结合共线性分析、BLAST比对等工具来确定首尾序列的相似性。如下示例通过调用blastn来完成。

##判断自环
#例如最长的序列单独取出得到一个 fasta 文件“long_scaffolds.fasta”
#我们再取其中 5' 端几百 bp 的碱基,得到“long_scaffolds.head.fasta”
 
#两个文件相互 blastn 核酸比对
#perl 脚本获取链接:https://pan.baidu.com/s/1-HkUh_C9JgYH9q-J2R7ZDA
perl blast_parser.pl long_scaffolds.head.fasta long_scaffolds.fasta
 
#输出结果和上述 blast 结果样式一致

 

如果最长的序列没能自环,那么我们就需要继续连接其他序列

首先看一下有没有合适的参考基因组可以做辅助。

BLAST比对结果中,给出了这些contigs/scaffolds序列最佳命中的参考线粒体基因组序列名称。我们可以从中找一条最相似的参考基因组,通过IDNCBIEMBL等数据库中下载它们,辅助我们确定这些contigs/scaffolds序列在基因组中的相对位置和方向(定位和定向),以继续往完整的环状线粒体基因组序列搭建。此外,参考基因组还能帮助我们确定自己线粒体基因组的最终长度范围。

6.png

##例如通过MUMmer 共线性分析定位 contigs/scaffolds 的顺序
#我们挑选出的 contigs/scaffolds 序列合并至一个 fasta 文件中,例如命名“select.fasta”
#参考基因组序列命名“ref1.fasta”(要单一物种的,不要混合物种的)
 
nucmer --mum -p mitochondria ref1.fasta select.fasta
delta-filter -m mitochondria.delta > mitochondria.filter
show-coords -T -r -l mitochondria.filter > mitochondria.1coords
mummerplot --postscript -p mitochondria mitochondria.delta
ps2pdf mitochondria.ps mitochondria.pdf

除了共线性分析,还有很多方法,大家可以自己试试,不再多说了。

但是有个问题,不同于叶绿体,动物线粒体的保守性没那么强,可能最相似的参考基因组和你的组装序列之间的差异仍然很大,这时候通过参考基因组去给我们的contigs/scaffolds序列定位定向就不是很妥当,这点需要大家特别注意!

 

如果有合适的参考基因组,确定了contigs/scaffolds序列间的顺序,我们下一步就可继续根据这些序列之间的重叠区,连接序列。期间该取反向互补的序列,注意反向互补。

如果没有找到合适的参考基因组,其实也不是很影响。首先动物线粒体本身不大,并且你挑选出的contigs/scaffolds序列也肯定不是很多,相互两两之间比较,查看是否具有重叠区,有的话就连接起来,期间该取反向互补的序列注意取反向互补,如此多试几次也差不多就能都连接起来了。

同样地,可以结合共线性分析、BLAST比对等工具来确定contigs/scaffolds序列间,末端之间是否存在有高可信度的一致区域,即重叠区。

 

4、填补gap


如果在第3步并未完全将线粒体基因组环起来,中间还存在gap,那么这一部分的内容将会是有用的。补gap肯定是最糟心的一步了,以下提到的方法主要基于个人经验。经验尚浅,大家若是有更好的方法,还请分享哈。

 

比较多个组装软件的拼接结果

在第1步组装时,建议使用多个软件作拼接。我们可知,不同的软件拼接效果肯定有所差异。

首先,运气好的话,某个软件的结果,其最长的那条属于线粒体的组装序列可以直接自环,就省事了。所以多试几种拼接工具是个挺不错的选择。

其次,对于出现gap的位置,可以结合多个软件的拼接序列,相互搭建,看不同软件所得序列之间是否能够相互填补空缺。

7.png

 

关于gap填补工具、手动延伸等其它方法

这些方法在小基因组的拼接中,很通用,我在前文“叶绿体组装”中有提到,可参考前文,因此这里就不再作更多叙述了。叶绿体基因组比线粒体要大得多,因此它适用的方法,动物线粒体更是可以。

 

其实下文“碱基校正”等过程也是通用的方法,也可参考前文,不过为了完整性,这里我还是把它们列在本文中了。

 

5、碱基校正


环化完成后,完整的线粒体基因组序列拼接好了,后续再加个碱基校正的过程。

Illumina的二代测序错误率约1%,也就是平均每100个碱基中会出现一个错误碱基,可能表现为“碱基的突变”,也可能表现为“碱基的插入/缺失”。在组装过程中,这种错误可能会被引入,使得基因组中的某些位点不太可靠(如,出现模糊碱基,或者极个别的N碱基)。

我们可以将测序clean reads重新与我们的组装结果基因组序列对齐,根据高深度的覆盖reads去识别、校正模糊的位点。

8.png

##例如,结合使用 BWA、SAMtools 实现测序 reads 和基因组的对齐过程,使用 Pilon 软件对碱基进行校正
#final.fasta 是最终环化后的结果
bwa index final.fasta
bwa mem -t 4 final.fasta test_mitochondria_cleandata_R1.fastq.gz test_mitochondria_cleandata_R2.fastq.gz | samtools view -@ 4 -bS -F 12 | samtools sort -@ 4 > bwa.sort.bam
samtools index bwa.sort.bam
pilon --genome final.fasta --frags bwa.sort.bam --output final.pilon1 --changes --threads 4

对于上述示例,Pilon校正后,得到两个文件,“final.pilon1.fasta”和“final.pilon1.changes”。若是changes文件大小为0(空文件),则其对应的fasta文件可作为最终结果;若不为0,则还需再基于新得到的fasta文件重新执行校正;以此类推,直到changes文件为0为止。动物线粒体,我不信校正不到0

事实上,由于Pilon也可以识别小段的N碱基序列,因此除了用于碱基校正外,或许也能用于填补小部分gap

 

6、最终检查


我们可以将原始的clean reads比对到最终的线粒体基因组上,得到bam文件,然后使用基因组浏览器(如IGVTablet等)打开bam文件以及基因组fasta文件,仔细查看测序reads在基因组序列中的覆盖情况,判断:

1reads覆盖是否均一,基因组中哪些区域reads覆盖深度较少,低深度的reads是否能有效支持这一段序列的拼接。

2)是否存在疑似错拼的位点,比方说某个位置两端的reads出现很大“断层”。

3)是否还有个别位置的碱基疑似由测序错误导致的类似于“SNP”、“InDel”情形。

4)对于重复序列较多的区域,reads是否能够有效跨越。

 

都确认无误了,再视情况调整下环状基因组序列的起点位置就可以了。可以参照近缘参考基因组的起点位置来调,或者根据基因组注释结果来。

 

7、关于动物线粒体功能注释


提到线粒体注释的话,动物线粒体的数据库还是挺丰富的(相较之下,植物线粒体连数据库都少的可怜)。尽管不同动物物种之间,线粒体基因组序列碱基组成差异很大,但是序列中主要结构、功能区类型却差不多是固定的。以高等动物线粒体为例,普遍都是由13个蛋白编码基因、2rRNA16tRNA1OL区和1OH区组成。当然,也会有例外的情形,不过也都差不多。注释起来也算是比较方便吧,同时注释结果也可为我们组装的基因组提供完整性参考:若缺少的结构区过多,那么组装过程中肯定出现了错误。

我在前文分享过两个可用于动物线粒体基因组注释的工具,MITOSGeSeq,希望有所帮助。当然还有很多方法,还请大家自己学习了。

 



http://blog.sciencenet.cn/blog-3406804-1199850.html

上一篇:KEGG自动注释工具KAAS简介及KO注释表的简单处理
下一篇:基因组重复序列检测工具RepeatMasker安装及使用

0

该博文允许注册用户评论 请点击登录 评论 (0 个评论)

数据加载中...
扫一扫,分享此博文

Archiver|手机版|科学网 ( 京ICP备07017567号-12 )

GMT+8, 2020-8-12 19:04

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部