qcoverage is a stand-alone Java application that calculates coverage statistics given a BAM file with mapped reads, and a GFF3 file with features of interest.
It can be run in 2 different modes:
Physical coverage includes both (a) the reads themselves, plus (b) the bases that lie between a pair of matched reads. So physical coverage includes all bases on a fragment - those at either end that were sequenced and can be seen in the reads PLUS those bases that can be inferred to be part of the sequenced fragment but that were not themselves sequenced. Physical coverage is only relevant for paired reads as for unpaired reads, sequence coverage is identical to physical coverage.
This link has more information on the differences between the two.
qcoverage requires java 7 and (ideally) a multi-core machine with at least 20GB of RAM.
Download the qcoverage tar file
Untar the tar file into a directory of your choice
You should see jar files for qcoverage and its dependencies:
[oholmes@minion0 qcoverage]$ tar xjvf qcoverage.tar.bz2 x antlr-3.2.jar x jopt-simple-3.2.jar x picard-1.110.jar x qbamfilter-1.0pre.jar x qcommon-0.1pre.jar x qpicard-0.1pre.jar X qio-0.1pre.jar x qcoverage-0.7pre.jar x sam-1.110.jar [oholmes@minion0 qcoverage]$
For sequence coverage:
java -jar qcoverage.jar -t seq --bam test1.bam --gff3 GRCh37.gff3 -o report.txt --log logfile
For physical coverage:
java -jar qcoverage.jar -t phys --bam test1.bam --gff3 GRCh37.gff3 -o report.txt --log logfile
Option Description ------ ----------- -V, -v, --version Print version info. --bai <BAI file> BAI file for the BAM file. --bam <BAM file> BAM file containing the reads. --gff3 <GFF3 file> The GFF3 file defining the features. -h, --help Shows this help message. --log File where log output will be directed (must have write permissions). --loglevel Logging level required, e.g. INFO, DEBUG. (Optional; if no parameter is specified, will default to INFO) -n <thread_count> Number of worker threads (yields n+1 total threads). (Optional) -o, --output <outputfile> Specifies the output file. --per-feature Perform per-feature coverage. (Optional) -q, --query <expression> The query string for selecting reads for coverage. (Optional) -t, --type <type> The type of coverage to perform (seq, sequence, phys, physical). --xml Output report in XML format. (Optional) --vcf Output report in VCF format. Needs the per-feature flag to also be set (Optional)
It is possible to have output from the --per-feature mode in either XML or VCF format. The XML format is very extensive with an XML table for each feature showing a full breakdown of how many bases of the feature were covered at what levels. For long regions with variable coverage, this table can run to hundreds of lines so an XML-format output file containing hundreds of thousands of features can run to tens of millions of lines - it can be as big as the BAM itself. In most cases, the level of detail provided by the XML output is overkill and the VCF format is sufficient. The VCF file contains a single line for each feature so it is almost identical in length to the GFF3 file that defines the features. Each line gives details about the feature including start and stop positions as well as a high level summary of the coverage. This information is sufficient to provide average coverage for the feature and to allow for the assessment of capture bait performance.
The VCF output looks like:
##fileformat=VCFv4.0 ##bam_file=/ABCD_1234/bamfile.bam ##gff_file=/SureSelect_All_Exon_50mb_filtered_baits_1-200_20110524_shoulders.gff3 ##FILTER=<ID=LowQual,Description="REQUIRED: QUAL < 50.0"> ##INFO=<ID=B,Number=.,Type=String,Description="Bait name"> ##INFO=<ID=BE,Number=.,Type=String,Description="Bait end position"> ##INFO=<ID=ZC,Number=.,Type=String,Description="bases with Zero Coverage"> ##INFO=<ID=NZC,Number=.,Type=String,Description="bases with Non Zero Coverage"> ##INFO=<ID=TOT,Number=.,Type=String,Description="Total number of sequenced bases"> #CHROM POS ID REF ALT QUAL FILTER INFO chr1 1 . . . . . B=fill;BE=14166;ZC=12240;NZC=1926;TOT=16840 chr1 14167 . . . . . B=bait_3_100;BE=14266;ZC=100;NZC=0;TOT=0 chr1 14267 . . . . . B=bait_2_100;BE=14366;ZC=100;NZC=0;TOT=0 chr1 14367 . . . . . B=bait_1_100;BE=14466;ZC=100;NZC=0;TOT=0 chr1 14467 . . . . . B=bait;BE=14587;ZC=64;NZC=57;TOT=445 chr1 14588 . . . . . B=bait_1_100;BE=14638;ZC=0;NZC=51;TOT=1455 chr1 14639 . . . . . B=bait;BE=14883;ZC=0;NZC=245;TOT=9763
If we look at one output line in more detail we see:
chr1 14467 . . . . . B=bait;BE=14587;ZC=64;NZC=57;TOT=445
which we interpret as:
To accelerate performance where hardware permits (e.g: multicore processors; multi-processor machines), pass the -n option to specify an appropriate number of worker threads.
For example, applying -n 15 to yield 15 worker threads plus 1 supervising thread:
java -Xmx20G -jar qcoverage.jar -n 15 -t seq --bam test1.bam --gff3 GRCh37_primary_chr13.gff3 -o report.txt --log logfile
Adequate node resources must be allocated to ensure speedup when using this technique on any cluster. It must also be noted that qcoverage uses more memory when in multithreaded mode.
qcoverage reuses the query engine underpinning qbamfilter, allowing the user to include only BAM reads satisfying a specified query expression. Like qbamfilter, this query expression is supplied with the -q,--query option.
Identical query-language semantics apply to qcoverage as for qbamfilter. A complete specification of the qbamfilter query language can be found on the qbamfilter page.
The default filter that is used by QCMG when running qcoverage is as follows:
sequence coverage
java -jar qcoverage.jar -t seq -q "and (flag_ReadFailsVendorQuality==false, flag_DuplicateRead==false, flag_ReadUnmapped==false, flag_NotprimaryAlignment==false)" -bam test1.bam --gff3 some.gff3 -o report.txt --log logfile
physical coverage
java -jar qcoverage.jar -t seq -q "and (flag_ReadFailsVendorQuality==false, flag_DuplicateRead==false, flag_ReadUnmapped==false, flag_NotprimaryAlignment==false)" -bam test1.bam --gff3 some.gff3 -o report.txt --log logfile
which are ensuring that only non-duplicate, mapped and primary alignment reads are considered for the coverage summaries.
Some other examples are:
Excluding reads above a defined ISIZE cutoff during physical coverage:
java -jar qcoverage.jar -q "ISIZE < 200" -t phys --bam test1.bam --gff3 some.gff3 -o report.txt --log logfile
Including only reads satisfying a ZP-quality of "AAA":
java -jar qcoverage.jar --query "ZP==AAA" -t seq --bam test1.bam --gff3 some.gff3 -o report.txt --log logfile
Excluding unpaired reads from physical coverage:
java -jar qcoverage.jar --query "flag_ReadPaired==true" -t phys --bam test1.bam --gff3 some.gff3 -o report.txt --log logfile