diff --git a/include/output_data.h b/include/output_data.h
index 10c1f39..1bfd0a4 100644
--- a/include/output_data.h
+++ b/include/output_data.h
@@ -122,6 +122,9 @@ class Output_BAM : public Output_FQ
     uint64_t num_reads_with_both_secondary_supplementary_alignment = ZeroDefault; // the number of long reads with both secondary and supplementary alignment.
     uint64_t forward_alignment = ZeroDefault;  // Total number of forward alignments
     uint64_t reverse_alignment = ZeroDefault;  // Total number of reverse alignments
+    int reads_with_mods = ZeroDefault; // Total number of reads with modification tags
+    int reads_with_mods_pos_strand = ZeroDefault; // Total number of reads with modification tags on the positive strand
+    int reads_with_mods_neg_strand = ZeroDefault; // Total number of reads with modification tags on the negative strand
 
     // Map of reads with supplementary alignments
     std::map<std::string, bool> reads_with_supplementary;
diff --git a/include/utils.h b/include/utils.h
new file mode 100644
index 0000000..9637f1e
--- /dev/null
+++ b/include/utils.h
@@ -0,0 +1,15 @@
+// Utility functions
+#ifndef UTILS_H
+#define UTILS_H
+
+/// @cond
+#include <string>
+/// @endcond
+
+// Print a message to stdout in a thread-safe manner
+void printMessage(std::string message);
+
+// Print an error message to stderr in a thread-safe manner
+void printError(std::string message);
+
+#endif // UTILS_H
diff --git a/src/hts_reader.cpp b/src/hts_reader.cpp
index dccec8c..180579a 100644
--- a/src/hts_reader.cpp
+++ b/src/hts_reader.cpp
@@ -14,6 +14,7 @@ Class for reading a set number of records from a BAM file. Used for multi-thread
 #include <algorithm>  // std::find
 
 #include "hts_reader.h"
+#include "utils.h"
 
 // HTSReader constructor
 HTSReader::HTSReader(const std::string & bam_file_name){
@@ -111,10 +112,10 @@ int HTSReader::readNextRecords(int batch_size, Output_BAM & output_data, std::mu
         // Follow here to get base modification tags:
         // https://github.com/samtools/htslib/blob/11205a9ba5e4fc39cc8bb9844d73db2a63fb8119/sam_mods.c
         // https://github.com/samtools/htslib/blob/11205a9ba5e4fc39cc8bb9844d73db2a63fb8119/htslib/sam.h#L2274
-        std::cout << "[TEST] Accessing base modification tags" << std::endl;
         hts_base_mod_state *state = hts_base_mod_state_alloc();
         if (bam_parse_basemod(record, state) >= 0) {
-            std::cout << "Base modification tags found" << std::endl;
+            printMessage("Base modification tags found");
+            // std::cout << "Base modification tags found" << std::endl;
             mod_tag_present = true;
 
             // Iterate over the state object to get the base modification tags
@@ -124,9 +125,18 @@ int HTSReader::readNextRecords(int batch_size, Output_BAM & output_data, std::mu
             int pos = 0;
             while ((n=bam_next_basemod(record, state, mods, 10, &pos)) > 0) {
                 for (int i = 0; i < n; i++) {
-                    std::cout << "Base modification at position " << pos << std::endl;
-                    std::cout << "Base modification type: " << mods[i].modified_base << std::endl;
-                    std::cout << "Base modification likelihood: " << mods[i].qual / 256.0 << std::endl;
+                    // Struct definition: https://github.com/samtools/htslib/blob/11205a9ba5e4fc39cc8bb9844d73db2a63fb8119/htslib/sam.h#L2226
+                    printMessage("Found base modification at position " + std::to_string(pos));
+                    printMessage("Modification type: " + std::string(1, mods[i].modified_base));
+                    printMessage("Canonical base: " + std::string(1, mods[i].canonical_base));
+                    printMessage("Likelihood: " + std::to_string(mods[i].qual / 256.0));
+                    printMessage("Strand: " + std::to_string(mods[i].strand));
+
+                    //
+                    // std::cout << "Base modification at position " << pos << std::endl;
+                    // std::cout << "Base modification type: " << mods[i].modified_base << std::endl;
+                    // std::cout << "Base modification likelihood: " << mods[i].qual / 256.0 << std::endl;
+                    // std::cout << "Base modification strand: " << mods[i].strand << std::endl;
                 }
             }
 
diff --git a/src/utils.cpp b/src/utils.cpp
new file mode 100644
index 0000000..4d27130
--- /dev/null
+++ b/src/utils.cpp
@@ -0,0 +1,26 @@
+#include "utils.h"
+
+/// @cond
+#include <stdio.h>
+#include <string>
+#include <iostream>
+#include <mutex>
+/// @endcond
+
+
+// Define print mutex
+std::mutex print_mtx;
+
+// Thread-safe print message function
+void printMessage(std::string message)
+{
+    std::lock_guard<std::mutex> lock(print_mtx);
+    std::cout << message << std::endl;
+}
+
+// Thread-safe print error function
+void printError(std::string message)
+{
+    std::lock_guard<std::mutex> lock(print_mtx);
+    std::cerr << message << std::endl;
+}