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