Adding total read count in BAM file to POPINS_SAMPLE_INFO file instead of command line.

src/crop_unmapped.h

@@ -352,6 +352,7 @@ findOtherReads(BamFileOut & matesStream,
 template<typename TAdapterTag>
 int
 crop_unmapped(double & avgCov,
+	unsigned & totalReads,
         Triple<CharString> & fastqFiles,
         CharString & matesBam,
         CharString const & mappingBam,
@@ -406,7 +407,7 @@ crop_unmapped(double & avgCov,
 
     // Iterate over the input file.
     BamAlignmentRecord record;
-    unsigned totalReads = 0;
+    totalReads = 0;
     unsigned long alignedBaseCount = 0;
     while (!atEnd(inStream))
     {
@@ -451,10 +452,6 @@ 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);
 
@@ -489,7 +486,8 @@ crop_unmapped(Triple<CharString> & fastqFiles,
         float as_factor)
 {
     double cov;
-    return crop_unmapped(cov, fastqFiles, matesBam, mappingBam, humanSeqs, tag, as_factor);
+    unsigned totalReads;
+    return crop_unmapped(cov, totalReads, fastqFiles, matesBam, mappingBam, humanSeqs, tag, as_factor);
 }
 
 #endif // #ifndef NOVINS_CROP_UNMAPPED_H_

src/popins2_crop_unmapped.h

@@ -86,17 +86,17 @@ int popins2_crop_unmapped(int argc, char const ** argv){
     // Crop unmapped reads and reads with unreliable mappings from the input bam file.
     if (options.adapters == "HiSeqX")
     {
-        if (crop_unmapped(info.avg_cov, fastqFiles, matesBam, options.mappingFile, options.humanSeqs, HiSeqXAdapters(), as_factor) != 0)
+        if (crop_unmapped(info.avg_cov, info.total_reads, fastqFiles, matesBam, options.mappingFile, options.humanSeqs, HiSeqXAdapters(), as_factor) != 0)
             return 7;
     }
     else if (options.adapters == "HiSeq")
     {
-        if (crop_unmapped(info.avg_cov, fastqFiles, matesBam, options.mappingFile, options.humanSeqs, HiSeqAdapters(), as_factor) != 0)
+        if (crop_unmapped(info.avg_cov, info.total_reads, fastqFiles, matesBam, options.mappingFile, options.humanSeqs, HiSeqAdapters(), as_factor) != 0)
             return 7;
     }
     else
     {
-        if (crop_unmapped(info.avg_cov, fastqFiles, matesBam, options.mappingFile, options.humanSeqs, NoAdapters(), as_factor) != 0)
+        if (crop_unmapped(info.avg_cov, info.total_reads, fastqFiles, matesBam, options.mappingFile, options.humanSeqs, NoAdapters(), as_factor) != 0)
             return 7;
     }
 

src/util.h

@@ -258,6 +258,7 @@ struct SampleInfo
     CharString sample_id;
     CharString bam_file;
     double avg_cov;
+    unsigned total_reads;
     unsigned read_len;
     CharString adapter_type;
 
@@ -310,6 +311,8 @@ readSampleInfo(SampleInfo & info, CharString & filename)
             info.bam_file = value;
         else if (field.compare("AVG_COV") == 0)
             info.avg_cov = lexicalCast<double>(value);
+	else if (field.compare("NUM_READS") == 0)
+	    lexicalCast<unsigned>(info.total_reads, value);
         else if (field.compare("READ_LEN") == 0)
             lexicalCast<unsigned>(info.read_len, value);
         else if (field.compare("ADAPTER_TYPE") == 0)
@@ -339,6 +342,7 @@ writeSampleInfo(SampleInfo & info, CharString & filename)
     stream << "SAMPLE_ID" << "\t" << info.sample_id << "\n";
     stream << "BAM_FILE" << "\t" << info.bam_file << "\n";
     stream << "AVG_COV" << "\t" << info.avg_cov << "\n";
+    stream << "NUM_READS" << "\t" << info.total_reads << "\n";
     stream << "READ_LEN" << "\t" << info.read_len << "\n";
     stream << "ADAPTER_TYPE" << "\t" << info.adapter_type << "\n";