Loading...
Searching...
No Matches
GreedyTupleAligner_test.cpp
/*
* Copyright (C) 2010 Regents of the University of Michigan
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#include "GreedyTupleAligner.h"
#include <vector>
#include <iostream>
#include <sstream>
#include <gtest/gtest.h>
// design a mock class for GenomeSequence
class MockGenomeSequence{
public:
MockGenomeSequence(std::string sequence) :
sequence(sequence) {};
char operator[] (int index) const {
if (index < 0 || index >= (int)sequence.size()) {
std::cerr << "exceeds boundary! at " << __FILE__ << ":" <<__LINE__ << std::endl;
return 'N';
}
return sequence[index];
}
char& operator[] (int index) {
if (index < 0 || index >= (int)sequence.size()) {
std::cerr << "exceeds boundary! at " << __FILE__ << ":" <<__LINE__ << std::endl;
return sequence[0];
}
return sequence[index];
}
int getNumberBases() const{
return sequence.size();
}
private:
std::string sequence;
};
void printRefQueryCigar(const std::vector<char>& prettyPrintReference,
const std::vector<char>& prettyPrintQuery,
CigarRoller& cs,
std::ostream& out){
out << " ref = ";
for(std::vector<char>::const_iterator i=prettyPrintReference.begin(); i<prettyPrintReference.end(); i++) out << *i;
out << std::endl;
out <<"query = ";
for(std::vector<char>::const_iterator i=prettyPrintQuery.begin(); i<prettyPrintQuery.end(); i++) out << *i;
out << std::endl;
out << "cigar = ";
out << cs.getString();
out << " ";
}
void runAlign(const char *query, const char *ref, CigarRoller& cs, int& matchPosition) {
Weight wt;
MockGenomeSequence reference (ref);
GreedyTupleAligner<const char *, MockGenomeSequence, int> ga(wt);
ga.Align(query, strlen(query), reference, 0, strlen(ref), cs, matchPosition);
}
void computePrettyOutput (std::vector<char>& prettyPrintReference,
const char* ref,
std::vector<char>& prettyPrintQuery,
const char* query,
const int matchPosition,
const CigarRoller& expectedCigar)
{
const char* pRef = ref;
const char* pQuery = query;
for (int index = 0; index < matchPosition; index++){
prettyPrintReference.push_back(*pRef++);
prettyPrintQuery.push_back(' ');
}
for (int i = 0; i< expectedCigar.size(); i++) {
switch( expectedCigar[i].operation) {
case CigarRoller::mismatch:
case CigarRoller::match:
for (unsigned int j = 0; j < expectedCigar[i].count; j++){
prettyPrintReference.push_back(*pRef++);
prettyPrintQuery.push_back(*pQuery++);
}
break;
case CigarRoller::del:
for (unsigned int j = 0; j < expectedCigar[i].count; j++){
prettyPrintReference.push_back(*pRef++);
prettyPrintQuery.push_back(' ');
}
break;
case CigarRoller::insert:
for (unsigned int j = 0; j < expectedCigar[i].count; j++){
prettyPrintReference.push_back(' ');
prettyPrintQuery.push_back(*pQuery++);
}
break;
default:
break;
}
} /// end for
while (*pRef !='\0')
prettyPrintReference.push_back(*pRef++);
while (*pQuery !='\0')
prettyPrintReference.push_back(*pQuery++);
}
bool verifyAlign(const char *query, const char *ref, const char *expectedCigarString, std::ostream& out = std::cout) {
out.seekp(std::ios_base::beg);
CigarRoller cs;
int matchPosition;
CigarRoller expectedCigar(expectedCigarString);
runAlign(query, ref, cs, matchPosition);
if (matchPosition < 0) {
fprintf(stderr, "No match in %s, %d \n", __FILE__, __LINE__);
return false;
}
std::vector<char> prettyPrintReference,prettyPrintQuery;
computePrettyOutput( prettyPrintReference, ref,
prettyPrintQuery, query,
matchPosition,
expectedCigar);
const char *str = cs.getString();
if((unsigned int)expectedCigar.getExpectedQueryBaseCount() != strlen(query)) {
printRefQueryCigar(prettyPrintReference, prettyPrintQuery, cs, out);
out << std::endl;
out << "Expected Cigar string length " << expectedCigar.getExpectedQueryBaseCount() << " does not match the length of the query " << strlen(query) << ". Please fix test case." << std::endl;
return false;
}
if((unsigned int)cs.getExpectedQueryBaseCount() != strlen(query)) {
printRefQueryCigar(prettyPrintReference, prettyPrintQuery, cs, out);
out << std::endl;
out << "Query Length of " << strlen(query) << " does not match computed cigar string length of " << cs.getExpectedQueryBaseCount() << std::endl;
return false;
}
if (strcmp(expectedCigarString, str) == 0) {
// printf("[Correct Answer = %s] \n", expectedCigarString) ;
return true;
} else {
printRefQueryCigar(prettyPrintReference, prettyPrintQuery, cs, out);
out << "[Correct Answer = " << expectedCigarString << "] --------------------- Wrong!" << std::endl;
return false;
}
return true;
}
TEST(GreedyTupleAlignerTest, AlignToShortReference) {
std::stringstream ss(std::stringstream::out);
#if 0
//exact align
EXPECT_TRUE( verifyAlign("12345", "123456789", "5M") ) << ss.str().c_str();
EXPECT_TRUE( verifyAlign("23456", "123456789", "5M") ) << ss.str().c_str();
//mismatch
EXPECT_TRUE( verifyAlign("123B567", "123456789", "7M") ) << ss.str().c_str();
EXPECT_TRUE( verifyAlign("234D678", "123456789", "7M") ) << ss.str().c_str();
// del
EXPECT_TRUE( verifyAlign("123467890","1234567890", "4M1D5M") ) << ss.str().c_str();
EXPECT_TRUE( verifyAlign("123467890","B1234567890B", "4M1D5M") ) << ss.str().c_str();
// ins
EXPECT_TRUE( verifyAlign("12345067890","1234567890", "5M1I5M") ) << ss.str().c_str();
EXPECT_TRUE( verifyAlign("12345067890","BBBB1234567890BBBB", "5M1I5M") ) << ss.str().c_str();
// soft clip
EXPECT_TRUE( verifyAlign("1234", "1235", "3M1S") ) << ss.str().c_str();
// The following will treat as two mismatches
// EXPECT_TRUE( verifyAlign("123456700", "123456789", "7M2S", ss) ) << ss.str().c_str();
#endif
EXPECT_TRUE( verifyAlign("1023456700", "123456789","1I7M2S") ) << ss.str().c_str();
}
TEST(GreedyTupleTestAligner, AlignToLongReference) {
std::stringstream ss(std::stringstream::out);
EXPECT_TRUE( verifyAlign("TTAGAATGCTATTGTGTTTGGAGATTTGAGGAAAGTGGGCGTGAAGACTTAGTGTTCATTTCCTCAACCTCTCTCTGTGTGAACATACGTCATCGGTCAGAAATTGGG"," ", "88M200D20M") ) << ss.str().c_str();
// This reads cigar string is the wrong length. By that I mean that the number
// of matches listed in the cigar should be the same as the length of
// the original read.
//
EXPECT_TRUE( verifyAlign("GTGAAACTCCATCTCAAAAATAAGTAAATAAATAAATACATACATAGGCACAGTGCAGTTGTTAGTCAGAATTAGGTCACACTGGATTAGGGTGAGTACTTAATGCAACAGGTCTGGGG"," ", "63M3D56M") ) << ss.str().c_str();
EXPECT_TRUE( verifyAlign("ATATTGTTTTTTTCAATGCATATCAAAACAATGTTTACAATATACTACAGCCTAAGTGTGCAATAGCATTATGTGTAGAAATGCACATACCATAATTAGTTTTTTTTTTTGAAAAAACT"," ", "109M249D10M") ) << ss.str().c_str();
EXPECT_TRUE( verifyAlign("CCAGACTATCTCAAGCAATCAACAGATTTAATGTAAGGAGTGTCAAAATCTGAATGATGCTTTTTGCAGAAATAGAAAATCCCTTTCTAATATTTTTATATTTTTGAC"," ", "99M200D9M") ) << ss.str().c_str();
// this is by no mean works, since our local aligner generates the output
// 54M200D34M47D8M7D2M1I5M4S , since local aligner prefer gapped alignment.
EXPECT_FALSE( verifyAlign("ATGAGGTCAGGAGATGGAGACCATCCTGGCTAACATGGTGAAACCCCATCTCTAAAAAAAGTGTAACAGAGGTGCATACTCAAAACTACAAAAGTCTCGTGAAAGGAA"," ", "54M200D54M") ) << ss.str().c_str();
}
#if 0
//
// In this test case, there is no indel that I am aware of, however, the
// aligner produces quite a long cigar string with many insertions
// and deletions. For now, I filter this case out later, but it would
// be nice if it would limit itself to one or a small number of indels.
//
EXPECT_TRUE( verifyAlign("GAATCAATACGCTCGGGATGCAGCGCCTAGCCGTTGGTTTGAGAATGGTTCTCTAGAGTTATCTTCACCCTCTACCTTGTGTGGCACTATTTCTTCTATGACCTTGAC"," ", "108M") ) << ss.str().c_str();
#endif
TEST(GreedyTupleTestAligner, EquivalentAlignment) {
std::stringstream ss(std::stringstream::out);
// The test cases here are different to what we expected,
// but hte alignment are both valid. See the example:
// ref: AAAT TGGGG (4M5D5M)
// ref: AAA TTGGGG (3M5D6M)
// read: AAAT CCCC TTGGGG
// We obtain 89M200D19M, the correct is 87M200D21M,
// the 2 matched bases shifted to the right
// if you check the output, you will see the 2 bases can be put before '200D' or after '200D'
// and they are essentially the same
EXPECT_FALSE( verifyAlign("TTTTCTTTTCAAAAATTTAAAAGTGACATACAAAATTATATGTGTATGTACAACAAAAGCTTAACTATAACACCTTGTTACATACTTTGGAATTGAAAGGCAGGAATG"," ", "87M200D21M") ) << ss.str().c_str();
// for the same reason as before
EXPECT_FALSE( verifyAlign("AACAGTTGAGAGGTACTAAAATTGAGTTTTCTTGAAAAATATATTTAATCTAAAGTACTGAAAATTTGGGGGAAAATGCTTAAGGTCATATTCCTTTTTTGAAAAGAT"," ", "80M200D28M") ) << ss.str().c_str();
}
CigarRoller
The purpose of this class is to provide accessors for setting, updating, modifying the CIGAR object....
Definition CigarRoller.h:67
CigarRoller::getString
const char * getString()
Get the string reprentation of the Cigar operations in this object, caller must delete the returned v...
Definition CigarRoller.cpp:272
Cigar::size
int size() const
Return the number of cigar operations.
Definition Cigar.h:364
Cigar::getExpectedQueryBaseCount
int getExpectedQueryBaseCount() const
Return the length of the read that corresponds to the current CIGAR string.
Definition Cigar.cpp:95
Cigar::del
@ del
deletion from the reference (the reference contains bases that have no corresponding base in the quer...
Definition Cigar.h:92
Cigar::mismatch
@ mismatch
mismatch operation. Associated with CIGAR Operation "M"
Definition Cigar.h:90
Cigar::match
@ match
match/mismatch operation. Associated with CIGAR Operation "M"
Definition Cigar.h:89
Cigar::insert
@ insert
insertion to the reference (the query sequence contains bases that have no corresponding base in the ...
Definition Cigar.h:91
GreedyTupleAligner
Definition GreedyTupleAligner.h:61
MockGenomeSequence
Definition GreedyTupleAligner_test.cpp:24
Weight
Weight includes various penalties(e.g.
Definition GreedyTupleAligner.h:41