宏基因组分析实战（2.1）-质量控制与评估（理论介绍）

1. 前言

回顾前文：宏基因组分析实战0：前言的内容，宏基因组测序数据分析的第一步是做质量控制（Quality Control, QC）。问题来了，为什么要做质控，需要做哪些质控？这个与二代测序的原理密不可分，首先说结论，我们需要做两类质控：

• 去除接头（adapter）

• 去除低质量序列

1.1 去除接头（adapter）

提到接头（adapter），不了解的伙伴肯定一脸懵，这是什么东西，通俗来讲就是一段DNA序列，用来干嘛呢，就是连接在被打成短片段的DNA待测序列的两端，以便于将这个片段固定在测序芯片（Flow cell）上并进行扩增测序。打个比方就是两个警察站在小偷两边和他手挽手，这样小偷就跑不掉了吧，乖乖去警察局（测序芯片）接受审讯（测序）吧。

细嗦一下接头，以老大哥illumina为例，大致分为三个部分：

1. P5/P7序列（红色部分）：主要是与测序芯片上的序列互补，用于固定待测序列；

2. index2/index2（黑色部分）：又称为barcode，就像我们生活中的条形码一样，用来打标签，避免混淆不同文库；

3. Rd1 SP/Rd2 SP(橙色/蓝色部分)：就是引物序列（Primer），用于结合待测序列

知道了这些之后，文库构建的细节也就大致清晰了：

1. 首先长的DNA片段被打断成小片段，然后将5’端磷酸化，3‘端加A尾（下图A）

2. 连接接头和DNA Insert，形成长Y接头（下图B）

3. 最后形成的就是完成后的文库了（下图C），后续就结合到测序芯片开始测序（下图D），便是经典的桥式PCR，边合成边测序（Sequencing by Synthesis, SBS）
   图片来源：

https://support.illumina.com/ko-kr/bulletins/2020/12/how-short-inserts-affect-sequencing-performance.html

有了这些知识后，我们就可以明白为什么要对下机数据去除接头了，因为接头是为了测序而额外增加的序列，不属于我们想要得到的目的序列，会对下游分析造成干扰，因此需要去除！

💡 这里的接头是指除了插入序列外的其他序列。也有将P5/P7序列定义为接头，区分了barcode和primer（如下文中的Trimmomatic）

1.2 去除低质量序列

这一类质控更好理解，其实就是根据测序质量来修剪序列，那么什么是测序质量呢，它是如何表征的呢，我们得先了解fastq文件的格式，这种文件每一条序列（sequence）的信息共分为四行，分别为

• 序列ID

• 碱基信息

• 标识符号“+”

• 碱基质量

以ERR1018205_1.fastq.gz为例，第一条序列的信息如下：

    @ERR1018205.1 FCD0JMKACXX:2:1101:1094:3418#TTAGGCAT/1
    CGTTGTTCACATATCCAATCTTCTGCAAAGGCTCGTAGTTCAGGATGTTGAAACCAAGCATGTTTTTCCAGCACAAATAGGATCGAAAACAGCCATTCAC
    +
    ???;:BDD>DBDDDA:EEEEI9AF+A<CB<3ABC8C7?:?9?DBBDDB8DDIICCC8(58)88=)7-@AC3)=;>.6;;6(..;>?>;?><=?(5:>>DA
    
    
    
      
      
      
      
      
    

我们关注的质量就在第四行 ，它是与第二行的碱基一一对应的，代表了每个碱基的测序质量。那有的同学可能就会问了，不是质量吗，为什么会是这一串奇奇怪怪的字符呢 ，该如何转换成我能看懂的质量表示形式呢。其实测序质量也是和二代测序的原理密不可分的，测序仪是通过捕捉带荧光标记碱基的荧光信号来实现对碱基的读取，根据捕捉到的信号强弱来判定读取碱基的可靠度（如下图），其结果用错误概率（error probability）来表示，进而将这个值通过取对数转化为质量得分（Quality Scores），就是现在比较通用的质量。其中具体转化公式如下：

其中P(A)就是指碱基A的错误概率。进而我们可以得到比较常用的质量得分和错误率的关系：

其中，如果这个碱基的质量得分为20，就表示错误概率为1%，换言之正确概率为99%；而若为30，则正确概率为99.9%。现在二代测序的测序质量基本都在30以上，所以精确度其实非常高。
图片来源：

https://support.illumina.com/content/dam/illumina-support/documents/documentation/system_documentation/hiseq2000/hiseq-2000-system-guide-15011190-03.pdf

可能有些伙伴还会问，这也和上面的字符串没关系呀 ，别急，关系这就来。由于质量分数的位数不固定，可能有个位数，也可能有十位数，不方便与碱基做对应。因此测序公司 还会将这个质量得分加一个质量体系（33或者64，现在更多为33），称为Phred ，然后根据ASCII码 进行转换，就得到测序数据第四行的结果啦

我们以第一个字符”?“来反推实践一下，首先下载一张ASCII表（如下图），其中Dec表示十进制数字，这里对应质量得分；Chr表示字符，这里对应fastq文件中第四行。通过查找，我们发现？对应的数字是63，由于这个是质量得分加33或64来的，只要减去就可以得到对应的质量得分。当然，这里64就明显不适用了，减完成负数了。也表明本次测序数据是使用33的质量体系（Phred33），减去33后得到质量得分为30，所以第一个碱基C的测序质量为30
图片来源：https://www.asciitable.com/

铺垫了这么多，相信大家也明白为什么要去除低质量的碱基（base）或者序列（sequence）了吧，因为低质量代表高错误率，会对后续的分析带来影响，所以必须去除！

一般来讲，我们把A、T、C、G称为碱基（base），把由碱基组成的一条片段称为序列（sequence）

💡 想了解更多关于fastq文件的知识，可以参考illumina官网的介绍~

本文由博客一文多发平台 OpenWrite 发布！