平时拿到数据后首先要进行质控环节，其中FastaQC软件的使用最为广泛，它可以为每一个样品生成一个html报告和一个‘zip’ 文件，zip解压之后生成fastqc_data.txt和summary.txt的文件，里面包含了测序样品的质量信息，但是如果有几十个甚至上百个数据时候，我们总不可能一一的打开每个html文件查看，因此Multiqc软件成为了解决这一问题的首选方案，可以基于结果zip文件进行多样品的整合分析，信息生成在一张报告上。然而哪里有问题，哪里就有R包，对于R语言爱好者，这里再提供一种可替代的方法--fastqcr包，好处就是不管有没有linux环境时候，都可以快速对数据进行质检，一起来学习下吧~~
主要功能
1) 安装和运行
- fastqc_install(): Unix环境下(MAC OSX和Linux)安装最新版的FastaQC
- fastqc(): R中运行FastQC
2) 汇总多样品FastQC报告
- qc <- qc_aggregate():合并多样品的FastQC报告至一个数据框
- summary(qc):针对qc_aggregate的结果生成汇总表
- qc_stats(): FastaQC报告的一般统计信息
3) 查看问题
- qc_fails(qc):显示失败的样本或模块。
- qc_warns(qc): 显示警告的样本或模块。
- qc_problems(qc): qc_fails()和qc_warns()的并集。 显示哪些样本或模块失败或警告。
4) 导入绘制Fastqc报告
- qc_read(): 读取FastQC数据到R种
- qc_plot(): 绘制FastQC数据
5）生成单样品和多样品QC报告
- qc_report():创建一个HTML文件，其中包含一个或多个文件的FastQC报告。 输入可以是包含多个FastQC报告的目录，也可以是单个样本FastQC报告。
安装加载
## 两种方式
# Cran
install.packages("fastqcr")
# Github
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/fastqcr")
library(fastqcr)
使用
fastqc函数，输入样品fq格式文件，输出样品的zip压缩包。
fastqc(fq.dir = "~/Documents/FASTQ", # FASTQ files directory
qc.dir = "~/Documents/FASTQC", # Results direcory
threads = 4 # Number of threads
)
Demo数据：在该包的目录文件下有现成的5个样品S1-S5的zip压缩包文件，这里直接作为示例数据进行演示
library(fastqcr)
# Demo QC dir
qc.dir <- system.file("fastqc_results", package = "fastqcr")
qc.dir
# 列出文件
list.files(qc.dir)
[1] "S1_fastqc.zip" "S2_fastqc.zip" "S3_fastqc.zip" "S4_fastqc.zip" "S5_fastqc.zip"
汇总FastQC报告
qc <- qc_aggregate(qc.dir)
qc
汇总信息
生成7列信息：
- sample： 样品名字
- module：fastqc模块
- status：每个样品的fastqc 模块特征
- tot.seq: 全部序列
- seq.length：测序长度
- pct.gc： GC含量
- pct.dup： 重复序列百分比
查看质控状态是“Warn”或者“FAIL”的样品及模块
library(dplyr)
qc %>%
select(sample, module, status) %>%
filter(status %in% c("WARN", "FAIL")) %>%
arrange(sample)
Warn和Fail信息
统计信息,查看每个模块各种状态上的样品个数
# Summary of qc
summary(qc)
# 生成一般统计信息
qc_stats(qc)
汇总信息
我们也可以针对模块特征去查询
# 哪些样品是失败的
qc_fails(qc, "module")
# 看那些样品是警告
qc_warns(qc, "module")
# 看哪些样品是有问题的 包含警告与失败
qc_problems(qc, "module")
# 另外加上compact = FALSE使表格不那么紧凑
qc_problems(qc, "module", compact = FALSE)
问题模块查询
生成HTML报告
直接全部样品生成报告，类似Multiqc功能。
qc.dir <- system.file("fastqc_results", package = "fastqcr")
qc.dir
# 生成报告
qc_report(qc.dir, result.file = "~/Desktop/multi-qc-result",
experiment = "Exome sequencing of colon cancer cell lines")
多样品汇总
我们也可以针对单样品生成报告,比如只生成S1样品，加上interpret= TRUE参数具有交互式。
qc.file <- system.file("fastqc_results", "S1_fastqc.zip", package = "fastqcr")
qc.file
# 生成报告
qc_report(qc.file, result.file = "one-sample-report-with-interpretation",
interpret = TRUE)
单样品
针对某个模块画图
我们都知道fastqc基本包括如下几个模块：
- “Summary”,
- “Basic Statistics”,
- “Per base sequence quality”,
- “Per sequence quality scores”,
- “Per base sequence content”,
- “Per sequence GC content”,
- “Per base N content”,
- “Sequence Length Distribution”,
- “Sequence Duplication Levels”,
- “Overrepresented sequences”,
- “Adapter Content”
若我们只想快速度看其中某个模块的情况，通过qc_plot(qc,"模块名称")即可出图
library(patchwork)
A1 <- qc_plot(qc.file, "Per sequence GC content")
A2 <- qc_plot(qc.file, "Per base sequence quality")
A3 <- qc_plot(qc.file, "Per sequence quality scores")
A4 <- qc_plot(qc.file, "Per base sequence content")
(A1+A2)/(A3+A4)
针对模块作图
总之该包的函数使用非常简单，花费一点点时间即可快速掌握，除了Multiqc软件，该方法其实也不失为一种备选方案。
网友评论