三代测序数据简单分析
三代测序数据简单分析
原创 saber universebiologygirl
image简单介绍:
三代测序技术读长较长,针对比较小的基因组像只有16kbp的人类线粒体DNA(mtDNA)是非常适用的,未来其应用应该会在测序错误率降低后得到显著提高。今天给大家介绍一下之前所做的mtDNA三代测序数据组装。因为也是初次接触数据组装操作,有不全面的地方请读者见谅,也可在留言区留言指正。
image(Wikipedia)
线粒体基因组(mitochondrial DNA,mtDNA)为环状多拷贝形式存在于线粒体中,不同mtDNA分子可能存在不同突变。
1,数据比对
首先拿到测序数据,如果已经有标准的参考基因组(例如mtDNA现在使用的是早期组装的英国人的一条序列rCRS作为参考),我们可以使用李恒编写的三代测序数据比对软件minimap2将原始数据GZ.fq比对到参考基因组reference.fa上。
这里简单介绍下这个软件,minimap2的安装以及基本的使用lh3已经在github上都作了说明(见(https://github.com/lh3/minimap2)):
安装
git clone https://github.com/lh3/minimap2
使用
# long sequences against a reference genome
值得注意的是Pacific Biosciences之前针对PacBio数据比对推荐使用的是Blasr,但是现在也建议使用minimap2:
Benchmarks show that pbmm2 outperforms BLASR in sequence identity, number of mapped bases, and especially runtime. pbmm2 is the official replacement for BLASR.
在它的github页面(https://github.com/PacificBiosciences/pbmm2),将minimap2整合成了pbmm2,其中可以将bam文件直接与参考基因组比对,用法非常简单友好:
A. Generate index file for reference and reuse it to align reads
这里我们使用minimap2将测序数据GZ.fq与参考基因组reference.fa比对:
minimap2 -ax map-ont reference.fa GZ.fq > GZ.sam
参数说明:
-a输出为sam格式: -a output in the SAM format (PAF by default)
这样可以检测测序数据中大概有多少是我们需要的序列、可以用来组装。
2,使用samtools查看比对之后的sam文件
samtools idxstats查看比对情况:
samtools idxstats file.sam
samtools view筛选数据,sam转bam等。
这里简单介绍一下samtools view的用法:
samtools view 用的较多的除了bam、sam文件之间的转换外,还有下面两个参数来筛选reads,根据序列比对之后的flags来筛选reads
-f INT only include reads with all of the FLAGs in INT present [0] ###只输出flags为对应值的序列
samtools flags 后面跟2进制的flags值,可以查看详细意义:
$ samtools flags 2304
3,数据组装
前面的步骤主要是为了检测数据,是检测mtDNA变异时的常见步骤。拿到测序数据,可以直接进行组装。这里使用的fastq文件比较小,如果文件较大,canu运算时间会非常久。
canu组装过程主要包括矫正、修剪、组装三个过程:
correction, trimming and assembly
每一步都可以单独进行,且对应很多参数,如果没有额外的参数需要调整,一般可以一条命令行执行所有操作,直接用canu组装原始数据:
canu -p GZ_asb1 -d ./result_dir genomeSize=16569 minReadLength=100 minOverlapLength=50 -nanopore-raw GZ.fq
参数说明:
-p <assembly-prefix> ###生成文件的前缀
组装完成后在result_dir目录下会得到组装后的序列文件
GZ_asb1.unitigs.fasta
4,组装后的序列优化
主要是nanopolish对组装后的数据GZ_asb1.unitigs.fasta进行优化(https://github.com/jts/nanopolish),MUMmer评价组装的好坏。
Nanopolish can calculate an improved consensus sequence for a draft genome assembly, detect base modifications, call SNPs and indels with respect to a reference genome and more.
首先使用nanopolish中的variants、vcf2fasta
Usage: nanopolish variants [OPTIONS] --reads reads.fa --bam alignments.bam --genome genome.fa
nanoplsh variants需要比对后的bam文件。将组装的后的序列作为参考基因组,将原始的fastq序列比对到其上:
minimap2 -ax map-ont -t 40 result_dir/GZ_asb1.unitigs.fasta GZ.fq |samtools sort -T tmp1 -o result_dir/GZ_asb1.sorted.bam -
(这里需要注意的是,之前我使用-r接的GZ.fq为fastq文件,好像也能运行,软件help中给出的最好是接fasta文件,所以这里可能需要转换一下,将fastq转为fasta文件) 使用参数说明如下:
-r, --reads=FILE the ONT reads are in fasta FILE ###原始的三代测序fastq文件
使用得到的vcf文件矫正组装的序列,得到新的序列:
nanopolish vcf2fasta --skip-checks -g GZ_asb1.unitigs.fasta GZ_asb1.polished.vcf > GZ_asb1_polished_genome.fa
注:GZ_asb1_polished_genome.fa为最后优化的组装结果。
MUMmer比较组装的fasta与参考基因组序列之间的差别
MUMmer3.23/dnadiff --prefix GZ_asb1_polished.dnadiff reference.fa GZ_asb1_polished_genome.fa
官方也给了一个例子,包括canu组装、优化等步骤(how to polish a genome assembly): https://nanopolish.readthedocs.io/en/latest/quickstart_consensus.html
image另外组装完成后,也可以用muscle、mafft这些多序列比对软件,或者blast,将组装的序列和参考序列进行比对(由于mtDNA只有一条,所以检测多序列比对结果,即可知道组装的好坏),这里小编自己组装的mtDNA序列,得到后,使用多序列比对与rCRS比对完成后,编写脚本检测的每个位置碱基差别,这里就不做过多介绍了。
结尾:因为也是第一次接触三代数据,以及组装数据,可能有些地方介绍的不够详细,如果大家还有什么不明白的可以在留言区指出,共同探讨,关于nanopore的组装今天就介绍到这里,希望能对大家有所帮助,请持续关注我们,谢谢!