From f4b4a665b5a9f1a72703230531649220ba6047b0 Mon Sep 17 00:00:00 2001
From: Paul Igor Costea <paul.igor.costea@embl.de>
Date: Mon, 12 Feb 2018 22:50:53 +0100
Subject: [PATCH] Refine the coverage counting to account for hard and soft
clipping.
---
src/qaTools/qaCompute.cpp | 33 ++++++++++++++++++++++++++-------
1 file changed, 26 insertions(+), 7 deletions(-)
diff --git a/src/qaTools/qaCompute.cpp b/src/qaTools/qaCompute.cpp
index e087283..aa77cae 100644
--- a/src/qaTools/qaCompute.cpp
+++ b/src/qaTools/qaCompute.cpp
@@ -526,13 +526,32 @@ int main(int argc, char *argv[])
++duplicates;
} else {
if (!userOpt.spanCov) {
- //All entries in SAM file are represented on the forward strand! (See specs of SAM format for details)
- ++entireChr[core->pos];
- ++usedReads;
- if ((uint32_t)(core->pos+core->l_qseq) >= chrSize)
- --entireChr[chrSize-1];
- else
- --entireChr[core->pos+core->l_qseq];
+ //All entries in SAM file are represented on the forward strand! (See specs of SAM format for details)
+ uint32_t* cigar = bam1_cigar(b);
+ uint32_t pp = core->pos+1;
+ int i = 0;
+ if(((*cigar) & BAM_CIGAR_MASK) == BAM_CSOFT_CLIP || ((*cigar) & BAM_CIGAR_MASK) == BAM_CHARD_CLIP) {//BWA is actively fucking with me.
+ //The start of the alignment is reported without this clip.
+ ++cigar; ++i;
+ }
+ while (i<core->n_cigar) {
+ ++i;
+ if (((*cigar) & BAM_CIGAR_MASK) != BAM_CMATCH) {
+ pp = pp + ((*cigar) >> BAM_CIGAR_SHIFT);
+ } else {
+ ++entireChr[pp];
+ pp = pp + ((*cigar) >> BAM_CIGAR_SHIFT);
+ --entireChr[pp];
+ }
+ ++cigar;
+ }
+
+ //++entireChr[core->pos];
+ ++usedReads;
+ //if ((uint32_t)(core->pos+core->l_qseq) >= chrSize)
+ // --entireChr[chrSize-1];
+ //else
+ // --entireChr[core->pos+core->l_qseq];
} else {
//Computing span coverage.
//Only consider first read in pair! and extend a bit to the end of the insert
--
GitLab