三代测序数据简单分析

三代测序数据简单分析

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