/*****************************************************************************
# 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<oligoTag>( 50,
"Assembly::>pChosenPrimersForAutofinishPCR_" );
}
pCP->setParametersForSequencingPrimersNotPCRPrimers( false );
pAssembly->figureOutContigOrderAndOrientation(); // this will set aRealContigs_
// find all contigEnd's that need multiplex pcr
RWTPtrOrderedVector<Contig> aContigEndsForMultiplexPCR;
RWTValOrderedVector<unsigned char> aContigEndsForMultiplexPCRWhichEnd;
// find all contigEnd pairs that need pairwise pcr
RWTPtrOrderedVector<contigEndPair> 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" );
}
}