include/search.h File Reference

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
void	search_protein1hit (struct PSSMatrix PSSMatrix, struct sequenceData *sequenceData, uint4 numSequences, uint4 tickFrequency)
void	search_protein2hit (struct PSSMatrix PSSMatrix, struct sequenceData *sequenceData, uint4 numSequences, uint4 tickFrequency)
void	search_nucleotide (struct PSSMatrix PSSMatrix, struct sequenceData *sequenceData, uint4 numSequences, uint4 tickFrequency)
void	search_nucleotide_longWord (struct PSSMatrix PSSMatrix, struct sequenceData *sequenceData, uint4 numSequences, uint4 tickFrequency)
void	search_nucleotide_largeTable (struct PSSMatrix PSSMatrix, struct sequenceData *sequenceData, uint4 numSequences, uint4 tickFrequency)
void	search_nucleotideIndex (unsigned char *address, unsigned char *endAddress, int4 tickFrequency, struct PSSMatrix PSSMatrix)
void	search_proteinSsearch (struct PSSMatrix PSSMatrix, struct sequenceData *sequenceData, uint4 numSequences, uint4 tickFrequency)
void	search_nucleotideSsearch (struct PSSMatrix PSSMatrix, struct sequenceData *sequenceData, uint4 numSequences, uint4 tickFrequency)

---

Function Documentation

void search_nucleotide	(	struct PSSMatrix	PSSMatrix,
		struct sequenceData *	sequenceData,
		uint4	numSequences,
		uint4	tickFrequency
	)

Definition at line 328 of file search.c.

References alignments_addUngappedExtension(), alignments_createNew(), alignments_currentAlignment, blast_numHits, blast_numTriggerExtensions, blast_numUngappedExtensions, PSSMatrix::bytePackedCodes, sequenceData::descriptionLength, sequenceData::descriptionStart, sequenceData::encodedLength, hitMatrix_diagonalMask, hitMatrix_furthest, hitMatrix_queryLength, int2, int4, nucleotideLookup, nucleotideLookup_additionalPositions, nucleotideLookup_bitLookup, nucleotideLookup_mask, parameters_outputType, parameters_tabular, parameters_verboseDloc, parameters_wordExtraLetters, parameters_wordSize, parameters_wordTableBytes, parameters_wordTableLetters, parameters_xml, PSSMatrix_packedLeftMatches, PSSMatrix_packedRightMatches, sequenceData::sequence, sequenceData::sequenceLength, uint4, ungappedExtension_bestScore, ungappedExtension_nucleotideExtend(), ungappedExtension_print(), and ungappedExtension_subjectEndReached.

Referenced by blast_search().

{
        uint4 sequenceCount = 0, codeword;
        uint4 descriptionStart = 0, descriptionLength = 0, encodedLength;
    uint4 numPackedBytes, numRemaining;
        unsigned char *subject, *sequenceEnd, prevByte, nextByte, rightMatches, leftMatches;
    unsigned char previousCode, *subjectPosition, *address;
        int4 subjectLength, subjectOffset, wordLengthMinusOne;
        struct ungappedExtension* ungappedExtension;
        int2 *queryOffsets, tableEntry, queryOffset, singleQueryOffset[2];
    int4 previousByteDistance;
    int4 diagonal;
    unsigned char** lastHit;

    wordLengthMinusOne = parameters_wordSize - 1;
        previousByteDistance = (parameters_wordTableLetters + 4);
        singleQueryOffset[1] = 0;

    while (sequenceCount < numSequences)
    {
        descriptionLength = sequenceData[sequenceCount].descriptionLength;
        descriptionStart = sequenceData[sequenceCount].descriptionStart;
        subjectLength = sequenceData[sequenceCount].sequenceLength;
        encodedLength = sequenceData[sequenceCount].encodedLength;
                address = subject = sequenceData[sequenceCount].sequence;

        // New sequence, new possible alignment
        alignments_currentAlignment = NULL;

        // Determine number of packed bytes and number of remaining letters
        numPackedBytes = subjectLength / 4;
                numRemaining = subjectLength % 4;

        // Only process sequence if at least as long as the word length
        if (subjectLength >= parameters_wordSize)
        {
            sequenceEnd = subject + numPackedBytes;

            // Read first char and advance
            previousCode = *address;
            address++;

            // Traverse until end of sequence
            while (address < sequenceEnd)
            {
                // Read next char and update codeword
                codeword = *address | (previousCode << 8);
                previousCode = *address;
                                tableEntry = nucleotideLookup[codeword];

                // Calculate subject offset
                subjectOffset = address - subject;

                #ifdef BITLOOKUP
                                if (nucleotideLookup_bitLookup[codeword >> 5] &
                    (1 << (codeword & nucleotideLookup_mask)))
                                #endif
                if (tableEntry)
                {
                        if (tableEntry > 0)
                    {
                        // Just one query position
                        singleQueryOffset[0] = queryOffset = tableEntry;
                                                queryOffsets = singleQueryOffset;
                    }
                    else
                    {
                        // Multiple query positions
                                                queryOffsets = nucleotideLookup_additionalPositions - tableEntry;
                        queryOffset = *queryOffsets;
                    }

                    subjectPosition = address + 1;
                    do
                    {
                        // Calculate number of matches to left and right of hit
                        nextByte = PSSMatrix.bytePackedCodes[queryOffset];
                        rightMatches = PSSMatrix_packedRightMatches[nextByte ^ *subjectPosition];
                        prevByte = PSSMatrix.bytePackedCodes[queryOffset - previousByteDistance];
                        leftMatches = PSSMatrix_packedLeftMatches[prevByte ^ *(address - parameters_wordTableBytes)];

                        if (rightMatches + leftMatches >= parameters_wordExtraLetters)
                        {
                            #ifndef NO_STAGE2
                            // Calculate subject offset
                            subjectOffset = subjectPosition - subject;

                            // Calculate the diagonal this hit is on
                            diagonal = (subjectOffset * 4 - queryOffset + hitMatrix_queryLength)
                                     & hitMatrix_diagonalMask;

                            // If we have not extended past this point on this diagonal
                            lastHit = hitMatrix_furthest + diagonal;

                            #ifdef VERBOSE
                            if (parameters_verboseDloc == descriptionStart)
                                printf("Hit %d,%d\n", queryOffset, subjectOffset * 4);
                            #endif

                            if (*lastHit < address)
                            {
                                // Perform ungapped extension
                                ungappedExtension
                                        = ungappedExtension_nucleotideExtend(queryOffset,
                                      subjectOffset, PSSMatrix, subject, subjectLength);

                                // Update furthest reached value for the diagonal
                                *lastHit = ungappedExtension_subjectEndReached;

                                #ifdef VERBOSE
                                if (parameters_verboseDloc == descriptionStart)
                                    printf("UngappedExtension %d,%d Score=%d\n", queryOffset,
                                    subjectOffset * 4, ungappedExtension_bestScore);
                                if (parameters_verboseDloc == descriptionStart && ungappedExtension)
                                        ungappedExtension_print(ungappedExtension);
                                #endif

                                // If extension scores high enough to trigger gapping
                                if (ungappedExtension)
                                {
                                    // Increment count of number of trigger extensions
                                    blast_numTriggerExtensions++;

                                    // Create new alignment if needed
                                    if (alignments_currentAlignment == NULL)
                                    {
                                        alignments_createNew(descriptionStart, descriptionLength, subject,
                                                             subjectLength, encodedLength);
                                    }

                                    // Add this extension to the alignment
                                    alignments_addUngappedExtension(ungappedExtension);
                                }

                                blast_numUngappedExtensions++;
                            }
                            #endif
                            blast_numHits++;
                        }

                        queryOffsets++;
                        queryOffset = *queryOffsets;
                                        }
                    while (queryOffset);
                }

                address++;
            }
        }

        sequenceCount++;

        // Every so often print status .
        if ((sequenceCount % tickFrequency == 0) &&
            parameters_outputType != parameters_xml && parameters_outputType != parameters_tabular)
        {
            #ifndef VERBOSE
            printf(".");
            fflush(stdout);
            #endif
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void search_nucleotide_largeTable	(	struct PSSMatrix	PSSMatrix,
		struct sequenceData *	sequenceData,
		uint4	numSequences,
		uint4	tickFrequency
	)

Definition at line 683 of file search.c.

References alignments_addUngappedExtension(), alignments_createNew(), alignments_currentAlignment, blast_numHits, blast_numTriggerExtensions, blast_numUngappedExtensions, PSSMatrix::bytePackedCodes, sequenceData::descriptionLength, sequenceData::descriptionStart, sequenceData::encodedLength, hitMatrix_diagonalMask, hitMatrix_furthest, hitMatrix_queryLength, int4, nucleotideLookup_additionalPositions_large, nucleotideLookup_bitLookup, nucleotideLookup_large, nucleotideLookup_mask, parameters_outputType, parameters_tabular, parameters_verboseDloc, parameters_wordExtraLetters, parameters_wordSize, parameters_wordTableBytes, parameters_wordTableLetters, parameters_xml, PSSMatrix_packedLeftMatches, PSSMatrix_packedRightMatches, sequenceData::sequence, sequenceData::sequenceLength, uint4, ungappedExtension_bestScore, ungappedExtension_nucleotideExtend(), ungappedExtension_print(), and ungappedExtension_subjectEndReached.

Referenced by blast_search().

{
        uint4 sequenceCount = 0, codeword;
        uint4 descriptionStart = 0, descriptionLength = 0, encodedLength;
    uint4 numPackedBytes, numRemaining;
        unsigned char *subject, *sequenceEnd, prevByte, nextByte, rightMatches, leftMatches;
    unsigned char previousCode, *subjectPosition, *address;
        int4 subjectLength, subjectOffset, wordLengthMinusOne;
        struct ungappedExtension* ungappedExtension;
        int4 *queryOffsets, tableEntry, queryOffset, singleQueryOffset[2];
    int4 previousByteDistance;
    int4 diagonal;
    unsigned char** lastHit;

    wordLengthMinusOne = parameters_wordSize - 1;
        previousByteDistance = (parameters_wordTableLetters + 4);
        singleQueryOffset[1] = 0;

    while (sequenceCount < numSequences)
    {
        descriptionLength = sequenceData[sequenceCount].descriptionLength;
        descriptionStart = sequenceData[sequenceCount].descriptionStart;
        subjectLength = sequenceData[sequenceCount].sequenceLength;
        encodedLength = sequenceData[sequenceCount].encodedLength;
                address = subject = sequenceData[sequenceCount].sequence;

        // New sequence, new possible alignment
        alignments_currentAlignment = NULL;

        // Determine number of packed bytes and number of remaining letters
        numPackedBytes = subjectLength / 4;
                numRemaining = subjectLength % 4;

        // Only process sequence if at least as long as the word length
        if (subjectLength >= parameters_wordSize)
        {
            sequenceEnd = subject + numPackedBytes;

            // Read first char and advance
            previousCode = *address;
            address++;

            // Traverse until end of sequence
            while (address < sequenceEnd)
            {
                // Read next char and update codeword
                codeword = *address | (previousCode << 8);
                previousCode = *address;
                                tableEntry = nucleotideLookup_large[codeword];

                // Calculate subject offset
                subjectOffset = address - subject;

                #ifdef BITLOOKUP
                                if (nucleotideLookup_bitLookup[codeword >> 5] &
                    (1 << (codeword & nucleotideLookup_mask)))
                                #endif
                if (tableEntry)
                {
                        if (tableEntry > 0)
                    {
                        // Just one query position
                        singleQueryOffset[0] = queryOffset = tableEntry;
                                                queryOffsets = singleQueryOffset;
                    }
                    else
                    {
                        // Multiple query positions
                                                queryOffsets = nucleotideLookup_additionalPositions_large - tableEntry;
                        queryOffset = *queryOffsets;
                    }

                    subjectPosition = address + 1;
                    do
                    {
                        // Calculate number of matches to left and right of hit
                        nextByte = PSSMatrix.bytePackedCodes[queryOffset];
                        rightMatches = PSSMatrix_packedRightMatches[nextByte ^ *subjectPosition];
                        prevByte = PSSMatrix.bytePackedCodes[queryOffset - previousByteDistance];
                        leftMatches = PSSMatrix_packedLeftMatches[prevByte ^ *(address - parameters_wordTableBytes)];

                        if (rightMatches + leftMatches >= parameters_wordExtraLetters)
                        {
                            #ifndef NO_STAGE2
                            // Calculate subject offset
                            subjectOffset = subjectPosition - subject;

                            // Calculate the diagonal this hit is on
                            diagonal = (subjectOffset * 4 - queryOffset + hitMatrix_queryLength)
                                     & hitMatrix_diagonalMask;

                            // If we have not extended past this point on this diagonal
                            lastHit = hitMatrix_furthest + diagonal;

                            #ifdef VERBOSE
                            if (parameters_verboseDloc == descriptionStart)
                                printf("Hit %d,%d\n", queryOffset, subjectOffset * 4);
                            #endif
                            if (*lastHit < address)
                            {
                                // Perform ungapped extension
                                ungappedExtension
                                        = ungappedExtension_nucleotideExtend(queryOffset,
                                      subjectOffset, PSSMatrix, subject, subjectLength);

                                // Update furthest reached value for the diagonal
                                *lastHit = ungappedExtension_subjectEndReached;

                                #ifdef VERBOSE
                                if (parameters_verboseDloc == descriptionStart)
                                    printf("UngappedExtension %d,%d Score=%d\n", queryOffset,
                                    subjectOffset * 4, ungappedExtension_bestScore);
                                if (parameters_verboseDloc == descriptionStart && ungappedExtension)
                                        ungappedExtension_print(ungappedExtension);
                                #endif

                                // If extension scores high enough to trigger gapping
                                if (ungappedExtension)
                                {
                                    // Increment count of number of trigger extensions
                                    blast_numTriggerExtensions++;

                                    // Create new alignment if needed
                                    if (alignments_currentAlignment == NULL)
                                    {
                                        alignments_createNew(descriptionStart, descriptionLength, subject,
                                                             subjectLength, encodedLength);
                                    }

                                    // Add this extension to the alignment
                                    alignments_addUngappedExtension(ungappedExtension);
                                }

                                blast_numUngappedExtensions++;
                            }
                            #endif
                            blast_numHits++;
                        }

                        queryOffsets++;
                        queryOffset = *queryOffsets;
                                        }
                    while (queryOffset);
                }

                address++;
            }
        }

        sequenceCount++;

        // Every so often print status .
        if ((sequenceCount % tickFrequency == 0) &&
            parameters_outputType != parameters_xml && parameters_outputType != parameters_tabular)
        {
            #ifndef VERBOSE
            printf(".");
            fflush(stdout);
            #endif
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void search_nucleotide_longWord	(	struct PSSMatrix	PSSMatrix,
		struct sequenceData *	sequenceData,
		uint4	numSequences,
		uint4	tickFrequency
	)

Definition at line 494 of file search.c.

References alignments_addUngappedExtension(), alignments_createNew(), alignments_currentAlignment, blast_numHits, blast_numTriggerExtensions, blast_numUngappedExtensions, PSSMatrix::bytePackedCodes, sequenceData::descriptionLength, sequenceData::descriptionStart, sequenceData::encodedLength, hitMatrix_diagonalMask, hitMatrix_furthest, hitMatrix_queryLength, int2, int4, nucleotideLookup, nucleotideLookup_additionalPositions, nucleotideLookup_bitLookup, nucleotideLookup_mask, parameters_outputType, parameters_tabular, parameters_verboseDloc, parameters_wordExtraBytes, parameters_wordExtraLetters, parameters_wordSize, parameters_wordTableBytes, parameters_wordTableLetters, parameters_xml, PSSMatrix_packedLeftMatches, PSSMatrix_packedRightMatches, sequenceData::sequence, sequenceData::sequenceLength, uint4, ungappedExtension_bestScore, ungappedExtension_nucleotideExtend(), ungappedExtension_print(), and ungappedExtension_subjectEndReached.

Referenced by blast_search().

{
        uint4 sequenceCount = 0, codeword;
        uint4 descriptionStart = 0, descriptionLength = 0, encodedLength;
    uint4 numPackedBytes, numRemaining;
        unsigned char *subject, *sequenceEnd, prevByte, nextByte, rightMatches, leftMatches;
    unsigned char previousCode, *subjectPosition, *address;
        int4 subjectLength, subjectOffset, wordLengthMinusOne;
        struct ungappedExtension* ungappedExtension;
        int2 *queryOffsets, tableEntry, queryOffset, singleQueryOffset[2];
    int4 previousByteDistance;
    int4 diagonal, extraBytesNeeded;
    unsigned char** lastHit;

    wordLengthMinusOne = parameters_wordSize - 1;
        previousByteDistance = parameters_wordTableLetters + parameters_wordExtraBytes * 4 + 4;
        singleQueryOffset[1] = 0;

    while (sequenceCount < numSequences)
    {
        descriptionLength = sequenceData[sequenceCount].descriptionLength;
        descriptionStart = sequenceData[sequenceCount].descriptionStart;
        subjectLength = sequenceData[sequenceCount].sequenceLength;
        encodedLength = sequenceData[sequenceCount].encodedLength;
                address = subject = sequenceData[sequenceCount].sequence;

        // Determine number of packed bytes and number of remaining letters
        numPackedBytes = subjectLength / 4;
                numRemaining = subjectLength % 4;

        // Only process sequence if at least as long as the word length
        if (subjectLength >= parameters_wordSize)
        {
            sequenceEnd = subject + numPackedBytes;

            // Read first char and advance
            previousCode = *address;
            address++;

            // Traverse until end of sequence
            while (address < sequenceEnd)
            {
                // Read next char and update codeword
                codeword = *address | (previousCode << 8);
                previousCode = *address;
                                tableEntry = nucleotideLookup[codeword];

                // Calculate subject offset
                subjectOffset = address - subject;

                #ifdef BITLOOKUP
                                if (nucleotideLookup_bitLookup[codeword >> 5] &
                    (1 << (codeword & nucleotideLookup_mask)))
                                #endif
                if (tableEntry)
                {
                        if (tableEntry > 0)
                    {
                        // Just one query position
                        singleQueryOffset[0] = queryOffset = tableEntry;
                                                queryOffsets = singleQueryOffset;
                    }
                    else
                    {
                        // Multiple query positions
                                                queryOffsets = nucleotideLookup_additionalPositions - tableEntry;
                        queryOffset = *queryOffsets;
                    }

                    subjectPosition = address + 1;
                    do
                    {
                                                extraBytesNeeded = parameters_wordExtraBytes;

                                                while (extraBytesNeeded)
                                                {
                                                        // Check for matching bytes to right
                            if (*subjectPosition != PSSMatrix.bytePackedCodes[queryOffset])
                                break;

//                                                      printf("Match %d,%d\n", queryOffset, (subjectPosition - subject)*4);

                            extraBytesNeeded--;
                            subjectPosition++;
                            subjectOffset++;
                            queryOffset+=4;
                                                }

                                                if (!extraBytesNeeded)
                                                {
                            // Calculate number of matches to left and right of hit
                            nextByte = PSSMatrix.bytePackedCodes[queryOffset];
                            rightMatches = PSSMatrix_packedRightMatches[nextByte ^ *(subjectPosition)];
                            prevByte = PSSMatrix.bytePackedCodes[queryOffset - previousByteDistance];
                            leftMatches = PSSMatrix_packedLeftMatches[prevByte ^ *(address - parameters_wordTableBytes)];

//                            printf("prev at %d,%d\n", queryOffset - previousByteDistance);
//                              printf("part matches=[%d,%d]\n", leftMatches, rightMatches);

                            if (rightMatches + leftMatches >= parameters_wordExtraLetters)
                            {
//                              printf("Hit! dloc=%d\n", descriptionStart);
                                #ifndef NO_STAGE2
                                // Calculate subject offset
                                subjectOffset = subjectPosition - subject;

                                // Calculate the diagonal this hit is on
                                diagonal = (subjectOffset * 4 - queryOffset + hitMatrix_queryLength)
                                        & hitMatrix_diagonalMask;

                                // If we have not extended past this point on this diagonal
                                lastHit = hitMatrix_furthest + diagonal;

                                #ifdef VERBOSE
                                if (parameters_verboseDloc == descriptionStart)
                                    printf("Hit %d,%d\n", queryOffset, subjectOffset * 4);
                                #endif
                                if (*lastHit < address)
                                {
                                    // Perform ungapped extension
                                    ungappedExtension
                                        = ungappedExtension_nucleotideExtend(queryOffset,
                                        subjectOffset, PSSMatrix, subject, subjectLength);

                                    // Update furthest reached value for the diagonal
                                    *lastHit = ungappedExtension_subjectEndReached;

                                    #ifdef VERBOSE
                                    if (parameters_verboseDloc == descriptionStart)
                                        printf("UngappedExtension %d,%d Score=%d\n", queryOffset,
                                        subjectOffset * 4, ungappedExtension_bestScore);
                                    if (parameters_verboseDloc == descriptionStart && ungappedExtension)
                                        ungappedExtension_print(ungappedExtension);
                                    #endif

                                    // If extension scores high enough to trigger gapping
                                    if (ungappedExtension)
                                    {
                                        // Increment count of number of trigger extensions
                                        blast_numTriggerExtensions++;

                                        // Create new alignment if needed
                                        if (alignments_currentAlignment == NULL)
                                        {
                                            alignments_createNew(descriptionStart, descriptionLength, subject,
                                                                subjectLength, encodedLength);
                                        }

                                        // Add this extension to the alignment
                                        alignments_addUngappedExtension(ungappedExtension);
                                    }

                                    blast_numUngappedExtensions++;
                                }
                                #endif
                                blast_numHits++;
                            }
                        }

                        queryOffsets++;
                        queryOffset = *queryOffsets;
                                        }
                    while (queryOffset);
                }

                address++;
            }
        }

        sequenceCount++;

        // Every so often print status .
        if ((sequenceCount % tickFrequency == 0) &&
            parameters_outputType != parameters_xml && parameters_outputType != parameters_tabular)
        {
            #ifndef VERBOSE
            printf(".");
            fflush(stdout);
            #endif
        }

        // New sequence, new possible alignment
        alignments_currentAlignment = NULL;
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void search_nucleotideIndex	(	unsigned char *	address,
		unsigned char *	endAddress,
		int4	tickFrequency,
		struct PSSMatrix	PSSMatrix
	)

void search_nucleotideSsearch	(	struct PSSMatrix	PSSMatrix,
		struct sequenceData *	sequenceData,
		uint4	numSequences,
		uint4	tickFrequency
	)

Definition at line 1105 of file search.c.

References alignments_addFinalAlignment(), alignments_addGappedExtension(), alignments_createNew(), alignments_currentAlignment, alignments_isFinalAlignment(), dpResults::best, dpResults::bestScore, blast_gappedNominalCutoff, sequenceData::descriptionLength, sequenceData::descriptionStart, sequenceData::encodedLength, encoding_byteUnpack(), encoding_insertWilds(), gappedExtension_score(), int4, gappedExtension::nominalScore, parameters_outputType, parameters_tabular, parameters_xml, gappedExtension::queryEnd, coordinate::queryOffset, sequenceData::sequence, sequenceData::sequenceLength, smithWatermanScoring_score(), smithWatermanScoring_scoreReverse(), smithWatermanTraceback_build(), gappedExtension::subjectEnd, coordinate::subjectOffset, and uint4.

Referenced by blast_search().

{
        uint4 sequenceCount = 0;
        uint4 descriptionStart = 0, descriptionLength = 0, encodedLength;
        unsigned char *subject, *unpackedSubject, *address;
        int4 subjectLength;
        struct gappedExtension* gappedExtension;
        struct dpResults dpResults, reverseDpResults;

    while (sequenceCount < numSequences)
    {
        descriptionLength = sequenceData[sequenceCount].descriptionLength;
        descriptionStart = sequenceData[sequenceCount].descriptionStart;
        subjectLength = sequenceData[sequenceCount].sequenceLength;
        encodedLength = sequenceData[sequenceCount].encodedLength;
                address = subject = sequenceData[sequenceCount].sequence;

        // New sequence, new possible alignment
        alignments_currentAlignment = NULL;

        // Unpack the subject
        unpackedSubject = encoding_byteUnpack(subject, subjectLength);

        // Re-insert wildcards
        encoding_insertWilds(unpackedSubject, subject + ((subjectLength + 3) / 4),
                             subject + encodedLength);

                // Perform SW score only
                dpResults = smithWatermanScoring_score(PSSMatrix, subjectLength, unpackedSubject);

        // If above e-value cutoff
                if (dpResults.bestScore >= blast_gappedNominalCutoff
            && alignments_isFinalAlignment(dpResults.bestScore))
                {
                // Perform SW alignment in the reverse direction, to find the start of the optimal alignment
                reverseDpResults = smithWatermanScoring_scoreReverse(PSSMatrix, subjectLength,
                               unpackedSubject, dpResults.best);

            // Collect traceback information and store alignment
                        gappedExtension = smithWatermanTraceback_build(PSSMatrix, subjectLength, unpackedSubject,
                                                           reverseDpResults.best, dpResults.best);

            if (reverseDpResults.bestScore != dpResults.bestScore ||
                dpResults.bestScore != gappedExtension->nominalScore)
            {
                fprintf(stderr, "Error: Forward and reverse Smith-Waterman alignment scores do not match\n");
                fprintf(stderr, "Forward Score=%d End=%d,%d\n", dpResults.bestScore,
                        dpResults.best.queryOffset, dpResults.best.subjectOffset);
                fprintf(stderr, "Reverse Score=%d End=%d,%d\n", reverseDpResults.bestScore,
                        reverseDpResults.best.queryOffset, reverseDpResults.best.subjectOffset);
                fprintf(stderr, "Traceback Score=%d End=%d,%d\n", gappedExtension->nominalScore,
                        gappedExtension->queryEnd, gappedExtension->subjectEnd);
//                exit(-1);
            }

            gappedExtension_score(gappedExtension);

                        alignments_createNew(descriptionStart, descriptionLength, unpackedSubject, subjectLength, 0);
                        alignments_addGappedExtension(alignments_currentAlignment, gappedExtension);
                        alignments_addFinalAlignment(dpResults.bestScore, alignments_currentAlignment);
                }
                else
        {
                free(unpackedSubject);
        }

        sequenceCount++;

        // Every so often print status .
        if ((sequenceCount % tickFrequency == 0) &&
            parameters_outputType != parameters_xml && parameters_outputType != parameters_tabular)
        {
            #ifndef VERBOSE
            printf(".");
            fflush(stdout);
            #endif
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void search_protein1hit	(	struct PSSMatrix	PSSMatrix,
		struct sequenceData *	sequenceData,
		uint4	numSequences,
		uint4	tickFrequency
	)

Definition at line 11 of file search.c.

References alignments_addUngappedExtension(), alignments_createNew(), alignments_currentAlignment, blast_numHits, blast_numTriggerExtensions, blast_numUngappedExtensions, constants_max_int2, sequenceData::descriptionLength, sequenceData::descriptionStart, sequenceData::encodedLength, hitMatrix_furthest, int4, PSSMatrix::matrix, group::nextGroups, group::nextWords, parameters_outputType, parameters_tabular, parameters_wordSize, parameters_xml, sequenceData::sequence, sequenceData::sequenceLength, uint2, uint4, ungappedExtension_oneHitExtend(), ungappedExtension_subjectEndReached, wordLookupDFA_additionalQueryPositions, wordLookupDFA_groups, and wordLookupDFA_numCodes.

Referenced by alignments_expandCluster(), and blast_search().

{
        uint4 sequenceCount = 0;
        uint4 descriptionStart = 0, descriptionLength = 0, encodedLength;
        unsigned char *subject, *sequenceEnd, *address;
        int4 subjectLength, subjectOffset, wordLengthMinusOne, count = 0;
        unsigned char currentWord, *currentBlock;
    struct group* currentGroup;
        uint2* queryOffsets, queryOffset;
        struct ungappedExtension* ungappedExtension;
    int4 diagonal;
    unsigned char** lastHit;

    wordLengthMinusOne = parameters_wordSize - 1;

    while (sequenceCount < numSequences)
    {
        descriptionLength = sequenceData[sequenceCount].descriptionLength;
        descriptionStart = sequenceData[sequenceCount].descriptionStart;
        subjectLength = sequenceData[sequenceCount].sequenceLength;
        encodedLength = sequenceData[sequenceCount].encodedLength;
                address = subject = sequenceData[sequenceCount].sequence;

        // New sequence, new possible alignment
        alignments_currentAlignment = NULL;

        // Only process sequence if at least as long as the word length
        if (subjectLength >= parameters_wordSize)
        {
            // Start at 000 state in Finite State Automata
            currentGroup = wordLookupDFA_groups;

            // Read first wordLength - 1 chars and advance
            count = wordLengthMinusOne;
            while (count > 0)
            {
                if (*address < wordLookupDFA_numCodes)
                    currentGroup = currentGroup->nextGroups + *address;
                else
                    currentGroup = currentGroup->nextGroups;

                address++;
                count--;
            }

            // Read the rest of the codes, using the Finite State Automata defined
            // by wordLookupDFA to get the query positions of int4erest
            sequenceEnd = subject + subjectLength;
            while (address < sequenceEnd)
            {
                currentBlock = currentGroup->nextWords;

                // If current code is a regular letter
                if (*address < wordLookupDFA_numCodes)
                {
                    // Use it
                    currentWord = currentBlock[*address];
                    currentGroup = currentGroup->nextGroups + *address;
                }
                else
                {
                    // Else check if we've reached end of the file
                    if (address >= sequenceEnd)
                        break;

                    // If not, we've read a wild code. Use first code instead
                    currentWord = *currentBlock;
                    currentGroup = currentGroup->nextGroups;
                }

                if (currentWord)
                {
                    // Calculate subject offset
                    subjectOffset = address - subject;

                    // At least one query position, stored at an extenal address
                    queryOffsets = ((uint2*)currentBlock) - currentWord;

                    // If the zero flag is stored at the first query position
                    if (!*queryOffsets)
                    {
                        // Go to an outside address for additional positions
                        queryOffsets = wordLookupDFA_additionalQueryPositions
                                     + (*(queryOffsets + 1) * constants_max_int2) + *(queryOffsets + 2);
                    }

                    do
                    {
                        queryOffset = *queryOffsets;

                        #ifndef NO_STAGE2
                        // Calculate the diagonal this hit is on
                        diagonal = subjectOffset - queryOffset;

                        // If we have not extended past this point on this diagonal
                        lastHit = hitMatrix_furthest + diagonal;
                        if (*lastHit < address)
                        {
                            // Increment tally number of extensions
                            blast_numUngappedExtensions++;

                            // If only one hit triggered this extension
                            ungappedExtension
                                = ungappedExtension_oneHitExtend(PSSMatrix.matrix + queryOffset,
                                                                 address, PSSMatrix, subject);

                            // Update furthest reached value for the diagonal
                            *lastHit = ungappedExtension_subjectEndReached;

                            // If extension scores high enough to trigger gapping
                            if (ungappedExtension)
                            {
                                // Increment count of number of trigger extensions
                                blast_numTriggerExtensions++;

                                // Create new alignment if needed
                                if (alignments_currentAlignment == NULL)
                                {
                                    // Create new alignment object
                                    alignments_createNew(descriptionStart, descriptionLength,
                                                         subject, subjectLength, encodedLength);
                                }

                                // Add this extension to the alignment
                                alignments_addUngappedExtension(ungappedExtension);
                            }
                        }
                        #endif

                        queryOffsets++; blast_numHits++;
                    }
                    while (*queryOffsets);
                }

                address++;
            }
        }

        sequenceCount++;

        // Every so often print status .
        if ((sequenceCount % tickFrequency == 0) &&
            parameters_outputType != parameters_xml && parameters_outputType != parameters_tabular)
        {
            #ifndef VERBOSE
            printf(".");
            fflush(stdout);
                        #endif
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void search_protein2hit	(	struct PSSMatrix	PSSMatrix,
		struct sequenceData *	sequenceData,
		uint4	numSequences,
		uint4	tickFrequency
	)

Definition at line 165 of file search.c.

References alignments_addUngappedExtension(), alignments_createNew(), alignments_currentAlignment, blast_numHits, blast_numTriggerExtensions, blast_numUngappedExtensions, constants_max_int2, sequenceData::descriptionLength, sequenceData::descriptionStart, sequenceData::encodedLength, hitMatrix_furthest, int4, PSSMatrix::matrix, group::nextGroups, group::nextWords, parameters_A, parameters_outputType, parameters_overlap, parameters_tabular, parameters_wordSize, parameters_xml, sequenceData::sequence, sequenceData::sequenceLength, uint2, uint4, ungappedExtension_extend(), ungappedExtension_subjectEndReached, wordLookupDFA_additionalQueryPositions, wordLookupDFA_groups, and wordLookupDFA_numCodes.

Referenced by alignments_expandCluster(), and blast_search().

{
        uint4 sequenceCount = 0;
        uint4 descriptionStart = 0, descriptionLength = 0, encodedLength;
        unsigned char *subject, *sequenceEnd;
        int4 subjectLength, subjectOffset, wordLengthMinusOne, count = 0;
        unsigned char currentWord, *currentBlock, *address;
    struct group* currentGroup;
        uint2* queryOffsets, queryOffset;
        struct ungappedExtension* ungappedExtension;
    int4 diagonal, distance;
    unsigned char** lastHit;

    wordLengthMinusOne = parameters_wordSize - 1;

    while (sequenceCount < numSequences)
    {
        descriptionLength = sequenceData[sequenceCount].descriptionLength;
        descriptionStart = sequenceData[sequenceCount].descriptionStart;
        subjectLength = sequenceData[sequenceCount].sequenceLength;
        encodedLength = sequenceData[sequenceCount].encodedLength;
                address = subject = sequenceData[sequenceCount].sequence;

        // New sequence, new possible alignment
        alignments_currentAlignment = NULL;

        // Only process sequence if at least as long as the word length
        if (subjectLength >= parameters_wordSize)
        {
            // Start at 000 state in Finite State Automata
            currentGroup = wordLookupDFA_groups;

            // Read first wordLength - 1 chars and advance
            count = wordLengthMinusOne;
            while (count > 0)
            {
                if (*address < wordLookupDFA_numCodes)
                    currentGroup = currentGroup->nextGroups + *address;
                else
                    currentGroup = currentGroup->nextGroups;

                address++;
                count--;
            }

            // Read the rest of the codes, using the Finite State Automata defined
            // by wordLookupDFA to get the query positions of int4erest
            sequenceEnd = subject + subjectLength;
            while (address < sequenceEnd)
            {
                currentBlock = currentGroup->nextWords;

                // If current code is a regular letter
                if (*address < wordLookupDFA_numCodes)
                {
                    // Use it
                    currentWord = currentBlock[*address];
                    currentGroup = currentGroup->nextGroups + *address;
                }
                else
                {
                    // Else check if we've reached end of the file
                    if (address >= sequenceEnd)
                        break;

                    // If not, we've read a wild code. Use first code instead
                    currentWord = currentBlock[0];
                    currentGroup = currentGroup->nextGroups;
                }

                if (currentWord)
                {
                    // Calculate subject offset
                    subjectOffset = address - subject;

                    // If at least one query position, stored at an extenal address
                    queryOffsets = ((uint2*)currentBlock) - currentWord;

                    // If the zero flag is stored at the first query position
                    if (!*queryOffsets)
                    {
                        // Go to an outside address for additional positions
                        queryOffsets = wordLookupDFA_additionalQueryPositions
                                     + (*(queryOffsets + 1) * constants_max_int2) + *(queryOffsets + 2);
                    }

                    do
                    {
                        queryOffset = *queryOffsets;

                        #ifndef NO_STAGE2
                        // Calculate the diagonal this hit is on
                        diagonal = subjectOffset - queryOffset;

                        // Calculate distance since last hit
                        lastHit = hitMatrix_furthest + diagonal;
                        distance = address - *lastHit;

                        if (distance >= parameters_A)
                        {
                            // Too far apart, update furthest
                            *lastHit = address;
                        }
                        else if (distance >= parameters_overlap)
                        {
                            // Not overlaping - extension triggered
                            // Increment tally number of extensions
                            blast_numUngappedExtensions++;

                            // Perform ungapped extension start between query/subject start/end
                            // and extending outwards in each direction
                            ungappedExtension
                                = ungappedExtension_extend(PSSMatrix.matrix + queryOffset,
                                                           address, *lastHit, PSSMatrix, subject);

                            // Update furthest reached value for the diagonal
                            *lastHit = ungappedExtension_subjectEndReached;

                            // If extension scores high enough to trigger gapping
                            if (ungappedExtension)
                            {
                                // Increment count of number of trigger extensions
                                blast_numTriggerExtensions++;

                                // Create new alignment if needed
                                if (alignments_currentAlignment == NULL)
                                {
                                    // Create new alignment object using subject with wilds
                                    alignments_createNew(descriptionStart, descriptionLength, subject,
                                                         subjectLength, encodedLength);
                                }

                                // Add this extension to the alignment
                                alignments_addUngappedExtension(ungappedExtension);
                            }
                        }
                        #endif

                        queryOffsets++; blast_numHits++;
                    }
                    while (*queryOffsets);
                }

                address++;
            }
        }

        sequenceCount++;

        // Every so often print status .
        if ((sequenceCount % tickFrequency == 0) &&
            parameters_outputType != parameters_xml && parameters_outputType != parameters_tabular)
        {
            #ifndef VERBOSE
            printf(".");
            fflush(stdout);
                        #endif
        }
    }
}

Here is the call graph for this function:

Here is the caller graph for this function:

void search_proteinSsearch	(	struct PSSMatrix	PSSMatrix,
		struct sequenceData *	sequenceData,
		uint4	numSequences,
		uint4	tickFrequency
	)

Definition at line 1032 of file search.c.

References alignments_addFinalAlignment(), alignments_addGappedExtension(), alignments_createNew(), alignments_currentAlignment, alignments_isFinalAlignment(), alignments_sortFinalAlignments(), dpResults::best, dpResults::bestScore, blast_gappedNominalCutoff, sequenceData::descriptionLength, sequenceData::descriptionStart, sequenceData::encodedLength, gappedExtension_score(), int4, gappedExtension::nominalScore, parameters_outputType, parameters_tabular, parameters_xml, sequenceData::sequence, sequenceData::sequenceLength, smithWatermanScoring_score(), smithWatermanScoring_scoreReverse(), smithWatermanTraceback_build(), and uint4.

Referenced by blast_search().

{
        uint4 sequenceCount = 0;
        uint4 descriptionStart = 0, descriptionLength = 0, encodedLength;
        unsigned char *subject, *address;
        int4 subjectLength;
        struct gappedExtension* gappedExtension;
        struct dpResults dpResults, reverseDpResults;

    while (sequenceCount < numSequences)
    {
        descriptionLength = sequenceData[sequenceCount].descriptionLength;
        descriptionStart = sequenceData[sequenceCount].descriptionStart;
        subjectLength = sequenceData[sequenceCount].sequenceLength;
        encodedLength = sequenceData[sequenceCount].encodedLength;
                address = subject = sequenceData[sequenceCount].sequence;

        // New sequence, new possible alignment
        alignments_currentAlignment = NULL;

        // Perform SW score only
                dpResults = smithWatermanScoring_score(PSSMatrix, subjectLength, subject);

//        printf("%d\n", dpResults.bestScore);

        // If above e-value cutoff
                if (dpResults.bestScore >= blast_gappedNominalCutoff
            && alignments_isFinalAlignment(dpResults.bestScore))
                {
            // Perform SW alignment in the reverse direction, to find the start of the optimal alignment
                reverseDpResults = smithWatermanScoring_scoreReverse(PSSMatrix, subjectLength, subject,
                                                                 dpResults.best);

            // Collect traceback information and store alignment
                        gappedExtension = smithWatermanTraceback_build(PSSMatrix, subjectLength, subject,
                                                           reverseDpResults.best, dpResults.best);

            if (reverseDpResults.bestScore != dpResults.bestScore ||
                dpResults.bestScore != gappedExtension->nominalScore)
            {
                fprintf(stderr, "Error: Forward and reverse Smith-Waterman alignment scores do not match\n");
                exit(-1);
            }

            gappedExtension_score(gappedExtension);

                        alignments_createNew(descriptionStart, descriptionLength, subject, subjectLength, encodedLength);
                        alignments_addGappedExtension(alignments_currentAlignment, gappedExtension);
                        alignments_addFinalAlignment(dpResults.bestScore, alignments_currentAlignment);
                }

        // Advance to next sequence
        address += encodedLength - 1;

        sequenceCount++;

        // Every so often print status .
        if ((sequenceCount % tickFrequency == 0) &&
            parameters_outputType != parameters_xml && parameters_outputType != parameters_tabular)
        {
            #ifndef VERBOSE
            printf(".");
            fflush(stdout);
                        #endif
        }
    }

        // Sort the final alignments by refined score
        alignments_sortFinalAlignments();
}

Here is the call graph for this function:

Here is the caller graph for this function:

---