工欲善其事必先利其器
QualiMap
QualiMap 是一款主要由Fernando Garcı ́a-Alcalde、Konstantin Okonechnikov 开发的用于评估高通量测序数据质量的工具。主要用于分析和可视化测序数据的质量指标。
其具有以下特性
- 全平台，QualiMap是用Java编写的，因此它无需编译可直接使用，适用于Linux、MacOS、Windows
- 用户界面友好，易于使用，同时也提供了命令行界面，方便自动化和集成到其他分析流程中
- 支持多种测序数据，包括但不限于RNA-Seq、ChIP-Seq、DNA-Seq等
- 软件提供了一系列质量指标，如比对率、覆盖率、GC含量等，帮助用户全面了解测序数据的质量状况
- 专门设计用于分析高通量测序数据，适用于大规模分析
发表文章
文1：Qualimap: evaluating next-generation sequencing alignment data（2012） 文2：Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data（2016） 期刊：Bioinformatics 作者&单位：Fernando Garcı ́a-Alcalde、Konstantin Okonechnikov [菲利普研究中心（西班牙）和 马克斯·普朗克感染生物学研究所（德国柏林）] DOI1：10.1093/bioinformatics/bts503 DOI2：10.1093/bioinformatics/btv566 官网：http://qualimap.conesalab.org/
简要用途
QualiMap的主要功能包括：
- 测序覆盖度分析：QualiMap可以计算测序覆盖度，包括平均覆盖度、覆盖度分布和覆盖度深度的统计信息。这有助于确定基因组或目标区域是否被均匀和充分地覆盖。
- 基因组比对质量评估：软件可以分析比对到参考基因组的读段（reads）的质量，包括比对的准确性、比对的多样性和潜在的错误。
- GC含量分析：QualiMap能够分析测序数据中的GC含量，并与预期的GC含量进行比较，以检测可能的偏差。
- 多态性和变异检测：软件可以帮助识别基因组中的多态性位点和变异，这对于遗传研究和变异分析非常重要。
- 可视化工具：QualiMap提供了丰富的图形和图表，使用户能够直观地查看和解释分析结果。
- 报告生成：用户可以生成包含所有重要统计数据和图形的综合报告，便于分享和进一步分析。
由于QualiMap提供了全面的质量评估工具，它在基因组学、转录组学和表观遗传学等领域的研究中非常有用。通过确保数据质量，研究人员可以更有信心地进行下游分析，如基因表达分析、变异检测和基因组注释。
如何安装
可以直接采用conda 安装管理软件，但是目前conda默认安装的版本是2.2.1
代码语言：javascript复制conda create -n wes ##先创建小环境，如果已经创建，无需再次执行
conda activate wes
codna install -y qualimap
mamba activate wes
mamba install qualimap
如果想使用最新版（2.3），也可以直接下载，解压，无需编译可直接使用
代码语言：javascript复制wget -c https://bitbucket.org/kokonech/qualimap/downloads/qualimap_v2.3.zip
unzip qualimap_v2.3.zip
cd qualimap_v2.3/
不过需要注意你的java版本：
代码语言：javascript复制$java -version
openjdk version "17.0.3-internal" 2022-04-19
OpenJDK Runtime Environment (build 17.0.3-internal 0-adhoc..src)
OpenJDK 64-Bit Server VM (build 17.0.3-internal 0-adhoc..src, mixed mode, sharing)
如果java版本过高，会报错：
这个错误信息表明Java虚拟机（JVM）无法识别MaxPermSize这个参数。这通常是因为你使用的Java版本是Java 8或更高版本，在这些版本中，MaxPermSize参数已经被移除，因为PermGen space 已经被元空间（Metaspace）所取代。
此时需要修改QualiMap使用的Java启动脚本或配置文件，移除-XX:MaxPermSize这个JVM参数
修改启动文件
正常启动
快速上手
bamqc —单个bam文件的QC统计
代码语言：javascript复制##最小化使用
qualimap bamqc -bam ERR089819.bam -outdir qualimap_results
-bam :必需参数，指定bam文件
-outdir : 可选参数，
rnaseq — RNA-Seq样本的bam文件的QC统计
代码语言：javascript复制qualimap rnaseq -bam kidney.bam -gtf human.64.gtf -outdir rnaseq_qc_results
-gtf : 必须参数，基因组注释文件
multi-bamqc —多样本的bam文件QC统计
代码语言：javascript复制qualimap multi-bamqc -d gh2ax_chip_seq.txt -outdir gh2ax_multibamqc
-d 必须参数：描述输入数据的配置文件；要求是2列或者3列的制表符分割文件，第一列是样品名，第二列是单个文件bamqc分析结果的路径（或者是样本bam文件路径，需要加 -r 参数），第三列是对应样本的分组
简单演示
代码语言：javascript复制##实例应用
qualimap bamqc -bam 231_0.bam -c -outformat PDF:HTML -outdir ~/qulimap_test/231_0_bam -nt 12 --java-mem-size=30G
-nt 线程数，默认是8
-outformat 设置输出文件格式
-c 在统计图表中绘制染色体边界
--java-mem-size 设置允许使用的内存大小，默认是1200M。不设置的话，很容易报错内存不够。（承担测试用的bam文件大小3G,报错 WARNING: out of memory!）
##实例
qualimap multi-bamqc -d test_seq.txt -outdir ~/qulimap_test/multi_out
$cat test_seq.txt
231_0 /home/project/mutiqc/231_0_bam_qc 0
231_12d /home/project/mutiqc/231_12d_bam_qc 12d
231_48h /home/project/mutiqc/231_48h_bam_qc 48h
qualimap rnaseq -bam LC231_0_1.sort.bam -gtf /home/jmzeng/rna/human/pipeline/gencode.v39.annotation.gtf -outdir ~/qulimap_test/rnaseq_qc_results --java-mem-size=4G
更多结果图表见：http://qualimap.conesalab.org/doc_html/samples.html#alignments
其余用法及参数
通用参数
代码语言：javascript复制-outdir ： 输出文件夹，输出HTML报告和原始数据
-outfile ： 输出文件，默认输出pdf报告文件
-outformat ：输出文件格式，可设置输出PDF或HTML格式，默认是HTML
bamqc
代码语言：javascript复制##生成PDF报告
qualimap bamqc -bam my_sample.bam -outdir my_sample_bamqc -outformat PDF
##分析特定区域
qualimap bamqc -bam my_sample.bam -outdir my_sample_bamqc -gff target_regions.gff
##链特异性分析
qualimap bamqc -bam my_sample.bam -outdir my_sample_bamqc -p strand-specific-forward
##可选参数
-cl： 设置每个bin的目标覆盖度直方图的上限，默认是50X
-dl： 设置重复率（duplication rate）上限，默认是50
-gd： 选择与基因组GC分布进行比较的物种，可选值为HUMAN或MOUSE。
-gff： 感兴趣区域的特征文件，格式可以是GFF/GTF或BED。
-hm ： 在插入缺失分析中考虑的同源多聚体的最小大小，默认是3。在分析过程中，如果一个同源多聚体的长度小于这个设定的大小，那么它所涉及的任何潜在插入缺失都不会被计入统计。可以帮助过滤掉较短的同源多聚体区域，因为它们可能对插入缺失的检测不太敏感，或者可能产生较多的假阳性结果。
-ip：激活重叠配对读段的检测。这样做的目的是识别那些重叠的读段对，并在计算平均覆盖度时适当地调整它们。具体来说，当检测到两个读段重叠时，它们共同覆盖的区域（overlap-region）只会被计算一次，而不是两次，这样可以得到一个更加准确的平均覆盖度（adapted mean coverage）【注意，双末端测序，在插入片段尺寸较小的情况下，配对末端读段可能会在很大比例上重叠。这意味着两个读段的测序结果可能会覆盖相同的基因组区域，导致这部分区域的覆盖度被过度估计。】
-nr：每个块分析的读段数量，默认是1000
-nw：使用的窗口数，默认是400
-oc：保存每个 base non-zero coverage 的文件。警告：对于大型基因组，预期会产生大文件。
-os：报告特征文件定义之外区域的信息（如果没有设置-gff选项，则忽略此选项）
-p：测序文库协议，可选 strand-specific-forward,strand-specific-reverse or non-strand-specific；默认是 non-strand-specific
-sd：激活此选项以从分析中跳过重复的比对。如果BAM文件中没有标记重复，则QualiMap会检测并可以选择跳过
-sdmode ：如果激活此选项，则跳过特定类型的重复比对
0 : only flagged duplicates 只有BAM文件中被标记为重复的比对结果会被跳过（默认项）；这意味着如果BAM文件中的比对结果已经被之前的工具（如Picard的MarkDuplicates或samtools的rmdup）标记为重复，那么Qualimap会忽略这些比对结果，并继续分析剩余的非重复比对结果。
1 : only estimated by Qualimap 只有QualiMap检测为重复的跳过，Qualimap会根据其内部的重复率估计方法来识别重复，并在分析中排除这些比对结果。
2 : both flagged and estimated 混合模式，既跳过BAM文件中已标记的重复比对结果，也跳过Qualimap自己检测到的重复比对结果
测序库的链特异性（Library strand specificity） 是指测序过程中能够保留RNA或DNA模板链方向信息的能力。这个信息对于正确解释测序结果非常重要，尤其是在转录组学（transcriptomics）和基因表达分析中。根据测序协议的不同，可以将测序库分为以下几种类型：
- 非链特异性（non-strand-specific）: 在这种协议中，测序读段不能提供关于其来源的模板链的方向信息。也就是说，读段可能来自模板的任一链，因此无法确定其确切的链方向。
- 正向链特异性（forward-stranded）: 在这种协议中，测序读段保留了与其来源的模板链相同的方向信息。也就是说，读段的方向与编码RNA或基因的方向一致。
- 反向链特异性（reverse-stranded）: 在这种协议中，测序读段保留了与其来源的模板链相反的方向信息。也就是说，读段的方向与编码RNA或基因的方向相反。
rnaseq
代码语言：javascript复制## 分析特定区域
qualimap rnaseq -bam my_sample.bam -gtf target_regions.gtf -outdir my_sample_rnaseq
## 分析双末端测序数据
qualimap rnaseq -bam my_sample.bam -gtf my_annotations.gtf -outdir my_sample_rnaseq -pe
## 检查链特异性
qualimap rnaseq -bam my_sample.bam -gtf my_annotations.gtf -outdir my_sample_rnaseq -p strand-specific-forward
##可选参数
-a：指定用于计数的算法
- uniquely-mapped-reads: 默认选项，只计算唯一映射到参考基因组的读段。
- proportional: 如果一个读段可以映射到多个位置，它的计数会按比例分配到这些位置
-npb：计算5'到3'偏差时，指定上游和下游核苷酸的数量，默认是100个
-ntb：计算5'到3'偏差时，指定要考虑的高表达转录本的数量，默认是1000个.
-oc ：指定输出计数结果的路径
-p：指定测序文库协议，同bamqc
-pe：如果设置此标志，表示实验是双末端测序，分析过程中将计算成对的片段（即一对配对读段）的数量，而不是单个读段的数量
-s：如果设置此标志，表示输入的文件已经按名称排序。如果未设置，将进行额外的按名称排序。这个参数只在配对末端分析时需要。
multi-bamqc
代码语言：javascript复制## 可选参数
-c : 仅在-r模式下使用。在图表中标记染色体边界
-gff: 仅在-r模式下使用。特征文件，包含感兴趣区域的信息，格式可以是GFF/GTF或BED
-hm ： 仅在-r模式下使用。在插入缺失(indel)分析中被考虑的同源多聚体的最小大小（默认为3）
-nr 仅在-r模式下使用。在一个块中分析的读段数量（默认是1000）。
-nw 仅在-r模式下使用。窗口数量（默认是400）
-r：输入为原始BAM文件。如果激活此选项，将首先对每个样本运行BAM QC过程，然后执行多样本分析
count — 转录组数据计数的统计，用于量化表达水平
代码语言：javascript复制qualimap counts -d GlcN_countsqc_input.txt -c -s mouse -outdir glcn_mice_counts
-d 必须参数：描述输入数据的配置文件；要求是4列的制表符分割文件。第一列是样品名，第二列是实验条件（ex:处理或未处理），第三列是样品计数数据的文件的路径；第四列是计数数据中包含计数值的列的索引（用于当所有样本的计数都包含在一个文件中，但需要统计不同样本列的情况）
##可选参数
-c ：条件比较。目前最多支持两个条件的比较。意味着你可以比较两组样本之间的表达量差异
-s ：使用给定物种的内置信息文件：HUMAN 或 MOUSE
-i ：包含基因GC含量、长度和类型信息的文件路径。这个文件用于提供额外的基因特征信息，以便在分析中使用
-k ：计数数量的阈值。这可以用来过滤掉低表达的基因，即只有当基因的表达量计数超过这个阈值时，它才会被包括在分析中
-R ：R脚本可执行文件的路径。默认情况下，可以从系统的 $PATH 环境变量中找到
clustering——表观遗传特征的聚类
代码语言：javascript复制qualimap clustering -sample clustering/hmeDIP.bam -control clustering/input.bam -regions annotations/transcripts.human.64.bed -outdir clustering_result
## 必须参数
-sample ：输入一个用逗号分隔的BAM文件的列表，处理组
-control ：输入一个用逗号分隔的BAM文件的列表，control组
-regions ： 设定区域文件的路径（BED或GFF格式的注释文件的路径）
## 可选参数
-b 用于设置bin（箱）的大小，默认为100
-c 设定用户想要将数据分成的组数。要求输入一个用逗号分隔的聚类大小列表。它涉及到数据分析中聚类方法的一个关键参数，用户可以通过逗号将多个数值分开，指定他们希望数据被划分成的不同组的数量。
-expr 指定实验的名称，用于生成报告或标识特定分析的实验
-f 设置片段的平滑长度，为了统一数据的长度，所有的测序读取（reads）都会被扩展或调整到这个指定的片段长度
-l 设置上游偏移，默认为2000
-name 输入一个用逗号分隔的重复名称列表，用于标识或区分在实验中重复的样本
-r 设置下游偏移，默认为500
-viz 设置可视化类型: heatmap or line
comp-couns——输入bam文件和注释文件，计算映射到每个区域reads的数量
代码语言：javascript复制qualimap comp-counts -bam kidney.bam -gtf ../annotations/human.64.gtf -out kidney.counts
## 可选参数
-a ：用于指定计数算法。可以选择使用的计数算法，有两种选择：“uniquely-mapped-reads”（默认）或者“proportional”
-id ：针对GTF文件的特定属性。用于指定GTF中用作特征ID的属性。具有相同ID的区域将作为同一特征的一部分进行汇总。默认为“gene_id”
-out：设置输出文件的路径
-p ：同ranseq
-pe：同ranseq
-s ：同ranseq
-type：针对GTF文件的特定属性。用于指定在计数时考虑的GTF的第三列的值。其他类型将被忽略。默认为“exon”
参考：
- http://qualimap.conesalab.org/doc_html/command_line.html
- http://qualimap.conesalab.org/doc_html/analysis.html#bamqc
- http://qualimap.conesalab.org/doc_html/tools.html#clustering
- https://zhuanlan.zhihu.com/p/613397485
