58.《Bioinformatics Data Skills》之
2021-08-24 本文已影响0人
DataScience
在分析基因组数据之前，我们首先需要关注测序的质量如何，回答下面两个问题：
- 测序技术错误如何分布？由何引起？
- 对下游分析产生如何影响？
前一节我们了解如何查看每个碱基的测序质量，然而面对数百万序列的测序，这样的观察便力不从心了。通过绘制碱基在read上的质量分布更为有效，最流行的测序质量评估工具为FastQC，输出结果包含各种图与测度。使用R的话可以使用类似功能的R包qrqc。
这里我们以qrqc为例展示碱基质量，进行如下实验：
- 分别使用
sickle与seqtk工具进行低质量碱基修剪 - 使用
qrqc工具进行修剪前后碱基质量得分展示，对比
本节数据下载地址：https://github.com/vsbuffalo/bds-files/tree/master/chapter-10-sequence
软件下载与安装参考官网：
sickle: https://github.com/najoshi/sickle
seqtk: https://github.com/lh3/seqtk
低质量碱基过滤
- 使用
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接结果文件名
- 使用
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%碱基质量范围
碱基质量的波动范围明显变小。
以上，其实进行碱基质量过滤非常简单（一行代码搞定），但是不能盲目地相信工具，所以需要通过图来展示工具的效果。谨慎的生物信息学分析通常都会不断调整参数，对比结果差别。
