肿瘤基因组分析教程：二、数据质控

本教程目的是对数据质控有初步理解，并熟悉以下工具的使用：FastQC，Skewer，FASTX-Toolkit。参考教程：Quality Control

在高通量测序过程中，很难避免引入测序错误，如碱基读取错误（base calling errors）、插入或缺失（small insertions/deletions），低质量碱基（poor quality reads）、引物与接头污染（primer/adapter contamination），以及最常见的碱基替换（substitution），错误发生频率约 0.5-2.0%，且多发生于3’区域。

因此，在数据分析之前，有必要对原始数据进行质量控制，如去除低质量碱基（trimming off low quality bases）、去除接头（cleaning up any sequencing adapters）、去除PCR重复（removing PCR duplicates），另外，在深度测序（> 15X）时，部分测序错误可以通过算法被修正。

一、质量值编码方案（Quality Value Encoding Schema）

质量值是测序中非常重要的参考标准，使用Phred编码方案，质量编码体现在FASTQ格式的第四行，如

@HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1
TTAATTGGTAAATAAATCTCCTAATAGCTTAGATNTTACCTTNNNNNNNNNNTAGTTTCTTGAGATTTGTTGGGGGAGACATTTTTGTGATTGCCTTGAT
+HWI-EAS209_0006_FC706VJ:5:58:5894:21141#ATCACG/1
efcfffffcfeefffcffffffddf`feed]`]_Ba_^__[YBBBBBBBBBBRTT\]][]dddd`ddd^dddadd^BBBBBBBBBBBBBBBBBBBBBBBB

碱基质量可以通过查表得出对应的Q-score，Q(A)=−10log10(P(A))，A表示错误概率，如下表所示:

二、FASTQ

现在我们的目录下有两个fastq文件，sam01_1.fastq & sam01_2.fastq，其实是来自同一个样本，双端测序生成的两个文件，我们只对其中一个检查质量

$ fastqc -f fastq sam01_1.fastq
Started analysis of sam01_1.fastq
Approx 5% complete for sam01_1.fastq
Approx 10% complete for sam01_1.fastq
Approx 15% complete for sam01_1.fastq
Approx 20% complete for sam01_1.fastq
Approx 25% complete for sam01_1.fastq
Approx 30% complete for sam01_1.fastq
Approx 35% complete for sam01_1.fastq
Approx 40% complete for sam01_1.fastq
Approx 45% complete for sam01_1.fastq
Approx 50% complete for sam01_1.fastq
Approx 55% complete for sam01_1.fastq
Approx 60% complete for sam01_1.fastq
Approx 65% complete for sam01_1.fastq
Approx 70% complete for sam01_1.fastq
Approx 75% complete for sam01_1.fastq
Approx 80% complete for sam01_1.fastq
Approx 85% complete for sam01_1.fastq
Approx 90% complete for sam01_1.fastq
Approx 95% complete for sam01_1.fastq
Analysis complete for sam01_1.fastq

# 运行完成，会生成两个文件：
sam01_1_fastqc.html # 在浏览器中查看
sam01_1_fastqc.zip # 存放结果的压缩文件

# 在浏览器中查看HTML文件，第一个表格是基本信息，可以看到平台、测序总长、读长，GC含量等信息

每个碱基的平均质量Q-score，也即前面提到的Phred值，这个测序结果质量相当好，平均质量在34以上，只有尾部稍差，也有30以上。

序列对应的质量分数，可以看到，绝大多数序列的Q值都很高，平均值约36

如果测序质量不好，则可以通过Read Trimming来去除低质量的读段。下面这张图显示的是读长，基本上是151 bp

可以看到有重复，后面会去掉

3’端有ILLUMINA的通用接头，红色的线

三、Read Trimming

$ conda install -c bioconda skewer # 安装skewer

# 得益于测序技术的进步，现在的测序数据大都很好，这一步可根据需要确定是否运行

$ skewer -t 4 -l 100 -q 30 -Q 30 -m pe -o trim sam01_1.fastq sam01_2.fastq
.--. .-.
: .--': :.-.
`. `. : `'.' .--. .-..-..-. .--. .--.
_`, :: . `.' '_.': `; `; :' '_.': ..'
`.__.':_;:_;`.__.'`.__.__.'`.__.':_;
skewer v0.2.2 [April 4, 2016] # 程序版本，上面这个图别出心裁，不过不仔细看还真看不出来
Parameters used: # 参数
-- 3' end adapter sequence (-x):        AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC # 接头序列
-- paired 3' end adapter sequence (-y): AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA
-- maximum error ratio allowed (-r):    0.100
-- maximum indel error ratio allowed (-d):      0.030
-- mean quality threshold (-Q):         30 # 平均质量
-- end quality threshold (-q):          30 # 末端质量
-- minimum read length allowed after trimming (-l):     100 # triming后允许最小长度
-- file format (-f):            Sanger/Illumina 1.8+ FASTQ (auto detected) # 文件格式
-- number of concurrent threads (-t):   4 # 线程数
Fri Nov 20 13:45:41 2020 >> started
|=================================================>| (100.00%)
Fri Nov 20 13:52:20 2020 >> done (399.163s) # 很快，7分钟
1582559 read pairs processed; of these: # 总配对reads
  16586 ( 1.05%) read pairs filtered out by quality control # 低质量的被过滤，只占1.05%
   2410 ( 0.15%) short read pairs filtered out after trimming by size control # 0.15%短读长，即不足100的reads被过滤
      0 ( 0.00%) empty read pairs filtered out after trimming by size control
1563563 (98.80%) read pairs available; of these: # 98.8%被保留
 449107 (28.72%) trimmed read pairs available after processing
1114456 (71.28%) untrimmed read pairs available after processing
log has been saved to "trim-trimmed.log".

# -t: number of threads to use
# -l: min length to keep after trimming
# -q: Quality threshold used for trimming at 3’ end
# -Q: mean quality threshold for a read
# -m: pair-end mode
# -o: Base name of output file # 输出文件

来看一下发生哪些变化，

-rwxrwxrwx 1 eric eric 525M Nov  3 11:53 sam01_1.fastq
-rwxrwxrwx 1 eric eric 525M Nov  3 11:53 sam01_2.fastq
-rwxrwxrwx 1 eric eric 514M Nov 20 13:52 trim-trimmed-pair1.fastq
-rwxrwxrwx 1 eric eric 513M Nov 20 13:52 trim-trimmed-pair2.fastq
-rwxrwxrwx 1 eric eric 2.3K Nov 20 13:52 trim-trimmed.log
# Trim之后，文件从525M/525M减小到514M/513M，生成一个日志文件

对TRIM后的文件查看质量

$ fastqc -f fastq trim-trimmed-pair1.fastq

观察一下，发生什么变化？首先，统计信息改变，总长度变为1563563（98.8%），跟skewer的日志中提供的信息一致，序列长度由151变为100-151，GC含量不变。

3’端的低质量碱基被过滤除掉

读长变为100-151

重复还在，可见不是这步被除掉的

接头序列被过滤除掉

假如，在建库过程中使用了特定的adaptor，你知道序列，但是没有被程序自动过滤除掉，则可以使用skewer来去掉接头，如接头序列是TGGAATTCTCGGGTGCCAAGGT,

skewer -x TGGAATTCTCGGGTGCCAAGGT -t 20 -l 10 -L 35 -q 30 adaptorQC.fastq.gz

# -x: adaptor sequence used
# -t: number of threads to use
# -l: min length to keep after trimming
# -L: Max length to keep after trimming, in this experiment we were expecting only small RNA fragments
# -Q: Quality threshold used for trimming at 3’ end. Use -m option to control the end you want to trim

另外，也可以锁定长度区间来FASTQ文件进行TRIM，这里不作介绍。