三代測序數(shù)據(jù)簡單分析
原創(chuàng) saber universebiologygirl
簡單介紹：
三代測序技術(shù)讀長較長，針對比較小的基因組像只有16kbp的人類線粒體DNA（mtDNA）是非常適用的，未來其應(yīng)用應(yīng)該會在測序錯誤率降低后得到顯著提高。今天給大家介紹一下之前所做的mtDNA三代測序數(shù)據(jù)組裝。因為也是初次接觸數(shù)據(jù)組裝操作，有不全面的地方請讀者見諒，也可在留言區(qū)留言指正。
（Wikipedia）
線粒體基因組（mitochondrial DNA，mtDNA）為環(huán)狀多拷貝形式存在于線粒體中，不同mtDNA分子可能存在不同突變。
1，數(shù)據(jù)比對
首先拿到測序數(shù)據(jù)，如果已經(jīng)有標(biāo)準(zhǔn)的參考基因組（例如mtDNA現(xiàn)在使用的是早期組裝的英國人的一條序列rCRS作為參考），我們可以使用李恒編寫的三代測序數(shù)據(jù)比對軟件minimap2將原始數(shù)據(jù)GZ.fq比對到參考基因組reference.fa上。
這里簡單介紹下這個軟件，minimap2的安裝以及基本的使用lh3已經(jīng)在github上都作了說明（見(https://github.com/lh3/minimap2)）：
安裝
git clone https://github.com/lh3/minimap2
使用
# long sequences against a reference genome
值得注意的是Pacific Biosciences之前針對PacBio數(shù)據(jù)比對推薦使用的是Blasr，但是現(xiàn)在也建議使用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將測序數(shù)據(jù)GZ.fq與參考基因組reference.fa比對:
minimap2 -ax map-ont reference.fa GZ.fq > GZ.sam
參數(shù)說明：
-a輸出為sam格式： -a output in the SAM format (PAF by default)
這樣可以檢測測序數(shù)據(jù)中大概有多少是我們需要的序列、可以用來組裝。
2，使用samtools查看比對之后的sam文件
samtools idxstats查看比對情況:
samtools idxstats file.sam
samtools view篩選數(shù)據(jù)，sam轉(zhuǎn)bam等。
這里簡單介紹一下samtools view的用法：
samtools view 用的較多的除了bam、sam文件之間的轉(zhuǎn)換外，還有下面兩個參數(shù)來篩選reads，根據(jù)序列比對之后的flags來篩選reads
-f INT only include reads with all of the FLAGs in INT present [0] ###只輸出flags為對應(yīng)值的序列
samtools flags 后面跟2進(jìn)制的flags值，可以查看詳細(xì)意義：
$ samtools flags 2304
3，數(shù)據(jù)組裝
前面的步驟主要是為了檢測數(shù)據(jù)，是檢測mtDNA變異時的常見步驟。拿到測序數(shù)據(jù)，可以直接進(jìn)行組裝。這里使用的fastq文件比較小，如果文件較大，canu運(yùn)算時間會非常久。
canu組裝過程主要包括矯正、修剪、組裝三個過程：
correction, trimming and assembly
每一步都可以單獨(dú)進(jìn)行，且對應(yīng)很多參數(shù)，如果沒有額外的參數(shù)需要調(diào)整，一般可以一條命令行執(zhí)行所有操作，直接用canu組裝原始數(shù)據(jù):
canu -p GZ_asb1 -d ./result_dir genomeSize=16569 minReadLength=100 minOverlapLength=50 -nanopore-raw GZ.fq
參數(shù)說明：
-p <assembly-prefix> ###生成文件的前綴
組裝完成后在result_dir目錄下會得到組裝后的序列文件
GZ_asb1.unitigs.fasta
4，組裝后的序列優(yōu)化
主要是nanopolish對組裝后的數(shù)據(jù)GZ_asb1.unitigs.fasta進(jìn)行優(yōu)化（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文件，好像也能運(yùn)行，軟件help中給出的最好是接fasta文件，所以這里可能需要轉(zhuǎn)換一下，將fastq轉(zhuǎn)為fasta文件） 使用參數(shù)說明如下：
-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為最后優(yōu)化的組裝結(jié)果。
MUMmer比較組裝的fasta與參考基因組序列之間的差別
MUMmer3.23/dnadiff --prefix GZ_asb1_polished.dnadiff reference.fa GZ_asb1_polished_genome.fa
官方也給了一個例子，包括canu組裝、優(yōu)化等步驟（how to polish a genome assembly）： https://nanopolish.readthedocs.io/en/latest/quickstart_consensus.html
另外組裝完成后，也可以用muscle、mafft這些多序列比對軟件，或者blast，將組裝的序列和參考序列進(jìn)行比對（由于mtDNA只有一條，所以檢測多序列比對結(jié)果，即可知道組裝的好壞），這里小編自己組裝的mtDNA序列，得到后，使用多序列比對與rCRS比對完成后，編寫腳本檢測的每個位置堿基差別，這里就不做過多介紹了。
結(jié)尾：因為也是第一次接觸三代數(shù)據(jù)，以及組裝數(shù)據(jù)，可能有些地方介紹的不夠詳細(xì)，如果大家還有什么不明白的可以在留言區(qū)指出，共同探討，關(guān)于nanopore的組裝今天就介紹到這里，希望能對大家有所幫助，請持續(xù)關(guān)注我們，謝謝！
