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