Changes

From Genome Analysis Wiki
Jump to navigationJump to search

Sam Library Usage Examples (view source)

Revision as of 09:07, 7 September 2011

1,908 bytes added ,  09:07, 7 September 2011
Remove unused code, add reference to BamUtil
Line 1: Line 1: −
This page contains examples on how to use the [[C++ Library: libbam|SAM/BAM Library]].
+
[[Category:libStatGen]]
 +
[[Category:libStatGen BAM]]
 +
 
 +
= Usage Examples =
 +
This page contains examples on how to use the [[C++ Library: libStatGen|SAM/BAM Library]].
 +
 
 +
See [[BamUtil]] for more examples of how it uses libStatGen to perform many operations.
    
== Simple Read/Write SAM/BAM files ==
 
== Simple Read/Write SAM/BAM files ==
Line 19: Line 25:  
   // Open the input file for reading.
 
   // Open the input file for reading.
 
   SamFile samIn;
 
   SamFile samIn;
−
   if(!samIn.OpenForRead(argv[1]))
+
   samIn.OpenForRead(argv[1]);
−
  {
  −
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
  −
      return(samIn.GetStatus());
  −
  }
      
   // Open the output file for writing.
 
   // Open the output file for writing.
 
   SamFile samOut;
 
   SamFile samOut;
−
   if(!samOut.OpenForWrite(argv[2]))
+
   samOut.OpenForWrite(argv[2]);
−
  {
  −
      fprintf(stderr, "%s\n", samOut.GetStatusMessage());
  −
      return(samOut.GetStatus());
  −
  }
      
   // Read the sam header.
 
   // Read the sam header.
 
   SamFileHeader samHeader;
 
   SamFileHeader samHeader;
−
   if(!samIn.ReadHeader(samHeader))
+
   samIn.ReadHeader(samHeader);
−
  {
  −
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
  −
      return(samIn.GetStatus());
  −
  }
      
   // Write the sam header.
 
   // Write the sam header.
−
   if(!samOut.WriteHeader(samHeader))
+
   samOut.WriteHeader(samHeader);
−
  {
  −
      fprintf(stderr, "%s\n", samOut.GetStatusMessage());
  −
      return(samOut.GetStatus());    
  −
  }
      
   SamRecord samRecord;
 
   SamRecord samRecord;
−
  SamValidationErrors samInvalidErrors;
      
     // Set returnStatus to success.  It will be changed
 
     // Set returnStatus to success.  It will be changed
Line 59: Line 48:  
     {
 
     {
 
         // Successfully read a record from the file, so write it.
 
         // Successfully read a record from the file, so write it.
−
         if(!samOut.WriteRecord(samHeader, samRecord))
+
         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();
   
     }
 
     }
   Line 81: Line 56:  
         samOut.GetCurrentRecordCount() << std::endl;
 
         samOut.GetCurrentRecordCount() << std::endl;
   −
     // Since the reads were successful, return the status based
+
     // Return success since a failure would have thrown
−
    // on the status of the writes.  If any failed, return
+
     // an exception.
−
     // their failure status.
   
     return(returnStatus);
 
     return(returnStatus);
 
  }
 
  }
Line 89: Line 63:       −
== Simple Read/Write Indexed BAM Files ==
+
== Read Only One Reference/Chromosome from Sorted/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.
 
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.
 
<source lang="cpp">
 
<source lang="cpp">
Line 101: Line 75:  
   // Open the input file for reading.
 
   // Open the input file for reading.
 
   SamFile samIn;
 
   SamFile samIn;
−
   if(!samIn.OpenForRead(inputFilename))
+
   samIn.OpenForRead(inputFilename);
−
  {
  −
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
  −
      return(samIn.GetStatus());
  −
  }
      
   // Open the bam index file for reading.
 
   // Open the bam index file for reading.
−
   if(!samIn.ReadBamIndex(indexFilename))
+
   samIn.ReadBamIndex(indexFilename);
−
  {
  −
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
  −
      return(samIn.GetStatus());
  −
  }
      
   // Open the output file for writing.
 
   // Open the output file for writing.
 
   SamFile samOut;
 
   SamFile samOut;
−
   if(!samOut.OpenForWrite(outputFilename))
+
   samOut.OpenForWrite(outputFilename);
−
  {
  −
      fprintf(stderr, "%s\n", samOut.GetStatusMessage());
  −
      return(samOut.GetStatus());
  −
  }
      
   // Read the sam header.
 
   // Read the sam header.
 
   SamFileHeader samHeader;
 
   SamFileHeader samHeader;
−
   if(!samIn.ReadHeader(samHeader))
+
   samIn.ReadHeader(samHeader);
−
  {
  −
      fprintf(stderr, "%s\n", samIn.GetStatusMessage());
  −
      return(samIn.GetStatus());
  −
  }
      
   // Write the sam header.
 
   // Write the sam header.
−
   if(!samOut.WriteHeader(samHeader))
+
   samOut.WriteHeader(samHeader);
−
  {
  −
      fprintf(stderr, "%s\n", samOut.GetStatusMessage());
  −
      return(samOut.GetStatus());    
  −
  }
      
   SamRecord samRecord;
 
   SamRecord samRecord;
−
  SamValidationErrors samInvalidErrors;
      
   // Loop through each Reference.
 
   // Loop through each Reference.
Line 150: Line 103:  
         numSectionRecords++;
 
         numSectionRecords++;
 
         // Successfully read a record from the file, so write it.
 
         // Successfully read a record from the file, so write it.
−
         if(!samOut.WriteRecord(samHeader, samRecord))
+
         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());
   
       }
 
       }
   Line 172: Line 115:  
}
 
}
 
</source>
 
</source>
 +
 +
 +
== Read a Specified Region from a Sorted & Indexed BAM File ==
 +
This example is found in pipeline/bam/WriteRegion.cpp and is part of the bam executable.
 +
 +
 +
<source lang="cpp">
 +
#include "SamFile.h"
 +
#include "Parameters.h"
 +
#include "BgzfFileType.h"
 +
 +
int writeRegion(int argc, char **argv)
 +
{
 +
    // Extract command line arguments.
 +
    static const int UNSPECIFIED_INT = -1;
 +
    String inFile = "";
 +
    String outFile = "";
 +
    String indexFile = "";
 +
    String refName = "";
 +
    int refID = UNSPECIFIED_INT;
 +
    int start = UNSPECIFIED_INT;
 +
    int end = UNSPECIFIED_INT;
 +
    bool noeof = false;
 +
 +
    ParameterList inputParameters;
 +
    BEGIN_LONG_PARAMETERS(longParameterList)
 +
        LONG_STRINGPARAMETER("in", &inFile)
 +
        LONG_STRINGPARAMETER("out", &outFile)
 +
        LONG_STRINGPARAMETER("bamIndex", &indexFile)
 +
        LONG_STRINGPARAMETER("refName", &refName)
 +
        LONG_INTPARAMETER("refID", &refID)
 +
        LONG_INTPARAMETER("start", &start)
 +
        LONG_INTPARAMETER("end", &end)
 +
        LONG_PARAMETER("noeof", &noeof)
 +
        END_LONG_PARAMETERS();
 +
 
 +
    inputParameters.Add(new LongParameters ("Input Parameters",
 +
                                            longParameterList));
 +
 +
    inputParameters.Read(argc-1, &(argv[1]));
 +
 +
    inputParameters.Status();
 +
 
 +
    // If no eof block is required for a bgzf file, set the bgzf file type to
 +
    // not look for it.
 +
    if(noeof)
 +
    {
 +
        // Set that the eof block is not required.
 +
        BgzfFileType::setRequireEofBlock(false);
 +
    }
 +
 +
    // Check to see if the in file was specified, if not, report an error.
 +
    if(inFile == "")
 +
    {
 +
        writeRegionUsage();
 +
        // mandatory argument was not specified.
 +
        std::cerr << "Missing mandatory argument: --in" << std::endl;
 +
        return(-1);
 +
    }
 +
    if(outFile == "")
 +
    {
 +
        writeRegionUsage();
 +
        // mandatory argument was not specified.
 +
        std::cerr << "Missing mandatory argument: --out" << std::endl;
 +
        return(-1);
 +
    }
 +
   
 +
    if(indexFile == "")
 +
    {
 +
        // In file was not specified, so set it to the in file
 +
        // + ".bai"
 +
        indexFile = inFile + ".bai";
 +
    }
 +
 +
    if(refID != -1 && refName.Length() != 0)
 +
    {
 +
        std::cerr << "Can't specify both refID and refName" << std::endl;
 +
        return(-1);
 +
    }
 +
 +
    SamFile samIn;
 +
    // Open the file for reading. 
 +
    samIn.OpenForRead(inFile);
 +
 +
    // If refName is set, use that.
 +
    if(refName.Length() != 0)
 +
    {
 +
        // Use Reference Name.
 +
        samIn.SetReadSection(refName.c_str(), start, end);
 +
    }
 +
    else
 +
    {
 +
        // Use Reference ID
 +
        samIn.SetReadSection(refID, start, end);
 +
    }
 +
 +
    // Open the output file for writing.
 +
    SamFile samOut;
 +
    samOut.OpenForWrite(outFile);
 +
 +
    // Open the bam index file for reading.
 +
    samIn.ReadBamIndex(indexFile);
 +
 +
    // Read & write the sam header.
 +
    SamFileHeader samHeader;
 +
    samIn.ReadHeader(samHeader);
 +
    samOut.WriteHeader(samHeader);
 +
 +
    // Read the sam records.
 +
    SamRecord samRecord;
 +
    // Track the status.
 +
    int numSectionRecords = 0;
 +
 +
    // Set returnStatus to success.
 +
    SamStatus::Status returnStatus = SamStatus::SUCCESS;
 +
 +
    // Keep reading records until they aren't anymore.
 +
    while(samIn.ReadRecord(samHeader, samRecord))
 +
    {
 +
        // Successfully read a record from the file, so write it.
 +
        samOut.WriteRecord(samHeader, samRecord);
 +
        ++numSectionRecords;
 +
    }
 +
 +
    std::cout << "Wrote " << outFile << " with " << numSectionRecords << " records.\n";
 +
    return(returnStatus);
 +
}
 +
</source>
 +
    
== Validate SAM/BAM Files Continuing Processing After Allowed Errors ==
 
== Validate SAM/BAM Files Continuing Processing After Allowed Errors ==
Line 190: Line 262:  
   }
 
   }
   −
   SamFile samIn;
+
   SamFile samIn(ErrorHandler::RETURN);
 
   // Open the file for reading.   
 
   // Open the file for reading.   
 
   if(!samIn.OpenForRead(inFile))
 
   if(!samIn.OpenForRead(inFile))
Retrieved from "http://genome.sph.umich.edu/wiki/Special:MobileDiff/1683...3695"

Navigation menu