纳米孔数据处理

背景 前面介绍了纳米孔测序的原理与碱基识别，本次带大家认识纳米孔测序数据的格式，以及怎么质控与处理。

一、原始测序数据 FAST5 文件

HDF（Hierarchical Data Format）是一种设计用于存储和组织大量数据的文件格式，是一种分级的数据文件，可以存储不同类型的图像和数码数据的文件格式，并且可以在不同类型的机器上传输，同时还有统一处理这种文件格式的函数库，最开始由美国国家超算中心研发，后来由一个非盈利组织 HDF Group 支持.HDF 支持多种商业及非商业的软件平台，包括MATLAB、Java、Python、R 和 Julia 等等，现在也提供了 Spark.其版本包括了 HDF4 和现在大量用的 HDF5。

关于 hdf5：https://www.neonscience.org/about-hdf5

nanopore 测序原始文件扩展名为 fast5，是一种经过修改的 HDF5 格式。fast5 文件由于包含内容较多，因此文件较大。通常是测序数据的 10-20 倍。

fastq 与 fast5 大小对应关系

fast5 文件层级结构示意图

HDFView 查看Fast5文件格式

二、利用 Guppy 进行 basecalling

2.1 Guppy 概述

Guppy 是 Oxford Nanopore 提供的碱基识别软件，是一款基于命令行的工具包，MinKNOW软件中内置了该工具，也可以单独安装。主要用于将 nanopore 各测序平台原始数据转换为碱基。Guppy 目前最为纳米孔官方生产环境的 basecalling 工具，可以替换之前的 Albacore工具，需要去 nanopore 官方网站注册下载。Guppy 升级很快，请选择最新版本下载使用。

guppy 主要包括三大功能：

1、碱基识别

Guppy 使用神经网络算法（RNN)将电信号数据解读为 DNA 或 RNA 的五个标准碱基，即腺嘌呤，鸟嘌呤，胞嘧啶，胸腺嘧啶和尿嘧啶。

2、拆分条码数据

添加条形码的每条序列，在序列开端和尾端均有 Oxford Nanopore Technologies 提供的条形码序列。Guppy 数据拆分功能基于识别出的碱基序列。

3、比对

用户可以提供 FASTA 或者 minimap2 index 的文件作为输入文件。比对过程基于 Oxford Nanopore Technologies 预设的参数，通过内置的 minimap2 将数据比对到参考序列。

2.2 Guppy 碱基识别

guppy 程序包主要包含三个程序，guppy_basecaller ，guppy_barcoder ，guppy_aligner 分别完成 basecalling，拆分 barcode 以及比对功能guppy_basecaller 也可以完成拆分barcode 功能，也可以先 basecalling，再使用 guppy_barcoder 进行拆分。

选择合适的操作系统版本，以及是否需要支持 GPU 计算，

利用 guppy 进行 basecalling 的操作并不难，主要输入原始测序文件夹，fast5 格式，配置文件，输出 fastq 格式文件或者转换完的 fast5，统计文件。

需要特别注意的是，一定要选择好正确的配置文件，包括芯片类型，建库方式，是否使用barcode 等信息。

basecalling
guppy_basecaller -i fast5_files/ -s output --config dna_r9.4.1_450bps_fast.cfg -r
--num_callers 4 --cpu_threads_per_caller 2 --fast5_out
合并数据
cat output/*.fastq >all.fastq
压缩数据
gzip all.fastq

必要参数：

-i 输入文件夹

-s 输出文件夹

--config 配置文件

-r 递归处理

可选参数：

--print_workflows 打印工作流程

--flowcell 芯片类型（FLO-MIN106/FLO-MIN107）

--fast5_out 指定输入 fast5 格式文件路径

--num_callers 并行处理，每个 num 产生一个 fastq 文件

--cpu_threads_per_caller 2 cpu 线程数

--compress_fastq 输出压缩格式

-q 每个文件最大条数，0 代表单独一个文件，给一个很大的值，都放在一起。

--num_callers 的数值跟生成的 fastq 文件数目一样

--records_per_fastq, per run ID or per caller。可以设成一个 caller

-x --device，GPU 计算，'auto', 或者 'cuda:<device_id>'

Guppy 结果文件：

1、FASTQ 文件

2、测序总结文档(Sequencing_summary.txt)

3、遥测信息文件(Sequence telemetry.js)一包含测序过程中的温度/电压/纳米孔状态等信息

4、Guppy 碱基识别软件日志(guppy_basecaller_1og-2019-**.**_**-**。*.log)

fastq_runid_4dd5f26403973742429ebeec317fec22caaed8da_0_0.fastq
fastq_runid_4dd5f26403973742429ebeec317fec22caaed8da_1_0.fastq
guppy_basecaller_log-2020-09-17_10-07-54.log
guppy_basecaller_log-2020-09-17_10-08-11.log
guppy_basecaller_log-2020-11-21_19-19-26.log
-rsequencing_summary.txt
sequencing_telemetry.js
workspace/

2.3 利用 Guppy 拆分条码数据

如果在建库时使用了 barcode，测序完可以通过 barcode 进行拆分样品。barcode 是在每一条序列两端添加固定序列，软件根据首位 barcode 序列不同进行拆分，类似分类彩球的过程。

由于存在连接 barcode 效率以及测序错误的存在，有些序列无法进行拆分，会被分配到unclassified 组。

均一化打分为 0-100 分的维度。大于 60 分以上默认为可信的拆分（阀值可调)

basecalling 同时拆分 barcode

guppy_basecaller -i fast5 -s fastq --config dna_r9.4.1_450bps_hac.cfg -r
--num_callers 4 --cpu_threads_per_caller 2 -x cuda:all --trim_barcodes
--barcode_kits SQK-RBK004 --min_score 10

basecalling 与拆分 barcode 分布完成

首先进行 basecalling
guppy_basecaller -i fast5 -s fastq --config dna_r9.4.1_450bps_hac.cfg -r
--num_callers 4 --cpu_threads_per_caller 2 -x cuda:all
拆分 barcode
guppy_barcoder -i fastq -s barcode --trim_barcodes --barcode_kits SQK-RBK004
--min_score 10

选项参数：

--barcode_kits ：barcodekit 名字

--min_score：barcode 比对的阈值，默认 60

--trim_barcodes：是否切除 barcode 序列，一般需要切除

输出文件:

1、目录 barcode{xx}一包含这个 barcode 号码下的.fastq 文件，每个文件包含 4000 条序列

2、目录 unclassified-包含所有分类打分低于 60 分的序列（阈值可调)

3、文件 barcoding_summary.txt-被鉴别出的 barcode，分值等，

4、文件 read_processor_Jog-YYYYMMDD_HHMMSS}.log-记录

所有文件的 barcode 分配情况以及完成情况

drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode01/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode02/
drwxr-xr-x. 2 meta bio 226 Oct 26 12:12 barcode03/
drwxr-xr-x. 2 meta bio 226 Oct 26 12:12 barcode04/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode05/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode06/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode07/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode08/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode09/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode10/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode11/
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 barcode12/
-rw-r--r--. 1 meta bio 11M Oct 26 12:12 barcoding_summary.txt
-rw-r--r--. 1 meta bio 11K Oct 26 12:12 read_processor_log-2020-10-26_12-12-25.log
drwxr-xr-x. 2 meta bio 82 Oct 26 12:12 unclassified/

三、利用 NanoPlot 进行数据质控

纳米孔测序与illumina、picbio测序的质控软件不同，不适用fastqc，会报错。主要还是Q值不同，介绍如下：

当前 nanopore 测序质量虽然有很大的改善，但准确性依然不及二代测序，例如 illumina 或者 BGIseq 等，前面介绍过，目前主流的 R9.4 芯片准确性在 92%左右。二代测序的质控标准质量值标准是 Q20，而纳米孔测序绝大部分碱基质量值都在 15 以下，质控标准为 Q7。

NanoPlot 可以用来对 nanopore 数据进行统计绘图，输入文件为 fastq 格式，绘图时需要调用 NanoStat 进行统计。NanoPlot 利用这些统计信息进行绘图，最终会生成一个网页格式文件，包括序列读长的直方图、序列读长与序列平均质量的散点图等。同时，该软件也可以对guppy 碱基识别后生成的统计文件 sequencing_summary.txt 进行绘图。

NanoPlot 绘制质控图

方法一：NanoPlot 质控
NanoPlot --fastq nanopore.fastq.gz -o nanoplot
方法二:对 guppy 生成的 sequencing_summary.txt 的绘图
NanoPlot --summary output/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

NanoPlot 质控图

还可以利用 filtlong 过滤数据

与二代测序不同，nanopre 测序为单端测序，每条数据长度不同，并且不包含通用 adapter接头，因此，数据过滤方式与 illumina 等二代测序完全不同，过滤数据主要是过滤掉一些长度过短的片段与整条测序质量比较低的片段，有些情况下如果包含 barcode，也需要对首尾进行切割。

filtlong 可以用于过滤、剪切长读长测序数据的软件。可以根据选定的质量值(Q 值）过滤，亦可根据选定的读长进行过滤。该软件同时对数据进行剪切，从读长的头或尾端截去特定数目的碱基。

https://github.com/rrwick/Filtlong

过滤数据

filtlong --min_length 2000 --min_mean_q 90 nanopore.fastq.gz | gzip >nanopore.filtlong.fq.gz
过滤完质控
NanoPlot --fastq clean.filtlong.fq.gz -o clean

--min_length ：最短长度

--min_mean_q：平均 Q 值

--keep_percent：保留最好数据的百分比

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