数据预处理综述

数据预处理综述

摘要：当今社会生物信息学已成为整个生命科学发展的重要组成部分，成为生命科学研究的前沿。随着测序技术的不断进步，获取基因序列的时间不断缩短，测序分析中的关键步骤之一的数据预处理也变得尤为重要。本文对基因测序的主要两种方法，数据预处理的概念及方法等方面进行了论述。随着技术的不断革新我们对生物信息学的掌握将更加深入更加灵活，数据预处理技术的要求也越来越高，它在功能基因的准确发现与识别、基因与蛋白质的表达与功能研究方面都将发挥关键的作用。

关键词：sanger测序法，Illumina，Sequencing by Synthesis ，FASTQC，Trimmomatic

1 主要的测序方法

重点描述sanger法和以Illumina/Solexa Genome Analyzer 的测序。

Sanger法是根据核苷酸在某一固定的点开始，随机在某一个特定的碱基处终止，并且在每个碱基后面进行荧光标记，产生以A、T、C、G结束的四组不同长度的一系列核苷酸，然后在尿素变性的PAGE胶上电泳进行检测，从而获得可见的DNA碱基序列。

原理:是利用一种DNA聚合酶来延伸结合在待定序列模板上的引物。直到掺入一种链终止核苷酸为止。每一次序列测定由一套四个单独的反应构成，每个反应含有所有四种脱氧核苷酸三磷酸(dNTP)，并混入限量的一种不同的双脱氧核苷三磷酸(ddNTP)。由于ddNTP缺乏延伸所需要的3-OH基团，使延长的寡聚核苷酸选择性地在G、A、T或C处终止。终止点由反应中相应的双脱氧而定。每一种dNTPs和ddNTPs的相对浓度可以调整，使反应得到一组长几百至几千碱基的链终止产物。它们具有共同的起始点，但终止在不同的的核苷酸上，可通过高分辨率变性凝胶电泳分离大小不同的片段，凝胶处理后可用X-光胶片放射自显影或非同位素标记进行检测。DNA的复制需要：DNA聚合酶，双链DNA模板，带有3'-OH末端的单链寡核苷酸引物，4种dNTP（dATP、dGTP、dTTP和dCTP）。聚合酶用模板作指导，不断地将dNTP加到引物的3'-OH末端，使引物延伸，合成出新的互补DNA链。如果加入一种特殊核苷酸，双脱氧核苷三磷酸（ddNTP），因它在脱氧核糖的3’位置缺少一个羟基，故不能同后续的dNTP形成磷酸二酯键。如，存在ddCTP、dCTP和三种其他的dNTP（其中一种为α-32P标记）的情况下，将引物、模板和DNA聚合酶一起保温，即可形成一种全部具有相同的5'-引物端和以ddC残基为3’端结尾的一系列长短不一片段的混合物。经变性聚丙烯酰胺凝胶电泳分离制得的放射性自显影区带图谱将为新合成的不同长度的DNA链中C的分布提供准确信息，从而将全部C的位置确定下来。类似的方法，在ddATP、ddGTP和ddTTP存在的条件下，可同时制得分别以ddA、ddG和ddT残基为3‘端结尾的三组长短不一的片段。将制得的四组混合物平行地点加在变性聚丙烯酰胺凝胶电泳板上进行电泳，每组制品中的各个组分将按其链长的不同得到分离，制得相应的放射性自显影图谱。从所得图谱即可直接读得DNA的碱基序列。与DNA复制不同的是sanger测序中的引物是单引物或者是单链。

第二代DNA序列测序技术（以Illumina/Solexa Genome Analyzer 测序为例）

核心思想：边合成边测序（Sequencing by Synthesis)，即通过捕捉新合成的末端的标记来确定DNA的序列

基本原理：Illumina/Solexa Genome Analyzer测序的基本原理是边合成边测序。在Sanger 等测序方法的基础上，通过技术创新，用不同颜色的荧光标记四种不同的dNTP，当DNA聚合酶合成互补链时，每添加一种dNTP就会释放出不同的荧光，根据捕捉的荧光信号并经过特定的计算机软件处理，从而获得待测DNA的序列信息。

操作流程：

1）测序文库的构建（Library Construction）:首先准备基因组DNA（虽然测序公司

要求样品量要达到200ng，但是Gnome Analyzer系统所需的样品量可低至100ng,能应用在很多样品有限的实验中），然后将DNA随机片段化成几百碱基或更短的小片段，并在两头加上特定的接头（Adaptor）。如果是转录组测序，则文库的构建要相对麻烦些，RNA片段化之后需反转成cDNA，然后加上接头，或者先将RNA反转成cDNA，然后再片段化并加上接头。片段的大小（Insert size）对于后面的数据分析有影响，可根据需要来选择。对于基因组测序来说，通常会选择几种不同的insert size，以便在组装（Assembly）的时候获得更多信息。

2）锚定桥接（Surface Attachment and Bridge Amplification）:Solexa测序的反应在叫做flow cell的玻璃管中进行，flow cell又被细分成8个Lane，每个Lane的内表面有无数的被固定的单链接头。上述步骤得到的带接头的DNA 片段变性成单链后与测序通道上的接头引物结合形成桥状结构以供后续的预扩增使用。

3）预扩增（Denaturation and Complete Amplification）:添加未标记的dNTP 和普通Taq 酶进行固相桥式PCR 扩增，单链桥型待测片段被扩增成为双链桥型片段。通过变性，释放出互补的单链，锚定到附近的固相表面。通过不断循环，将会在Flow cell 的固相表面上获得上百万条成簇分布的双链待测片段。

4）单碱基延伸测序（Single Base Extension and Sequencing）:在测序的flow cell中加入四种荧光标记的dNTP 、DNA聚合酶以及接头引物进行扩增，在每一个测序簇延伸互补链时，每加入一个被荧光标记的dNTP就能释放出相对应的荧光，测序仪通过捕获荧光信号，并通过计算机软件将光信号转化为测序峰，从而获得待测片段的序列信息。从荧光信号获取待测片段的序列信息的过程叫做Base Calling，Illumina公司Base Calling所用的软件是Illumina’s Genome Analyzer Sequencing Control Software andPipeline Analysis Software。读长会受到多个引起信号衰减的因素所影响，如荧光标记的不完全切割。随着读长的增加，错误率也会随之上升。

5）数据分析（Data Analyzing）:这一步严格来讲不能算作测序操作流程的一部分，但是只有通过这一步前面的工作才显得有意义。测序得到的原始数据是长度只有几十个碱基的序列，要通过生物信息学工具将这些短的序列组装成长的Contigs甚至是整个基因组的框架，或者把这些序列比对到已有的基因组或者相近物种基因组序列上，并进一步分析得到有生物学意义的结果。

2 数据预处理的步骤及方法：

1)Fastqc

当二代测序的原始数据拿到手之后，第一步要做的就是看一看原始reads的质量。

常用的工具就是fastqc我们在服务器上用命令行来运行fastqc：

fastqc [-o output dir] [--(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 ..

seqfileN

-o用来指定输出文件的所在目录，注意是不能自动新建目录的。输出的结果是.zip 文件，默认自动解压缩，命令里加上--noextract则不解压缩。-f用来强制指定输入文件格式，默认会自动检测。-c用来指定一个contaminant文件，fastqc会把overrepresented sequences往这个contaminant文件里搜索。contaminant文件的格式是"Name\tSequences"，#开头的行是注释。加上-q 会进入沉默模式，即不出现下面的提示：

Started analysis of target.fq

Approx 5% complete for target.fq

Approx 10% complete for target.fq

如果输入的fastq文件名是target.fq，fastqc的输出的压缩文件将是target.fq_fastqc.zip。

解压后，查看html格式的结果报告。结果分为如下几项：

结果分为绿色的"PASS"，黄色的"WARN"和红色的"FAIL"。其中各项的意义如下：

1 Basic statistics

2Per base sequence quality：quality就是Fred值，-10*log10(p)，p为测错的概率。所以一条reads某位置出错概率为0.01时，其quality就是20。横轴代表位置，纵轴quality。红色表

示中位数，黄色是25%-75%区间，触须是10%-90%区间，蓝线是平均数。

若任一位置的下四分位数低于10或中位数低于25，报"WARN"；若任一位置的下四分位数低于5或中位数低于20，报"FAIL".

3 Per Sequence Quality Scores：每条reads的quality的均值的分布，横轴为quality，纵轴是reads数目。当出现上图的情况时，我们就会知道有一部分reads具有比较差的质量。

当峰值小于27（错误率0.2%）时报"WARN"，当峰值小于20（错误率1%）时报"FAIL"。

4 Per Base Sequence Content:对所有reads的每一个位置，统计ATCG四种碱基（正常情况）

的分布横轴为位置，纵轴为百分比。正常情况下四种碱基的出现频率应该是接近的，而且没有位置差异。因此好的样本中四条线应该平行且接近。当部分位置碱基的比例出现bias 时，即四条线在某些位置纷乱交织，往往提示我们有overrepresented sequence的污染。当所有位置的碱基比例一致的表现出bias时，即四条线平行但分开，往往代表文库有bias (建库过程或本身特点)，或者是测序中的系统误差。当任一位置的A/T比例与G/C比例相差超过10%，报"WARN"；当任一位置的A/T比例与G/C比例相差超过20%，报"FAIL"。

5 Per Base GC Content:对所有reads的每个位置，统计GC含量。如果建库足够均匀，reads 的每个位置应当是没有差异的，所以GC含量的线应当平行于X轴，反映样品（基因组、转录组等）的GC含量。当部分位置GC含量出现bias时，往往提示我们有overrepresented sequence的污染。当所有位置的GC含量一致的表现出bias时，往往代表文库有bias (建库过程或本身特点)，或者是测序中的系统误差。当任一位置的GC含量偏离均值的5%时，报"WARN"；当任一位置的GC含量偏离均值的10%时，报"FAIL"。

6 Per Sequence GC Content：统计reads的平均GC含量的分布。红线是实际情况，蓝线是理论分布（正态分布，均值不一定在50%，而是由平均GC含量推断的）。曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差（overrepresented reads）。形状接近

正态但偏离理论分布的情况提示我们可能有系统偏差。偏离理论分布的reads超过15%时，报"WARN"；偏离理论分布的reads超过30%时，报"FAIL"。

7 Per Base N Content：当测序仪器不能辨别某条reads的某个位置到底是什么碱基时，就会产生“N”。对所有reads的每个位置，统计N的比率。正常情况下N的比例是很小的，所以图上常常看到一条直线，但放大Y轴之后会发现还是有N的存在，这不算问题。当Y轴在0%-100%的范围内也能看到“鼓包”时，说明测序系统出了问题。当任意位置的N的比例超过5%，报"WARN"；当任意位置的N的比例超过20%，报"FAIL"。

8 Sequence Length Distribution :reads长度的分布。当reads长度不一致时报"WARN"；当有长

度为0的read时报“FAIL”。

9 Duplicate Sequences:统计序列完全一样的reads的频率。测序深度越高，越容易产生一定程度的duplication，这是正常的现象，但如果duplication的程度很高，就提示我们可能有bias的存在（如建库过程中的PCR duplication）。横坐标是duplication的次数，纵坐标是duplicated reads的数目，以unique reads的总数作为100%。上图的情况中，相当于unique reads数目～20%的reads是观察到两个重复的，～7%是观察到三次重复的，依此类推。

可以想象，如果原始数据很大（事实往往如此），做这样的统计将非常慢，所以fastqc中用fq数据的前200,000条reads统计其在全部数据中的重复情况。重复数目大于等于10的reads 被合并统计，这也是为什么我们看到上图的最右侧略有上扬。大于75bp的reads只取50bp （不知道怎么选的）进行比较。但由于reads越长越不容易完全相同（由测序错误导致），所以其重复程度仍有可能被低估。当非unique的reads占总数的比例大于20%时，报"WARN"；当非unique的reads占总数的比例大于50%时，报"FAIL“。

10 Overrepresented Sequences :如果有某个序列大量出现，就叫做over-represented。fastqc

的标准是占全部reads的0.1%以上。和上面的duplicate analysis一样，为了计算方便，只取了fq数据的前200,000条reads进行统计，所以有可能over-represented reads不在里面。而且大于75bp的reads也是只取50bp。如果命令行中加入了-c contaminant file，出现的over-represented sequence会从contaminant_file里面找匹配的hit（至少20bp且最多一个mismatch），可以给我们一些线索。

当发现超过总reads数0.1%的reads时报”WARN“，当发现超过总reads数1%的reads时报”FAIL“。

11 Overrepresented Kmers:如果某k个bp的短序列在reads中大量出现，其频率高于统计期望的话，fastqc将其记为over-represented k-mer。默认的k = 5，可以用-k --kmers选项来调节，范围是2-10。出现频率总体上3倍于期望或是在某位置上5倍于期望的k-mer被认为是over-represented。fastqc除了列出所有over-represented k-mers，还会把前6个的per base distribution画出来。当有出现频率总体上3倍于期望或是在某位置上5倍于期望的k-mer 时，报”WARN“；当有出现频率在某位置上10倍于期望的k-mer时报"FAIL"。

2）2222222

2）Trimmomatic

它是一个针对Illumina高通量测序的reads trim的工具。即能够针对paired-end 也能弄single ended.它能够利用FASTQ文件（phred + 33 或者是phred + 64 碱基质量格式，取决于Illumina 测序的机器）.对于single-ended，一个输入文件和一个输出文件，加上参数。对于paired-end 数据，两个输入文件，4个输出文件，分别为2个是'paired'，2个是'unpaired'（一个为forward 的,一个为reverse的）。Trimmomatic用两种策略来去除adapter: Palindrome and Simple simple trimming是利用每一个adapter序列去跟reads匹配，如果匹配上，就删除read的这部分序列。

Palindrome trimming 是在adapter序列中reading through一个短的片段。在这个方法中，

相应的adapter序列在硅片上的连接被连接到了reads的开始，形成了adapter+read序列，正链和反链比对，如果比对显示read-through，正链剪切掉adapter，反链去掉（由于它没有包含新的数据）Trimmomatic功能强大包括如下功能：

ILLUMINACLIP: Cut adapter and other illumina-specific sequences from the read. SLIDINGWINDOW: Perform a sliding window trimming, cutting once the average quality within the window falls below a threshold.

LEADING: Cut bases off the start of a read, if below a threshold quality

TRAILING: Cut bases off the end of a read, if below a threshold quality

CROP: Cut the read to a specified length

HEADCROP: Cut the specified number of bases from the start of the read

MINLEN: Drop the read if it is below a specified length

TOPHRED33: Convert quality scores to Phred-33将碱基质量转换为pred33格式

TOPHRED64: Convert quality scores to Phred-64 将碱基质量转换为pred64格式AVGQUAL:丢掉质量低于这个值的reads