nanopore 分析流程

作者: wo_monic | 来源:发表于2020-10-16 20:35 被阅读0次

单分子纳米孔测序

https://poretools.readthedocs.io/en/latest/index.html#

三代测序技术
 测序技术
 2019年8月综述纳米孔测序
 液体活检的纳米孔测序：肺癌患者无细胞DNA拷贝数变异的分析
 纳米孔测序及其临床应用

nanopor技术专栏
技术专栏详细介绍了从下载数据到软件分析的全流程的过程，提供了示例代码和示例数据。

一个人的全基因分析流程

使用的软件

使用HDFView 查看Fast5文件格式。
https://www.hdfgroup.org/downloads/hdfview/
使用ont_fast5_api软件对fast5文件进行拆分与合并

ont_fast5_api 网址

https://github.com/nanoporetech/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.软件安装

比对 minimap2 github
拼接 miniasm
gfatools
ngmlr
graphmap
拼接last
MarginAlign
过滤 filtlong
拼接 canu github
拼接 Smartdenovo
拼接 wtdbg2
拼接 NextDenovo
纠错、校正NextPolish
拼接 flye
基因组装配的质量评估工具 quast
比对 mummer 4
nanopack
pip install nanopack

基因组可视化 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格式展示纠错的位点。