/***************************************************************************** # Copyright (C) 1994-2008 by David Gordon. # All rights reserved. # # This software is part of a beta-test version of the Consed/Autofinish # package. It should not be redistributed or # used for any commercial purpose, including commercially funded # sequencing, without written permission from the author and the # University of Washington. # # This software is provided ``AS IS'' and any express or implied # warranties, including, but not limited to, the implied warranties of # merchantability and fitness for a particular purpose, are disclaimed. # In no event shall the authors or the University of Washington be # liable for any direct, indirect, incidental, special, exemplary, or # consequential damages (including, but not limited to, procurement of # substitute goods or services; loss of use, data, or profits; or # business interruption) however caused and on any theory of liability, # whether in contract, strict liability, or tort (including negligence # or otherwise) arising in any way out of the use of this software, even # if advised of the possibility of such damage. # # Building Consed from source is error prone and not simple which is # why I provide executables. Due to time limitations I cannot # provide any assistance in building Consed. Even if you do not # modify the source, you may introduce errors due to using a # different version of the compiler, a different version of motif, # different versions of other libraries than I used, etc. For this # reason, if you discover Consed bugs, I can only offer help with # those bugs if you first reproduce those bugs with an executable # provided by me--not an executable you have built. # # Modifying Consed is also difficult. Although Consed is modular, # some modules are used by many other modules. Thus making a change # in one place can have unforeseen effects on many other features. # It may takes months for you to notice these other side-effects # which may not seen connected at all. It is not feasable for me to # provide help with modifying Consed sources because of the # potentially huge amount of time involved. # #*****************************************************************************/ #include "createPrimerCandidates.h" #include "findPrimerMatchesElsewhere.h" #include "findPrimerSelfMatches.h" #include "whichPrimersAreAcceptable.h" #include "sortAcceptablePrimers.h" #include "sortUnacceptablePrimers.h" #include "displayPrimers.h" #include "assembly.h" #include "consed.h" #include "calculatePrimerMeltingTemperatures.h" #include #include "findPrimerMatchesAgainstSequencesInFile.h" #include "copyPrimerMatchesFromShorterToLongerPrimers.h" #include "consedParameters.h" #include "checkPrimersForMononucleotideRepeats.h" #include "findTemplatesForPrimers.h" #include "checkPrimersForACGT.h" #include "countsOfWhyPrimersAreNotAcceptable.h" #include "popupErrorMessage.h" #include "textbox.h" #ifdef NO_POUND_POUND_MACROS #define messageIfChanged( PARAMETER ) \ if ( pOriginal->PARAMETER != \ pCurrent->PARAMETER ) { \ soParametersChanged += #PARAMETER " was changed from "; \ soParametersChanged += RWCString( (long) pOriginal->PARAMETER ); \ soParametersChanged += " to "; \ soParametersChanged += RWCString( (long) pCurrent->PARAMETER ); \ soParametersChanged += "\n"; \ bPrimerPickingParameterChanged = true; \ } #else #define messageIfChanged( PARAMETER ) \ if ( pOriginal-> ##PARAMETER != \ pCurrent-> ##PARAMETER ) { \ soParametersChanged += #PARAMETER " was changed from "; \ soParametersChanged += RWCString( (long) pOriginal->##PARAMETER ); \ soParametersChanged += " to "; \ soParametersChanged += RWCString( (long) pCurrent->##PARAMETER ); \ soParametersChanged += "\n"; \ bPrimerPickingParameterChanged = true; \ } #endif static bool bDidPrimerPickingParametersChange( consedParameters* pOriginal, consedParameters* pCurrent, RWCString& soParametersChanged ) { bool bPrimerPickingParameterChanged = false; messageIfChanged( nPrimersMaxMatchElsewhereScoreToUse_ ); messageIfChanged( nPrimersMinQuality_ ); messageIfChanged( nPrimersMaxSelfMatchScore_ ); messageIfChanged( nPrimersMaxLengthOfMononucleotideRepeat_ ); messageIfChanged( nPrimersMinMeltingTempToUse_ ); messageIfChanged( nPrimersMaxMeltingTempToUse_ ); messageIfChanged( nPrimersWindowSizeInLookingToUse_ ); messageIfChanged( bPrimersOKToChoosePrimersInSingleSubcloneRegion_ ); messageIfChanged( bPrimersOKToChoosePrimersWhereHighQualityDiscrepancies_ ); messageIfChanged( bPrimersOKToChoosePrimersWhereUnalignedHighQualityRegion_ ); return( bPrimerPickingParameterChanged ); } void pickPrimers2( Contig* pContigOfCursor, const int nPaddedConsPosOfCursor, const bool bForwardNotReversePrimer, const bool bCloneNotSubcloneTemplate, ContigWin* pContigWin, countsOfWhyPrimersAreNotAcceptable* pCounts, primerType*& pPrimerArray, int& nNumberOfPrimers ) { createPrimerCandidates( pContigOfCursor, nPaddedConsPosOfCursor, bForwardNotReversePrimer, &pPrimerArray, nNumberOfPrimers ); cout << "finding primer matches elsewhere" << endl; // this screen eliminate by far the most // Do it first and then we'll just have 10% or so to apply the // other screens to. findPrimerMatchesElsewhere( pPrimerArray, nNumberOfPrimers, bForwardNotReversePrimer, bCloneNotSubcloneTemplate ); if (consedParameters::pGetConsedParameters()->bPrimersScreenForVector_ ) { cout << "find primer matches against sequences in file" << endl; findPrimerMatchesAgainstSequencesInFile( pPrimerArray, nNumberOfPrimers, bCloneNotSubcloneTemplate ); } copyPrimerMatchesFromShorterToLongerPrimers( pPrimerArray, nNumberOfPrimers, bForwardNotReversePrimer); cout << "calculating primer melting temperatures" << endl; calculatePrimerMeltingTemperatures( pPrimerArray, nNumberOfPrimers ); cout << "finding if primer sticks to itself" << endl; findPrimerSelfMatches( pPrimerArray, nNumberOfPrimers ); int nNumberOfAcceptablePrimers; whichPrimersAreAcceptable( pPrimerArray, nNumberOfPrimers, nNumberOfAcceptablePrimers ); cout << "checking for mononucleotide repeats" << endl; checkPrimersForMononucleotideRepeats( pPrimerArray, nNumberOfPrimers, nNumberOfAcceptablePrimers ); cout << "checking for just acgt in primers" << endl; checkPrimersForACGT( pPrimerArray, nNumberOfPrimers, nNumberOfAcceptablePrimers ); if ( !bCloneNotSubcloneTemplate && consedParameters::pGetConsedParameters()->bPrimersPickTemplatesForPrimers_ ) { cout << "finding templates..."; cout.flush(); pContigOfCursor->setPaddedPositionsArray(); findTemplatesForPrimers( pPrimerArray, nNumberOfPrimers, nNumberOfAcceptablePrimers, bForwardNotReversePrimer ); cout << "done" << endl; cout.flush(); } if ( nNumberOfAcceptablePrimers > 0 ) sortAcceptablePrimers( pPrimerArray, nNumberOfPrimers, bForwardNotReversePrimer ); pCounts->accountForAllPrimers( pPrimerArray, nNumberOfPrimers ); } void pickPrimers( Contig* pContigOfCursor, const int nPaddedConsPosOfCursor, const bool bForwardNotReversePrimer, const bool bCloneNotSubcloneTemplate, ContigWin* pContigWin ) { Assembly* pAssembly = ConsEd::pGetAssembly(); cout << "about to set padded to unpadded conversion table..."; cout.flush(); pAssembly->setContigMatchTablesInAllContigs(); cout << "done" << endl; cout << "about to set unpadded to padded conversion table..."; cout.flush(); pCP->setParametersForSequencingPrimersNotPCRPrimers( true ); // needed by nPaddedIndexFast in displayPrimers.cpp if ( bCloneNotSubcloneTemplate ) pAssembly->setAllPaddedPositionsArrays(); else pContigOfCursor->setPaddedPositionsArray(); cout << "done" << endl; cout.flush(); if ( !bCloneNotSubcloneTemplate && consedParameters::pGetConsedParameters()->bPrimersPickTemplatesForPrimers_ ) { cout << "creating list of subclone templates..."; cout.flush(); pAssembly->setContigTemplateArrays(); cout << "done" << endl; cout.flush(); } // save a copy of the original primer-picking parameters. consedParameters* pOriginalPrimerPickingParameters = new consedParameters(); pOriginalPrimerPickingParameters->copyPrimerPickingParametersToSelf( pCP ); countsOfWhyPrimersAreNotAcceptable* pCounts = new countsOfWhyPrimersAreNotAcceptable(); primerType* pPrimerArray = NULL; // will be returned by pickPrimers2 int nNumberOfPrimers; const int nExitedLoopBecauseTriedTooManyTimes = 1; const int nExitedLoopBecauseFoundPrimers = 2; const int nExitedLoopBecauseCouldNotFindAnyUnacceptablePrimers = 3; const int nStartedLoop = 4; const int nExitedLoopBecauseCouldNotFigureOutWhatToDo = 5; int nReasonExitedLoop = nStartedLoop; const int nMaxTimesTried = 40; int nTimesTried = 0; while( 1 ) { ++nTimesTried; pickPrimers2( pContigOfCursor, nPaddedConsPosOfCursor, bForwardNotReversePrimer, bCloneNotSubcloneTemplate, pContigWin, pCounts, pPrimerArray, nNumberOfPrimers ); if ( pCounts->nAcceptablePrimers_ > 0 ) { nReasonExitedLoop = nExitedLoopBecauseFoundPrimers; break; } if ( nTimesTried >= nMaxTimesTried ) { nReasonExitedLoop = nExitedLoopBecauseTriedTooManyTimes; break; } if ( pCounts->nUnacceptablePrimers_ == 0 ) { // this is the case in which no acceptable or unacceptable primers // could be found if ( pCP->nPrimersWindowSizeInLookingToUse_ < 700 ) { pCP->nPrimersWindowSizeInLookingToUse_ = 700; printf( "changed PrimersWindowSizeInLooking to %d\n", pCP->nPrimersWindowSizeInLookingToUse_ ); continue; } else { nReasonExitedLoop = nExitedLoopBecauseCouldNotFindAnyUnacceptablePrimers; break; } } // if reached here, we need to relax the parameters further // But what parameters? First let's follow Jeremy Schmutz' // recommendations: // // matchelsewhere, 18-20-22 the second would be primer quality. 30 > 20 > 15. // //Also you cna turn off don't pick where HQDs. and the other paramter i // can't rememeber right now ( oh yea pick in single subclone region). // //jeremy // Jeremy Schmutz suggests the following procedure: // 1) increase matchElsewhere to 18 // 2) increase matchElsewhere to 20 // *) decrease quality to 20 // *) allow primers in single subclone regions // *) allow primers where hqd's are // *) allow primers where there are unaligned high quality regions // 3) increase matchElsewhere to 22 // *) change melting temperatures, if this would help // *) increase matchElsewhere to 25 // Kelsi Scott follows the following procedure: // 1) increases length to 15-30 or 15-35 // 2) increases window of melting temp to 55-65 // 3) allows primers in single subclone, hqd, and unaligned high quality // 4) increases MatchElsewhere to 20 or 25 // 5) increases mononucleotide repeats to 5 if ( pCP->nPrimersMaxMatchElsewhereScoreToUse_ < 18 ) { printf( "changing pCP->nPrimersMaxMatchElsewhereScoreToUse_ from %d to 18\n", pCP->nPrimersMaxMatchElsewhereScoreToUse_ ); pCP->nPrimersMaxMatchElsewhereScoreToUse_ = 18; } else if ( pCP->nPrimersMaxMatchElsewhereScoreToUse_ < 20 ) { printf( "changing pCP->nPrimersMaxMatchElsewhereScoreToUse_ from %d to 20\n", pCP->nPrimersMaxMatchElsewhereScoreToUse_ ); pCP->nPrimersMaxMatchElsewhereScoreToUse_ = 20; } else if ( pCP->nPrimersMinQuality_ > 25 ) { printf( "changing pCP->nPrimersMinQuality_ from %d to 25\n", pCP->nPrimersMinQuality_ ); pCP->nPrimersMinQuality_ = 25; } else if ( pCP->nPrimersMinQuality_ > 20 ) { printf( "changing pCP->nPrimersMinQuality_ from %d to 20\n", pCP->nPrimersMinQuality_ ); pCP->nPrimersMinQuality_ = 20; } else if ( !pCP->bPrimersOKToChoosePrimersInSingleSubcloneRegion_ ) { printf( "changed pCP->bPrimersOKToChoosePrimersInSingleSubcloneRegion_ from false to true\n" ); pCP->bPrimersOKToChoosePrimersInSingleSubcloneRegion_ = true; } else if ( !pCP->bPrimersOKToChoosePrimersWhereHighQualityDiscrepancies_ ) { printf( "changing pCP->bPrimersOKToChoosePrimersWhereHighQualityDiscrepancies_ from false to true\n" ); pCP->bPrimersOKToChoosePrimersWhereHighQualityDiscrepancies_ = true; } else if ( !pCP->bPrimersOKToChoosePrimersWhereUnalignedHighQualityRegion_ ) { printf( "changingpCP->bPrimersOKToChoosePrimersWhereUnalignedHighQualityRegion_ from false to true\n" ); pCP->bPrimersOKToChoosePrimersWhereUnalignedHighQualityRegion_ = true; } else if ( pCP->nPrimersMaxMatchElsewhereScoreToUse_ < 22 ) { printf( "changing pCP->nPrimersMaxMatchElsewhereScoreToUse_ from %d to 22\n", pCP->nPrimersMaxMatchElsewhereScoreToUse_ ); pCP->nPrimersMaxMatchElsewhereScoreToUse_ = 22; } else { // see if we can widen the melting temperature range, if // that would help float fUnacceptablePrimers = pCounts->nUnacceptablePrimers_; bool bChangedSomething = false; if ( (float) pCounts->nBAD_PRIMER_MELTING_TEMPERATURE_TOO_LOW_ / fUnacceptablePrimers > 0.0 ) { if ( ABS( pCP->nPrimersMinMeltingTempToUse_ - pOriginalPrimerPickingParameters->nPrimersMinMeltingTempToUse_ ) < 10 ) { pCP->nPrimersMinMeltingTempToUse_ -= 5; printf( "changed PrimersMinMeltingTemp to %d\n", pCP->nPrimersMinMeltingTempToUse_ ); bChangedSomething = true; } } if ( (float) pCounts->nBAD_PRIMER_MELTING_TEMPERATURE_TOO_HIGH_ / fUnacceptablePrimers > 0.0 ) { if ( ABS( pCP->nPrimersMaxMeltingTempToUse_ - pOriginalPrimerPickingParameters->nPrimersMaxMeltingTempToUse_ ) < 10 ) { pCP->nPrimersMaxMeltingTempToUse_ += 5; printf( "changed PrimersMaxMeltingTemp to %d\n", pCP->nPrimersMaxMeltingTempToUse_ ); bChangedSomething = true; } } // don't try anything else until done with trying widening // melting temperature if ( !bChangedSomething ) { if ( pCP->nPrimersMaxMatchElsewhereScoreToUse_ < 25 ) { pCP->nPrimersMaxMatchElsewhereScoreToUse_ = 25; printf( "changed PrimersMaxMatchElsewhereScore to %d\n", pCP->nPrimersMaxMatchElsewhereScoreToUse_ ); bChangedSomething = true; } } if ( !bChangedSomething ) { if ( (float) pCounts->nBAD_PRIMER_MONONUCLEOTIDE_REPEAT_ / fUnacceptablePrimers > 0.0 && pCP->nPrimersMaxLengthOfMononucleotideRepeat_ < 5 ) { pCP->nPrimersMaxLengthOfMononucleotideRepeat_ = 5; printf( "changed PrimersMaxLengthOfMononucleotideRepeat to %d\n", pCP->nPrimersMaxLengthOfMononucleotideRepeat_ ); bChangedSomething = true; } } if (!bChangedSomething ) { // the problem was none of the above. What could it be? popupErrorMessage( "could not figure out what to relax" ); nReasonExitedLoop = nExitedLoopBecauseCouldNotFigureOutWhatToDo; break; } } } // while( 1 ) if ( nReasonExitedLoop == nExitedLoopBecauseFoundPrimers ) { // case of found primers. // Display them. (This assumes that the acceptable primers // are at the beginning of pPrimerArray which is done by // sortAcceptablePrimers in pickPrimers2) displayPrimers( pPrimerArray, pCounts->nAcceptablePrimers_, bForwardNotReversePrimer, bCloneNotSubcloneTemplate, pContigWin ); // Did we have to relax the parameters // to find any? RWCString soParametersChanged; if ( bDidPrimerPickingParametersChange( pOriginalPrimerPickingParameters, pCP, soParametersChanged ) ) { TextBox* pTextBox = new TextBox( "Warning: Had to Relax Parameters", 20 ); RWCString soWarningMessage = "Consed had to relax the following parameters in order to find any primers at all:\n"; soWarningMessage += soParametersChanged; soWarningMessage += "\nThus Consed cannot guarantee that any of these primers will work.\nIf you believe that some other primer will work,\nConsed can tell you what is wrong with that primer\nby doing the following: In the Aligned Reads Window, \npoint to the \"Misc\" menu, hold down the left mouse \nbutton and release on \"Check Primer\". Enter the \nleft and right consensus positions of the primer, \ncheck which strand, and whether the primer is to use \nsubclone templates or the whole clone as a template. \nConsed will tell you all that is wrong with that primer."; pTextBox->append( soWarningMessage ); pTextBox->makeVisible(); } } else if ( nReasonExitedLoop == nExitedLoopBecauseTriedTooManyTimes || nReasonExitedLoop == nExitedLoopBecauseCouldNotFigureOutWhatToDo ) { // could not find any parameters, despite relaxing parameters RWCString soWhyArePrimersUnacceptable; pCounts->getCountsInString( soWhyArePrimersUnacceptable ); RWCString soErrorMessage; // consider case of couldn't find an acceptable template if ( pCounts->nUnacceptablePrimers_ > 0 && (float) pCounts->nBAD_PRIMER_HAS_NO_TEMPLATE_ / (float) pCounts->nUnacceptablePrimers_ > .8 ) { soErrorMessage += "Since most of the primers are unacceptable \nbecause no acceptable templates could be \nfound for them, perhaps you need to try \nwhole clone primers.\n\n"; } soErrorMessage += soWhyArePrimersUnacceptable; RWCString soParametersChanged; if ( bDidPrimerPickingParametersChange( pOriginalPrimerPickingParameters, pCP, soParametersChanged ) ) soErrorMessage += soParametersChanged; TextBox* pTextBox = new TextBox( "Why No Acceptable Primers", 20 ); pTextBox->append( soErrorMessage ); pTextBox->makeVisible(); } else if ( nReasonExitedLoop == nExitedLoopBecauseCouldNotFindAnyUnacceptablePrimers ) { RWCString soWhyArePrimersUnacceptable; pCounts->getCountsInString( soWhyArePrimersUnacceptable ); RWCString soErrorMessage = soWhyArePrimersUnacceptable; TextBox* pTextBox = new TextBox( "Why No Acceptable Primers", 20 ); pTextBox->append( "Failed to find any primers, even unacceptable ones.\nCould it be that you have set the length of the primers too short or that the lengths are incompatible with the melting temperature?" ); pTextBox->append( soErrorMessage ); pTextBox->makeVisible(); } else assert( false ); // restore original primer picking parameters pCP->copyPrimerPickingParametersToSelf( pOriginalPrimerPickingParameters ); }