src/conversion/HapConverter.cpp

// $Id: HapConverter.cpp,v 1.24 2011/03/07 06:08:48 bobgian Exp $

/*
  Copyright 2002 Patrick Colacurcio, Peter Beerli, Mary Kuhner,
  Jon Yamato and Joseph Felsenstein

  This software is distributed free of charge for non-commercial use
  and is copyrighted.  Of course, we do not guarantee that the software
  works, and are not responsible for any damage you may cause or have.
*/

#include <cassert>

#include "Converter_HapConverter.h"
#include "Converter_DataSourceException.h"
#include "Converter_ParserUtil.h"
#include "random.h"
#include "constants.h"

//------------------------------------------------------------------------------------

HapConverter::HapConverter(RegionDS& region, Random& random) :
    ConverterIf(), m_region(region), m_random(random),
    m_hapnames(region.GetAllSeqNames())
{
    // m_hapnames should end up with the same order of sequences as
    // they existed in the user file, don't know if this works!

} // HapConverter::ctor

//------------------------------------------------------------------------------------

HapConverter::~HapConverter()
{
} // HapConverter::dtor

//------------------------------------------------------------------------------------

vector<long> HapConverter::ParsePhaseInfo(ifstream& input) const
// EWFIX.P3.BUG.351 -- this doesn't handle the case where you have no
// phase known/unknown info, but you're just trying to assign
// samples to individuals
{
    skipWhiteSpace(input);
    string theline;
    getLine(input, theline);

    istringstream linestream(theline);
    string somestring;
    vector<long> sites;
    getNumber(linestream,somestring);

    while(!somestring.empty())
    {
        long position = FlagCheck(somestring,"marker position");
        sites.push_back(position);
        skipWhiteSpace(linestream);
        somestring.erase();
        getNumber(linestream,somestring);
    }

    return sites;

} // HapConverter::ParsePhaseInfo

//------------------------------------------------------------------------------------

vector<string> HapConverter::PopHapNames(long ntopop)
{
    vector<string>::iterator start = m_hapnames.begin(),
        stop = m_hapnames.begin()+ntopop;
    vector<string> hnames(start,stop);
    m_hapnames.erase(start,stop);

    return hnames;

} // HapConverter::PopHapNames

//------------------------------------------------------------------------------------

string HapConverter::PopTilDelimiter(istream& input, const char& dlm) const
{
    char ch(input.get());
    string buffer;

    while((ch != EOF) && (ch != dlm))
    {
        buffer += ch;
        ch = input.get();
    }

    if (ch == EOF) input.putback(ch);

    return buffer;

} // HapConverter::PopTilDelimiter

//------------------------------------------------------------------------------------

IndDSVec HapConverter::ParseFirstLine(istringstream& firstline, const string& filename, ifstream& filestr)
{
    string somestring;
    IndDSVec individuals;

    if (getNumber(firstline,somestring))
    {
        m_nindividuals = FlagCheck(somestring, string("number of individuals"));
    }
    else
    {
        string errormsg = "The file " + filename + " is incorrectedly ";
        errormsg += "formated.\nThe first token must be an integer ";
        errormsg += "equal to the number of individuals in the file.";
        throw FileFormatError(errormsg);
    }

    //  Its possible that there will be an adjacency key.
    somestring.erase();
    if (getWord(firstline, somestring))
    {
        if (CaselessStrCmp(somestring,"adjacent"))
        {
            long nhaps;
            // adjacent sequences are haplotypes
            somestring.erase();
            if (getNumber(firstline,somestring))
                nhaps = FlagCheck(somestring,string("number of haplotypes"));
            else
            {
                string errormsg = "The file " + filename + " is ";
                errormsg += "incorrectly formated.\n";
                errormsg += "The converter expected to find the";
                errormsg += " global number of haplotypes for each";
                errormsg += " individual.";
                throw FileFormatError(errormsg);
            }

            // now get phase information for all individuals
            long nsites = m_region.getNmarkers();
            vector<long> psites(nsites);
            somestring.erase();
            if (getWord(firstline,somestring))
            {
                if (CaselessStrCmp(somestring,"all"))
                {
                    for(long site=0; site < nsites; ++site)
                        psites[site] = site;
                }
                else
                {
                    string errormsg = "The file " + filename + " is ";
                    errormsg += "incorrectly formated.\nThe unknown ";
                    errormsg += "keyword, " + somestring + ", was ";
                    errormsg += "encountered.";
                    throw FileFormatError(errormsg);
                }
            }
            else
            {
                psites = ParsePhaseInfo(filestr);
            }

            unsigned long totalhaps = m_nindividuals*nhaps;
            if (m_hapnames.size() != totalhaps)
            {
                string errormsg = "The first line of " + filename;
                errormsg += " specifies " + ToString(m_nindividuals);
                errormsg += " individuals\nwith " + ToString(nhaps);
                errormsg += " haplotypes each (requiring a total of ";
                errormsg += ToString(totalhaps) + " haplotypes).\nThe";
                errormsg += " genetic data provides only ";
                errormsg += ToString(m_hapnames.size()) + " haplotypes.";
                throw FileFormatError(errormsg);
            }

            // now setup the individuals
            for(int ind = 0; ind < m_nindividuals; ++ind)
            {
                string randomname = m_random.Name();
                IndividualDS individual(randomname);
                vector<string> hnames(PopHapNames(nhaps));
                individual.AddHap(hnames);
                individual.SetUnknownPhase(psites);
                individuals.push_back(individual);
            }
        }
        else
        {
            string errormsg = "The file " + filename + " is ";
            errormsg += "incorrectly formated.\nThe unknown ";
            errormsg += "keyword, " + somestring + ", was ";
            errormsg += "encountered.";
            throw FileFormatError(errormsg);
        }
    }

    return individuals;

} // HapConverter::ParseFirstLine

//------------------------------------------------------------------------------------

IndDSVec HapConverter::ReadHapInfo(const string& filename)
{
    ifstream haplotypefile(filename.c_str(),ios::in);

    if (!haplotypefile)
    {
        string errormsg = "Couldn't find the file: " + filename;
        throw FileFormatError (errormsg);
    }

    string line;

    if (!getLine(haplotypefile,line))
    {
        string errormsg = "The file " + filename + " appears to";
        errormsg += " be empty.";
        throw FileFormatError (errormsg);
    }

    istringstream linestream(line);

    IndDSVec individuals(ParseFirstLine(linestream,filename,haplotypefile));

    if (!individuals.empty()) return individuals;

    // read in the haplotype name delimiter character;
    line.erase();
    getLine(haplotypefile,line);
    const char hapname_dlm(line[0]);

    // EWFIX.P3 -- need to insist that hapname_dlm is not a digit

    // Now start parsing the individuals
    for(int ind = 0; ind < m_nindividuals; ++ind)
    {
        line.erase();
        if (getLine(haplotypefile,line))
        {
            istringstream linestr(line);
            string somestring;
            getName(linestr,somestring);
            IndividualDS individual(somestring);

            somestring.erase();
            getNumber(linestr,somestring);
            long int nhapnames = FlagCheck(somestring,string("number of haplotypes"));

            for(int hname = 0; hname < nhapnames; ++hname)
            {
                somestring = PopTilDelimiter(linestr,hapname_dlm);
                StripLeadingSpaces(somestring);
                StripTrailingSpaces(somestring);
                individual.AddHap(somestring);
            }

            vector<long> sites(ParsePhaseInfo(haplotypefile));
            individual.SetUnknownPhase(sites);

            individuals.push_back(individual);
        }
        else
        {
            string errormsg = "The haplotypefile is incorrectedly formated.\n";
            errormsg += "The converter expected to find ";
            errormsg += ToString(m_nindividuals) + " individuals ";
            errormsg += "but only found " + indexToKey(ind);
            throw FileFormatError(errormsg);
        }
    }

    return individuals;

} // HapConverter::ReadHapInfo

//------------------------------------------------------------------------------------

void HapConverter::ReplaceIndividualsWith(IndDSVec& individuals)
{

    vector<string> popnames(individuals.size());
    unsigned long whichind;
    for(whichind = 0; whichind < individuals.size(); ++whichind)
    {
        // can't do population assignment in FailToAdd due to identical
        // sequence names possibly present
        if (m_region.FailToAdd(individuals[whichind],popnames[whichind]))
        {
            string errormsg = "The haplotypes, ";
            errormsg += individuals[whichind].GetHapNamesForPrint();
            errormsg += " were not part of the data set for the";
            errormsg += " region, " + m_region.getName();
            throw FileFormatError(errormsg);
        }
    }

    for(whichind = 0; whichind < individuals.size(); ++whichind)
    {
        m_region.AddIndividual(individuals[whichind],popnames[whichind]);
    }

} // HapConverter::ReplaceIndividualsWith

//------------------------------------------------------------------------------------

void HapConverter::addConvertedLamarcDS(LamarcDS&)
{
    // DEBUG debug -- currently a no-op function!
    // lamarc.mergeTo(m_lamarc);
    assert(false);  // this code never ought to be called
} // HapConverter::addConvertedLamarcDS

//____________________________________________________________________________________