src/genome.cpp

/*****************************************************************************
* This file is provided under the Creative Commons Attribution 3.0 license.
*
* You are free to share, copy, distribute, transmit, or adapt this work
* PROVIDED THAT you attribute the work to the authors listed below.
* For more information, please see the following web page:
* http://creativecommons.org/licenses/by/3.0/
*
* This file is a component of the Sleipnir library for functional genomics,
* authored by:
* Curtis Huttenhower (chuttenh@princeton.edu)
* Mark Schroeder
* Maria D. Chikina
* Olga G. Troyanskaya (ogt@princeton.edu, primary contact)
*
* If you use this library, the included executable tools, or any related
* code in your work, please cite the following publication:
* Curtis Huttenhower, Mark Schroeder, Maria D. Chikina, and
* Olga G. Troyanskaya.
* "The Sleipnir library for computational functional genomics"
*****************************************************************************/
#include "stdafx.h"
#include "genome.h"
#include "meta.h"
#include "annotation.h"

namespace Sleipnir {

CGene::CGene( const std::string& strID ) : CGeneImpl(strID) { }

CGeneImpl::CGeneImpl( const std::string& strName ) : m_strName(strName),
    m_iOntologies(0), m_apOntologies(NULL), m_apveciAnnotations(NULL),
    m_iSynonyms(0), m_astrSynonyms(NULL), m_fRNA(false), m_fDubious(false) { }

CGeneImpl::~CGeneImpl( ) {
    size_t  i;

    if( m_iOntologies ) {
        delete[] m_apOntologies;
        for( i = 0; i < m_iOntologies; ++i )
            delete m_apveciAnnotations[ i ];
        delete[] m_apveciAnnotations; }
    if( m_iSynonyms )
        delete[] m_astrSynonyms; }

CGeneImpl& CGeneImpl::operator=( const CGeneImpl& Gene ) {
    size_t  i, j;

    m_strName = Gene.m_strName;
    if( m_iSynonyms = Gene.m_iSynonyms ) {
        m_astrSynonyms = new string[ m_iSynonyms ];
        for( i = 0; i < m_iSynonyms; ++i )
            m_astrSynonyms[ i ] = Gene.m_astrSynonyms[ i ]; }
    if( m_iOntologies = Gene.m_iOntologies ) {
        m_apOntologies = new IOntology const*[ m_iOntologies ];
        m_apveciAnnotations = new vector<size_t>*[ m_iOntologies ];
        for( i = 0; i < m_iOntologies; ++i ) {
            m_apOntologies[ i ] = Gene.m_apOntologies[ i ];
            m_apveciAnnotations[ i ] = new vector<size_t>( );
            m_apveciAnnotations[ i ]->resize( Gene.m_apveciAnnotations[ i ]->size( ) );
            for( j = 0; j < m_apveciAnnotations[ i ]->size( ); ++j )
                (*m_apveciAnnotations[ i ])[ j ] = (*Gene.m_apveciAnnotations[ i ])[ j ]; } }

    return *this; }

bool CGene::AddAnnotation( const IOntology* pOntology, size_t iTerm ) {
    size_t  i, iOnto;

    for( iOnto = 0; iOnto < m_iOntologies; ++iOnto )
        if( m_apOntologies[ iOnto ] == pOntology )
            break;
    if( iOnto >= m_iOntologies )
        IncrementOntologies( pOntology );

    for( i = 0; i < m_apveciAnnotations[ iOnto ]->size( ); ++i )
        if( (*m_apveciAnnotations[ iOnto ])[ i ] == iTerm )
            return false;
    m_apveciAnnotations[ iOnto ]->push_back( iTerm );
    return true; }

void CGeneImpl::IncrementOntologies( const IOntology* pOntology ) {
    const IOntology**   apOntologies;
    vector<size_t>**    apveciAnnotations;

    apOntologies = new IOntology const*[ m_iOntologies + 1 ];
    if( m_apOntologies ) {
        memcpy( apOntologies, m_apOntologies, m_iOntologies * sizeof(*m_apOntologies) );
        delete[] m_apOntologies; }
    apOntologies[ m_iOntologies ] = pOntology;
    m_apOntologies = apOntologies;

    apveciAnnotations = new vector<size_t>*[ m_iOntologies + 1 ];
    if( m_apveciAnnotations ) {
        memcpy( apveciAnnotations, m_apveciAnnotations, m_iOntologies *
            sizeof(*m_apveciAnnotations) );
        delete[] m_apveciAnnotations; }
    apveciAnnotations[ m_iOntologies++ ] = new vector<size_t>;
    m_apveciAnnotations = apveciAnnotations; }

bool CGene::IsAnnotated( const IOntology* pOntology ) const {
    size_t  i;

    for( i = 0; i < m_iOntologies; ++i )
        if( pOntology == m_apOntologies[ i ] )
            return true;

    return false; }

bool CGene::IsAnnotated( const IOntology* pOntology, size_t iTerm ) const {
    size_t  i, j;

    for( i = 0; i < m_iOntologies; ++i )
        if( pOntology == m_apOntologies[ i ] ) {
            for( j = 0; j < m_apveciAnnotations[ i ]->size( ); ++j )
                if( iTerm == (*m_apveciAnnotations[ i ])[ j ] )
                    return true;
            break; }

    return false; }

bool CGene::AddSynonym( const std::string& strName ) {
    size_t  i;
    string* astrSynonyms;

    if( !strName.length( ) )
        g_CatSleipnir( ).warn( "CGene::AddSynonym( %s ) adding null synonym to %s",
            strName.c_str( ), m_strName.c_str( ) );
    if( strName == m_strName )
        return false;
    for( i = 0; i < m_iSynonyms; ++i )
        if( strName == m_astrSynonyms[ i ] )
            return false;

    astrSynonyms = new string[ ++m_iSynonyms ];
    for( i = 0; ( i + 1 ) < m_iSynonyms; ++i )
        astrSynonyms[ i ] = m_astrSynonyms[ i ];
    astrSynonyms[ i ] = strName;

    if( m_astrSynonyms )
        delete[] m_astrSynonyms;
    m_astrSynonyms = astrSynonyms;
    return true; }

bool CGene::SetWeight(float weight){
    m_weight = weight;
    return true;
};



const char  CGenomeImpl::c_szDubious[]  = "Dubious";
const char  CGenomeImpl::c_szORF[]      = "ORF";
const char* CGenomeImpl::c_aszRNA[]     = { "ncRNA", "rRNA", "snRNA", "snoRNA", "tRNA",
    NULL };

CGenomeImpl::~CGenomeImpl( ) {
    size_t  i;

    for( i = 0; i < m_vecpGenes.size( ); ++i )
        delete m_vecpGenes[ i ]; }

bool CGenome::Open( std::istream& istmFeatures ) {
    static const size_t c_iBuf  = 4096;
    char            szBuf[ c_iBuf ];
    vector<string>  vecstrLine, vecstrNames;
    size_t          i;

    while( istmFeatures.peek( ) != EOF ) {
        istmFeatures.getline( szBuf, c_iBuf - 1 );
        vecstrLine.clear( );
        CMeta::Tokenize( szBuf, vecstrLine );
        if( vecstrLine.size( ) < 16 )
            return false;

        if( vecstrLine[ 1 ] != c_szORF ) {
            for( i = 0; c_aszRNA[ i ]; ++i )
                if( vecstrLine[ 1 ] == c_aszRNA[ i ] )
                    break;
            if( !c_aszRNA[ i ] )
                continue; }

        {
            CGene&  Gene    = AddGene( vecstrLine[ 3 ] );

//  1   Type        ORF, CDS, RNA, etc.
//  2   Qualifier   Dubious, Verified, etc.
//  3   ORF
//  4   Name
//  5   Aliases
//  15  DESC
            Gene.SetRNA( vecstrLine[ 1 ] != c_szORF );
            Gene.SetDubious( vecstrLine[ 2 ] == c_szDubious );
            if( vecstrLine[ 4 ].length( ) )
                AddSynonym( Gene, vecstrLine[ 4 ] );
            if( vecstrLine[ 5 ].length( ) ) {
                vecstrNames.clear( );
                CMeta::Tokenize( vecstrLine[ 5 ].c_str( ), vecstrNames, "|" );
                for( i = 0; i < vecstrNames.size( ); ++i )
                    AddSynonym( Gene, vecstrNames[ i ] ); }
            Gene.SetGloss( vecstrLine[ 15 ] );
        } }

    return true; }

bool CGenome::Open( const std::vector<std::string>& vecstrGenes ) {
    size_t  i;

    for( i = 0; i < vecstrGenes.size( ); ++i )
        AddGene( vecstrGenes[ i ] );

    return true; }

bool CGenome::Open( const char* szFile, std::vector<CGenes*>& vecpGenes ) {
    ifstream    ifsm;

    if( !szFile )
        return Open( cin, vecpGenes );

    ifsm.open( szFile );
    return ( ifsm.is_open( ) ? Open( ifsm, vecpGenes ) : false ); }

bool CGenome::Open( std::istream& istmGenes, std::vector<CGenes*>& vecpGenes ) {
    char            szLine[ c_iBufferSize ];
    string          strLine;
    vector<string>  vecstrLine;
    CGenes*         pGenes;

    while( istmGenes.peek( ) != EOF ) {
        istmGenes.getline( szLine, c_iBufferSize - 1 );
        szLine[ c_iBufferSize - 1 ] = 0;
        if( !( strLine = CMeta::Trim( szLine ) ).length( ) )
            continue;
        vecstrLine.clear( );
        CMeta::Tokenize( strLine.c_str( ), vecstrLine );
        pGenes = new CGenes( *this );
        if( !pGenes->Open( vecstrLine ) ) {
            delete pGenes;
            g_CatSleipnir( ).error( "CGenome::Open( ) could not open line: %s", szLine );
            return false; }
        vecpGenes.push_back( pGenes ); }

    return true; }

CGene& CGenome::AddGene( const std::string& strID ) {
    TMapStrI::const_iterator    iterGene;
    CGene*                      pGene;

    if( ( iterGene = m_mapGenes.find( strID ) ) != m_mapGenes.end( ) )
        return GetGene( iterGene->second );

    pGene = new CGene( strID );
    m_vecpGenes.push_back( pGene );
    m_mapGenes[ strID ] = m_vecpGenes.size( ) - 1;
    return *pGene; }

size_t CGenome::FindGene( const std::string& strGene ) const {
    size_t  i, j, iRet;

    if( ( iRet = GetGene( strGene ) ) != -1 )
        return iRet;

    for( i = 0; i < GetGenes( ); ++i )
        for( j = 0; j < GetGene( i ).GetSynonyms( ); ++j )
            if( strGene == GetGene( i ).GetSynonym( j ) )
                return i;

    return -1; }

vector<string> CGenome::GetGeneNames( ) const {
    vector<string>  vecstrRet;
    size_t          i;

    vecstrRet.resize( m_vecpGenes.size( ) );
    for( i = 0; i < vecstrRet.size( ); ++i )
        vecstrRet[ i ] = m_vecpGenes[ i ]->GetName( );

    return vecstrRet; }

size_t CGenome::CountGenes( const IOntology* pOntology ) const {
    size_t  i, j, iRet;

    for( iRet = i = 0; i < m_vecpGenes.size( ); ++i )
        for( j = 0; j < m_vecpGenes[ i ]->GetOntologies( ); ++j )
            if( pOntology == m_vecpGenes[ i ]->GetOntology( j ) ) {
                iRet++;
                break; }

    return iRet; }

bool CGenome::AddSynonym( CGene& Gene, const std::string& strName ) {

    if( ( strName != Gene.GetName( ) ) && ( Gene.AddSynonym( strName ) ) ) {
        m_mapGenes[ strName ] = m_mapGenes[ Gene.GetName( ) ];
        return true; }

    return false; }

bool CGenes::Open( const char* szFile, CGenome& Genome, std::vector<std::string>& vecstrNames, std::vector<CGenes*>& vecpGenes ) {
    ifstream        ifsm;
    vector<char>    veccBuffer;
    istream*        pistm;
    bool            fRet;

    if( szFile ) {
        ifsm.open( szFile );
        pistm = &ifsm; }
    else
        pistm = &cin;
    if( !ifsm.is_open( ) ) {
        g_CatSleipnir( ).error( "CGenes::Open( %s ) could not open file", szFile ? szFile : "stdin" );
        return false; }

    veccBuffer.resize( CFile::GetBufferSize( ) );
    fRet = false;
    while( !pistm->eof( ) ) {
        vector<string>  vecstrLine;

        pistm->getline( &veccBuffer[0], veccBuffer.size( ) - 1 );
        CMeta::Tokenize( &veccBuffer[0], vecstrLine );
        if( vecstrLine.empty( ) )
            continue;
        vecstrNames.push_back( vecstrLine[0] );
        vecstrLine.erase( vecstrLine.begin( ) );
        vecpGenes.push_back( new CGenes( Genome ) );
        if( !( fRet = vecpGenes.back( )->Open( vecstrLine ) ) )
            break; }
    if( szFile )
        ifsm.close( );

    return fRet; }

CGenes::CGenes( CGenome& Genome ) : CGenesImpl( Genome ) { }

CGenesImpl::CGenesImpl( CGenome& Genome ) : m_Genome(Genome),isWeighted(false) { }

bool CGenes::Open( std::istream& istm, bool fCreate ) {
    static const size_t c_iBuffer   = 1024;
    char    szBuf[ c_iBuffer ];
    CGene*  pGene;
    size_t  i, iGene;
    char*   pc;

    if( istm.rdstate( ) != ios_base::goodbit )
        return false;

    m_mapGenes.clear( );
    m_vecpGenes.clear( );
    while( istm.peek( ) != EOF ) {
        istm.getline( szBuf, c_iBuffer - 1 );
        if( pc = strchr( szBuf, '\t' ) )
            *pc = 0;
        if( !szBuf[ 0 ] || ( szBuf[ 0 ] == c_cComment ) )
            continue;
        if( fCreate )
            pGene = &m_Genome.AddGene( szBuf );
        else {
            if( ( iGene = m_Genome.FindGene( szBuf ) ) == -1 ) {
                g_CatSleipnir( ).warn( "CGenes::Open( %d ) unknown gene: %s", fCreate, szBuf );
                continue; }
            pGene = &m_Genome.GetGene( iGene ); }
        for( i = 0; i < m_vecpGenes.size( ); ++i )
            if( m_vecpGenes[ i ] == pGene )
                break;
        if( i != m_vecpGenes.size( ) )
            continue;
        m_mapGenes[ pGene->GetName( ) ] = m_vecpGenes.size( );
        m_vecpGenes.push_back( pGene ); }
    isWeighted = false;
    return true; }

bool CGenes::OpenWeighted( std::istream& istm, bool fCreate ) {
    static const size_t c_iBuffer   = 1024;
    char    szBuf[ c_iBuffer ];
    CGene*  pGene;
    size_t  i, iGene;
    char*   pc;
    vector<string> vecstrTokens;

    if( istm.rdstate( ) != ios_base::goodbit )
        return false;

    m_mapGenes.clear( );
    m_vecpGenes.clear( );
    while( istm.peek( ) != EOF ) {
        istm.getline( szBuf, c_iBuffer - 1 );
        //if( pc = strchr( szBuf, '\t' ) )
        //  *pc = 0;
        if( !szBuf[ 0 ] || ( szBuf[ 0 ] == c_cComment ) )
            continue;
        szBuf[c_iBuffer - 1] = 0;
        vecstrTokens.clear();
        CMeta::Tokenize(szBuf, vecstrTokens);
        if (vecstrTokens.empty())
            continue;
        if (vecstrTokens.size() != 2) {
            //cerr << "Illegal label line (" << vecstrTokens.size() << "): "
                //  << szBuf << endl;
            return false;
        }

        if( fCreate )
            pGene = &m_Genome.AddGene( vecstrTokens[0] );
        else {
            if( ( iGene = m_Genome.FindGene( vecstrTokens[0] ) ) == -1 ) {
                g_CatSleipnir( ).warn( "CGenes::Open( %d ) unknown gene: %s", fCreate, vecstrTokens[0].c_str() );
                continue; }
            pGene = &m_Genome.GetGene( iGene ); }
        pGene->SetWeight(atof(vecstrTokens[1].c_str()));
        for( i = 0; i < m_vecpGenes.size( ); ++i )
            if( m_vecpGenes[ i ] == pGene )
                break;
        if( i != m_vecpGenes.size( ) )
            continue;
        m_mapGenes[ pGene->GetName( ) ] = m_vecpGenes.size( );
        m_vecpGenes.push_back( pGene ); }
    isWeighted = true;
    return true; }


size_t CGenes::CountAnnotations( const IOntology* pOntology, size_t iTerm, bool fRecursive,
    const CGenes* pBackground ) const {
    size_t  i, iRet;

    for( iRet = i = 0; i < m_vecpGenes.size( ); ++i )
        if( ( !pBackground || pBackground->IsGene( m_vecpGenes[ i ]->GetName( ) ) ) &&
            pOntology->IsAnnotated( iTerm, *m_vecpGenes[ i ], fRecursive ) )
            iRet++;

    return iRet; }

bool CGenes::Open( const std::vector<std::string>& vecstrGenes, bool fCreate ) {
    size_t  i, iGene;
    CGene*  pGene;

    m_mapGenes.clear( );
    m_vecpGenes.clear( );
    for( i = 0; i < vecstrGenes.size( ); ++i ) {
        if( !fCreate && ( ( iGene = m_Genome.FindGene( vecstrGenes[ i ] ) ) == -1 ) )
            continue;

        pGene = fCreate ? &m_Genome.AddGene( vecstrGenes[ i ] ) : &m_Genome.GetGene( iGene );
        m_mapGenes[ vecstrGenes[ i ] ] = m_vecpGenes.size( );
        m_vecpGenes.push_back( pGene ); }
    isWeighted = false;
    return true; }

void CGenes::Filter( const CGenes& GenesExclude ) {
    size_t  i, j, iSize;

    iSize = m_vecpGenes.size( );
    for( i = 0; i < GenesExclude.GetGenes( ); ++i )
        for( j = 0; j < iSize; ++j )
            if( m_vecpGenes[ j ] == &GenesExclude.GetGene( i ) ) {
                m_mapGenes.erase( m_vecpGenes[ j ]->GetName( ) );
                m_vecpGenes[ j ] = m_vecpGenes[ m_vecpGenes.size( ) - 1 ];
                iSize--;
                break; }
    m_vecpGenes.resize( iSize ); }

vector<string> CGenes::GetGeneNames( ) const {
    vector<string>  vecstrRet;
    size_t          i;

    vecstrRet.resize( m_vecpGenes.size( ) );
    for( i = 0; i < vecstrRet.size( ); ++i )
        vecstrRet[ i ] = m_vecpGenes[ i ]->GetName( );

    return vecstrRet; }

}