Changes

From Genome Analysis Wiki
Jump to navigationJump to search

Sam Library Usage Examples (view source)

Revision as of 10:42, 2 August 2010

4,961 bytes added ,  10:42, 2 August 2010
no edit summary
Line 90: Line 90:       −
== 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 173: Line 173:  
}
 
}
 
</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. 
 +
    if(!samIn.OpenForRead(inFile))
 +
    {
 +
        fprintf(stderr, "%s\n", samIn.GetStatusMessage());
 +
        return(samIn.GetStatus());
 +
    }
 +
 +
    // If refName is set, use that.
 +
    if(refName.Length() != 0)
 +
    {
 +
        // Use Reference Name.
 +
        if(!samIn.SetReadSection(refName.c_str(), start, end))
 +
        {
 +
            fprintf(stderr, "%s\n", samIn.GetStatusMessage());
 +
            return(samIn.GetStatus());
 +
        }
 +
    }
 +
    else
 +
    {
 +
        // Use Reference ID
 +
        if(!samIn.SetReadSection(refID, start, end))
 +
        {
 +
            fprintf(stderr, "%s\n", samIn.GetStatusMessage());
 +
            return(samIn.GetStatus());
 +
        }
 +
    }
 +
 +
    // Open the output file for writing.
 +
    SamFile samOut;
 +
    if(!samOut.OpenForWrite(outFile))
 +
    {
 +
        fprintf(stderr, "%s\n", samOut.GetStatusMessage());
 +
        return(samOut.GetStatus());
 +
    }
 +
 +
    // Open the bam index file for reading.
 +
    if(!samIn.ReadBamIndex(indexFile))
 +
    {
 +
        fprintf(stderr, "%s\n", samIn.GetStatusMessage());
 +
        return(samIn.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());   
 +
    }
 +
 +
    // Read the sam records.
 +
    SamRecord samRecord;
 +
    // Track the status.
 +
    int numSectionRecords = 0;
 +
 +
    // 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 they aren't anymore.
 +
    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();
 +
        }
 +
        ++numSectionRecords;
 +
    }
 +
 +
    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 << "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 ==
Retrieved from "http://genome.sph.umich.edu/wiki/Special:MobileDiff/1892"

Navigation menu