3,045
edits
Changes
From Genome Analysis Wiki
Jump to navigationJump to searchno edit summary
|}
|}
=== Usage Example ===
=== Usage Examples ===
==== 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.
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.
==== 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.
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">
#include "SamFile.h"
#include "SamValidation.h"
int ReadIndexedBam(const char* inputFilename,
int ReadIndexedBam(const char* inputFilename,
const char* outputFilename,
const char* outputFilename,
return(0);
return(0);
}
</source>
==== 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.
<source>
#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);
}
}
</source>
</source>
