生信基础知识复习（二）Fastq与质控

豆豆写于18.11.28
花花昨天将我们“生信星球”的开篇语修改了，改的很棒，特此表扬！🐶
上次说了关于测序的一些个人认为的重点，这次接着来关于得到数据的下一步质控的有关知识点【这个不是全面的，只是我在第二次学习中认为比较重要的】
生信基础知识复习（一）测序

最常见的数据存储格式-fastq

不管是Illumina、Ion torrent、Pacbio都是以fastq文件进行存储，简单说几个重要点：

• 经常见到xxx.fastq.gz这是fastq的压缩文件，是二进制文件，直接打开是不对的，想要查看有两种办法：第一种：先解压缩，利用gunzip xxx.fastq.gz ，然后用less；第二种，利用zless直接查看

• fastq文件扩展名一般是fq或者fastq，以@开头（而fasta是以>开头），四行一组，分为ID信息（包括Flowcell、lane、tile、横纵坐标、index、reads编号等）、碱基序列、一个+号、质量值【我们需要的是第2、4行】

• 关于碱基与质量对应关系：我们得到了序列信息，想要知道测序质量如何，就需要进行质量转换，但是目前采用的质量体系并不是说碱基质量是20，就在fastq中对应表示为20，而是需要加上33或者64后再转换成对应的ASCII值
  这是因为： 在ASCII表中，小于32的ASCII值都是一些形状符号，32又是一个空值，所以最小也只能从33开始。
  

  

  早期的测序仪是采用Phred 64质量体系，ASCII 64是@符号，正好对应了fastq的开头符；现在基本上都转为了Phred 33的质量体系

• 那么如何快速区分Phred 64 或者33呢？

  如果放眼望去，基本上都是大写字母，那么基本可以判断时Phred 64的数据；如果包含一些#/%/@等比较怪异的符号，那么一般是Phred 33的数据

• 质量值和错误率怎么联系的呢？

  比如质量值P计算得到结果是20，那么错误率E就是百分之一，即Q20；同理P为30，即Q30，错误率为千分之一（表示该位点碱基被测1000次，出错几率为1次）

fastq质控

指标一：碱基含量分布

碱基含量就是测序数据中，ATCG的分布情况。由于测序的随机性，理论上应该A=T，G=C；另外测序数据的GC含量应该和基因组中的碱基含量一致（对于全基因组来讲），也就是说，测序是把所有cluster的第一个碱基拿出来统计质量，而这些cluster也覆盖到了全部的基因组

可能出现碱基含量分布异常：比如波动明显，但这也不能表示测序就是错误的，因为有可能测序样本是混合的或者测序一次不能达到饱和（测序一次不能覆盖到基因组的全部碱基，例如RNA-seq）；另外仅仅看一条reads的fastqc结果可能也会在Per base sequence content部分报警，因为单独一条序列不一定就准确匹配上全基因组的碱基比例

指标二：碱基质量分布

对于单个碱基，用Q数值描述；对于整体碱基，用Q20/Q30等百分比描述

Q20百分比就是：质量值大于或等于20的碱基占全部碱基的比例

一般尾部的数据质量较差是正常的

fastq处理

关于数据过滤，不要求将fastqc结果全部调成绿色/PASS，只要不影响后面分析就好；另外有时参数的设置拿不准可以先用一个大概的值，后面结果不满意可以重新过滤数据

我们知道，测序时需要在插入片段（也就是我们要测的片段）两端加很多东西：index、primer、接头，这些东西一般是提前去除的，但也不排除特殊情况。另外，还有一些判断标准

低质量数据

如果检测变异时发现某个位点发生了变异，那么是不是因为测序错误造成的“假阳性”导致的呢？因此，低质量数据一定要提前去除，否则后续分析会很麻烦。

那么怎么定义低质量呢？低于Q20的碱基比率高于30%就是低质量reads

另外，PE数据中有一条不满足要求，另外一条也要去掉

Adapter接头

一般在给到我们raw data时，基本已经去除干净，但是不排除我们自己的序列长度过短，小于PCR扩增循环（cycle）数，导致测序测多了（测完了插入片段还继续向后测，测到了adapter）；另外还可能由于adapter之间没加插入片段就自己连接，得到的测序结果就只有adapter

因此，为了确保我们的数据是真实的插入片段数据，还需要根据给定的adapter序列进行比对，如果发现了3‘端含有adapter的reads（出现频率并不高），就要把整条reads去掉

N碱基

假如100bp的reads中出现了10个以上的N碱基，即N碱基数量占reads总长10%以上（N的位置连续或者不连续），就要去除整条reads。

Duplication reads

定义：两对PE reads，其中它们的reads1完全相同，并且reads2完全相同

这里就有一个问题：这里的PE reads只是DNA片段的两头100多bp，还有好几百bp没有测，那么直接认为这两个DNA都是一致的吗？

其实是可以的，因为两个数据的两端各100多bp完全一致，这概率已经很小了，既然出现了这种情况，那么可以认为，在构建文库时将整条DNA打断后产生的两个DNA小片段是完全一致的。

另外，这也侧面反映了打断过程随机性不好，也就是说同一个位置被反复测序很多次，一般在大片段文库中存在较高的duplication。

举个例子，整条长100bp的DNA，以10bp为单位打断的随机性是不是比以50bp为单位打断的随机性要高？这50bp就相当于构建大片度文库的片段，随着文库片段越来越大，重复也是随之增长

但是，需要注意的是，在RNA-seq中，duplication是不能去除的，因为数据本身就存在基因表达丰度的差别（这不就是我们想要的么？），这种duplication并不是由于打断不随机造成的；如果去除了RNA-seq中的duplication，许多基因的丰度信息就会缺失