src/api/internal/SamHeaderValidator_p.cpp

   1 // ***************************************************************************
   2 // SamHeaderValidator.cpp (c) 2010 Derek Barnett
   3 // Marth Lab, Department of Biology, Boston College
   4 // All rights reserved.
   5 // ---------------------------------------------------------------------------
   6 // Last modified: 13 January 2011 (DB)
   7 // ---------------------------------------------------------------------------
   8 // Provides functionality for validating SamHeader data
   9 // ***************************************************************************
  10
  11 #include <api/SamConstants.h>
  12 #include <api/SamHeader.h>
  13 #include <api/internal/SamHeaderValidator_p.h>
  14 #include <api/internal/SamHeaderVersion_p.h>
  15 using namespace BamTools;
  16 using namespace BamTools::Internal;
  17
  18 #include <iostream>
  19 #include <set>
  20 #include <sstream>
  21 using namespace std;
  22
  23 // ------------------------------------------------------------------------
  24 // Allow validation rules to vary, as needed, between SAM header versions
  25 //
  26 // use SAM_VERSION_X_Y to tag important changes
  27 //
  28 // Together, they will allow for comparisons like:
  29 // if ( m_version < SAM_VERSION_2_0 ) {
  30 //     // use some older rule
  31 // else
  32 //     // use rule introduced with version 2.0
  33
  34 static const SamHeaderVersion SAM_VERSION_1_0 = SamHeaderVersion(1,0);
  35 static const SamHeaderVersion SAM_VERSION_1_3 = SamHeaderVersion(1,3);
  36
  37 // TODO: This functionality is currently unused.
  38 //       Make validation "version-aware."
  39 //
  40 // ------------------------------------------------------------------------
  41
  42 const string SamHeaderValidator::ERROR_PREFIX = "ERROR: ";
  43 const string SamHeaderValidator::WARN_PREFIX  = "WARNING: ";
  44 const string SamHeaderValidator::NEWLINE      = "\n";
  45
  46 SamHeaderValidator::SamHeaderValidator(const SamHeader& header)
  47     : m_header(header)
  48 { }
  49
  50 SamHeaderValidator::~SamHeaderValidator(void) { }
  51
  52 bool SamHeaderValidator::Validate(bool verbose) {
  53
  54     // validate header components
  55     bool isValid = true;
  56     isValid &= ValidateMetadata();
  57     isValid &= ValidateSequenceDictionary();
  58     isValid &= ValidateReadGroupDictionary();
  59     isValid &= ValidateProgramData();
  60
  61     // report errors if desired
  62     if ( verbose ) {
  63         PrintErrorMessages();
  64         PrintWarningMessages();
  65     }
  66
  67     // return validation status
  68     return isValid;
  69 }
  70
  71 bool SamHeaderValidator::ValidateMetadata(void) {
  72     bool isValid = true;
  73     isValid &= ValidateVersion();
  74     isValid &= ValidateSortOrder();
  75     isValid &= ValidateGroupOrder();
  76     return isValid;
  77 }
  78
  79 bool SamHeaderValidator::ValidateVersion(void) {
  80
  81     const string& version = m_header.Version;
  82
  83     // warn if version not present
  84     if ( version.empty() ) {
  85         AddWarning("Version (VN) missing. Not required, but strongly recommended");
  86         return true;
  87     }
  88
  89     // invalid if version does not contain a period
  90     const size_t periodFound = version.find(Constants::SAM_PERIOD);
  91     if ( periodFound == string::npos ) {
  92         AddError("Invalid version (VN) format: " + version);
  93         return false;
  94     }
  95
  96     // invalid if major version is empty or contains non-digits
  97     const string majorVersion = version.substr(0, periodFound);
  98     if ( majorVersion.empty() || !ContainsOnlyDigits(majorVersion) ) {
  99         AddError("Invalid version (VN) format: " + version);
 100         return false;
 101     }
 102
 103     // invalid if major version is empty or contains non-digits
 104     const string minorVersion = version.substr(periodFound + 1);
 105     if ( minorVersion.empty() || !ContainsOnlyDigits(minorVersion) ) {
 106         AddError("Invalid version (VN) format: " + version);
 107         return false;
 108     }
 109
 110     // TODO: check if version is not just syntactically OK,
 111     // but is also a valid SAM version ( 1.0 .. CURRENT )
 112
 113     // all checked out this far, then version is OK
 114     return true;
 115 }
 116
 117 // assumes non-empty input string
 118 bool SamHeaderValidator::ContainsOnlyDigits(const string& s) {
 119     const size_t nonDigitPosition = s.find_first_not_of(Constants::SAM_DIGITS);
 120     return ( nonDigitPosition == string::npos ) ;
 121 }
 122
 123 bool SamHeaderValidator::ValidateSortOrder(void) {
 124
 125     const string& sortOrder = m_header.SortOrder;
 126
 127     // warn if sort order not present
 128     if ( sortOrder.empty() ) {
 129         AddWarning("Sort order (SO) missing. Not required, but strongly recommended");
 130         return true;
 131     }
 132
 133     // if sort order is valid keyword
 134     if ( sortOrder == Constants::SAM_HD_SORTORDER_COORDINATE ||
 135          sortOrder == Constants::SAM_HD_SORTORDER_QUERYNAME  ||
 136          sortOrder == Constants::SAM_HD_SORTORDER_UNSORTED
 137        )
 138     { return true; }
 139
 140     // otherwise
 141     AddError("Invalid sort order (SO): " + sortOrder);
 142     return false;
 143 }
 144
 145 bool SamHeaderValidator::ValidateGroupOrder(void) {
 146
 147     const string& groupOrder = m_header.GroupOrder;
 148
 149     // if no group order, no problem, just return OK
 150     if ( groupOrder.empty() ) return true;
 151
 152     // if group order is valid keyword
 153     if ( groupOrder == Constants::SAM_HD_GROUPORDER_NONE  ||
 154          groupOrder == Constants::SAM_HD_GROUPORDER_QUERY ||
 155          groupOrder == Constants::SAM_HD_GROUPORDER_REFERENCE
 156        )
 157     { return true; }
 158
 159     // otherwise
 160     AddError("Invalid group order (GO): " + groupOrder);
 161     return false;
 162 }
 163
 164 bool SamHeaderValidator::ValidateSequenceDictionary(void) {
 165
 166     // TODO: warn/error if no sequences ?
 167
 168     bool isValid = true;
 169
 170     // check for unique sequence names
 171     isValid &= ContainsUniqueSequenceNames();
 172
 173     // iterate over sequences
 174     const SamSequenceDictionary& sequences = m_header.Sequences;
 175     SamSequenceConstIterator seqIter = sequences.ConstBegin();
 176     SamSequenceConstIterator seqEnd  = sequences.ConstEnd();
 177     for ( ; seqIter != seqEnd; ++seqIter ) {
 178         const SamSequence& seq = (*seqIter);
 179         isValid &= ValidateSequence(seq);
 180     }
 181
 182     // return validation state
 183     return isValid;
 184 }
 185
 186 bool SamHeaderValidator::ContainsUniqueSequenceNames(void) {
 187
 188     bool isValid = true;
 189     set<string> sequenceNames;
 190     set<string>::iterator nameIter;
 191
 192     // iterate over sequences
 193     const SamSequenceDictionary& sequences = m_header.Sequences;
 194     SamSequenceConstIterator seqIter = sequences.ConstBegin();
 195     SamSequenceConstIterator seqEnd  = sequences.ConstEnd();
 196     for ( ; seqIter != seqEnd; ++seqIter ) {
 197         const SamSequence& seq = (*seqIter);
 198         const string& name = seq.Name;
 199
 200         // lookup sequence name
 201         nameIter = sequenceNames.find(name);
 202
 203         // error if found (duplicate entry)
 204         if ( nameIter != sequenceNames.end() ) {
 205             AddError("Sequence name (SN): " + name + " is not unique");
 206             isValid = false;
 207         }
 208
 209         // otherwise ok, store name
 210         sequenceNames.insert(name);
 211     }
 212
 213     // return validation state
 214     return isValid;
 215 }
 216
 217 bool SamHeaderValidator::ValidateSequence(const SamSequence& seq) {
 218     bool isValid = true;
 219     isValid &= CheckNameFormat(seq.Name);
 220     isValid &= CheckLengthInRange(seq.Length);
 221     return isValid;
 222 }
 223
 224 bool SamHeaderValidator::CheckNameFormat(const string& name) {
 225
 226     // invalid if name is empty
 227     if ( name.empty() ) {
 228         AddError("Sequence entry (@SQ) is missing SN tag");
 229         return false;
 230     }
 231
 232     // invalid if first character is a reserved char
 233     const char firstChar = name.at(0);
 234     if ( firstChar == Constants::SAM_EQUAL || firstChar == Constants::SAM_STAR ) {
 235         AddError("Invalid sequence name (SN): " + name);
 236         return false;
 237     }
 238     // otherwise OK
 239     return true;
 240 }
 241
 242 bool SamHeaderValidator::CheckLengthInRange(const string& length) {
 243
 244     // invalid if empty
 245     if ( length.empty() ) {
 246         AddError("Sequence entry (@SQ) is missing LN tag");
 247         return false;
 248     }
 249
 250     // convert string length to numeric
 251     stringstream lengthStream(length);
 252     unsigned int sequenceLength;
 253     lengthStream >> sequenceLength;
 254
 255     // invalid if length outside accepted range
 256     if ( sequenceLength < Constants::SAM_SQ_LENGTH_MIN || sequenceLength > Constants::SAM_SQ_LENGTH_MAX ) {
 257         AddError("Sequence length (LN): " + length + " out of range");
 258         return false;
 259     }
 260
 261     // otherwise OK
 262     return true;
 263 }
 264
 265 bool SamHeaderValidator::ValidateReadGroupDictionary(void) {
 266
 267     // TODO: warn/error if no read groups ?
 268
 269     bool isValid = true;
 270
 271     // check for unique read group IDs & platform units
 272     isValid &= ContainsUniqueIDsAndPlatformUnits();
 273
 274     // iterate over read groups
 275     const SamReadGroupDictionary& readGroups = m_header.ReadGroups;
 276     SamReadGroupConstIterator rgIter = readGroups.ConstBegin();
 277     SamReadGroupConstIterator rgEnd  = readGroups.ConstEnd();
 278     for ( ; rgIter != rgEnd; ++rgIter ) {
 279         const SamReadGroup& rg = (*rgIter);
 280         isValid &= ValidateReadGroup(rg);
 281     }
 282
 283     // return validation state
 284     return isValid;
 285 }
 286
 287 bool SamHeaderValidator::ContainsUniqueIDsAndPlatformUnits(void) {
 288
 289     bool isValid = true;
 290     set<string> readGroupIds;
 291     set<string> platformUnits;
 292     set<string>::iterator idIter;
 293     set<string>::iterator puIter;
 294
 295     // iterate over sequences
 296     const SamReadGroupDictionary& readGroups = m_header.ReadGroups;
 297     SamReadGroupConstIterator rgIter = readGroups.ConstBegin();
 298     SamReadGroupConstIterator rgEnd  = readGroups.ConstEnd();
 299     for ( ; rgIter != rgEnd; ++rgIter ) {
 300         const SamReadGroup& rg = (*rgIter);
 301
 302         // --------------------------------
 303         // check for unique ID
 304
 305         // lookup read group ID
 306         const string& id = rg.ID;
 307         idIter = readGroupIds.find(id);
 308
 309         // error if found (duplicate entry)
 310         if ( idIter != readGroupIds.end() ) {
 311             AddError("Read group ID (ID): " + id + " is not unique");
 312             isValid = false;
 313         }
 314
 315         // otherwise ok, store id
 316         readGroupIds.insert(id);
 317
 318         // --------------------------------
 319         // check for unique platform unit
 320
 321         // lookup platform unit
 322         const string& pu = rg.PlatformUnit;
 323         puIter = platformUnits.find(pu);
 324
 325         // error if found (duplicate entry)
 326         if ( puIter != platformUnits.end() ) {
 327             AddError("Platform unit (PU): " + pu + " is not unique");
 328             isValid = false;
 329         }
 330
 331         // otherwise ok, store platform unit
 332         platformUnits.insert(pu);
 333     }
 334
 335     // return validation state
 336     return isValid;
 337 }
 338
 339 bool SamHeaderValidator::ValidateReadGroup(const SamReadGroup& rg) {
 340     bool isValid = true;
 341     isValid &= CheckReadGroupID(rg.ID);
 342     isValid &= CheckSequencingTechnology(rg.SequencingTechnology);
 343     return isValid;
 344 }
 345
 346 bool SamHeaderValidator::CheckReadGroupID(const string& id) {
 347
 348     // invalid if empty
 349     if ( id.empty() ) {
 350         AddError("Read group entry (@RG) is missing ID tag");
 351         return false;
 352     }
 353
 354     // otherwise OK
 355     return true;
 356 }
 357
 358 bool SamHeaderValidator::CheckSequencingTechnology(const string& technology) {
 359
 360     // if no technology provided, no problem, just return OK
 361     if ( technology.empty() ) return true;
 362
 363     // if technology is valid keyword
 364     if ( Is454(technology)      ||
 365          IsHelicos(technology)  ||
 366          IsIllumina(technology) ||
 367          IsPacBio(technology)   ||
 368          IsSolid(technology)
 369        )
 370     { return true; }
 371
 372     // otherwise
 373     AddError("Invalid read group sequencing platform (PL): " + technology);
 374     return false;
 375 }
 376
 377 bool SamHeaderValidator::Is454(const string& technology) {
 378     return ( technology == Constants::SAM_RG_SEQTECHNOLOGY_454 ||
 379              technology == Constants::SAM_RG_SEQTECHNOLOGY_LS454_LOWER ||
 380              technology == Constants::SAM_RG_SEQTECHNOLOGY_LS454_UPPER
 381            );
 382 }
 383
 384 bool SamHeaderValidator::IsHelicos(const string& technology) {
 385     return ( technology == Constants::SAM_RG_SEQTECHNOLOGY_HELICOS_LOWER ||
 386              technology == Constants::SAM_RG_SEQTECHNOLOGY_HELICOS_UPPER
 387            );
 388 }
 389
 390 bool SamHeaderValidator::IsIllumina(const string& technology) {
 391     return ( technology == Constants::SAM_RG_SEQTECHNOLOGY_ILLUMINA_LOWER ||
 392              technology == Constants::SAM_RG_SEQTECHNOLOGY_ILLUMINA_UPPER
 393            );
 394 }
 395
 396 bool SamHeaderValidator::IsPacBio(const string& technology) {
 397     return ( technology == Constants::SAM_RG_SEQTECHNOLOGY_PACBIO_LOWER ||
 398              technology == Constants::SAM_RG_SEQTECHNOLOGY_PACBIO_UPPER
 399            );
 400 }
 401
 402 bool SamHeaderValidator::IsSolid(const string& technology) {
 403     return ( technology == Constants::SAM_RG_SEQTECHNOLOGY_SOLID_LOWER ||
 404              technology == Constants::SAM_RG_SEQTECHNOLOGY_SOLID_UPPER
 405            );
 406 }
 407
 408 bool SamHeaderValidator::ValidateProgramData(void) {
 409     bool isValid = true;
 410     isValid &= ContainsUniqueProgramIds();
 411     isValid &= ValidatePreviousProgramIds();
 412     return isValid;
 413 }
 414
 415 bool SamHeaderValidator::ContainsUniqueProgramIds(void) {
 416     bool isValid = true;
 417     // TODO: once we have ability to handle multiple @PG entries,
 418     // check here for duplicate ID's
 419     // but for now, just return true
 420     return isValid;
 421 }
 422
 423 bool SamHeaderValidator::ValidatePreviousProgramIds(void) {
 424     bool isValid = true;
 425     // TODO: check that PP entries are valid later, after we get multiple @PG-entry handling
 426     // just return true for now
 427     return isValid;
 428 }
 429 void SamHeaderValidator::AddError(const string& message) {
 430     m_errorMessages.push_back(ERROR_PREFIX + message + NEWLINE);
 431 }
 432
 433 void SamHeaderValidator::AddWarning(const string& message) {
 434     m_warningMessages.push_back(WARN_PREFIX + message + NEWLINE);
 435 }
 436
 437 void SamHeaderValidator::PrintErrorMessages(void) {
 438
 439     // skip if no error messages
 440     if ( m_errorMessages.empty() ) return;
 441
 442     // print error header line
 443     cerr << "* SAM header has " << m_errorMessages.size() << " errors:" << endl;
 444
 445     // print each error message
 446     vector<string>::const_iterator errorIter = m_errorMessages.begin();
 447     vector<string>::const_iterator errorEnd  = m_errorMessages.end();
 448     for ( ; errorIter != errorEnd; ++errorIter )
 449         cerr << (*errorIter);
 450 }
 451
 452 void SamHeaderValidator::PrintWarningMessages(void) {
 453
 454     // skip if no warning messages
 455     if ( m_warningMessages.empty() ) return;
 456
 457     // print warning header line
 458     cerr << "* SAM header has " << m_warningMessages.size() << " warnings:" << endl;
 459
 460     // print each warning message
 461     vector<string>::const_iterator warnIter = m_warningMessages.begin();
 462     vector<string>::const_iterator warnEnd  = m_warningMessages.end();
 463     for ( ; warnIter != warnEnd; ++warnIter )
 464         cerr << (*warnIter);
 465 }