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: 20 July 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 } // namespace BamTools
  30
  31 // ---------------------------------------------
  32 // RandomSettings implementation
  33
  34 struct RandomTool::RandomSettings {
  35
  36     // flags
  37     bool HasAlignmentCount;
  38     bool HasInput;
  39     bool HasOutput;
  40     bool HasRegion;
  41
  42     // parameters
  43     unsigned int AlignmentCount;
  44     vector<string> InputFiles;
  45     string OutputFilename;
  46     string Region;
  47
  48     // constructor
  49     RandomSettings(void)
  50         : HasAlignmentCount(false)
  51         , HasInput(false)
  52         , HasOutput(false)
  53         , HasRegion(false)
  54         , AlignmentCount(RANDOM_MAX_ALIGNMENT_COUNT)
  55     { }
  56 };
  57
  58 // ---------------------------------------------
  59 // RandomTool implementation
  60
  61 RandomTool::RandomTool(void)
  62     : AbstractTool()
  63     , m_settings(new RandomSettings)
  64 {
  65     // set program details
  66     Options::SetProgramInfo("bamtools random", "grab a random subset of alignments", "[-in <filename> -in <filename> ...] [-out <filename>] [-region <REGION>]");
  67
  68     // set up options
  69     OptionGroup* IO_Opts = Options::CreateOptionGroup("Input & Output");
  70     Options::AddValueOption("-in",  "BAM filename", "the input BAM file",  "", m_settings->HasInput,  m_settings->InputFiles,     IO_Opts, Options::StandardIn());
  71     Options::AddValueOption("-out", "BAM filename", "the output BAM file", "", m_settings->HasOutput, m_settings->OutputFilename, IO_Opts, Options::StandardOut());
  72
  73     OptionGroup* FilterOpts = Options::CreateOptionGroup("Filters");
  74     Options::AddValueOption("-n", "count", "number of alignments to grab.  Note - no duplicate checking is performed (currently)", "", m_settings->HasAlignmentCount, m_settings->AlignmentCount, FilterOpts, RANDOM_MAX_ALIGNMENT_COUNT);
  75     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);
  76 }
  77
  78 RandomTool::~RandomTool(void) {
  79     delete m_settings;
  80     m_settings = 0;
  81 }
  82
  83 int RandomTool::Help(void) {
  84     Options::DisplayHelp();
  85     return 0;
  86 }
  87
  88 int RandomTool::Run(int argc, char* argv[]) {
  89
  90     // TODO: Handle BAM input WITHOUT index files.
  91
  92     // parse command line arguments
  93     Options::Parse(argc, argv, 1);
  94
  95     // set to default input if none provided
  96     if ( !m_settings->HasInput )
  97         m_settings->InputFiles.push_back(Options::StandardIn());
  98
  99     // open our BAM reader
 100     BamMultiReader reader;
 101     reader.Open(m_settings->InputFiles);
 102     string headerText    = reader.GetHeaderText();
 103     RefVector references = reader.GetReferenceData();
 104
 105     // check that reference data is available, used for generating random jumps
 106     if ( references.empty() ) {
 107         cerr << "No reference data available... quitting." << endl;
 108         reader.Close();
 109         return 1;
 110     }
 111
 112     // see if user specified a REGION
 113     BamRegion region;
 114     if ( m_settings->HasRegion ) {
 115         if ( Utilities::ParseRegionString(m_settings->Region, reader, region) )
 116             reader.SetRegion(region);
 117     }
 118
 119     // open out BAM writer
 120     BamWriter writer;
 121     writer.Open(m_settings->OutputFilename, headerText, references);
 122
 123     // seed our random number generator
 124     srand (time(NULL) );
 125
 126     // grab random alignments
 127     BamAlignment al;
 128     unsigned int i = 0;
 129     while ( i < m_settings->AlignmentCount ) {
 130
 131         int randomRefId    = 0;
 132         int randomPosition = 0;
 133
 134         // use REGION constraints to generate random refId & position
 135         if ( m_settings->HasRegion ) {
 136
 137             int lowestRefId  = region.LeftRefID;
 138             int highestRefId = region.RightRefID;
 139             int rangeRefId   = (highestRefId - lowestRefId) + 1;
 140             randomRefId = lowestRefId + (int)(rangeRefId * (double)(rand()/((double)RAND_MAX + 1)));
 141
 142             int lowestPosition  = ( (randomRefId == region.LeftRefID)  ? region.LeftPosition  : 0 );
 143             int highestPosition = ( (randomRefId == region.RightRefID) ? region.RightPosition : references.at(randomRefId).RefLength - 1 );
 144             int rangePosition = (highestPosition - lowestPosition) + 1;
 145             randomPosition = lowestPosition + (int)(rangePosition * (double)(rand()/((double)RAND_MAX + 1)));
 146         }
 147
 148         // otherwise generate 'normal' random refId & position
 149         else {
 150
 151             // generate random refId
 152             int lowestRefId = 0;
 153             int highestRefId = references.size() - 1;
 154             int rangeRefId = (highestRefId - lowestRefId) + 1;
 155             randomRefId = lowestRefId + (int)(rangeRefId * (double)(rand()/((double)RAND_MAX + 1)));
 156
 157             // generate random position
 158             int lowestPosition = 0;
 159             int highestPosition = references.at(randomRefId).RefLength - 1;
 160             int rangePosition = (highestPosition - lowestPosition) + 1;
 161             randomPosition = lowestPosition + (int)(rangePosition * (double)(rand()/((double)RAND_MAX + 1)));
 162         }
 163
 164         // if jump & read successful, save alignment
 165         if ( reader.Jump(randomRefId, randomPosition) ) {
 166             while ( reader.GetNextAlignmentCore(al) ) {
 167                 if ( al.RefID == randomRefId && al.Position >= randomPosition ) {
 168                     writer.SaveAlignment(al);
 169                     ++i;
 170                     break;
 171                 }
 172             }
 173         }
 174     }
 175
 176     // close reader & writer
 177     reader.Close();
 178     writer.Close();
 179     return 0;
 180 }