博文

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

已有 18089 次阅读 2017-8-10 12:59 |个人分类:RNA-seq|系统分类:科研笔记| 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的时候有些信息可能有误，如本人用hisat2做reads 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相差较大（基于此可以联合rnaseq和bamqc使用），主要包括以下内容：

基本比对情况汇总表

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

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数据。