src/api/internal/SamHeaderValidator_p.cpp

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