TestEquals.cpp

00001 /*
00002  *  Copyright (C) 2010  Regents of the University of Michigan
00003  *
00004  *   This program is free software: you can redistribute it and/or modify
00005  *   it under the terms of the GNU General Public License as published by
00006  *   the Free Software Foundation, either version 3 of the License, or
00007  *   (at your option) any later version.
00008  *
00009  *   This program is distributed in the hope that it will be useful,
00010  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *   GNU General Public License for more details.
00013  *
00014  *   You should have received a copy of the GNU General Public License
00015  *   along with this program.  If not, see <http://www.gnu.org/licenses/>.
00016  */
00017 
00018 #include "TestEquals.h"
00019 #include <assert.h>
00020 
00021 
00022 
00023 void testSeqEquals()
00024 {
00025     // Call generic test which since the sam and bam are identical, should
00026     // contain the same results.
00027     EqualsTest::testEq(EqualsTest::SAM);
00028     EqualsTest::testEq(EqualsTest::BAM);
00029 }
00030 
00031 
00032 const char* EqualsTest::READ_NAMES[] = 
00033     {"01:====", "02:===X", "03:==X=", "04:==XX", "05:=X==", "06:=X=X",
00034      "07:=XX=", "08:=XXX", "09:X===", "10:X==X", "11:X=X=", "12:X=XX",
00035      "13:XX==", "14:XX=X", "15:XXX=", "16:XXXX", "Read:GGCCTA;Ref:CCTA",
00036      "Read:CCTA;Ref:CCTA", "Read:CCGTxxxC;Ref:CCxTAACC", 
00037      "Read:CCxxAC;Ref:CCTAACC"};
00038 
00039 const char* EqualsTest::READ_SEQS_BASES[] = 
00040     {"CCTA", "CCTT", "CCAA", "CCAT", "CTTA", "CTTT", "CTAA", "CTAT", "TCTA", "TCTT", "TCAA", "TCAT", "TTTA", "TTTT", "TTAA", "TTAT", "GGCCTA",
00041      "CCTA", "CCGTC", "CCAC"};
00042 
00043 const char* EqualsTest::READ_SEQS_EQUALS[] = 
00044     {"====", "===T", "==A=", "==AT", "=T==", "=T=T", "=TA=", "=TAT", "T===", "T==T", "T=A=", "T=AT", "TT==", "TT=T", "TTA=", "TTAT", "GG====",
00045      "====", "==G==", "===="};
00046 
00047 const char* EqualsTest::READ_SEQS_MIXED[] = 
00048     {"C===", "=C=T", "==AA", "==AT", "=TTA", "CT=T", "=TAA", "=TAT", "T=TA", "TC=T", "TCA=", "TCAT", "TT=A", "TT=T", "TTA=", "TTAT", "GGC=T=",
00049      "C=T=", "C=GT=", "C=A="};
00050 
00051 const char* EqualsTest::expectedReferenceName;
00052 const char* EqualsTest::expectedMateReferenceName;
00053 const char* EqualsTest::expectedMateReferenceNameOrEqual;
00054 const char* EqualsTest::expectedCigar;
00055 const char* EqualsTest::expectedQuality;
00056 
00057 std::vector<unsigned int> EqualsTest::expectedCigarHex;
00058 
00059 int EqualsTest::expected0BasedAlignmentEnd;
00060 int EqualsTest::expected1BasedAlignmentEnd;
00061 int EqualsTest::expectedAlignmentLength;
00062 int EqualsTest::expected0BasedUnclippedStart;
00063 int EqualsTest::expected1BasedUnclippedStart;
00064 int EqualsTest::expected0BasedUnclippedEnd;
00065 int EqualsTest::expected1BasedUnclippedEnd;
00066 bamRecordStruct EqualsTest::expectedRecord;
00067 
00068 void EqualsTest::testEq(FileType inputType)
00069 {
00070     reset();
00071     SamFile inSam;
00072 
00073     std::string outputBase = "results/out";
00074 
00075     if(inputType == SAM)
00076     {
00077         assert(inSam.OpenForRead("testFiles/testEq.sam"));
00078         outputBase += "SamEq";
00079     }
00080     else
00081     {
00082         assert(inSam.OpenForRead("testFiles/testEq.bam"));
00083         outputBase += "BamEq";
00084     }
00085 
00086    // Read the SAM Header.
00087     SamFileHeader samHeader;
00088     assert(inSam.ReadHeader(samHeader));
00089 
00090     std::string outputName = outputBase + "Bases.sam";
00091     SamFile outBasesSam( outputName.c_str(), SamFile::WRITE);
00092     outputName = outputBase + "Equals.sam";
00093     SamFile outEqualsSam(outputName.c_str(), SamFile::WRITE);
00094     outputName = outputBase + "Orig.sam";
00095     SamFile outOrigSam(  outputName.c_str(), SamFile::WRITE);
00096     outputName = outputBase + "Bases.bam";
00097     SamFile outBasesBam( outputName.c_str(), SamFile::WRITE);
00098     outputName = outputBase + "Equals.bam";
00099     SamFile outEqualsBam(outputName.c_str(), SamFile::WRITE);
00100     outputName = outputBase + "Orig.bam";
00101     SamFile outOrigBam(  outputName.c_str(), SamFile::WRITE);
00102     assert(outBasesSam.WriteHeader(samHeader));
00103     assert(outEqualsSam.WriteHeader(samHeader));
00104     assert(outOrigSam.WriteHeader(samHeader));
00105     assert(outBasesBam.WriteHeader(samHeader));
00106     assert(outEqualsBam.WriteHeader(samHeader));
00107     assert(outOrigBam.WriteHeader(samHeader));
00108 
00109     outBasesSam.SetWriteSequenceTranslation(SamRecord::BASES);
00110     outEqualsSam.SetWriteSequenceTranslation(SamRecord::EQUAL);
00111     outOrigSam.SetWriteSequenceTranslation(SamRecord::NONE);
00112     outBasesBam.SetWriteSequenceTranslation(SamRecord::BASES);
00113     outEqualsBam.SetWriteSequenceTranslation(SamRecord::EQUAL);
00114     outOrigBam.SetWriteSequenceTranslation(SamRecord::NONE);
00115 
00116     GenomeSequence reference("testFiles/chr1_partial.fa");
00117 
00118     inSam.SetReference(&reference);
00119 
00120     SamRecord samRecord;
00121 
00122     // The set of 16 variations are repeated 3 times: once with all charcters
00123     // 1) Matches have the actual bases in them.
00124     // 2) Matches have '='
00125     // 3) Matches are mixed between bases and '='
00126     // Since Sequences are 4 characters long, there are 16 variations
00127     // of match/mismatch.
00128     for(int j = 0; j < 16; j++)
00129     {
00130         assert(inSam.ReadRecord(samHeader, samRecord) == true);
00131         validateEqRead(samRecord, j, READ_SEQS_BASES[j]);
00132         assert(outBasesSam.WriteRecord(samHeader, samRecord));
00133         assert(outEqualsSam.WriteRecord(samHeader, samRecord));
00134         assert(outOrigSam.WriteRecord(samHeader, samRecord));
00135         assert(outBasesBam.WriteRecord(samHeader, samRecord));
00136         assert(outEqualsBam.WriteRecord(samHeader, samRecord));
00137         assert(outOrigBam.WriteRecord(samHeader, samRecord));
00138     }
00139     for(int j = 0; j < 16; j++)
00140     {
00141         assert(inSam.ReadRecord(samHeader, samRecord) == true);
00142         validateEqRead(samRecord, j, READ_SEQS_EQUALS[j]);
00143         assert(outBasesSam.WriteRecord(samHeader, samRecord));
00144         assert(outEqualsSam.WriteRecord(samHeader, samRecord));
00145         assert(outOrigSam.WriteRecord(samHeader, samRecord));
00146         assert(outBasesBam.WriteRecord(samHeader, samRecord));
00147         assert(outEqualsBam.WriteRecord(samHeader, samRecord));
00148         assert(outOrigBam.WriteRecord(samHeader, samRecord));
00149     }
00150     for(int j = 0; j < 16; j++)
00151     {
00152         assert(inSam.ReadRecord(samHeader, samRecord) == true);
00153         validateEqRead(samRecord, j, READ_SEQS_MIXED[j]);
00154         assert(outBasesSam.WriteRecord(samHeader, samRecord));
00155         assert(outEqualsSam.WriteRecord(samHeader, samRecord));
00156         assert(outOrigSam.WriteRecord(samHeader, samRecord));
00157         assert(outBasesBam.WriteRecord(samHeader, samRecord));
00158         assert(outEqualsBam.WriteRecord(samHeader, samRecord));
00159         assert(outOrigBam.WriteRecord(samHeader, samRecord));
00160     }
00161 
00162     expectedCigar = "2S4M";
00163     expectedCigarHex.clear();
00164     expectedCigarHex.push_back(0x24);
00165     expectedCigarHex.push_back(0x40);
00166     expected0BasedUnclippedStart = expectedRecord.myPosition-2;
00167     expected1BasedUnclippedStart = expected0BasedUnclippedStart + 1;
00168     expectedRecord.myBlockSize = 70;
00169     expectedRecord.myReadNameLength = 21;
00170     expectedRecord.myCigarLength = 2;
00171     expectedRecord.myReadLength = 6;
00172     expectedQuality = "??I00?";
00173     assert(inSam.ReadRecord(samHeader, samRecord) == true);
00174     validateEqRead(samRecord, 16, READ_SEQS_MIXED[16]);
00175         assert(outBasesSam.WriteRecord(samHeader, samRecord));
00176         assert(outEqualsSam.WriteRecord(samHeader, samRecord));
00177         assert(outOrigSam.WriteRecord(samHeader, samRecord));
00178         assert(outBasesBam.WriteRecord(samHeader, samRecord));
00179         assert(outEqualsBam.WriteRecord(samHeader, samRecord));
00180         assert(outOrigBam.WriteRecord(samHeader, samRecord));
00181 
00182     expectedCigar = "4M4H";
00183     expectedCigarHex.clear();
00184     expectedCigarHex.push_back(0x40);
00185     expectedCigarHex.push_back(0x45);
00186     expected0BasedUnclippedStart = expectedRecord.myPosition;
00187     expected1BasedUnclippedStart = expected0BasedUnclippedStart + 1;
00188     expected0BasedUnclippedEnd = expectedRecord.myPosition + 7;
00189     expected1BasedUnclippedEnd = expected0BasedUnclippedEnd + 1;
00190     expectedRecord.myBlockSize = 65;
00191     expectedRecord.myReadNameLength = 19;
00192     expectedRecord.myCigarLength = 2;
00193     expectedRecord.myReadLength = 4;
00194     expectedQuality = "I00?";
00195     assert(inSam.ReadRecord(samHeader, samRecord) == true);
00196     validateEqRead(samRecord, 17, READ_SEQS_MIXED[17]);
00197         assert(outBasesSam.WriteRecord(samHeader, samRecord));
00198         assert(outEqualsSam.WriteRecord(samHeader, samRecord));
00199         assert(outOrigSam.WriteRecord(samHeader, samRecord));
00200         assert(outBasesBam.WriteRecord(samHeader, samRecord));
00201         assert(outEqualsBam.WriteRecord(samHeader, samRecord));
00202         assert(outOrigBam.WriteRecord(samHeader, samRecord));
00203 
00204     expectedCigar = "1M1P1M1I1M3D1M";
00205     expectedCigarHex.clear();
00206     expectedCigarHex.push_back(0x10);
00207     expectedCigarHex.push_back(0x16);
00208     expectedCigarHex.push_back(0x10);
00209     expectedCigarHex.push_back(0x11);
00210     expectedCigarHex.push_back(0x10);
00211     expectedCigarHex.push_back(0x32);
00212     expectedCigarHex.push_back(0x10);
00213     expected0BasedAlignmentEnd = expectedRecord.myPosition + 6;
00214     expected1BasedAlignmentEnd = expected0BasedAlignmentEnd + 1;
00215     expectedAlignmentLength = 7;
00216     expected0BasedUnclippedStart = expectedRecord.myPosition;
00217     expected1BasedUnclippedStart = expected0BasedUnclippedStart + 1;
00218     expected0BasedUnclippedEnd = expected0BasedAlignmentEnd;
00219     expected1BasedUnclippedEnd = expected0BasedUnclippedEnd + 1;
00220     expectedRecord.myBlockSize = 95;
00221     expectedRecord.myReadNameLength = 27;
00222     expectedRecord.myCigarLength = 7;
00223     expectedRecord.myReadLength = 5;
00224     expectedQuality = "I00??";
00225     assert(inSam.ReadRecord(samHeader, samRecord) == true);
00226     validateEqRead(samRecord, 18, READ_SEQS_MIXED[18]);
00227         assert(outBasesSam.WriteRecord(samHeader, samRecord));
00228         assert(outEqualsSam.WriteRecord(samHeader, samRecord));
00229         assert(outOrigSam.WriteRecord(samHeader, samRecord));
00230         assert(outBasesBam.WriteRecord(samHeader, samRecord));
00231         assert(outEqualsBam.WriteRecord(samHeader, samRecord));
00232         assert(outOrigBam.WriteRecord(samHeader, samRecord));
00233 
00234     expectedCigar = "2M2N2M";
00235     expectedCigarHex.clear();
00236     expectedCigarHex.push_back(0x20);
00237     expectedCigarHex.push_back(0x23);
00238     expectedCigarHex.push_back(0x20);
00239     expected0BasedAlignmentEnd = expectedRecord.myPosition + 5;
00240     expected1BasedAlignmentEnd = expected0BasedAlignmentEnd + 1;
00241     expectedAlignmentLength = 6;
00242     expected0BasedUnclippedStart = expectedRecord.myPosition;
00243     expected1BasedUnclippedStart = expected0BasedUnclippedStart + 1;
00244     expected0BasedUnclippedEnd = expected0BasedAlignmentEnd;
00245     expected1BasedUnclippedEnd = expected0BasedUnclippedEnd + 1;
00246     expectedRecord.myBlockSize = 74;
00247     expectedRecord.myReadNameLength = 24;
00248     expectedRecord.myCigarLength = 3;
00249     expectedRecord.myReadLength = 4;
00250     expectedQuality = "I00?";
00251     assert(inSam.ReadRecord(samHeader, samRecord) == true);
00252     validateEqRead(samRecord, 19, READ_SEQS_MIXED[19]);
00253         assert(outBasesSam.WriteRecord(samHeader, samRecord));
00254         assert(outEqualsSam.WriteRecord(samHeader, samRecord));
00255         assert(outOrigSam.WriteRecord(samHeader, samRecord));
00256         assert(outBasesBam.WriteRecord(samHeader, samRecord));
00257         assert(outEqualsBam.WriteRecord(samHeader, samRecord));
00258         assert(outOrigBam.WriteRecord(samHeader, samRecord));
00259 
00260 }
00261 
00262 
00263 void EqualsTest::reset()
00264 {
00265     expectedReferenceName = "1";
00266     expectedMateReferenceName = "1";
00267     expectedMateReferenceNameOrEqual = "=";
00268     expectedCigar = "4M";
00269     expectedQuality = "I00?";
00270 
00271     // The First cigar is 4M which is 4 << 4 | 0 = 0x40 = 64
00272     expectedCigarHex.clear();
00273     expectedCigarHex.push_back(0x40);
00274 
00275     expectedRecord.myBlockSize = 50;
00276     expectedRecord.myReferenceID = 0;
00277     expectedRecord.myPosition = 10010;
00278     expectedRecord.myReadNameLength = 8;
00279     expectedRecord.myMapQuality = 0;
00280     expectedRecord.myBin = 4681;
00281     expectedRecord.myCigarLength = 1;
00282     expectedRecord.myFlag = 73;
00283     expectedRecord.myReadLength = 4;
00284     expectedRecord.myMateReferenceID = 0;
00285     expectedRecord.myMatePosition = 10008;
00286     expectedRecord.myInsertSize = 0;
00287 
00288     expected0BasedAlignmentEnd = 10013;
00289     expected1BasedAlignmentEnd = expected0BasedAlignmentEnd + 1;
00290     expectedAlignmentLength = 4;
00291     expected0BasedUnclippedStart = expectedRecord.myPosition;
00292     expected1BasedUnclippedStart = expected0BasedUnclippedStart + 1;
00293     expected0BasedUnclippedEnd = expected0BasedAlignmentEnd;
00294     expected1BasedUnclippedEnd = expected1BasedAlignmentEnd;
00295 }
00296 
00297 void EqualsTest::validateEqRead(SamRecord& samRecord, 
00298                                 int readIndex,
00299                                 const char* actualExpectedSequence)
00300 {
00301     char tag[3];
00302     char type;
00303     void* value;
00304 
00305     //////////////////////////////////////////
00306     // Validate Record 1
00307     // Check the alignment end
00308     assert(samRecord.get0BasedAlignmentEnd() == expected0BasedAlignmentEnd);
00309     assert(samRecord.get1BasedAlignmentEnd() == expected1BasedAlignmentEnd);
00310     assert(samRecord.getAlignmentLength() == expectedAlignmentLength);
00311     assert(samRecord.get0BasedUnclippedStart() == expected0BasedUnclippedStart);
00312     assert(samRecord.get1BasedUnclippedStart() == expected1BasedUnclippedStart);
00313     assert(samRecord.get0BasedUnclippedEnd() == expected0BasedUnclippedEnd);
00314     assert(samRecord.get1BasedUnclippedEnd() == expected1BasedUnclippedEnd);
00315 
00316     // Check the accessors.
00317     assert(samRecord.getBlockSize() == expectedRecord.myBlockSize);
00318     assert(samRecord.getReferenceID() == expectedRecord.myReferenceID);
00319     assert(strcmp(samRecord.getReferenceName(), expectedReferenceName) == 0);
00320     assert(samRecord.get1BasedPosition() == expectedRecord.myPosition + 1);
00321     assert(samRecord.get0BasedPosition() == expectedRecord.myPosition);
00322     assert(samRecord.getReadNameLength() == 
00323            expectedRecord.myReadNameLength);
00324     assert(samRecord.getMapQuality() == expectedRecord.myMapQuality);
00325     assert(samRecord.getBin() == expectedRecord.myBin);
00326     assert(samRecord.getCigarLength() == expectedRecord.myCigarLength);
00327     assert(samRecord.getFlag() == expectedRecord.myFlag);
00328     assert(samRecord.getReadLength() == expectedRecord.myReadLength);
00329     assert(samRecord.getMateReferenceID() == 
00330            expectedRecord.myMateReferenceID);
00331     assert(strcmp(samRecord.getMateReferenceName(),
00332                   expectedMateReferenceName) == 0);
00333     assert(strcmp(samRecord.getMateReferenceNameOrEqual(),
00334                   expectedMateReferenceNameOrEqual) == 0);
00335     assert(samRecord.get1BasedMatePosition() ==
00336            expectedRecord.myMatePosition + 1);
00337     assert(samRecord.get0BasedMatePosition() == 
00338            expectedRecord.myMatePosition);
00339     assert(samRecord.getInsertSize() == expectedRecord.myInsertSize);
00340     assert(strcmp(samRecord.getReadName(), READ_NAMES[readIndex]) == 0);
00341     assert(strcmp(samRecord.getCigar(), expectedCigar) == 0);
00342     samRecord.setSequenceTranslation(SamRecord::BASES);
00343     assert(strcmp(samRecord.getSequence(), READ_SEQS_BASES[readIndex]) == 0);
00344     assert(strcmp(samRecord.getQuality(), expectedQuality) == 0);
00345 
00346     assert(samRecord.getSequence(0) == READ_SEQS_BASES[readIndex][0]);
00347     assert(samRecord.getQuality(0) == expectedQuality[0]);
00348     assert(samRecord.getSequence(1)== READ_SEQS_BASES[readIndex][1]);
00349     assert(samRecord.getQuality(1) == expectedQuality[1]);
00350     assert(samRecord.getSequence(2) == READ_SEQS_BASES[readIndex][2]);
00351     assert(samRecord.getQuality(2) == expectedQuality[2]);
00352     assert(samRecord.getSequence(3) == READ_SEQS_BASES[readIndex][3]);
00353     assert(samRecord.getQuality(3) == expectedQuality[3]);
00354 
00355     assert(strcmp(samRecord.getSequence(SamRecord::EQUAL), 
00356                   READ_SEQS_EQUALS[readIndex]) == 0);
00357     assert(samRecord.getSequence(0, SamRecord::EQUAL) == 
00358            READ_SEQS_EQUALS[readIndex][0]);
00359     assert(samRecord.getQuality(0) == expectedQuality[0]);
00360     assert(samRecord.getSequence(1, SamRecord::EQUAL) == 
00361            READ_SEQS_EQUALS[readIndex][1]);
00362     assert(samRecord.getQuality(1) == expectedQuality[1]);
00363     assert(samRecord.getSequence(2, SamRecord::EQUAL) == 
00364            READ_SEQS_EQUALS[readIndex][2]);
00365     assert(samRecord.getQuality(2) == expectedQuality[2]);
00366     assert(samRecord.getSequence(3, SamRecord::EQUAL) == 
00367            READ_SEQS_EQUALS[readIndex][3]);
00368     assert(samRecord.getQuality(3) == expectedQuality[3]);
00369 
00370     assert(strcmp(samRecord.getSequence(SamRecord::NONE), 
00371                   actualExpectedSequence) == 0);
00372     assert(samRecord.getSequence(0, SamRecord::NONE) == 
00373            actualExpectedSequence[0]);
00374     assert(samRecord.getQuality(0) == expectedQuality[0]);
00375     assert(samRecord.getSequence(1, SamRecord::NONE) == 
00376            actualExpectedSequence[1]);
00377     assert(samRecord.getQuality(1) == expectedQuality[1]);
00378     assert(samRecord.getSequence(2, SamRecord::NONE) == 
00379            actualExpectedSequence[2]);
00380     assert(samRecord.getQuality(2) == expectedQuality[2]);
00381     assert(samRecord.getSequence(3, SamRecord::NONE) == 
00382            actualExpectedSequence[3]);
00383     assert(samRecord.getQuality(3) == expectedQuality[3]);
00384 
00385     samRecord.setSequenceTranslation(SamRecord::NONE);
00386     assert(strcmp(samRecord.getSequence(), 
00387                   actualExpectedSequence) == 0);
00388     assert(samRecord.getSequence(0) == 
00389            actualExpectedSequence[0]);
00390     assert(samRecord.getQuality(0) == expectedQuality[0]);
00391     assert(samRecord.getSequence(1) == 
00392            actualExpectedSequence[1]);
00393     assert(samRecord.getQuality(1) == expectedQuality[1]);
00394     assert(samRecord.getSequence(2) == 
00395            actualExpectedSequence[2]);
00396     assert(samRecord.getQuality(2) == expectedQuality[2]);
00397     assert(samRecord.getSequence(3) == 
00398            actualExpectedSequence[3]);
00399     assert(samRecord.getQuality(3) == expectedQuality[3]);
00400 
00401     // No tags, should return false.
00402     
00403     assert(samRecord.getNextSamTag(tag, type, &value) == false);
00404 
00405     // Get the record ptr.   
00406     samRecord.setSequenceTranslation(SamRecord::BASES);
00407     validateEqReadBuffer(samRecord, READ_SEQS_BASES[readIndex]);
00408     samRecord.setSequenceTranslation(SamRecord::NONE);
00409     validateEqReadBuffer(samRecord, actualExpectedSequence);
00410     samRecord.setSequenceTranslation(SamRecord::EQUAL);
00411     validateEqReadBuffer(samRecord, READ_SEQS_EQUALS[readIndex]);
00412 }
00413 
00414 
00415 void EqualsTest::validateEqReadBuffer(SamRecord& samRecord, 
00416                                       const char* expectedSequence)
00417 {
00418     const bamRecordStruct* bufferPtr;
00419     unsigned char* varPtr;
00420 
00421     bufferPtr = (const bamRecordStruct*)samRecord.getRecordBuffer();
00422     // Validate the buffers match.
00423     assert(bufferPtr->myBlockSize == expectedRecord.myBlockSize);
00424     assert(bufferPtr->myReferenceID == expectedRecord.myReferenceID);
00425     assert(bufferPtr->myPosition == expectedRecord.myPosition);
00426     assert(bufferPtr->myReadNameLength == expectedRecord.myReadNameLength);
00427     assert(bufferPtr->myMapQuality == expectedRecord.myMapQuality);
00428     assert(bufferPtr->myBin == expectedRecord.myBin);
00429     assert(bufferPtr->myCigarLength == expectedRecord.myCigarLength);
00430     assert(bufferPtr->myFlag == expectedRecord.myFlag);
00431     assert(bufferPtr->myReadLength == expectedRecord.myReadLength);
00432     assert(bufferPtr->myMateReferenceID == 
00433            expectedRecord.myMateReferenceID);
00434     assert(bufferPtr->myMatePosition == expectedRecord.myMatePosition);
00435     assert(bufferPtr->myInsertSize == expectedRecord.myInsertSize);
00436 
00437     // Validate the variable length fields in the buffer.
00438     // Set the pointer to the start of the variable fields.
00439     varPtr = (unsigned char*)(&(bufferPtr->myData[0]));
00440 
00441     // Validate the readname.
00442     for(int i = 0; i < expectedRecord.myReadNameLength; i++)
00443     {
00444         assert(*varPtr == samRecord.getReadName()[i]);
00445         varPtr++;
00446     }
00447 
00448     // Validate the cigar.
00449     for(int i = 0; i < expectedRecord.myCigarLength; i++)
00450     {
00451         assert(*(unsigned int*)varPtr == expectedCigarHex[i]);
00452         // Increment the varptr the size of an int.
00453         varPtr += 4;
00454     }
00455 
00456     // Validate the sequence.
00457     int expectedSeqHex = 0;
00458     for(int i = 0; i < expectedRecord.myReadLength; i++)
00459     {
00460         int hexChar;
00461         switch(expectedSequence[i])
00462         {
00463             case '=':
00464                 hexChar = 0x0;
00465                 break;
00466             case 'A':
00467             case 'a':
00468                 hexChar = 0x1;
00469                 break;
00470             case 'C':
00471             case 'c':
00472                 hexChar = 0x2;
00473                 break;
00474             case 'G':
00475             case 'g':
00476                 hexChar = 0x4;
00477                 break;
00478             case 'T':
00479             case 't':
00480                 hexChar = 0x8;
00481                 break;
00482             case 'N':
00483             case 'n':
00484                 hexChar = 0xF;
00485                 break;
00486         }
00487         if(i%2 == 0)
00488         {
00489             expectedSeqHex = hexChar << 4;
00490         }
00491         else
00492         {
00493             expectedSeqHex |= hexChar;
00494             assert(*varPtr == expectedSeqHex);
00495             varPtr++;         
00496         }
00497     }
00498     if((expectedRecord.myReadLength%2) != 0)
00499     {
00500        // Odd number of sequences, so need to check the last one.
00501         assert(*varPtr == expectedSeqHex);
00502         varPtr++;
00503     }
00504 
00505     // Validate the Quality
00506     for(int i = 0; i < expectedRecord.myReadLength; i++)
00507     {
00508         assert(*varPtr == samRecord.getQuality()[i] - 33);
00509         varPtr++;
00510     }
00511 
00512 }
00513 
00514 
Generated on Thu Dec 9 12:22:13 2010 for StatGen Software by  doxygen 1.6.3