三代测序nanopore下机数据分析实战

  现在三代测序越来越便宜了，很多小伙伴都对三代数据分析有着需求和好奇，这次直接一次性把三代流程给到大家，快上手操练起来吧。

1、软件安装

略（需要的话请留言或者加群）

2、质控（2选1）

（1）nanoQC

mkdir nanoQC#检测测序质量nanoQC  输入数据路径.fastq.gz -o nanoQC#统计质量信息NanoStat --fastq  输入数据路径.sra.fastq.gz --outdir statreports

（2）NanoPlot

#输出质量图谱NanoPlot --fastq 输入数据路径.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

3、过滤（2选1）

这部分的参数选择应该参考上一步的统计结果

（1）Nanofilt

#nanofilt无法识别压缩文件，需要先解压gunzip -c 输入数据路径.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

注：输出结果需要用”>”导入文件中，否则只是输出到屏幕上

（2）Filtlong

#fitlong无需解压直接使用，此外还可以使用二代测序数据作为参考进行过滤filtlong --min_length 1000 --min_mean_q 7 输入数据路径.fastq.gz |  gzip >clean.filtlong.fq.gz

参数：         
--min_length ：最短长度

--min_mean_q：平均Q值

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

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

注：输出结果需要用”>”导入文件中，否则只是输出到屏幕上

4、比对

（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

优点主要是：将nanopore测序数据比对到参考序列上就是整个nanopore数据分析的核心，有多种比对功能，除了支持nanopre数据之外，还支持pacbio数据。比对模式可以是reads与reads之前比对，reads与基因组比对，基因组与基因组比对，以及短序列与基因组的比对，不同的比对有不同的作用

参数：

-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文件中不输出碱基质量

5、排序

#samtools处理samtools sort -@ 8 -o bam -o s0137.sorted.bam s1037.samsamtools index s0137.sorted.bamsamtools faidx ref.fq

6、拼接 

这个部分的软件五花八门，以后考虑出个测评，此处简单举两个栗子

（1）Canu

canu  -p 输出前缀 -d 输出结果目录  genomeSize=5g maxThreads=96 -nanopore-raw 原始的nanopore输入数据 >  canu.log

注意：建议使用这个软件时，数据先纠错后，再组装。

优点主要是：cano的安装特别容易，使用的最普遍，会自动进行纠错、修剪、组装。

（2）Unicycler

unicycler -1 二代数据_1.fastq.gz -2 二代数据_2.fastq.gz -l 三代数据.fastq.gz -o output_dir

优点主要是：可以同时做二代和三代组装外，对于三代组装的操作也更为简单，软件内部即可完成组装、纠错、成环等一系列步骤。现已获得上千引用。

7.拼接质量查看*

（1）Quast（此步骤可选）

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

注：该工具常用于比较多个组装软件的拼接质量

8.组装结果优化

（1）Medaka

medaka_consensus -i 原始reads数据.fastq.gz -d 组装后的assembly数据.fasta -o medaka_result -m r941_min_high_g360 -v medaka.vcf -t 24 > medaka.log

参数：

        -i 输入测序 reads

-d 需要纠错的拼接结果

        -o 输出结果文件

-m 芯片类型,需要选好

-t 线程数

（2）Pilon（该软件建议与minimap2联合用）

使用二代数据对三代数据polish

#对拼接后结果建立索引mkdir pilonln -s 拼接后数据路径/consensus.fasta assembly.fasta
#对bam进行处理samtools sort -@ 4 -O bam -o miniasm.sorted.bam miniasm.samsamtools index miniasm.sorted.bam
#利用pilon进行纠错java -Xmx32G -jar 软件路径/pilon-1.23.jar --genome assembly.fasta --fix all --changes --bam miniasm.sorted.bam --output all_pillon --outdir all_pillon --threads 24 --vcf 2> pilon.log

参数：

--frags: 表示输入的是1kb以内的paired-end文库，

--jumps表示 大于1k以上的mate pair文库,

--bam输入的bam文件

-vcf: 输出一个vcf文件，包含每个碱基的信息

--fix: Pilon将会处理的内容，基本上选snps和indels就够了

--variant: 启发式的变异检测，等价于--vcf

--fix all,breaks, 如果是polish不要使用该选项

--minmq: 用于Pilon堆叠的read最低比对质量，默认是0。

注意，pilon可以重复运行进行多次纠错

只需将pilon结果再建立索引进行输入就行

bwa index pilon.fasta

剩余代码与之前一致，只需要更改文件名即可

参考：

https://www.jianshu.com/p/1053bf88f210

https://github.com/lh3/minimap2#general

https://lh3.github.io/minimap2/minimap2.html#10

https://doi.org/10.1101/169557

https://link.zhihu.com/?target=https%3A//github.com/marbl/canu

https://link.zhihu.com/?target=https%3A//canu.readthedocs.io/en/latest/tutorial.html

https://cloud.tencent.com/developer/article/2139549

随着苗苗班的粉丝越来越多，有好几个小伙伴都在问我有没有培训或者课程。想着我们的粉丝确实很多是刚刚到生信领域的小苗苗，这个需求我之前一直倒是没注意到。我初步打算第一课就开简单实用的全基因组分析实战，大概就讲周六周日两个全天（讲不完就再加时间）。主要内容就是手把手带着大家从测序原理，环境配置，到软件安装，数据实战。不过我也不知道有多少人会报名，想征求一下大家的意见，有需求的就麻烦举个手（留言或者加微信！！），如果能凑够一个小班就开课。或者大家有什么别的想听的主题，或者意见建议都欢迎提出哦~

(1)参加培训、（2）商务合作、(3)付费分析、(4)进交流群

请加以下微信并备注目的。

这里是苗苗班，希望和大家共同交流进步，一起成长为参天有根树😊