0基础学习基因组三代测序组装（1）——数据质控

当第二代高通量测序技术进入成熟阶段后，读长过短、PCR扩增带来的偏向性等问题开始日益凸显；作为基因组学上新的转折点，以PacBio单分子实时测序技术及纳米孔单分子测序技术为首的第三代高通量测序技术（Third-generation Sequencing）开始进入科研应用，从单分子水平上对DNA分子的实时测序，成功解决了二代测序几大困扰：极端 GC含量区域覆盖度低、高度重复区域无法较好地拼装、大片段变异难以准确检测、不能直接检测碱基修饰等问题。

ONT（Oxford Nanopore Technologies）牛津纳米孔测序技术作为第三代单分子实时测序技术，其原理是基于高分子膜两侧电压和其中的蛋白质纳米孔，当单分子DNA从纳米孔通过时，会引起孔两侧电位差来实现信号检测， 而ATCG四种碱基的带电性质不一样，因此利用电信号的差异就能检测出通过的碱基类型，从而实现测序。

Nanopore商业化平台有三个：MinION、GridION及PromethION。本系列笔记的三代测序数据来源于PromethION测序平台测序的一个cell，PromethION测序仪拥有48个流动槽，每个流动槽拥有3000个纳米孔通道（总计144000个），适用于大样本量的高通量快速测序。

basecalling

在ONT的测序平台中，将通过纳米孔的DNA或RNA链产生的电位信号转化为相应的碱基序列的过程，称为basecalling。Nanopore测序的下机数据的原始数据格式为包含原始测序电信号的fast5格式，官方有提供工具Guppy进行basecalling，以mean_qscore_template的数值大于等于7为标准（也就是测序质量大于7的reads）得到原始测序数据，这样得到的basecalling数据为fastq格式（.fastq或者.fq结尾），所以我拿到的就是已经basecalling后的结果。

• fast5: 原始电信号文件，以.fast5为文件结尾。此文件既有测序得到的序列信息，还有甲基化修饰信息（甲基化位点电信号会不一样）。
• fastq: fast5文件转换而来，四行一个单位，序列和碱基质量一一对应。

basecalling的同时还可以一起拆分barcode条码序列，这里我没用到guppy这个软件，了解一下就行。经过basecalling后，文件会分为fail和pass两部分，pass部分就是满足Q值>7的序列（二代测序质控标准是Q20，这里的三代测序质控标准是Q7，准确性不及二代测序）。

还有一个summary.txt文件，这是一个测序汇总文件，结构如下：

第一眼看上去很乱，几个重要的列含义如下：

1. filename_fastq fasq文件名
2. filename_fast5 fast5文件名
3. read_id 每条read对应的id号
4. run_id 这一次运行产生的id号，一个flowcell通常为一个run
5. channel mux 该条read在哪个channel测的
6. start_time 这条read测序起始时间
7. duration 这条read测序经过时间
8. passes_filtering Q值大于7为TRUE否则为FALSE
9. sequence_length_template read长度，三代测序数据过滤的指标之一
10. mean_qscore_template 非常重要的指标，每一个read的平均Q值
11. 有关barcode的都是标签序列相关参数，因为不同样品接头会添加不同的标签序列，混测的时候根据标签序列与样品的对应关系来区分不同样品。

返回的数据是guppy处理过的，也就是raw reads，接下来质控的过程就需要自己动手了。

nanoplot质控

先说明下为什么要用这个工具，三代测序的数据读长比二代测序长很多，而且每条序列的长度都是不一样的。不能用之前转录组数据分析中的fastq工具，会报错，因此使用nanoplot工具来生成质检报告，同样也是会生成各种html文件方便浏览结果。

先创建一个nanoplot专用的环境，下载nanoplot，之后的质控过程都在这个环境下进行。

1
2
3
4
5
6
7
8

# 创建环境并下载nanoplot
conda create -n nanoplot -y -c bioconda nanoplot
# 激活环境
. activate nanoplot
# 生成质检报告。可以用pass.fq文件，也可以直接用summery.txt文件。-o参数后面是输出文件夹名称。
NanoPlot --summary summary.txt --loglength -o summary-plots-log-transformed
NanoPlot -t 4 --fastq pass.fq.gz --plots hex dot -o nanoplot_out
# 详细参数设置可以在NanoPlot --help中查看

运行结束之后会生成summary-plots-log-transformed这个文件夹，我们可以用xftp工具查看里面的html结果文件，也可以挑取一些数据做数据质量统计表。

放一张原始测序数据读长分布图示意一下：

点击这里查看用summary.txt生成的质控报告

点击这里查看用pass.fq.gz生成的质控报告（推荐用这个）

两种方法生成的质控报告略微有点差别，因为summary.txt文件中记录了所有序列，可以看到有部分序列质量在Q5-Q7之间；而pass.fg.gz生成的质控报告中，所有序列的质量都在Q7及以上。后续以分析pass.fq.gz文件生成的质控报告为准，对这个文件序列的长度进行过滤。

如果不需要图，只需要知道有多少条reads、reads平均长度、N50、N90这些数据做表格的话，还有一个比较实用的perl脚本，怎么使用就不赘述了，源代码放底下参考。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

#/usr/bin/perl -w
use strict;
use warnings;

### *fastq.gz 数据统计 N50 N90 num_seqs sum_len min_len avg_len max_len 
### usage: perl stat.fastq.gz.N50.N90.pl *.fastq.gz

my $fastq_gz = $ARGV[0];
open(IN,"gzip -dc $fastq_gz|") or die ("can not open $fastq_gz\n");
open(OUT,">$fastq_gz.stat");
print OUT "name num_seqs    sum_len min_len avg_len max_len N50 N90\n";
print OUT "$fastq_gz\t";

my ($len,$total,$num_seqs,$min_len,$max_len)=(0,0,0,0,0);
my @length_list;

while(<IN>){
    my $title = $_;
    my $seq = <IN>;
    my $add = <IN>;
    my $quality = <IN>;
    $seq =~ s/\r|\n|\r\n//mg;
    $len = length($seq);
    if($len>0){
        $total += $len;
        push @length_list,$len;
    }
    if($min_len == 0){$min_len = $len;}elsif($min_len > $len){$min_len = $len;}
    if($max_len == 0){$max_len = $len;}elsif($max_len < $len){$max_len = $len;} 
    $len=0;
    $num_seqs++;
}

my $avg_len = $total/$num_seqs;
print OUT "$num_seqs\t";
print OUT "$total\t";
print OUT "$min_len\t";
print OUT "$avg_len\t";
print OUT "$max_len\t";

@length_list=sort{$b<=>$a} @length_list;

my ($count,$half)=(0,0);
for (my $j=0;$j<@length_list;$j++){
    $count+=$length_list[$j];
    if (($count>=$total/2)&&($half==0)){
        print OUT "$length_list[$j]\t";
        $half=$length_list[$j]
    }elsif ($count>=$total*0.9){
        print OUT "$length_list[$j]\t\n";
        exit;
    }
}

filtlong过滤数据

前面说过，二代测序是双端测序，三代测序是单端测序，两者过滤数据的要求不同。三代测序主要是过滤长度过短的序列和测序质量较低的序列。在basecalling中我们过滤了Q值小于7的序列，现在还要过滤read长度小于1000bp的序列。过滤后的序列可以直接用于后续的组装。

1
2
3
4

# 安装filtlong软件
conda install -y filtlong
# 设置序列最短为1000bp，压缩结果文件到新文件中
filtlong --min_length 1000 pass.fq.gz | gzip > clean_filter.fq.gz

可以看到过滤了一部分数据，用过滤后的数据再跑一次NanoPlot

1

NanoPlot -t 4 --fastq clean_filter.fq.gz --plots hex dot -o filt_nanoplot_out

测序数据读长分布如下，可以看到已经没有1kb以下的reads了：

点击这里查看过滤后的质控报告

至此质控过滤流程结束，我们可以做一个下机数据质控统计表：