二倍体线虫基因组组装实战（1）——数据质控

​ 最近拿到了一个线虫基因组（秀丽隐杆线虫，Caenorhabditis elegans）的二代HIFI测序数据和三代ONT测序数据。网上查询出来线虫应该是个二倍体，打算拿这个数据整体从0基础开始进行多倍体的一个基因组组装。

1.技术背景

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

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

HiFi测序数据（High Fidelity Sequencing）是基于PacBio的SMRT（Single Molecule, Real-Time）测序技术的一种高质量长读长测序数据。单分子实时测序：HiFi测序技术通过对单个DNA分子进行实时测序，使用高保真（HiFi）算法来提高测序的准确性。多重读取：在进行HiFi测序时，同一分子可以被多次读取，这种重复测序的方式大大提高了准确性。通过对相同DNA片段的多次读取，HiFi数据能有效纠正由于测序过程中的错误，从而产生高保真的序列数据。HiFi测序的读长通常在10-20 kb之间，适合对长重复序列的测序和组装。

HiFi测序的准确性通常低于1%的错误率，这使得它在进行变异检测（如SNP和插入/缺失）时非常可靠。通过使用ONT长读长数据和HIFI数据，可以尽可能的保证基因组组装的质量。

2.ONT数据质控

basecalling

在ONT（Oxford Nanopore Technologies）测序平台上，basecalling 是将通过纳米孔的 DNA 或 RNA 链产生的电位信号转化为碱基序列的过程。其测序原始数据格式为 fast5，包含了原始电信号信息。通过官方工具 Guppy 进行 basecalling 时，常以 mean_qscore_template ≥ 7（Q值大于7）作为质量标准，筛选出较高质量的测序数据。经过 basecalling，输出数据转换为 fastq 格式，每条序列由四行组成，包含序列和相应碱基质量。

• fast5 文件：包含原始电信号，记录了测序的序列信息以及甲基化修饰信息（因甲基化会影响电信号）。
• fastq 文件：从 fast5 转换而来，常以 .fastq 或 .fq 结尾，包含筛选后的碱基序列和质量信息。

在 basecalling 时还可选择拆分条码（barcode）序列，但此流程未使用 Guppy 的该功能。basecalling 结果会分为 fail 和 pass 两部分，满足 Q值大于7的序列放入 pass 文件夹（相比二代测序的 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质控

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

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

#进入目标目录
cd qc.Cae_elegans
# 创建环境并下载nanoplot
conda create -n nanoplot -y -c bioconda nanoplot
# 激活环境
. activate nanoplot
# 生成质检报告。可以用.fq文件，也可以直接用summery.txt文件。-o参数后面是输出文件夹名称。
NanoPlot --fastq ./ont.Cae_elegans.fq --loglength --legacy hex --plots dot -o ./qc.ont.raw.report -t 80
# 详细参数设置可以在NanoPlot --help中查看

常用的 NanoPlot 命令功能(详细可阅读官方文档):

1. 基本质量分析：
   • 使用 --fastq 或 --bam 指定 FASTQ 或 BAM 文件。
   • 结合 --loglength 对读长进行对数变换，查看数据的分布特征。
   • 常见输出：读长、质量分布图。
2. 生成详细统计信息：
   • 添加 --N50、--qualmean 等选项计算和显示 N50、平均质量等指标。
   • 使用 --onlystats 仅输出统计信息，而不生成可视化图表。
3. 多种图表类型组合：
   • 使用 --plots 指定所需的图表类型（如 hex、dot、cumulative、box 等），可组合生成不同的图表以便分析。
4. 数据过滤：
   • 设置 --minlength 和 --minqual 过滤低质量或短读长的数据，以便生成更有代表性的统计结果。
   • –minlength 1000：仅包含读长大于1000 bp的序列。
   • –minqual 7：仅包含质量值大于7的序列。
5. 生成交互式 HTML 报告：
   • 使用 --interactive 生成交互式 HTML 报告，方便在浏览器中查看数据细节。
6. 图像定制：
   • 使用 --figsize 设置输出图像的尺寸，以适应特定的展示需求。

运行结束之后会生成ont.Cae_elegans_combined_plots这个文件夹，我们可以查看里面的html结果文件，也可以挑取一些数据做数据质量统计表。

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

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

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

如果不需要图，只需要知道有多少条reads、reads平均长度、N50、N90这些数据做表格的话，还有一个比较实用的perl脚本，以供参考。

#/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的序列。过滤后的序列可以直接用于后续的组装。

过滤掉 1000 bp 以下的短读长可以带来以下几个好处：

1. 提高组装质量：短读长（<1000 bp）往往难以与长读长片段对齐或组装，可能导致错配或组装碎片化。过滤掉这些短片段有助于提升组装结果的连贯性和准确性。
2. 减少噪音：在 ONT 数据中，短片段的质量通常不如长片段稳定，且可能来自于测序过程中的噪音或污染。去除短片段可以减少背景噪声，提升数据质量。
3. 优化计算资源：短读长会增加下游分析的计算量，但对完整组装或结构变异检测的贡献较少。去除这些片段有助于降低计算资源的消耗，使分析更高效。
4. 适用于特定分析需求：在结构变异分析或大型基因组 de novo 组装中，长读长（>1000 bp）的数据能更好地覆盖复杂区域，短读长的贡献有限。对于这些研究需求，保留长读长数据能够更专注于有用信息。

#进入目标目录
cd ..
cd Cae_elegans_data

# 安装filtlong软件
conda install -y filtlong

# 设置序列最短为1000bp，结果文件到新文件中
filtlong --min_length 1000 ont.Cae_elegans.fq > ont.filter1000-.fq

用过滤后的数据再跑一次NanoPlot

#用过滤后的数据再跑一次NanoPlot,详细参数设置可以在NanoPlot --help中查看 
. activate nanoplot
cd ../qc.Cae_elegans/qc.ont.clean.filter_1000_reoprt

NanoPlot -t 120 --fastq ../Cae_elegans_data/ont.filter1000-.fq --loglength --legacy hex --plots dot -o qc.ont.filter1000-.reoprt

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