58.《Bioinformatics Data Skills》之

在分析基因组数据之前，我们首先需要关注测序的质量如何，回答下面两个问题：

• 测序技术错误如何分布？由何引起？
• 对下游分析产生如何影响？

前一节我们了解如何查看每个碱基的测序质量，然而面对数百万序列的测序，这样的观察便力不从心了。通过绘制碱基在read上的质量分布更为有效，最流行的测序质量评估工具为FastQC，输出结果包含各种图与测度。使用R的话可以使用类似功能的R包qrqc。

这里我们以qrqc为例展示碱基质量，进行如下实验：

1. 分别使用sickle与seqtk工具进行低质量碱基修剪
2. 使用qrqc工具进行修剪前后碱基质量得分展示，对比

本节数据下载地址：https://github.com/vsbuffalo/bds-files/tree/master/chapter-10-sequence
软件下载与安装参考官网：
sickle: https://github.com/najoshi/sickle
seqtk: https://github.com/lh3/seqtk

低质量碱基过滤

1. 使用sickle进行低质量碱基过滤：

$ sickle se -f untreated1_chr4.fq -t sanger -o untreated1_chr4_sickle.fq
SE input file: untreated1_chr4.fq
Total FastQ records: 204355
FastQ records kept: 203121
FastQ records discarded: 1234

其中se代表单末端（single pair），-f接文件名，-t接质量类型，-o接结果文件名

2. 使用seqtk子命令trimfq进行低质量碱基结果过滤：

$ seqtk trimfq untreated1_chr4.fq > untreated1_chr4_trimfq.fq

过滤前后碱基质量对比

> library(qrqc)
> library(dplyr)
> library(ggplot2)

使用R包qrqc readSeqFile函数导入过滤前后的测序文件（只关注碱基质量，忽略其它特征）：

> fqfiles <- c(none="untreated1_chr4.fq", sickle="untreated1_chr4_sickle.fq", trimfq="untreated1_chr4_trimfq.fq")
> seq_info <- lapply(fqfiles, readSeqFile, hash = F, kmer = F)

采用getQual函数计算各个文件包含的碱基质量，并整合为一个data.frame：

> sqs_qual <- lapply(seq_info, getQual) %>% bind_rows(.id = "trimmer")
> head(sqs_qual)
  trimmer position ymin alt.lower    lower   middle    upper alt.upper ymax
1    none        1    0  31.18736 32.10936 32.40624 32.70312  32.88125   33
2    none        2    0  28.74042 31.50799 32.55130 33.24792  33.69917   34
3    none        3    0  28.08877 31.15024 32.41778 33.09209  33.63684   34
4    none        4    0  27.88972 31.12117 32.41488 33.09941  33.63976   34
5    none        5    0  28.28632 31.20898 32.45740 33.16739  33.66696   34
6    none        6    0  27.74141 30.94298 32.36516 33.05458  33.62183   34
      mean
1 32.33774
2 32.03997
3 31.78172
4 31.75728
5 31.91736
6 31.70643

使用ggplot包绘制不同read位置碱基平均质量分布（图1）：

> ggplot(sqs_qual, aes(x = position, y = mean, linetype = trimmer)) +
    geom_line() +
    labs(y = "quality (mean)") +
    theme_bw()  

图1 没有过滤和分别使用sickle，trimfq工具过滤碱基前后平均质量分布

可以看出过滤之后碱基平均质量有所提升，不过随着read长度增加碱基质量明显下降。

使用qualPlot函数展示碱基质量在10%~90%之间的分布（图2）：

> qualPlot(seq_info, quartile.color=NULL, mean.color=NULL) +
    scale_y_continuous("quality (sanger)") +
    theme_bw()

图2 10%~90%碱基质量范围

碱基质量的波动范围明显变小。

以上，其实进行碱基质量过滤非常简单（一行代码搞定），但是不能盲目地相信工具，所以需要通过图来展示工具的效果。谨慎的生物信息学分析通常都会不断调整参数，对比结果差别。