ReferenceSequence.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 "assert.h"
#include "ctype.h"
#include "stdio.h"
#include "Error.h"


#include "Generic.h"
#include "ReferenceSequence.h"

#include <algorithm>
#include <istream>
#include <fstream>
#include <sstream>
#include <stdexcept>

//
// Given a buffer with a fasta format contents, count
// the number of chromsomes in it and return that value.
//
bool getFastaStats(const char *fastaData, size_t fastaDataSize, uint32_t &chromosomeCount, uint64_t &baseCount)
{
    chromosomeCount = 0;
    baseCount = 0;
    bool atLineStart = true;

    //
    // loop over the fasta file, essentially matching for the
    // pattern '^>.*$' and counting them.
    //
    for (size_t fastaIndex = 0; fastaIndex < fastaDataSize; fastaIndex++)
    {
        switch (fastaData[fastaIndex])
        {
            case '\n':
            case '\r':
                atLineStart = true;
                break;
            case '>':
            {
                if (!atLineStart) break;
                chromosomeCount++;
                //
                // eat the rest of the line
                //
                while (fastaIndex < fastaDataSize && fastaData[fastaIndex]!='\n' && fastaData[fastaIndex]!='\r')
                {
                    fastaIndex++;
                }
                break;
            }
            default:
                baseCount++;
                atLineStart = false;
                break;
        }

    }
    return false;
}

#if 0
// turn this into a template on read/quality/etc...
int GenomeSequence::debugPrintReadValidation(
    std::string &read,
    std::string &quality,
    char   direction,
    genomeIndex_t   readLocation,
    int sumQuality,
    int mismatchCount,
    bool recurse
) 
{
    int validateSumQ = 0;
    int validateMismatchCount = 0;
    int rc = 0;
    std::string genomeData;

    for (uint32_t i=0; i<read.size(); i++)
    {
        if (tolower(read[i]) != tolower((*this)[readLocation + i]))
        {
            validateSumQ += quality[i] - '!';
            // XXX no longer valid:
            if (direction=='F' ? i<24 : (i >= (read.size() - 24))) validateMismatchCount++;
            genomeData.push_back(tolower((*this)[readLocation + i]));
        }
        else
        {
            genomeData.push_back(toupper((*this)[readLocation + i]));
        }
    }
    assert(validateSumQ>=0);
    if (validateSumQ != sumQuality && validateMismatchCount == mismatchCount)
    {
        printf("SUMQ: Original Genome: %s  test read: %s : actual sumQ = %d, test sumQ = %d\n",
               genomeData.c_str(),
               read.c_str(),
               validateSumQ,
               sumQuality
              );
        rc++;
    }
    else if (validateSumQ == sumQuality && validateMismatchCount != mismatchCount)
    {
        printf("MISM: Original Genome: %s  test read: %s : actual mismatch %d test mismatches %d\n",
               genomeData.c_str(),
               read.c_str(),
               validateMismatchCount,
               mismatchCount
              );
        rc++;
    }
    else if (validateSumQ != sumQuality && validateMismatchCount != mismatchCount)
    {
        printf("BOTH: Original Genome: %s  test read: %s : actual sumQ = %d, test sumQ = %d, actual mismatch %d test mismatches %d\n",
               genomeData.c_str(),
               read.c_str(),
               validateSumQ,
               sumQuality,
               validateMismatchCount,
               mismatchCount
              );
        rc++;
    }

    if (recurse && abs(validateMismatchCount - mismatchCount) > (int) read.size()/2)
    {
        printf("large mismatch difference, trying reverse strand: ");
        std::string reverseRead = read;
        std::string reverseQuality = quality;
        getReverseRead(reverseRead);
        reverse(reverseQuality.begin(), reverseQuality.end());
        rc = debugPrintReadValidation(reverseRead, reverseQuality, readLocation, sumQuality, mismatchCount, false);
    }
    return rc;
}
#endif