From 572f3abf2181f99f59ef346e0923f83b2a3a6377 Mon Sep 17 00:00:00 2001
From: Birte Kehr <keb35159@rhskl26.uni-regensburg.de>
Date: Wed, 6 Dec 2023 11:28:35 +0100
Subject: [PATCH] Adding total read count in BAM file to command line output.
---
src/crop_unmapped.h | 8 ++++++++
1 file changed, 8 insertions(+)
diff --git a/src/crop_unmapped.h b/src/crop_unmapped.h
index eea8f3e..8d4a68a 100644
--- a/src/crop_unmapped.h
+++ b/src/crop_unmapped.h
@@ -406,6 +406,7 @@ crop_unmapped(double & avgCov,
// Iterate over the input file.
BamAlignmentRecord record;
+ unsigned totalReads = 0;
unsigned long alignedBaseCount = 0;
while (!atEnd(inStream))
{
@@ -416,6 +417,9 @@ crop_unmapped(double & avgCov,
if (hasFlagDuplicate(record) or hasFlagSecondary(record) or
hasFlagQCNoPass(record) or hasFlagSupplementary(record)) continue;
+ // Count 'interesting' reads..
+ totalReads += 1;
+
if (!hasFlagUnmapped(record))
alignedBaseCount += length(record.seq);
@@ -447,6 +451,10 @@ crop_unmapped(double & avgCov,
avgCov = (double)alignedBaseCount / (double)genomeLength;
std::ostringstream msg;
+ msg << "BAM file lists alignments for " << totalReads << " reads after duplicate removal and filtering on QC pass flag.";
+ printStatus(msg);
+
+ msg.str("");
msg << "Map of low quality mates has " << otherReads.size() << " records.";
printStatus(msg);
--
GitLab