luria的个人博客分享 http://blog.sciencenet.cn/u/luria

博文

Qualimap2.多样本高通量数据高级QC工具

已有 4883 次阅读 2017-8-10 12:59 |个人分类:RNA-seq|系统分类:科研笔记|关键词:bam,,sam,,QC,,RNAseq,,Qualimap,,归一化图,,reads| Sam, bam, RNAseq, Qualimap, 归一化图

今天推荐一款叫做Qualimap2的软件,它可以对reads mapping的结果进行统计,如果是RNAseq数据,甚至还可能算出counts并进行简单的后续分析。此外,输入的注释文件可以是gff/gtf/bed,相比于其它同类软件使用起来更加方便。 ==========================================================================

Qualimap2的安装:

从官网http://qualimap.bioinfo.cipf.es/doc_html/intro.html下载Qualimap2(当前最新版为v2.2.1),因为它是用java写的,可直接运行。但是程序调用了一些R包,需先做如下安装:

R #进入R环境

install.packages(“optparse”)

在弹出的镜像选择列表中选China(Beijing) [注意:不是China(Beijing)[https]!]。可直接安装。

source(http://bioconductor.org/biocLite.R)

biocLite(c(“NOISeq”, “Repitools”, “Rsamtools”, “GenomicFeatures”, “rtracklayer))

注:也可以使用其自带的安装程序安装,但是自动安装通常会报错,因此不建议。自动安装代码如下:

cd qualimap_v2.2.1

Rscript scripts/installDependencies.r

==========================================================================

Qualimap2的使用:

图形界面版的没有什么好说的,直接运行qualimap即可.

因为处理的数据量通常会比较大(多),一般会在linux下批量运行,因此本博文重点关注命令行操作代码。

1. bamqc方法

bam文件进行QC

./qualimap bamqc -bam xx.srt.bam -outdir bamqc_result -outformat PDF:HTML

# -bam指定bam文件路径,-outdir指定输出的文件夹,-outformat指定输出的报告格式,可选PDF/HTML/PDF:HTML,这里建议PDF:HTML,表示同时生成两种格式的报告

# 不同软件进行比对时可能对reads的处理不同(例如是否保留未比对上的reads,以及multi-reads的处理等)。因此在使Qualimap2的时候有些信息可能有误,如本人用hisat2reads mapping,在Qualimap2统计ATCGN时,N超过100%,这显然不可能,这可能是算法上存在bug

不过绝大多部结果还是可信的,这里仅列出几项最重要的,大家可以自行查看生成的网页版报告(下同):

比对结果信息汇总表

Insert Size统计图 (点击可查看大图,下同)

从统计图中可获取到文库约为350bp

如果有reference对应的注释文件则

./qualimap bamqc -bam xx.srt.bam -gff yy.gtf -oc count.matrix -outdir bamqc -outformat PDF:HTML

# -bam指定bam文件路径;-gff指定参考基因组注释信息,输入的格式可以是GFF/GTF/BED-oc保存有read覆盖位点的深度信息于此文件中(没有覆盖的位点不写出);-outdir指定输出文件夹;-outformat指定生成报告格式。

生成的报告中会在上述的基础上增加一个区域内信息汇总表如下:

区域内信息汇总表

2. rnaseq方法

RNA比对数据进行QC

需要提供bam文件以及注释文件用于生成基因表达均一化

./qualimap rnaseq -bam xx.srt.bam -gtf yy.gtf -oc count.matrix -outdir rnaseq_result -pe -outformat PDF:HTML

# 其它选项与上述bamqc相同,仅增加-pe表示是PE测序的

# 虽然代码相似但运行完生成的结果文件与bamqc相差较大(基于此可以联合rnaseqbamqc使用)主要包括以下内容:

基本比对情况汇总表


比对到的区域注释,同时还会给出对应的图


5’ 3’bias情况图


新转录本占比情况

3. multi-bamqc方法

Qualimap2还可以对多个bam文件进行批量统计先构建一个mapfile1.txt文件如下:

各列用tab分开,第一列是样品名,第二列是样品对应的bam文件所在位置(如上图,也可以输入bamqc方法的结果文件夹,建议直接用bam文件作为输入),第三列是分组信息如果各样本独立,每个Group都不一样;如果是一组,写什么文字无关紧要,同一组文字相同即可。

./qualimap multi-bamqc -r -d mapfile1.txt -gff xx.gtf -outdir multi_bamqc -outformat PDF:HTML

# multi-bamqc可以利用bamqc方法的结果也可以直接输入bam文件,流程化完成,如果输入的是bam文件,则需指定-r选项。-d指定上述mapfile1.txt文件路径。其它选项与bamqc相同。

# 运行完生成的报告中,统计了一些内容如下(以下未列出的可查看生成的网页版报告):

各样本的整体统计表


各物种比对情况PCA


比对重复率情况

各样品Insert Size统计图

在网页版的输出结果中,每个图都可以通过右键查看和保存。

4. comp-counts方法

可以利用Qualimap2生成counts matrix,具体代码如下:

awk ‘{print \$1,\$2}’ mapfile1.txt | while read a b; do ./qualimap comp-counts -bam \$b -gtf xx.gtf -out \$a.count -pe; done

# 运行完成后生成一系列count文件,head查看结果如下:

5. counts方法

使用counts方法前需生成样品信息文件,这里取名为mapfile2.txt,软件说明中写到,到目前最新版,仅能对两组数据进行比较,如果有多组会报错。mapfile2.txt的格式如下:

第一列是样品名,第二列是组别(具体用什么文字无关紧要,能区分样品组别关系即可),第三列是counts文件的路径,可以是multi-bamqc生成的counts文件(也可以是其它软件生成的,但为了衔接流畅,建议用multi-bamqc生成counts),第四列指定在counts文件中哪一列记录count数据。

./qualimap counts -c -d mapfile2.txt -outdir counts_result -outformat PDF:HTML

# -c表示两组间进行比对;-d指定上述mapfile2.txt文件的路径,其它选项同bamqc

# 完成后在counts_result文件夹中生成一些报告,其中ComparisonReport.html是按组分析的结果

两组counts分布盒形图

低表达敏感性分析

GlobalReport.html是所有样品画到同一张同中,不分组。报告中,有个Options表,其中可以看出Counts threshold5,表明程序默认的counts阈值为5,这个值可以通过-k参数调节。

counts的密度分布图

注:横坐标为log2(counts)

样本间相关情况散点图

各样本测序饱和度分析


各样品的表达量盒形图



http://blog.sciencenet.cn/blog-2970729-1070459.html

上一篇:Alpha多样性指数之simpson指数
下一篇:BioConda --生信工作者的福音

0

该博文允许注册用户评论 请点击登录 评论 (2 个评论)

数据加载中...

Archiver|手机版|科学网 ( 京ICP备14006957 )

GMT+8, 2019-6-25 15:41

Powered by ScienceNet.cn

Copyright © 2007- 中国科学报社

返回顶部