专栏首页生信补给站R|fastqcr QC数据处理

R|fastqcr QC数据处理

FastQC是一款较常用的高通量数据质控软件,每个样本会得到一个zip和html的结果文件,查看略有不便。

fastqcr R包可以整理并分析多样不的zip格式的结果报告,当然也可以直接做fastQC分析。

一 加载安装R包

install.packages("fastqcr")
library(fastqcr)

二 R跑FastQC程序

fastqc(fq.dir = "~/Documents/FASTQ", # FASTQ files directory
      qc.dir = "~/Documents/FASTQC", # Results direcory
      threads = 4                    # Number of threads
      )

三 Fastcq结果整理绘图

3.1 读入fastqc的zip结果文件

qc.dir <- system.file("fastqc_results", package = "fastqcr")
qc.dir
# [1] "C:/Users/MJ/Documents/R/win-library/3.3/fastqcr/fastqc_results"

查看文件夹中文件

list.files(qc.dir)
#[1] "A1.R2.clean_fastqc.zip" "A1_R1.clean_fastqc.zip" "A2.R1.clean_fastqc.zip"
#[4] "A2.R2.clean_fastqc.zip"

3.2 整体统计结果 Aggregate & summary

1)Aggregate函数,Sample的基本统计;qc_stats()类似

qc <- qc_aggregate(qc.dir)
head(qc,3)
 sample                      module        status tot.seq seq.length pct.gc pct.dup
 <chr>                       <chr>         <chr>  <chr>   <chr>       <dbl>   <dbl>
1 170197A-Ti-1.R1.clean.fq.gz Basic Statis~ PASS   364589~ 150           51.    18.0
2 170197A-Ti-1.R1.clean.fq.gz Per base seq~ PASS   364589~ 150           51.    18.0
3 170197A-Ti-1.R1.clean.fq.gz Per tile seq~ WARN   364589~ 150           51.    18.0

其中:

sample: sample names
module: fastqc modules
status: fastqc module status for each sample
tot.seq: total sequences (i.e.: the number of reads)
seq.length: sequence length
pct.gc: percentage of GC content
pct.dup: percentage of duplicate reads

2)summary函数,Module的基本统计

head(summary(qc),3)

module nb_samples nb_fail nb_pass nb_warn failed warned

  <chr>                 <dbl>   <dbl>   <dbl>   <dbl> <chr>                   <chr>
1 Adapter Content          4.      0.      4.      0. NA                      NA    
2 Basic Statistics         4.      0.      4.      0. NA                      NA    
3 Kmer Content             4.      4.      0.      0. 170197A-Ti-1.R1.clean.~ NA        

其中:

module: fastqc modules
nb_samples: the number of samples tested
nb_pass, nb_fail, nb_warn: the number of samples that passed, failed and warned, respectively.
failed, warned: the name of samples that failed and warned, respectively.

3)dplyr函数

aggregated完可以使用dplyr得到重点关注的信息:如WARN和FALL样本:

library(dplyr)
qc %>%select(sample, module, status) %>% filter(status %in% c("WARN", "FAIL")) %>%arrange(sample)  sample                       module                    status
  <chr>                        <chr>                     <chr>
1 170197A-Ti-1.R1.clean.fq.gz  Per tile sequence quality WARN  
2 170197A-Ti-1.R1.clean.fq.gz  Per sequence GC content   WARN  
3 170197A-Ti-1.R1.clean.fq.gz  Kmer Content              FAIL  

附:dplyr简单介绍 dplyr

3.3 异常信息统计

1) Module异常:qc_warns & qc_problems函数

qc_warns(qc, "module") # See which module warned in the most samples
 module                    nb_problems sample                                      
 <chr>                           <int> <chr>                                      
1 Per sequence GC content             4 170197A-Ti-1.R1.clean.fq.gz, 170197A-Ti-1.R~
2 Per tile sequence quality           2 170197A-Ti-1.R1.clean.fq.gz, 170197A-Ti-10.~

查看特定模块的异常

qc_problems(qc, "module",  name = "Per sequence GC content")
 module                  nb_problems sample                       status
 <chr>                         <int> <chr>                        <chr>
1 Per sequence GC content           4 170197A-Ti-1.R1.clean.fq.gz  WARN  
2 Per sequence GC content           4 170197A-Ti-1.R2.clean.fq.gz  WARN  
3 Per sequence GC content           4 170197A-Ti-10.R1.clean.fq.gz WARN  
4 Per sequence GC content           4 170197A-Ti-10.R2.clean.fq.gz WARN  

2)Sample 异常 : qc_fails & qc_problems 函数

qc_problems(qc, "sample", compact = FALSE) #fail 和 warn的 ,compact展示方式
sample                       nb_problems module                    status
  <chr>                              <int> <chr>                     <chr>
1 170197A-Ti-1.R1.clean.fq.gz            3 Kmer Content              FAIL  
2 170197A-Ti-1.R1.clean.fq.gz            3 Per tile sequence quality WARN  
3 170197A-Ti-1.R1.clean.fq.gz            3 Per sequence GC content   WARN  
4 170197A-Ti-1.R2.clean.fq.gz            3 Per tile sequence quality FAIL  
5 170197A-Ti-1.R2.clean.fq.gz            3 Kmer Content              FAIL  
6 170197A-Ti-1.R2.clean.fq.gz            3 Per sequence GC content   WARN  

四 绘图及生成报告

4.1 导入单个样本fastq结果

qc.file <- system.file("fastqc_results", "170197A-Ti-1.R1.clean_fastqc.zip", package = "fastqcr")
qc <- qc_read(qc.file) #读入qc的结果,后续绘图使用
names(qc) #查看包含的元素

4.2 绘制样本基本图形

qc_plot(qc, "Per sequence GC content")
qc_plot(qc, "Per base sequence quality")
qc_plot(qc, "Per sequence quality scores")
qc_plot(qc, "Per base sequence content")
qc_plot(qc, "Sequence duplication levels")

4.3 生成单样本网页报告

1)图形有解释信息的报告

qc_report(qc.file, result.file = "one-sample-report-with-interpretation",interpret = TRUE)

在工作目录下生成html的网页报告,图形含有解释信息(如下),表格可交互。

2)图形无解释信息的报告

qc_report(qc.file, result.file = "one-sample-report",interpret = FALSE)

4.4 多样本网页报告

qc.dir <- system.file("fastqc_results", package = "fastqcr")
qc_report(qc.dir, result.file = "~/Desktop/multi-qc-result",experiment = "Exome sequencing")

五 参考资料

fastqcr: An R Package Facilitating Quality Controls of Sequencing Data for Large Numbers of Samples

本文分享自微信公众号 - 生信补给站(Bioinfo_R_Python),作者:MJ

原文出处及转载信息见文内详细说明,如有侵权,请联系 yunjia_community@tencent.com 删除。

原始发表时间:2018-04-26

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

我来说两句

0 条评论
登录 后参与评论

相关文章

  • TCGA数据挖掘 | Xena - TCGA数据下载

    TCGA (The Cancer Genome Atlas)作为目前超常用的癌症基因信息的数据库,有多种肿瘤的表达谱数据,变异信息(mutation,copy ...

    西游东行
  • 数据分析|R-异常值处理

    前面介绍了拿到脏数据后,对缺失值的处理数据分析|R-缺失值处理,今天分享一下另一种脏数据-异常值的处理。

    西游东行
  • ggplot2|详解八大基本绘图要素

    ggplot2是由Hadley Wickham创建的一个十分强大的可视化R包。按照ggplot2的绘图理念,Plot(图)= data(数据集)+ Aesthe...

    西游东行
  • Python 面向对象 高阶-描述符与设计模式#学习猿地

    > 当一个类中,包含了三个魔术方法(`__get__,__set__,__delete__`)之一,或者全部时,那么这个类就称为描述符类

    学习猿地
  • Python 面向对象 高阶-描述符与设计模式#学习猿地

    > 当一个类中,包含了三个魔术方法(`__get__,__set__,__delete__`)之一,或者全部时,那么这个类就称为描述符类

    学习猿地
  • shell技巧分享(三)

    songleo
  • Android编程设计模式之单例模式实例详解

    本文实例讲述了Android编程设计模式之单例模式。分享给大家供大家参考,具体如下:

    砸漏
  • SAP 脚本录制与回放功能

    SAP系统的脚本录制功能,支持VB Script,可以将屏幕操作记录下来,转换成VB Script代码,VB Script代码编辑修改后可用在office软件...

    用户5495712
  • 编程小白 | 每日一练(75)

    这道理放在编程上也一并受用。在编程方面有着天赋异禀的人毕竟是少数,我们大多数人想要从编程小白进阶到高手,需要经历的是日积月累的学习,那么如何学习呢?当然是每天都...

    C语言入门到精通
  • JavaScript异步编程:Generator与Async

    从Promise开始,JavaScript就在引入新功能,来帮助更简单的方法来处理异步编程,帮助我们远离回调地狱。 Promise是下边要讲的Generator...

    贾顺名

扫码关注云+社区

领取腾讯云代金券