C++ Genome::GetChrIdFromContigId方法代码示例

//.........这里部分代码省略.........
                    for ( i = 0 ; i < tag ; ++i )
                        newNParts.push_back( nparts[i] ) ;
                    if ( newNParts.size() > 1 )
                        connects.push_back( newNParts ) ;

                    newNParts.clear() ;
                    for ( ; i < size ; ++i  )
                        newNParts.push_back( nparts[i] ) ;
                    if ( newNParts.size() > 1 )
                        connects.push_back( newNParts ) ;

                    tag = 0 ;
                }
            }
        }
        fclose( rascafFile ) ;
    }
    if ( contigLevel == false )
    {
        genome.SetIsOpen( contigLevel ) ;
    }

    // Build the graph
    int contigCnt = genome.GetContigCount() ;
    int edgeCnt = 0 ;
    int csize = connects.size() ;

    for ( i = 0 ; i < csize ; ++i )
        edgeCnt += connects[i].size() ;

    ContigGraph contigGraph( contigCnt, contigCnt + edgeCnt ) ;
    for ( i = 0 ; i < contigCnt - 1 ; ++i )
    {
        if ( genome.GetChrIdFromContigId( i ) == genome.GetChrIdFromContigId( i + 1 ) )
        {
            contigGraph.AddEdge( i, 1, i + 1, 0 ) ;
        }
    }
    for ( i = 0 ; i < csize ; ++i )
    {
        std::vector<struct _part> &parts = connects[i] ;
        int size = parts.size() ;
        for ( int j = 0 ; j < size - 1 ; ++j )
        {
            struct _part &a = parts[j] ;
            struct _part &b = parts[j + 1] ;

            // Two dummy nodes for each contig. Left is 0, right is 1
            int dummyU = 0 ;
            int dummyV = 0 ;
            if ( a.strand == '+' )
                dummyU = 1 ;
            if ( b.strand == '-' )
                dummyV = 1 ;
            contigGraph.AddEdge( a.contigId, dummyU, b.contigId, dummyV, true ) ;
        }
    }

    // Check the cycles in the contig graph. This may introduces when combining different rascaf outputs.
    int *visitTime = new int[contigCnt] ;
    struct _pair *neighbors = new struct _pair[ MAX_NEIGHBOR ] ;

    bool *isInCycle = new bool[contigCnt] ;
    std::vector<int> cycleNodes ;
    memset( visitTime, -1, sizeof( int ) * contigCnt ) ;
    memset( isInCycle, false, sizeof( bool ) * contigCnt ) ;

// Search one dangling paths
void SearchDangling( int u, int inDummy, bool *used, int time, int *visitTime, ContigGraph &contigGraph, bool add, int chosenNodes[], int chosenDummyNodes[], int &chosenCnt, Genome &genome )
{
    //if ( u == 10 )
    //	printf( "hi\n" ) ;
    if ( visitTime[u] == time )
        return ;
    visitTime[u] = time ;
    if ( add ) // Only one path will be added.
    {
        chosenNodes[ chosenCnt ] = u ;
        chosenDummyNodes[ chosenCnt ] = inDummy ;
        ++chosenCnt ;
    }
    //if ( u == 10246 )
    //	printf( "hi2 %d\n", chosenCnt ) ;
    struct _pair *neighbors = new struct _pair[ MAX_NEIGHBOR ] ;
    int ncnt ;
    int i ;
    ncnt = contigGraph.GetNeighbors( u, 1 - inDummy, neighbors, MAX_NEIGHBOR ) ;
    bool newAdd = true ;
    //if ( u == 10 )
    //	printf( "hi3 %d %d\n", ncnt, 1 - inDummy ) ;

    // Firstly, bias towards a connection onto another scaffold
    for ( i = 0 ; i < ncnt ; ++i )
    {
        if ( used[ neighbors[i].a ] )
            continue ;
        if ( genome.GetChrIdFromContigId( u ) != genome.GetChrIdFromContigId( neighbors[i].a ) )
        {
            SearchDangling( neighbors[i].a, neighbors[i].b, used, time, visitTime,
                            contigGraph, newAdd, chosenNodes, chosenDummyNodes, chosenCnt, genome ) ;
            newAdd = false ;
        }
    }

    // Lastly, bias towards the closest contig in the raw assembly, since we know the connections are on the same scaffold
    int direction ;
    if ( inDummy == 1 )
        direction = -1 ;
    else
        direction = 1 ;

    int min = genome.GetContigCount() + 1 ;
    int mintag = -1 ;

    for ( i = 0 ; i < ncnt ; ++i )
    {
        //if ( u == 10 )
        //	printf( "%d\n", i ) ;
        if ( used[ neighbors[i].a ] )
            continue ;
        //if ( u == 34674 && neighbors[i].a == 144159 )
        //	printf( "hi\n" ) ;
        //SearchDangling( neighbors[i].a, neighbors[i].b, used, time, visitTime, contigGraph, newAdd, chosenNodes, chosenDummyNodes, chosenCnt, genome ) ;
        //newAdd = false ;
        int tmp = direction * ( neighbors[i].a - u ) ;
        if ( tmp < min )
        {
            min = tmp ;
            mintag = i ;
        }
    }

    if ( mintag != -1 )
    {
        SearchDangling( neighbors[ mintag ].a, neighbors[ mintag ].b, used, time, visitTime,
                        contigGraph, newAdd, chosenNodes, chosenDummyNodes, chosenCnt, genome ) ;
        newAdd = false ;
    }
    delete[] neighbors ;
}