/***************************************************************************** # 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 "autoFinishPCR.h" #include "subcloneTTemplate.h" #include "contigEndPair.h" #include "assembly.h" #include "rwtptrorderedvector.h" #include "rwtvalorderedvector.h" #include "contig.h" #include "pickCombinationsOfPCRPrimers.h" #include "consed.h" #include "consedParameters.h" #include "fileDefines.h" #include "contigEnd.h" #include "regionToAmplify.h" void autoFinishPCR() { Assembly* pAssembly = ConsEd::pGetAssembly(); if ( pCP->bWillDoExperiments_ && pCP->bAutoFinishTagOligosWhenDoExperiments_ ) { pAssembly->pChosenPrimersForAutofinishPCR_ = new RWTPtrOrderedVector( 50, "Assembly::>pChosenPrimersForAutofinishPCR_" ); } pCP->setParametersForSequencingPrimersNotPCRPrimers( false ); pAssembly->figureOutContigOrderAndOrientation(); // this will set aRealContigs_ // find all contigEnd's that need multiplex pcr RWTPtrOrderedVector aContigEndsForMultiplexPCR; RWTValOrderedVector aContigEndsForMultiplexPCRWhichEnd; // find all contigEnd pairs that need pairwise pcr RWTPtrOrderedVector aContigEndPairsNeedingPairwisePCR; int nContig; for( nContig = 0; nContig < pAssembly->aRealContigs_.length(); ++nContig ) { Contig* pContig = pAssembly->aRealContigs_[ nContig ]; pContig->determineIfTagsPreventExtendingContig(); pContig->determineIfTagsPreventPCR(); for( int nWhichEndI = 0; nWhichEndI <= 1; ++nWhichEndI ) { int nWhichEnd = ( nWhichEndI == 0 ? nLeftGap : nRightGap ); if ( nWhichEnd == nLeftGap && pContig->bTagsSayDoNotExtendToLeft_ ) { fprintf( pAO, "not considering pcr for left end of %s since tags say to not extend it\n", pContig->soGetName().data() ); continue; } if ( nWhichEnd == nRightGap && pContig->bTagsSayDoNotExtendToRight_ ) { fprintf( pAO, "not considering pcr for right end of %s since tags say to not extend it\n", pContig->soGetName().data() ); continue; } if ( nWhichEnd == nLeftGap && pContig->bTagsSayDoNotDoPCRWithLeftEnd_ ) { fprintf( pAO, "not considering pcr for left end of %s since tags say to not do pcr with it\n", pContig->soGetName().data() ); continue; } if ( nWhichEnd == nRightGap && pContig->bTagsSayDoNotDoPCRWithRightEnd_ ) { fprintf( pAO, "not considering pcr for right end of %s since tags say to not do pcr with it\n", pContig->soGetName().data() ); continue; } if ( nWhichEnd == nLeftGap && pContig->bReadsWillExtendLeftEnd_ && pCP->bAutoFinishDoNotDoPCRIfEndIsExtendedByReads_ ) { fprintf( pAO, "not considering pcr for left end of %s since autofinish already found walks and consed.autoFinishDoNotDoPCRIfEndIsExtendedByReads is true\n", pContig->soGetName().data() ); continue; } if ( nWhichEnd == nRightGap && pContig->bReadsWillExtendRightEnd_ && pCP->bAutoFinishDoNotDoPCRIfEndIsExtendedByReads_ ) { fprintf( pAO, "not considering pcr for right end of %s since autofinish already found walks and consed.autoFinishDoNotDoPCRIfEndIsExtendedByReads is true\n", pContig->soGetName().data() ); continue; } fprintf( pAO, "considering pcr for %s end of contig %s...\n", szWhichGap( nWhichEnd ), pContig->soGetName().data() ); if ( pContig->bIsThisEndTheEndOfTheClone( nWhichEnd ) ) continue; if ( !pContig->pContig_[ nWhichEnd ] ) { // there is not sufficient linking fwd/rev pair information // so this end must be part of the multiplex pcr fprintf( pAO, " not sufficient linkage for this contig-end\n" ); aContigEndsForMultiplexPCR.insert( pContig ); aContigEndsForMultiplexPCRWhichEnd.insert( nWhichEnd ); } else { // there is sufficient linking fwd/rev pair information. // Let's see how many of these linking templates are // available for finishing. contigEndPair* pCEP = pContig->pCEP_[ nWhichEnd ]; if ( !pCEP ) { // this can happen as follows: // -------- ----------- // A B (lots of links) // -------- -- // A C (a few links) // so C is inserted between A and B, but there // are no links between C and B. // In such a case, shall we do PCR? Sure. Contig* pOtherContig = pContig->pContig_[ nWhichEnd ]; int nEndOfOtherContig = pContig->nWhichEnd_[ nWhichEnd ]; fprintf( pAO, "contigEndPair between %s of %s and %s of %s has no templates, so doing PCR\n", szWhichGap( nWhichEnd ), pContig->soGetName().data(), szWhichGap( nEndOfOtherContig ), pOtherContig->soGetName().data() ); contigEndPair* pCEP = new contigEndPair( pContig, pOtherContig, nWhichEnd, nEndOfOtherContig ); if ( !aContigEndPairsNeedingPairwisePCR.bContains( pCEP ) ) { aContigEndPairsNeedingPairwisePCR.insert( pCEP ); } continue; } if ( !pCEP->pArrayOfTemplates_ ) { // don't think this ever happens... if ( !pCEP->bUserDefined_ ) { fprintf( pAO, "contigEndPair between %s and %s has no templates, but is not user defined. How can this be?\n" ); continue; } // Probably *no* linking templates. if ( aContigEndPairsNeedingPairwisePCR.index( pCEP ) == RW_NPOS ) aContigEndPairsNeedingPairwisePCR.insert( pCEP ); continue; } int nNumberOfAvailableTemplates = 0; for( int nSub = 0; nSub < pCEP->pArrayOfTemplates_->length(); ++nSub ) { subcloneTTemplate* pSub = (*(pCEP->pArrayOfTemplates_))[ nSub ]; fprintf( pAO, " linking template %s is available: %s\n", pSub->soGetName().data(), ( pSub->bOKToUseThisTemplate_ ? "yes" : "no" ) ); if ( pSub->bOKToUseThisTemplate_ ) ++nNumberOfAvailableTemplates; } if ( nNumberOfAvailableTemplates >= pCP->nAutoFinishDoNotDoPCRIfThisManyAvailableGapSpanningTemplates_ ) { fprintf( pAO, " %d available gap-spanning templates so not doing pcr with this contig-end\n", nNumberOfAvailableTemplates ); continue; } // do pairwise pcr if ( aContigEndPairsNeedingPairwisePCR.index( pCEP ) == RW_NPOS ) aContigEndPairsNeedingPairwisePCR.insert( pCEP ); } } // for( int nWhichEndI } // for( nContig = 0 // do pairwise pcr int nContigEndPair; for( nContigEndPair = 0; nContigEndPair < aContigEndPairsNeedingPairwisePCR.length(); ++nContigEndPair ) { contigEndPair* pCEP = aContigEndPairsNeedingPairwisePCR[ nContigEndPair ]; fprintf( pAO, "\nPCR FOR ORIENTED CONTIGS{\n" ); fprintf( pAO, "doing pairwise pcr for %s of %s and %s of %s\n", szWhichGap( pCEP->nWhichEnd_[0] ), pCEP->pContig_[0]->soGetName().data(), szWhichGap( pCEP->nWhichEnd_[1] ), pCEP->pContig_[1]->soGetName().data() ); if ( pCEP->bGapSizeSet_ && pCEP->nGapSize_ > pCP->nAutoFinishDoNotDoOrientedPCRIfGapSizeLargerThanThis_ ) { fprintf( pAO, "not doing PCR for this pair since gap size = %d\n", pCEP->nGapSize_ ); fprintf( pAO, "which is greater than maximum of %d\n", pCP->nAutoFinishDoNotDoOrientedPCRIfGapSizeLargerThanThis_ ); } else { regionToAmplify regionGap; regionGap.autoFinishPCRToCloseAGap( pCEP->pContig_[0], pCEP->nWhichEnd_[0], pCEP->pContig_[1], pCEP->nWhichEnd_[1] ); } fprintf( pAO, "}PCR FOR ORIENTED CONTIGS\n\n" ); fflush( pAO ); } if ( aContigEndsForMultiplexPCR.length() != 0 && pCP->bAutoFinishAllowUnorientedPCRReactions_ ) { fprintf( pAO, "\nPCR FOR UNORIENTED CONTIGS{\n" ); pickCombinationsOfPCRPrimers pickComb; pickComb.suggestPCRIfNecessaryForClosingGaps( &aContigEndsForMultiplexPCR, &aContigEndsForMultiplexPCRWhichEnd ); fprintf( pAO, "}PCR FOR UNORIENTED CONTIGS\n\n" ); } else { fprintf( pAO, "\nNot doing pcr for unoriented contigs because\n" ); if ( aContigEndsForMultiplexPCR.length() == 0 ) fprintf( pAO, "there are no unoriented contig ends\n\n" ); else if ( !pCP->bAutoFinishAllowUnorientedPCRReactions_ ) fprintf( pAO, "consed.autoFinishAllowUnorientedPCRReactions is false\n\n" ); } }