Difference between revisions of "C++ Class: SamFile"

From Genome Analysis Wiki
Jump to navigationJump to search
 
(26 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 +
[[Category:C++]]
 +
[[Category:libStatGen]]
 +
[[Category:libStatGen BAM]]
 +
 
== Reading/Writing SAM/BAM Files In Your Program ==
 
== Reading/Writing SAM/BAM Files In Your Program ==
 
The '''SamFile''' class allows a user to easily read/write a SAM/BAM file.
 
The '''SamFile''' class allows a user to easily read/write a SAM/BAM file.
Line 6: Line 10:
 
'''Future Enhancements:''' Add the ability to read alignments that match a given start, end position for a specific reference sequence.  
 
'''Future Enhancements:''' Add the ability to read alignments that match a given start, end position for a specific reference sequence.  
  
 +
This class is part of [[C++ Library: libStatGen|C++ Library: libStatGen]].
 +
 +
=== Class Documentation ===
  
−
=== Class Methods ===
+
See: http://csg.sph.umich.edu//mktrost/doxygen/current/classSamFile.html
  
−
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
+
== Child Classes ==
−
|-style="background: #f2f2f2; text-align: center;"  '''SamFile Class Methods'''
+
=== SamFileReader ===
−
! Method Name !!  Description
+
http://csg.sph.umich.edu//mktrost/doxygen/current/classSamFileReader.html
−
|-
 
−
| <code>bool SamFile::IsEOF()</code>
 
−
| bool: true if the end of file has been reached, false if not.
 
−
Be careful using this method when you are only reading a specific section since you may reach the end of your section without hitting the end of the file
 
  
−
|-
+
=== SamFileWriter ===
−
| <code>bool SamFile::OpenForRead(const char* filename)</code>
+
http://csg.sph.umich.edu//mktrost/doxygen/current/classSamFileWriter.html
−
| Opens the specified file for reading.
 
−
Determines if it is a BAM/SAM file by reading the beginning of the file.
 
−
Returns true if successfully opened reading, false if not.
 
−
|-
 
−
| <code>bool SamFile::OpenForWrite(const char * filename)</code>
 
−
| bool: true if successfully opened, false if not.
 
−
Opens as BAM file if the specified filename ends in .bam. Otherwise it is opened as a SAM file.
 
−
Returns true if successfully opened for writing, false if not.
 
−
|-
 
−
| <code>bool SamFile::ReadBamIndex(const char * filename)</code>
 
−
| bool: true if the bam index file was successfully read, false if not.
 
−
Reads the specified bam index file.  It must be read prior to setting a read section, for seeking and reading portions of a bam file.
 
−
|-
 
−
| <code>void SamFile::Close()</code>
 
−
| Close the file if there is one open.
 
−
|-
 
−
| <code>bool SamFile::ReadHeader(SamFileHeader& header)</code>
 
−
| Reads the header section from the file and stores it in the passed in header.
 
−
Returns true if successfully read, false if not.
 
−
|-
 
−
| <code>bool SamFile::WriteHeader(SamFileHeader& header)</code>
 
−
| Writes the specified header into the file.
 
−
Returns true if successfully written, false if not.
 
−
|- 
 
−
| <code>bool SamFile::ReadRecord(SamFileHeader& header, SamRecord& record)</code>
 
−
| Reads the next record from the file and stores it in the passed in record.
 
  
−
If it is an indexed BAM file and SetReadSection was called, only alignments in the section specified by <code>SetReadSection</code> are read.  If they all have already been read, this method returns false.
+
== Statistics ==
 +
=== Statistic Generation ===
  
−
Validates that the record is sorted according to the value set by <code>setSortedValidation</code>.  No sorting validation is done if specified to be unsorted, or <code>setSortedValidation</code> was never called.
+
The following statistics can be optionally recorded when reading a SamFile by specifying <code>SamFile::GenerateStatistics()</code> and displayed with <code>SamFile::PrintStatistics()</code>
  
−
Returns false if the record was not successfully read or not properly sorted.  Returns true if successfully read and properly sorted.
+
The statistics only reflect alignments that were successfully read from the BAM file. Alignments that failed to parse from the file are not reflected in the statistics, but alignments that are invalid for other reasons may show up in the statistics.
−
|-
 
−
| <code>bool SamFile::WriteRecord(SamFileHeader& header, SamRecord& record)</code>
 
−
| Writes the specified record into the file.
 
−
Validates that the record is sorted according to the value set by <code>setSortedValidation</code>.  No sorting validation is done if specified to be unsorted, or <code>setSortedValidation</code> was never called.  Returns false and does not write the record if the record was not properly sorted.
 
  
−
Returns false if the record was not properly sorted or not successfully written.  Returns true if properly sorted and successfully written.
+
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
 +
|-style="background: #f2f2f2; text-align: center;"
 +
|+ '''Read Counts'''
 +
! Statistic !! Description
 +
|-
 +
|TotalReads
 +
| Total number of alignments that were successfully read from the file.
 
|-
 
|-
−
| <code>void SamFile::setSortedValidation(SortedType sortType)<\code>
+
|MappedReads
−
| Set the flag to validate that the file is sorted as it is read/written.  Must be called after the file has been opened.
+
| Total number of alignments that were successfully read from the file with FLAG bit 0x004 set to 0 (not unmapped).
−
sortType specifies the type of sort to be checked for.
 
 
|-
 
|-
−
| <code>uint32_t SamFile::GetCurrentRecordCount()</code>
+
|PairedReads
−
| Return the number of records that have been read/written so far.
+
| Total number of alignments that were successfully read from the file with FLAG bit 0x001 set to 1 (paired).
 
|-
 
|-
−
| <code>SamStatus::Status SamFile::GetFailure()</code>
+
|ProperPair
−
| Get the type of failure that occurred on a method failure.
+
| Total number of alignments that were successfully read from the file with FLAG bits 0x001 set to 1 (paired) AND 0x002 (proper pair).
 
|-
 
|-
−
| <code>bool SamFile::SetReadSection(int32_t refID)</code>
+
|DuplicateReads
−
| Tell the class which reference ID should be read from the BAM file.  This is the index into the BAM Index list of reference information: 0 - #references.  The records for that reference id will be retrieved on each <code>ReadRecord</code> call.  When all records have been retrieved for the specified reference id, <code>ReadRecord</code> will return false until a new read section is set.
+
| Total number of alignments that were successfully read from the file with FLAG bit 0x400 set to 1 (PCR or optical duplicate).
−
Pass in -1 in order to read the section of the bam file not associated with any reference ID.
+
|-
−
Returns true if the read section was successfully set, false if not.  False is returned if the BAM Index File has not yet been read or if a BAM file is not open for reading.
+
|QCFailureReads
 +
| Total number of alignments that were successfully read from the file with FLAG bit 0x200 set to 1 (failed quality checks).
 
|}
 
|}
  
−
=== Class Enums ===
 
 
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
 
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
−
|-style="background: #f2f2f2; text-align: center;" '''SamFile Class Methods'''
+
|-style="background: #f2f2f2; text-align: center;"
−
! colspan="2"| enum SortedType
+
|- '''Read Percentages'''
 +
! Statistic !! Description
 
|-
 
|-
−
! Enum Value !!  Description
+
|MappingRate(%)
 +
| 100 * MappedReads/TotalReads
 
|-
 
|-
−
| UNSORTED
+
|PairedReads(%)
−
| Do not do any sorting validation.
+
| 100 * PairedReads/TotalReads
 
|-
 
|-
−
| FLAG
+
|ProperPair(%)
−
| Validate that the file is sorted by the type specified in the SO Tag of the file's header.
+
| 100 * ProperPair/TotalReads
 
|-
 
|-
−
| COORDINATE
+
|DupRate(%)
−
| Validate that the file is sorted by Coordinate.
+
| 100 * DuplicateReads/TotalReads
 
|-
 
|-
−
| QUERY_NAME
+
|QCFailRate(%)
−
| Validate that the file is sorted by Query Name.
+
| 100 * QCFailureReads/TotalReads
 
|}
 
|}
  
−
=== Usage Example ===
+
{| style="margin: 1em 1em 1em 0; background-color: #f9f9f9; border: 1px #aaa solid; border-collapse: collapse;" border="1"
−
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.
+
|-style="background: #f2f2f2; text-align: center;"
−
 
+
|- '''Base Counts'''
−
<source lang="cpp">
+
! Statistic !! Description
−
int main(int argc, char ** argv)
+
|-
−
{
+
|TotalBases
−
  if(argc != 3)
+
| Sum of the SEQ lengths for all alignments that were successfully read from the file.
−
  {
+
NOTE: Includes bases that are 'N'.
−
      printf("./bam <inputFile> <outputFile.sam/bam>\n");
+
|-
−
      exit(-1);
+
|BasesInMappedReads
−
  }
+
| Sum of the SEQ lengths for all alignments that were successfully read from the file with FLAG bit 0x004 set to 0 (not unmapped).
−
 
+
NOTE: Includes bases that are 'N'.
−
 
+
|}
−
  SamFile samIn;
 
−
     
 
−
  samIn.OpenForRead(argv[1]);
 
−
 
 
−
  SamFile samOut;
 
−
 
 
−
  samOut.OpenForWrite(argv[2]);
 
−
 
 
−
  // Read the sam header.
 
−
  SamFileHeader samHeader;
 
−
  samIn.ReadHeader(samHeader);
 
−
 
 
−
  samOut.WriteHeader(samHeader);
 
−
 
 
−
  // Read the first sam record.
 
−
  SamRecord samRecord;
 
−
 
 
−
  // Track if any of the records are invalid. Initialize to 0, it will
 
−
  // be set to -1 on an invalid record.
 
−
  int status = 0;
 
−
 
 
−
  // Keep reading records until the end of the file is reached.
 
−
  int numValidRecords = 0;
 
−
  while(!samIn.IsEOF())
 
−
  {
 
−
      if(samIn.ReadRecord(samHeader, samRecord) == true)
 
−
      {
 
−
        // Successfully read a record from the file, so check to see
 
−
        // if it is valid.
 
−
        if(samRecord.isValid(samHeader))
 
−
        {
 
−
            //  It is valid, so write it.
 
−
            numValidRecords++;
 
−
            samOut.WriteRecord(samHeader, samRecord);
 
−
        }
 
−
        else
 
−
        {
 
−
            // There is at least one invalid record, return -1.
 
−
            status = -1;
 
−
        }
 
−
      }
 
−
      else
 
−
      {
 
−
        // Failed to read the record.  Check to see if it failed due to
 
−
        // there being no more records, which is not a failure.  Any
 
−
        // other failure reason counts as a failure.
 
−
        if(samIn.GetFailure() != SamStatus::NO_MORE_RECS)
 
−
        {
 
−
            // Read failed, invalid.
 
−
            status = -1;
 
−
        }
 
−
      }
 
−
  }
 
−
  std::cout << "Number of records read = " <<
 
−
      samIn.GetCurrentRecordCount() << std::endl;
 
−
  std::cout << "Number of records written = " <<
 
−
      samOut.GetCurrentRecordCount() << std::endl;
 
−
  std::cout << "Number of valid records = " << numValidRecords << std::endl;
 
−
 
 
−
  return(status);
 
−
}
 
−
</source>
 
−
 
 
−
 
 
−
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">
 
−
int ReadIndexedBam(const char* inputFilename,
 
−
                  const char* outputFilename,
 
−
                  const char* indexFilename)
 
−
{
 
−
  SamFile samIn;
 
−
     
 
−
  samIn.OpenForRead(inputFilename);
 
−
  samIn.ReadBamIndex(indexFilename);
 
−
 
 
−
  SamFile samOut;
 
  
−
  samOut.OpenForWrite(outputFilename);
+
NOTE: If the TotalReads is greater than 10^6, then the Read Counts and Base Counts specify the total counts divided by 10^6. This is indicated in the output with a (e6) appended to the field name.
  
−
  // Read the sam header.
+
==== Example Statistics Output ====
−
  SamFileHeader samHeader;
+
<pre>
−
  samIn.ReadHeader(samHeader);
+
TotalReads(e6) 18.90
−
  // Write the sam header.
+
MappedReads(e6) 14.77
−
  samOut.WriteHeader(samHeader);
+
PairedReads(e6) 18.90
 +
ProperPair(e6) 11.28
 +
DuplicateReads(e6) 0.00
 +
QCFailureReads(e6) 0.00
  
−
  SamRecord samRecord;
+
MappingRate(%) 78.17
−
 
+
PairedReads(%) 100.00
−
  int numValidRecords = 0;
+
ProperPair(%) 59.68
−
  int numRecords = 0;
+
DupRate(%) 0.00
 +
QCFailRate(%) 0.00
  
−
  // Loop through each Reference.
+
TotalBases(e6) 699.30
−
  for(int i = -1; i < 23; i++)
+
BasesInMappedReads(e6) 546.67
−
  {
+
</pre>
−
      int numSectionRecords = 0;
 
−
      samIn.SetReadSection(i);
 
−
      // Keep reading records until they aren't read successfully.
 
−
      while(samIn.ReadRecord(samHeader, samRecord) == true)
 
−
      {
 
−
        numSectionRecords++;
 
−
        numRecords++;
 
−
        // Successfully read a record from the file, so check to see
 
−
        // if it is valid.
 
−
        if(samRecord.isValid(samHeader))
 
−
        {
 
−
            // It is valid, so write it.
 
−
            numValidRecords++;
 
−
            samOut.WriteRecord(samHeader, samRecord);
 
−
        }
 
−
      }
 
−
      std::cout << "Reference ID " << i << " has " << numSectionRecords
 
−
                << " records" << std::endl;
 
−
  }
 
−
     
 
−
  std::cout << "Number of records = " << numRecords << std::endl;
 
−
  std::cout << "Number of valid records = " << numValidRecords << std::endl;
 
  
−
  return(0);
+
== Usage Examples ==
−
}
+
[[Sam Library Usage Examples]]
−
</source>
 

Latest revision as of 11:03, 2 February 2017


Reading/Writing SAM/BAM Files In Your Program

The SamFile class allows a user to easily read/write a SAM/BAM file.

The SamFile class contains additional functionality that allows a user to read specific sections of sorted & indexed BAM files. In order take advantage of this capability, the index file must be read prior to setting the read section. This logic saves the time of having to read the entire file and takes advantage of the seeking capability of BGZF files.

Future Enhancements: Add the ability to read alignments that match a given start, end position for a specific reference sequence.

This class is part of C++ Library: libStatGen.

Class Documentation

See: http://csg.sph.umich.edu//mktrost/doxygen/current/classSamFile.html

Child Classes

SamFileReader

http://csg.sph.umich.edu//mktrost/doxygen/current/classSamFileReader.html

SamFileWriter

http://csg.sph.umich.edu//mktrost/doxygen/current/classSamFileWriter.html

Statistics

Statistic Generation

The following statistics can be optionally recorded when reading a SamFile by specifying SamFile::GenerateStatistics() and displayed with SamFile::PrintStatistics()

The statistics only reflect alignments that were successfully read from the BAM file. Alignments that failed to parse from the file are not reflected in the statistics, but alignments that are invalid for other reasons may show up in the statistics.

Read Counts
Statistic Description
TotalReads Total number of alignments that were successfully read from the file.
MappedReads Total number of alignments that were successfully read from the file with FLAG bit 0x004 set to 0 (not unmapped).
PairedReads Total number of alignments that were successfully read from the file with FLAG bit 0x001 set to 1 (paired).
ProperPair Total number of alignments that were successfully read from the file with FLAG bits 0x001 set to 1 (paired) AND 0x002 (proper pair).
DuplicateReads Total number of alignments that were successfully read from the file with FLAG bit 0x400 set to 1 (PCR or optical duplicate).
QCFailureReads Total number of alignments that were successfully read from the file with FLAG bit 0x200 set to 1 (failed quality checks).
Statistic Description
MappingRate(%) 100 * MappedReads/TotalReads
PairedReads(%) 100 * PairedReads/TotalReads
ProperPair(%) 100 * ProperPair/TotalReads
DupRate(%) 100 * DuplicateReads/TotalReads
QCFailRate(%) 100 * QCFailureReads/TotalReads
Statistic Description
TotalBases Sum of the SEQ lengths for all alignments that were successfully read from the file.

NOTE: Includes bases that are 'N'.

BasesInMappedReads Sum of the SEQ lengths for all alignments that were successfully read from the file with FLAG bit 0x004 set to 0 (not unmapped).

NOTE: Includes bases that are 'N'.

NOTE: If the TotalReads is greater than 10^6, then the Read Counts and Base Counts specify the total counts divided by 10^6. This is indicated in the output with a (e6) appended to the field name.

Example Statistics Output

TotalReads(e6)	18.90
MappedReads(e6)	14.77
PairedReads(e6)	18.90
ProperPair(e6)	11.28
DuplicateReads(e6)	0.00
QCFailureReads(e6)	0.00

MappingRate(%)	78.17
PairedReads(%)	100.00
ProperPair(%)	59.68
DupRate(%)	0.00
QCFailRate(%)	0.00

TotalBases(e6)	699.30
BasesInMappedReads(e6)	546.67

Usage Examples

Sam Library Usage Examples

Retrieved from "http://genome.sph.umich.edu/w/index.php?title=C%2B%2B_Class:_SamFile&oldid=14607"