本文转载自知乎 nanopor技术专栏，修改了其中部分代码。

单分子纳米孔测序

使用的软件

使用HDFView 查看Fast5文件格式。

https://www.hdfgroup.org/downloads/hdfview/

使用ont_fast5_api软件对fast5文件进行拆分与合并

ont_fast5_api 网址

安装

直接使用pip安装

pip install ont-fast5-api

自行编译安装

git clone https://github.com/nanoporetech/ont_fast5_api

cd ont_fast5_api

python3 setup.py install

使用案例

多条合并为一个文件

single_to_multi_fast5 -i fast5_files/ -s multi -n 4000 --recursive

一个文件拆分为多条

multi_to_single_fast5 -i multi/ -s single --recursive -t 6

1.碱基识别工具

DeepNano

poretools 官网

poretools安装

2.序列比对工具

GraphMap

MarginAlign

3.从头组装工具

nanocorrect，首先利用graph- based，greedy partial order aligner方法进行纠错，然后利用Celera Assembler将纠错后的reads进行组装，最后利用nanopolish对组装结果进行进一步提升

4、单核苷酸变异检测工具

PoreSeq 低通量下，仍然有高准确率

Nanopolish

MarginAlign中的marginCaller模块

5、共有序列的测序(consensus sequencing)方法

实战

1.软件安装

基因组可视化 win+unix tablet

软件给出官方地址，简单的自行安装。比较复杂的会写出详细安装过程

2. 质控

mkdir nanoQC

#检测测序质量

nanoQC ../2.rawdata/minion/all..fastq.gz -o nanoQC

#统计质量信息

NanoStat --fastq ../2.rawdata/minion/all.sra.fastq.gz --outdir statreports

2.1 质控

NanoPlot --fastq ../2.rawdata/minion/all..fastq.gz -t 16 --maxlength 40000 --plots hex dot pauvre kde -o nanoplot

Nanoplot --summary sequencing_summary.txt --loglength -o summary

选项参数：

-t：线程数目

-o, --outdir：输出结果目录

-p, --prefix：输出结果前缀

--color：点的颜色

--N50 表示在序列读长的直方图中显示N50的标识

--title：标题

--downsample ：在输入文件中随机抽取n条序列进行处理

--minlength：忽略nbp以下的reads

-- fastq：输入fastq格式文件

-f：图片类型

--plots：绘图类型，kde,hex,dot,pauvre

2.2 过滤(nanofilt和filtlong二选一即可)

使用nanofilt过滤，nanofilt无法识别压缩文件，需要先解压。

gunzip -c ../2.rawdata/minion/all.fastq.gz | NanoFilt -q 7 -l 1000 --headcrop 50 --tailcrop 50| gzip > clean.NanoFilt.fastq.gz

选项参数：

-l ,--length :过滤掉小于此长度的序列

--maxlength ：过滤掉超过此长度的序列

-q , --quality ：过滤掉低质量序列

--minGC：过滤掉GC含量小于此百分比的序列

--maxGC：过滤掉GC含量大于此百分比的序列

--headcrop：从头部切掉n bp

--tailcrop：尾部切掉n bp

结果解读

输入一个压缩格式的fastq文件，结果软件处理，已经将长度小于1000，同时整条序列平均质量值小于Q7的过滤掉，同时将每条序列首位各剪掉50bp。这些剩余的序列则为“clean data”。可以使用NanoPlot重新进行一下质控。

使用filtlong过滤

filtlong --min_length 1000 --min_mean_q 7 ../2.rawdata/minion/all.fastq.gz | gzip >clean.filtlong.fq.gz

filtlong还可以以二代的数据为参考进行过滤。

选项参数：

--min_length ：最短长度

--min_mean_q：平均Q值

--keep_percent：保留最好数据的百分比，80%，直接写80，不能写0.8

--target_bases：保留多少数据，单位为bp

2.3 比对

2.3.1 使用minimap2进行比对

构建基因组的索引文件

minimap2 -d ref.mmi ref.fa ＃索引

比对

minimap2 -a ref.mmi reads.fq > alignment.sam ＃对齐

#等价于

minimap2 -ax map-ont ref.mmi reads.fq >alignment.sam

常用选项参数

主要分成五大类，索引(Indexing)，回帖(Mapping)，比对(Alignment),输入/输出(Input/Output),预设值(Preset)。

-x ：非常中要的一个选项，软件预测的一些值，针对不同的数据选择不同的值

map-pb/map-ont: pb或者ont数据与参考序列比对；

ava-pb/ava-ont: 寻找pd数据或者ont数据之间的overlap关系；

asm5/asm10/asm20: 拼接结果与参考序列进行比对，适合~0.1/1/5% 序列分歧度；

splice: 长reads的切割比对

sr: 短reads比对

-d :创建索引文件名

-a ：指定输出格式为sa格式，默认为PAF

-Q ：sam文件中不输出碱基质量

-R ：reads Group信息，与bwa比对中的-R一致

-t：线程数，默认为3

2.3.2

2.4 sam转bam,排序

#samtools处理

samtools sort -@ 8 -o bam -o s0137.sorted.bam s1037.sam

samtools index s0137.sorted.bam

samtools faidx ref.fq

2.5tablet可视化基因组

tablet支持Windows，linux.

下面是Windows上的操作

2.5.1 将准备好的文件，至少四个，参考序列fasta格式及索引fai文件，比对并排序后的bam文件及索引bai文件；

ref.fq

ref.fai

ref.gff

s1037.sorted.bam

s1037.sorted.bam.bai

2.5.2、打开tablet，首先导入bam文件，然后导入fasta文件，索引文件无需导入，但必须与对应文件在同一目录下，(tablet对中文支持不好,路径不要有中文)；

2.5.3 、可以通过鼠标调整显示模式；

2.5.4、tablet还支持导入gff格式结果，方便查看固定区域

2.6 重新拼接基因组(多个软件分别拼接后，再从中选择最优的)

2.6.1 canu

软件安装一定要安装github的版本。

canu -d canu -p canu genomeSize=3g maxThreads=96 -nanopore-raw ../4.filter/clean.filtlong.fq.gz >canu.log

选项参数：

-p：输出前缀

-d：输出结果目录

-nanopore-raw：输入的为没有纠错过得nanopore数据

-num_threads：CPU线程数

genomeSize：设置预估的基因组大小，这用于让Canu估计测序深度；

maxThreads：设置运行的最大线程数；

rawErrorRate：用来设置两个未纠错read之间最大期望差异碱基数；

correctedErrorRate：则是设置纠错后read之间最大期望差异碱基数，这个参数需要在 组装 时多次调整；

minReadLength：表示只使用大于阈值的序列

minOverlapLength：表示Overlap的最小长度。提高minReadLength可以提高运行速度，增加minOverlapLength可以降低假阳性的overlap。

总结

1、有纠错步骤；

2、部分基因组拼接效果比较好；

3、默认会占用所有CPU，非常耗时，慎重使用；

4、有些数据无法拼接出结果。

2.6.2 miniasm拼接，gfatools转换格式

使用miniasm拼接首先需要使用minim2将测序数据进行自身比对，查找共有区域，生成paf格式文件。注意使用minimap2比对的时候一定要正确设置好-x选项，nanopore拼接需要使用ava-ont选项。然后使用miniasm进行拼接，miniasm拼接不会直接生成fasta序列，而是会生成gfa格式文件。最后使用gfatools提出拼接好的序列。

minimap2 -x ava-ont -t 12 ../4.filter/clean.filtlong.fq.gz ../4.filter/clean.filtlong.fq.gz | gzip -1 > reads.paf.gz

miniasm -f ../4.filter/clean.filtlong.fq.gz reads.paf.gz >reads.gfa

gfatools gfa2fa reads.gfa >miniasm.fa

注意：第一步，是用minimap2对测序的文件自身进行比对。和前面的比对到基因组上不一样。

总结

1、首先利用minimap2比对，运行速度非常快；

2、软件不能直接生成fasta格式，需要使用gfatools生成；

3、没有纠错过程，碱基错误率较高，需要后续进行纠错；

4、拼出来的基因组可能比真实基因组要小，有时候小很多，需要注意。

2.6.3 Smartdenovo拼接(阮珏 出品)

# 生成脚本

smartdenovo.pl -c 1 -t 8 -J 500 -k 16 -p smartdenovo ../4.filter/clean.filtlong.fq.gz >smartdenovo.mak

#运行脚本

make -f smartdenovo.mak

选项参数：

-p STR 输出结果前缀，默认wtasm

-e STR overlap方法，zmo或者dmo

-t INT 线程数，默认为8

-k INT overlap连接kmer大小

-J INT 最小read长度

-c INT 是否生成一致性序列

2.6.4 wtdbg2拼接 (阮珏 出品)

快速运行模式，直接使用wtdbg2.pl脚本

wtdbg2.pl -t 96 -x ont -g 3g -o wtdbg2 reads.fa.gz

选项参数：

-t 线程数量

-x 数据模式(nanopore的数据选择ont)

-g 基因组大小(根据实际物种决定)

-o 输出目录

总结：

1、对于nanopore数据，wtdbg2可能会产生比真实基因组小的集合。

2、当输入fasta和fastq格式的多个文件时，请先输入fastq，然后再输入fasta。否则，程序将无法在fastq中找到'>'，并在一次读取后附加所有fastq。

2.6.5 NextDenovo拼接基因组

支持py2和py3

先用pip安装依赖包Psutil和Drmaa (Only required by running under non-local system)

#生成一个配置文件，可以直接从软件安装目录下拷贝

cp /ifs1/Software/biosoft/NextDenovo/doc/run.cfg ./

#生成一个reads列表，名为input.fofn

ls ../../filter/clean.filtlong.fq.gz > input.fofn

#将input.fofn添加到配置文件run.cfg的input_fofn关键字后面，默认不用修改

input_fofn = ./input.fofn

#运行软件

nohup time nextDenovo run.cfg &

拼接结束之后，最终结果隐藏的比较深，在以下目录中，可以配合nextpolish进行纠错，得到更优的结果。

01_rundir/03.ctg_graph/01.ctg_graph.sh.work/ctg_graph00/nextgraph.assembly.contig.fasta

2.6.6 flye拼接基因组

flye --nano-raw ../4.filter/clean.filtlong.fq.gz -g 3g -o output -t 48 -i 3

常用选项参数：(不习惯长选项的可以直接使用短选项)

--pacbio-raw ：输入原始pacbio数据

--pacbio-corr ：输入纠错后的pacbio数据

--nano-raw：输入原始nanopore数据

--nano-corr ：输入纠错后的nanopore数据

--genome-size：预估基因组大小，用于评估覆盖深度

--out-dir：输出结果文件路径

--threads：cpu线程数据

--iterations：纠错迭代次数

--min-overlap：最小overlap连接大小

--meta: 拼接宏基因组数据

--plasmids: 拼接质粒数据

输出结果

最后结果目录中有三个文件比较重要。

1、assembly.fasta ：最终拼接得到的基因组序列，fasta格式。

2、assembly_graph.{gfa|gv} ：拼接过程中用到的repeat graph。

3、assembly_info.txt：拼接结果统计信息，也可以自己单独使用seqkit工具统计。

2.7用quast评估多个拼接结果

quast评估案例

软件适合以下应用场景：

1、得到不同软件拼接的基因组序列，想比较一下哪个软件拼接效果更好；

2、同一软件，比较不同选项参数拼接的结果，哪个更好；

3、比较拼接得到的基因组，与参考序列的相似性。

quast.py -r ref.fa canu.fa miniasm.fa wtdbg2.cns.fa smartdenovo.fa -o quast

-o --output-dir 输出结果目录

-r 参考序列文件

-g --genes 参考序列基因坐标，一般BED或者GFF格式文件

-m --min-contig 最小contig长度，也就是小于这个阈值的不参与计算

-t --threads 使用线程数目，默认使用四分之一的cpu

--help 列出全部选项参数，大部分情况下，默认这些选项即可

2.8 组装结果优化(polish)

2.8.0 组装优化前后对比

使用Mummer软件包中dnadiff工具将纠错前后的序列进行基因组的比对。dnadiff会直接给出一个统计报告，里面会列出两条序列之间差异的部分，然后我们可以使用mummerplot工具对统计进行进行粗略的可视化。

#拼接结果优化前后比较

dnadiff ../before.fasta after.fasta

mummerplot --filter --png -p all out.delta

gnuplot all.gp

-i 输入测序reads

-d 需要纠错的拼接结果

-o 输出结果文件

-m 芯片类型

-t 并行计算

2.8.1 利用medaka组装结果纠错

安装方式 pip install medaka

medaka_consensus -i ../4.filter/clean.filtlong.fq.gz -d ../6.quast/miniasm.fa -o medaka_output -m r941_min_high -t 4 >medaka.log

-i 输入测序reads

-d 需要纠错的拼接结果

-o 输出结果文件

-m 芯片类型

-t 并行计算

最终生成的 consensus.fasta则为纠错后的结果。

calls_to_draft.bam

calls_to_draft.bam.bai

consensus.fasta

consensus.fasta.split/

consensus_probs.hdf

优化前后比较

dnadiff ../6.quast/miniasm.fa medaka/consensus.fasta

mummerplot --filter --png -p all out.delta

gnuplot all.gp

2.8.2 pilon纠错

下载wget https://github.com/broadinstitute/pilon/releases/download/v1.23/pilon-1.23.jar，直接是二进制文件。

#minimap2比对

minimap2 -d miniasm.fa.min miniasm.fa

minimap2 -ax map-ont miniasm.fa.min ../4.filter/clean.filtlong.fq.gz >miniasm.sam

#对bam进行处理

samtools sort -@ 4 -O bam -o miniasm.sorted.bam miniasm.sam

samtools index miniasm.sorted.bam

#利用pilon进行纠错

java -Xmx32G -jar /share/softwares/pilon/pilon-1.23.jar --genome miniasm.fa --fix all --changes --bam miniasm.sorted.bam --output all_pillon --outdir all_pillon --threads 6 --vcf 2> pilon.log

参数详情

常用选项参数：

--genome提供输入参考基因组

--frags 表示输入 < 1kb 的文库BAM --jumps 输入 > 1kb 的文库BAM

-unpaired 输入非配对的BAM。

--output表示输出的前缀

--outdir表示输出文件夹

--changes 列出发生变化的部分,以FASTA形式保存

--vcf 以VCF形式保存。

--fix 声明对参考基因组做哪方面的改进，包括“snps”,”indels”,”gaps”,”local”, 默认是”all”

利用二代数据进行纠错

#对拼接结果建立索引

bwa index draft.fasta

#illumina比对排序建索引

bwa mem -t 2 draft.fasta ../0.rawdata/SRX5341420_1.fastq.gz ../0.rawdata/SRX5341420_2.fastq.gz | samtools view - -Sb | samtools sort - -@ 14 -o illumina.sorted.bam

samtools index illumina.sorted.bam

#利用pilon进行纠错

java -Xmx32G -jar /share/softwares/pilon/pilon-1.23.jar --genome draft.fasta --fix all --changes --frag --genome draft.fasta --fix all --changes --frags illumina.sorted.bam --output assembly_pillon --outdir assembly_pillon --threads 6 --vcf 2> pillon.log

结果解析

可能是因为java写的程序的缘故(不太确定)，pilon的运行速度比较慢。最终会生成*

*_pillon.changes：展示哪些位点经过了纠错；

*_pillon.fasta ：最终纠错得到的结果；

*_pillon.vcf ：vcf格式展示纠错的位点。

2.8.3 利用racon纠错