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 // ---------------------------------------------------------------------------
   5 // Last modified: 7 April 2011 (DB)
   6 // ---------------------------------------------------------------------------
   7 // Grab a random subset of alignments (testing tool)
   8 // ***************************************************************************
   9
  10 #include "bamtools_random.h"
  11
  12 #include <api/BamMultiReader.h>
  13 #include <api/BamWriter.h>
  14 #include <utils/bamtools_options.h>
  15 #include <utils/bamtools_utilities.h>
  16 using namespace BamTools;
  17
  18 #include <ctime>
  19 #include <cstdlib>
  20 #include <iostream>
  21 #include <string>
  22 #include <vector>
  23 using namespace std;
  24
  25 namespace BamTools {
  26
  27 // define constants
  28 const unsigned int RANDOM_MAX_ALIGNMENT_COUNT = 10000;
  29
  30 // utility methods for RandomTool
  31 int getRandomInt(const int& lowerBound, const int& upperBound) {
  32     const int range = (upperBound - lowerBound) + 1;
  33     return ( lowerBound + (int)(range * (double)rand()/((double)RAND_MAX + 1)) );
  34 }
  35
  36 } // namespace BamTools
  37
  38 // ---------------------------------------------
  39 // RandomSettings implementation
  40
  41 struct RandomTool::RandomSettings {
  42
  43     // flags
  44     bool HasAlignmentCount;
  45     bool HasInput;
  46     bool HasOutput;
  47     bool HasRegion;
  48     bool IsForceCompression;
  49
  50     // parameters
  51     unsigned int AlignmentCount;
  52     vector<string> InputFiles;
  53     string OutputFilename;
  54     string Region;
  55
  56     // constructor
  57     RandomSettings(void)
  58         : HasAlignmentCount(false)
  59         , HasInput(false)
  60         , HasOutput(false)
  61         , HasRegion(false)
  62         , IsForceCompression(false)
  63         , AlignmentCount(RANDOM_MAX_ALIGNMENT_COUNT)
  64         , OutputFilename(Options::StandardOut())
  65     { }
  66 };
  67
  68 // ---------------------------------------------
  69 // RandomToolPrivate implementation
  70
  71 struct RandomTool::RandomToolPrivate {
  72
  73     // ctor & dtor
  74     public:
  75         RandomToolPrivate(RandomTool::RandomSettings* settings)
  76             : m_settings(settings)
  77         { }
  78
  79         ~RandomToolPrivate(void) { }
  80
  81     // interface
  82     public:
  83         bool Run(void);
  84
  85     // data members
  86     private:
  87         RandomTool::RandomSettings* m_settings;
  88 };
  89
  90 bool RandomTool::RandomToolPrivate::Run(void) {
  91
  92     // set to default stdin if no input files provided
  93     if ( !m_settings->HasInput )
  94         m_settings->InputFiles.push_back(Options::StandardIn());
  95
  96     // open our reader
  97     BamMultiReader reader;
  98     if ( !reader.Open(m_settings->InputFiles) ) {
  99         cerr << "bamtools random ERROR: could not open input BAM file(s)... Aborting." << endl;
 100         return false;
 101     }
 102
 103     // look up index files for all BAM files
 104     reader.LocateIndexes();
 105
 106     // make sure index data is available
 107     if ( !reader.HasIndexes() ) {
 108         cerr << "bamtools random ERROR: could not load index data for all input BAM file(s)... Aborting." << endl;
 109         reader.Close();
 110         return false;
 111     }
 112
 113     // get BamReader metadata
 114     const string headerText = reader.GetHeaderText();
 115     const RefVector references = reader.GetReferenceData();
 116     if ( references.empty() ) {
 117         cerr << "bamtools random ERROR: no reference data available... Aborting." << endl;
 118         reader.Close();
 119         return false;
 120     }
 121
 122     // determine compression mode for BamWriter
 123     bool writeUncompressed = ( m_settings->OutputFilename == Options::StandardOut() &&
 124                               !m_settings->IsForceCompression );
 125     BamWriter::CompressionMode compressionMode = BamWriter::Compressed;
 126     if ( writeUncompressed ) compressionMode = BamWriter::Uncompressed;
 127
 128     // open BamWriter
 129     BamWriter writer;
 130     writer.SetCompressionMode(compressionMode);
 131     if ( !writer.Open(m_settings->OutputFilename, headerText, references) ) {
 132         cerr << "bamtools random ERROR: could not open " << m_settings->OutputFilename
 133              << " for writing... Aborting." << endl;
 134         reader.Close();
 135         return false;
 136     }
 137
 138     // if user specified a REGION constraint, attempt to parse REGION string
 139     BamRegion region;
 140     if ( m_settings->HasRegion && !Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
 141         cerr << "bamtools random ERROR: could not parse REGION: " << m_settings->Region << endl;
 142         cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
 143              << endl;
 144         reader.Close();
 145         writer.Close();
 146         return false;
 147     }
 148
 149     // seed our random number generator
 150     srand( time(NULL) );
 151
 152     // grab random alignments
 153     BamAlignment al;
 154     unsigned int i = 0;
 155     while ( i < m_settings->AlignmentCount ) {
 156
 157         int randomRefId    = 0;
 158         int randomPosition = 0;
 159
 160         // use REGION constraints to select random refId & position
 161         if ( m_settings->HasRegion ) {
 162
 163             // select a random refId
 164             randomRefId = getRandomInt(region.LeftRefID, region.RightRefID);
 165
 166             // select a random position based on randomRefId
 167             const int lowerBoundPosition = ( (randomRefId == region.LeftRefID)
 168                                              ? region.LeftPosition
 169                                              : 0 );
 170             const int upperBoundPosition = ( (randomRefId == region.RightRefID)
 171                                              ? region.RightPosition
 172                                              : (references.at(randomRefId).RefLength - 1) );
 173             randomPosition = getRandomInt(lowerBoundPosition, upperBoundPosition);
 174         }
 175
 176         // otherwise select from all possible random refId & position
 177         else {
 178
 179             // select random refId
 180             randomRefId = getRandomInt(0, (int)references.size() - 1);
 181
 182             // select random position based on randomRefId
 183             const int lowerBoundPosition = 0;
 184             const int upperBoundPosition = references.at(randomRefId).RefLength - 1;
 185             randomPosition = getRandomInt(lowerBoundPosition, upperBoundPosition);
 186         }
 187
 188         // if jump & read successful, save first alignment that overlaps random refId & position
 189         if ( reader.Jump(randomRefId, randomPosition) ) {
 190             while ( reader.GetNextAlignmentCore(al) ) {
 191                 if ( al.RefID == randomRefId && al.Position >= randomPosition ) {
 192                     writer.SaveAlignment(al);
 193                     ++i;
 194                     break;
 195                 }
 196             }
 197         }
 198     }
 199
 200     // cleanup & exit
 201     reader.Close();
 202     writer.Close();
 203     return true;
 204 }
 205
 206 // ---------------------------------------------
 207 // RandomTool implementation
 208
 209 RandomTool::RandomTool(void)
 210     : AbstractTool()
 211     , m_settings(new RandomSettings)
 212     , m_impl(0)
 213 {
 214     // set program details
 215     Options::SetProgramInfo("bamtools random", "grab a random subset of alignments", "[-in <filename> -in <filename> ...] [-out <filename>] [-forceCompression] [-n] [-region <REGION>]");
 216
 217     // set up options
 218     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
 219     Options::AddValueOption("-in",  "BAM filename", "the input BAM file",  "", m_settings->HasInput,  m_settings->InputFiles,     IO_Opts, Options::StandardIn());
 220     Options::AddValueOption("-out", "BAM filename", "the output BAM file", "", m_settings->HasOutput, m_settings->OutputFilename, IO_Opts, Options::StandardOut());
 221     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);
 222     Options::AddValueOption("-region", "REGION", "only pull random alignments from within this genomic region. Index file is recommended for better performance, and is used automatically if it exists. See \'bamtools help index\' for more details on creating one", "", m_settings->HasRegion, m_settings->Region, IO_Opts);
 223
 224     OptionGroup* SettingsOpts = Options::CreateOptionGroup("Settings");
 225     Options::AddValueOption("-n", "count", "number of alignments to grab. Note - no duplicate checking is performed", "", m_settings->HasAlignmentCount, m_settings->AlignmentCount, SettingsOpts, RANDOM_MAX_ALIGNMENT_COUNT);
 226 }
 227
 228 RandomTool::~RandomTool(void) {
 229
 230     delete m_settings;
 231     m_settings = 0;
 232
 233     delete m_impl;
 234     m_impl = 0;
 235 }
 236
 237 int RandomTool::Help(void) {
 238     Options::DisplayHelp();
 239     return 0;
 240 }
 241
 242 int RandomTool::Run(int argc, char* argv[]) {
 243
 244     // parse command line arguments
 245     Options::Parse(argc, argv, 1);
 246
 247     // initialize RandomTool with settings
 248     m_impl = new RandomToolPrivate(m_settings);
 249
 250     // run RandomTool, return success/fail
 251     if ( m_impl->Run() )
 252         return 0;
 253     else
 254         return 1;
 255 }