src/toolkit/bamtools_random.cpp

   1 // ***************************************************************************
   2 // bamtools_random.cpp (c) 2010 Derek Barnett, Erik Garrison
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 3 September 2010 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Grab a random subset of alignments.
   9 // ***************************************************************************
  10
  11 #include <ctime>
  12 #include <cstdlib>
  13 #include <iostream>
  14 #include <string>
  15 #include <vector>
  16 #include "bamtools_random.h"
  17 #include "bamtools_options.h"
  18 #include "bamtools_utilities.h"
  19 #include "BamMultiReader.h"
  20 #include "BamWriter.h"
  21 using namespace std;
  22 using namespace BamTools;
  23
  24 namespace BamTools {
  25
  26 // define constants
  27 const unsigned int RANDOM_MAX_ALIGNMENT_COUNT = 10000;
  28
  29 // utility methods for RandomTool
  30 const int getRandomInt(const int& lowerBound, const int& upperBound) {
  31     const int range = (upperBound - lowerBound) + 1;
  32     return ( lowerBound + (int)(range * (double)rand()/((double)RAND_MAX + 1)) );
  33 }
  34
  35 } // namespace BamTools
  36
  37 // ---------------------------------------------
  38 // RandomSettings implementation
  39
  40 struct RandomTool::RandomSettings {
  41
  42     // flags
  43     bool HasAlignmentCount;
  44     bool HasInput;
  45     bool HasOutput;
  46     bool HasRegion;
  47     bool IsForceCompression;
  48
  49     // parameters
  50     unsigned int AlignmentCount;
  51     vector<string> InputFiles;
  52     string OutputFilename;
  53     string Region;
  54
  55     // constructor
  56     RandomSettings(void)
  57         : HasAlignmentCount(false)
  58         , HasInput(false)
  59         , HasOutput(false)
  60         , HasRegion(false)
  61         , IsForceCompression(false)
  62         , AlignmentCount(RANDOM_MAX_ALIGNMENT_COUNT)
  63         , OutputFilename(Options::StandardOut())
  64     { }
  65 };
  66
  67 // ---------------------------------------------
  68 // RandomTool implementation
  69
  70 RandomTool::RandomTool(void)
  71     : AbstractTool()
  72     , m_settings(new RandomSettings)
  73 {
  74     // set program details
  75     Options::SetProgramInfo("bamtools random", "grab a random subset of alignments", "[-in <filename> -in <filename> ...] [-out <filename>] [-forceCompression] [-n] [-region <REGION>]");
  76
  77     // set up options
  78     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
  79     Options::AddValueOption("-in",  "BAM filename", "the input BAM file",  "", m_settings->HasInput,  m_settings->InputFiles,     IO_Opts, Options::StandardIn());
  80     Options::AddValueOption("-out", "BAM filename", "the output BAM file", "", m_settings->HasOutput, m_settings->OutputFilename, IO_Opts, Options::StandardOut());
  81     Options::AddOption("-forceCompression", "if results are sent to stdout (like when piping to another tool), default behavior is to leave output uncompressed. Use this flag to override and force compression", m_settings->IsForceCompression, IO_Opts);
  82
  83     OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
  84     Options::AddValueOption("-n", "count", "number of alignments to grab.  Note - no duplicate checking is performed", "", m_settings->HasAlignmentCount, m_settings->AlignmentCount, FilterOpts, RANDOM_MAX_ALIGNMENT_COUNT);
  85     Options::AddValueOption("-region", "REGION", "limit source of random alignment subset to a particular genomic region. Index file is recommended for better performance, and is read automatically if it exists as <filename>.bai or <filename>.bti. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, FilterOpts);
  86 }
  87
  88 RandomTool::~RandomTool(void) {
  89     delete m_settings;
  90     m_settings = 0;
  91 }
  92
  93 int RandomTool::Help(void) {
  94     Options::DisplayHelp();
  95     return 0;
  96 }
  97
  98 int RandomTool::Run(int argc, char* argv[]) {
  99
 100     // parse command line arguments
 101     Options::Parse(argc, argv, 1);
 102
 103     // set to default stdin if no input files provided
 104     if ( !m_settings->HasInput )
 105         m_settings->InputFiles.push_back(Options::StandardIn());
 106
 107     // open our reader
 108     BamMultiReader reader;
 109     if ( !reader.Open(m_settings->InputFiles) ) {
 110         cerr << "ERROR: Could not open input BAM file(s)." << endl;
 111         return 1;
 112     }
 113
 114     // make sure index data is available
 115     if ( !reader.IsIndexLoaded() ) {
 116         cerr << "ERROR: Could not load index data for all input BAM file(s)." << endl;
 117         cerr << "\'bamtools random\' requires valid index files to provide efficient performance." << endl;
 118         reader.Close();
 119         return 1;
 120
 121     }
 122
 123     // get BamReader metadata
 124     const string headerText = reader.GetHeaderText();
 125     const RefVector references = reader.GetReferenceData();
 126     if ( references.empty() ) {
 127         cerr << "ERROR: No reference data available - required to perform random access throughtout input file(s)." << endl;
 128         reader.Close();
 129         return 1;
 130     }
 131
 132     // open our writer
 133     BamWriter writer;
 134     bool writeUncompressed = ( m_settings->OutputFilename == Options::StandardOut() && !m_settings->IsForceCompression );
 135     if ( !writer.Open(m_settings->OutputFilename, headerText, references, writeUncompressed) ) {
 136         cerr << "ERROR: Could not open BamWriter." << endl;
 137         reader.Close();
 138         return 1;
 139     }
 140
 141     // if user specified a REGION constraint, attempt to parse REGION string
 142     BamRegion region;
 143     if ( m_settings->HasRegion && !Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 144         cerr << "ERROR: Could not parse REGION: " << m_settings->Region << endl;
 145         cerr << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
 146         reader.Close();
 147         writer.Close();
 148         return 1;
 149     }
 150
 151     // seed our random number generator
 152     srand( time(NULL) );
 153
 154     // grab random alignments
 155     BamAlignment al;
 156     unsigned int i = 0;
 157     while ( i < m_settings->AlignmentCount ) {
 158
 159         int randomRefId    = 0;
 160         int randomPosition = 0;
 161
 162         // use REGION constraints to select random refId & position
 163         if ( m_settings->HasRegion ) {
 164
 165             // select a random refId
 166             randomRefId = getRandomInt(region.LeftRefID, region.RightRefID);
 167
 168             // select a random position based on randomRefId
 169             const int lowerBoundPosition = ( (randomRefId == region.LeftRefID)  ? region.LeftPosition  : 0 );
 170             const int upperBoundPosition = ( (randomRefId == region.RightRefID) ? region.RightPosition : (references.at(randomRefId).RefLength - 1) );
 171             randomPosition = getRandomInt(lowerBoundPosition, upperBoundPosition);
 172         }
 173
 174         // otherwise select from all possible random refId & position
 175         else {
 176
 177             // select random refId
 178             randomRefId = getRandomInt(0, (int)references.size() - 1);
 179
 180             // select random position based on randomRefId
 181             const int lowerBoundPosition = 0;
 182             const int upperBoundPosition = references.at(randomRefId).RefLength - 1;
 183             randomPosition = getRandomInt(lowerBoundPosition, upperBoundPosition);
 184         }
 185
 186         // if jump & read successful, save first alignment that overlaps random refId & position
 187         if ( reader.Jump(randomRefId, randomPosition) ) {
 188             while ( reader.GetNextAlignmentCore(al) ) {
 189                 if ( al.RefID == randomRefId && al.Position >= randomPosition ) {
 190                     writer.SaveAlignment(al);
 191                     ++i;
 192                     break;
 193                 }
 194             }
 195         }
 196     }
 197
 198     // cleanup & exit
 199     reader.Close();
 200     writer.Close();
 201     return 0;
 202 }