生信

三代测序数据简单分析

2020-06-02  本文已影响0人  天明豆豆

三代测序数据简单分析

原创 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的组装今天就介绍到这里,希望能对大家有所帮助,请持续关注我们,谢谢!

上一篇下一篇

猜你喜欢

热点阅读