Sam Library Usage Examples

Usage Examples

This page contains examples on how to use the SAM/BAM Library.

Simple Read/Write SAM/BAM files

The following example reads in a sam/bam file and writes it out as a sam/bam file. The file format of the input sam/bam is determined by the SamFile class based on reading the type from the file. The file format of the output sam/bam file is determined by the SamFile class based on the extension of the output file. A ".bam" extension indicates a BAM file. All other extensions indicate SAM files.

#include "SamFile.h"
#include "SamValidation.h"

int main(int argc, char ** argv)
{
   // Check for the appropriate number of arguments.
   if(argc != 3)
   {
      printf("./bam <inputFile> <outputFile.sam/bam>\n");
      exit(-1);
   }

   // Open the input file for reading.
   SamFile samIn;
   if(!samIn.OpenForRead(argv[1]))
   {
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
      return(samIn.GetStatus());
   }

   // Open the output file for writing.
   SamFile samOut;
   if(!samOut.OpenForWrite(argv[2]))
   {
      fprintf(stderr, "%s\n", samOut.GetStatusMessage());
      return(samOut.GetStatus());
   }

   // Read the sam header.
   SamFileHeader samHeader;
   if(!samIn.ReadHeader(samHeader))
   {
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
      return(samIn.GetStatus());
   }

   // Write the sam header.
   if(!samOut.WriteHeader(samHeader))
   {
      fprintf(stderr, "%s\n", samOut.GetStatusMessage());
      return(samOut.GetStatus());     
   }

   SamRecord samRecord;
   SamValidationErrors samInvalidErrors;

    // Set returnStatus to success.  It will be changed
    // to the failure reason if any of the writes fail.
    SamStatus::Status returnStatus = SamStatus::SUCCESS;

    // Keep reading records until ReadRecord returns false.
    while(samIn.ReadRecord(samHeader, samRecord))
    {
        // Successfully read a record from the file, so write it.
        if(!samOut.WriteRecord(samHeader, samRecord))
        {
            // Failed to write a record.
            fprintf(stderr, "%s\n", samOut.GetStatusMessage());
            returnStatus = samOut.GetStatus();
        }
    }

    if(samIn.GetStatus() != SamStatus::NO_MORE_RECS)
    {
        // Failed to read a record.
        fprintf(stderr, "%s\n", samIn.GetStatusMessage());
        // Set the return status to the reason why
        // the read failed.
        returnStatus = samIn.GetStatus();
    }

    std::cout << std::endl << "Number of records read = " << 
        samIn.GetCurrentRecordCount() << std::endl;
    std::cout << "Number of records written = " << 
        samOut.GetCurrentRecordCount() << std::endl;

    // Since the reads were successful, return the status based
    // on the status of the writes.  If any failed, return
    // their failure status.
    return(returnStatus);
 }

Simple Read/Write Indexed BAM Files

This example reads in the inputFilename bam file and writes it back out section by section to the specified outputFilename, starting with section -1. It also prints a count of the number of records in each section.

#include "SamFile.h"
#include "SamValidation.h"

int ReadIndexedBam(const char* inputFilename,
                   const char* outputFilename,
                   const char* indexFilename)
{
   // Open the input file for reading.
   SamFile samIn;
   if(!samIn.OpenForRead(inputFilename))
   {
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
      return(samIn.GetStatus());
   }

   // Open the bam index file for reading.
   if(!samIn.ReadBamIndex(indexFilename))
   {
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
      return(samIn.GetStatus());
   }

   // Open the output file for writing.
   SamFile samOut;
   if(!samOut.OpenForWrite(outputFilename))
   {
      fprintf(stderr, "%s\n", samOut.GetStatusMessage());
      return(samOut.GetStatus());
   }

   // Read the sam header.
   SamFileHeader samHeader;
   if(!samIn.ReadHeader(samHeader))
   {
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
      return(samIn.GetStatus());
   }

   // Write the sam header.
   if(!samOut.WriteHeader(samHeader))
   {
      fprintf(stderr, "%s\n", samOut.GetStatusMessage());
      return(samOut.GetStatus());     
   }

   SamRecord samRecord;
   SamValidationErrors samInvalidErrors;

   // Loop through each Reference.
   for(int i = -1; i < 23; i++)
   {
      int numSectionRecords = 0;
      samIn.SetReadSection(i);
      // Keep reading records until they aren't more.
      while(samIn.ReadRecord(samHeader, samRecord))
      {
         numSectionRecords++;
         // Successfully read a record from the file, so write it.
         if(!samOut.WriteRecord(samHeader, samRecord))
         {
            // Failed to write a record.
            fprintf(stderr, "%s\n", samOut.GetStatusMessage());
         }
      }

      if(samIn.GetStatus() != SamStatus::NO_MORE_RECS)
      {
         // Failed to read a record.
         fprintf(stderr, "%s\n", samIn.GetStatusMessage());
      }

      std::cout << "Reference ID " << i << " has " << numSectionRecords 
                << " records" << std::endl;
   }
   
   std::cout << "Number of records = " << samIn.GetCurrentRecordCount() << std::endl;
   
   return(0);
}

Validate SAM/BAM Files Continuing Processing After Allowed Errors

This example demonstrates calling validation on a SAM/BAM file. It also demonstrates a more complicated read loop that continues until a configurable number of errors are hit.

#include "SamFile.h"
#include "SamValidation.h"

int Validate(String inFile, int quitAfterErrorNum, int maxReportedErrors, SamFile::SortedType sortType)
{
   // Check to see if the in file was specified.  If it wasn't this is an error.
   if(inFile == "")
   {
      // In file was not specified but it is mandatory.
      std::cerr << "--in is a mandatory argument for validate, "
                << "but was not specified" << std::endl;
      return(-1);
   }

   SamFile samIn;
   // Open the file for reading.   
   if(!samIn.OpenForRead(inFile))
   {
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
      return(samIn.GetStatus());
   }

   // Set the sorting validation type.
   samIn.setSortedValidation(sortType);

   // Read the sam header.
   SamFileHeader samHeader;
   if(!samIn.ReadHeader(samHeader))
   {
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
      return(samIn.GetStatus());
   }

   // Read the sam records.
   SamRecord samRecord;

   // Track the status.
   SamStatus::Status status = SamStatus::SUCCESS;

   // Keep reading records until the end of the file is reached.
   int numValidRecords = 0;
   int numInvalidRecords = 0;
   int numErrorRecords = 0;
   int numRecords = 0;
   int numReportedErrors = 0;
   int totalErrorRecords = 0;

   SamValidationErrors invalidSamErrors;

   // Keep reading records from the file until SamFile::ReadRecord
   // indicates to stop (returns false).
   while( (samIn.ReadRecord(samHeader, samRecord))||
          (SamStatus::isContinuableStatus(samIn.GetStatus()) &&
           ( (quitAfterErrorNum < 0) ||
             (numErrorRecords < quitAfterErrorNum) )) )
   {
      ++numRecords;
      if(samIn.GetStatus() == SamStatus::SUCCESS)
      {
         // Successfully set the record, so check to see if it is valid.
         // Clear any errors in the list.
         invalidSamErrors.clear();
         if(!SamValidator::isValid(samHeader, samRecord, invalidSamErrors))
         {
            // The record is not valid.
            ++numInvalidRecords;
            ++totalErrorRecords;
            if(numReportedErrors < maxReportedErrors)
            {
               std::cout << "Record " << numRecords << std::endl
                         << invalidSamErrors << std::endl;
               numReportedErrors += invalidSamErrors.numErrors();
               status = SamStatus::INVALID;
            }
         }
         else
         {
            // Valid record, so increment the counter.
            ++numValidRecords;
         }
      }
      else
      {
         // Error reading the record.
         ++numErrorRecords;
         if(numReportedErrors < maxReportedErrors)
         {
            // report error.
            std::cout << "Record " << numRecords << std::endl
                      << samIn.GetStatusMessage() << std::endl
                      << std::endl;
            ++numReportedErrors;
            ++totalErrorRecords;
            status = samIn.GetStatus();
         }
      }
   }

   if(samIn.GetStatus() != SamStatus::NO_MORE_RECS)
   {
      ++numErrorRecords;
      // The last read call had a failure, so report it.
      if(numReportedErrors < maxReportedErrors)
      {
         std::cout << "Record " << numRecords << ": ";
         std::cout << std::endl << samIn.GetStatusMessage() << std::endl;
      }
      status = samIn.GetStatus();
   }

   printf("\nNumber of records read = %d\n", numRecords);
   printf("Number of valid records = %d\n", numValidRecords);
   printf("Returning: %d (%s)\n", status, SamStatus::getStatusString(status));
   return(status);
}