/*****************************************************************************
# 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 "contig.h"
#include "readTypes.h"
#include <math.h>
#include "consedParameters.h"
#include "numutil.h"
#include "bPrimerHasMatchesElsewhere.h"
#include "createPrimerCandidatesWith3PrimeEnd.h"
#include "bPrimerHasMatchesAgainstSequencesInFile.h"
#include "whichPrimersHaveAcceptableMeltingTemp.h"
#include "whichPrimersHaveAcceptableSelfMatch.h"
#include "pFindShortestAcceptablePrimer.h"
#include <stdlib.h>
#include "ucGetReadType.h"
#include "bIsTopStrandRead.h"
#include <iostream.h>
#include <iomanip.h>
#include "experiment.h"
#include "guiapp.h"
#include "consed.h"
#include "please_wait.h"
#include "paddedFromUnpadded.h"
#include "checkPrimersForMononucleotideRepeats.h"
#include "whyIsPrimerNotAcceptableTypes.h"
#include <stdio.h>
#include "mbtPtrVector.h"
#include "fileDefines.h"
#include "dErrorRateFromQuality.h"
#include "dGetErrorsFixedAtBase.h"
#include "nGetAmountOfOverlap.h"
#include "abs.h"
#include "nMAX_QUALITY_LEVEL2.h"
#include "multiplyDistributions.h"
#include "multiplyDistributions2.h"
#include "multiplyDistributions3.h"
#include "multiplyDistributionWithMinimum.h"
#include "makeDistributionAllOnes.h"
#include "assignDistributions.h"
#include "dumpDistribution.h"
#include "library.h"
#include "bigArrayOfLittleArraysOfDistributions.h"
#include "checkPrimersForACGT.h"
void Contig :: setUpForCreatingExperimentsList() {
// prepare for creating primers
// this is done already in autoFinish.cpp and in guiFindFinishingReads
// by setContigMatchTablesInAllContigs
// setContigMatchTables();
setTargetNumberOfErrorsForThisContig();
createFakeSubcloneTTemplatesForWholeCloneReads();
// now this is set by Assembly::setAllUnpaddedErrorProbablities
// setUnpaddedErrorProbabilities();
setExpectedNumberOfErrorsForThisContig();
setUnpaddedReadsThatCouldImproveQualities();
createCumulativeErrorArray();
if ( pCP->bAutoFinishSuggestSpecialChemistryForRunsAndStops_ ) {
determineLocationsOfRunsAndStops();
}
pCandidateExperiments_ = new mbtPtrOrderedVector<experiment>;
pChosenExperiments_ = new mbtPtrOrderedVector<experiment>;
pUniversalPrimerExperimentsForLowQualityRegions_ =
new mbtPtrOrderedVector<experiment>;
}
void Contig :: batchFindFinishingReads() {
fprintf( pAO, "\n\nCONTIG: %s ( with %d reads and length %d )\n",
(char*) soGetName().data(),
nGetNumberOfFragsInContig(),
nGetUnpaddedSeqLength() );
setUpForCreatingExperimentsList();
fprintf( pAO, "Estimated number of errors in consensus (not including doNotFinish tags): %.2f\n",
dExpectedNumberOfErrorsInConsensus_ );
fprintf( pAO, "Estimated number of errors in contig (includes gaps on ends): %.2f\n",
dExpectedNumberOfErrorsInExtendedContig_ );
fprintf( pAO, "Target number of errors for contig: %.2f\n\n",
dTargetNumberOfErrors_ );
if ( bIsThisEndTheEndOfTheClone( nLeftGap ) ) {
fprintf( pAO, "Left end of contig %s is one end of the clone\n",
soGetName().data() );
}
if ( bIsThisEndTheEndOfTheClone( nRightGap ) ) {
fprintf( pAO, "Right end of contig %s is one end of the clone\n",
soGetName().data() );
}
// this will examine fwd/rev pair information, and possibly
// suggest minilibraries or suggest reads to flank the gaps
dealWithContigEnds();
if ( pCP->bAutoFinishEmulate9_66Behavior_ ) {
if ( consedParameters::pGetConsedParameters()->bAutoFinishCoverLowConsensusQualityRegions_
||
consedParameters::pGetConsedParameters()->bAutoFinishCloseGaps_ ) {
createExperimentCandidatesListFor9_66();
chooseExperimentsToCoverLowQualityRegionsAndGapsFor9_66();
}
goto skipFor9_66Behavior;
}
if ( consedParameters::pGetConsedParameters()->bAutoFinishCoverLowConsensusQualityRegions_
||
consedParameters::pGetConsedParameters()->bAutoFinishCloseGaps_ ) {
saveStateOfConsensus();
if ( pCP->bAutoFinishNearGapsSuggestEachMissingReadOfReadPairs_ ) {
if ( pCP->bAutoFinishAllowDeNovoUniversalPrimerSubcloneReads_ ) {
fprintf( pAO,
"\n\nChoosing de novo universal primer reads to try to close gaps\n\n" );
nearGapsSuggestEachMissingReadOfReadPairs();
}
else {
fprintf( pAO,
"\n\nYou have specified:\nconsed.autoFinishAllowDeNovoUniversalPrimerSubcloneReads: false\nconsed.autoFinishNearGapsSuggestEachMissingReadOfReadPairs: true\nso autofinish will *not* suggest each missing read of read pairs\n" );
}
}
fprintf( pAO,
"\n\nChoosing universal primer reads to cover low quality regions and close gaps...\n\n" );
createUniversalPrimerExperimentCandidatesList();
chooseExperimentsToCoverLowQualityRegionsAndGaps( true );
// universal primer phase 1
if ( pCP->dAutoFinishRedundancy_ > 1.0 ) {
restoreStateOfConsensus();
recalculateCurrentExpectedNumberOfErrorsForThisContig(
dExpectedNumberOfErrorsInConsensus_,
dExpectedNumberOfErrorsInExtendedContig_ );
recalculateErrorsFixedOfExperimentCandidates();
fprintf( pAO, "\n\nChoosing redundant univeral primer reads to cover low quality regions and close gaps...\n\n" );
chooseExperimentsToCoverLowQualityRegionsAndGaps( true );
// universal primer phase 2
restoreStateOfConsensus();
recalculateConsensusUsingAllChosenUniversalPrimerExperiments();
recalculateCurrentExpectedNumberOfErrorsForThisContig(
dExpectedNumberOfErrorsInConsensus_,
dExpectedNumberOfErrorsInExtendedContig_ );
clearExperimentCandidatesList();
}
createCustomPrimerExperimentCandidatesList();
fprintf( pAO, "\n\nChoosing custom primer reads to cover low quality regions and close gaps...\n\n" );
chooseExperimentsToCoverLowQualityRegionsAndGaps( false );
// const bool bUniversalPrimerNotCustomPrimerReads )
if ( pCP->bAutoFinishCloseGaps_ )
warnUserIfCouldNotExtendContig();
}
skipFor9_66Behavior:
if ( consedParameters::pGetConsedParameters()->bAutoFinishCoverSingleSubcloneRegions_ ) {
fprintf( pAO, "\n\nChoosing custom primer reads to cover single subclone regions for contig %s\n\n",
(char*) soGetName().data() );
coverSingleSubcloneRegions();
}
fprintf( pAO, "Contig %s had %d suggested experiments.\n\n\n",
(char*) soGetName().data(),
pChosenExperiments_->length() );
}
void Contig :: chooseExperimentsToCoverLowQualityRegionsAndGaps(
const bool bUniversalPrimerNotCustomPrimerReads ) {
// Possible problem: if the consensus is high quality, but you want
// autofinish to extend the consensus, it might not do so.
// (DG, Feb 2000)
if ( dExpectedNumberOfErrorsInExtendedContig_ <= dTargetNumberOfErrors_ ) {
fprintf( pAO, "Estimated number of errors in extended contig %.2f is already less than the maximum %.2f for this contig\n",
dExpectedNumberOfErrorsInExtendedContig_,
dTargetNumberOfErrors_ );
return;
}
chooseExperimentsToCoverLowQualityRegionsAndGaps2(
bUniversalPrimerNotCustomPrimerReads );
}
void Contig :: warnUserIfCouldNotExtendContig() {
if ( !bIsThisEndTheEndOfTheClone( nLeftGap ) &&
bDoTagsAllowExtendingThisEnd( nLeftGap ) ) {
if ( nUnpaddedHighQualitySegmentStart_ ==
nUpdatedUnpaddedHighQualitySegmentStart_ )
fprintf( pAO, "WARNING: COULD NOT EXTEND CONTIG TO LEFT--MIGHT NEED PCR\n" );
else {
ConsEd::pGetAssembly()->bAnyReadsToCloseGaps_ = true;
bReadsWillExtendLeftEnd_ = true;
}
}
if ( !bIsThisEndTheEndOfTheClone( nRightGap ) &&
bDoTagsAllowExtendingThisEnd( nRightGap ) ) {
if ( nUnpaddedHighQualitySegmentEnd_ ==
nUpdatedUnpaddedHighQualitySegmentEnd_ )
fprintf( pAO, "WARNING: COULD NOT EXTEND CONTIG TO RIGHT--MIGHT NEED PCR\n" );
else {
ConsEd::pGetAssembly()->bAnyReadsToCloseGaps_ = true;
bReadsWillExtendRightEnd_ = true;
}
}
}
void Contig :: chooseExperimentsToCoverLowQualityRegionsAndGaps2(
const bool bUniversalPrimerNotCustomPrimerReads) {
if ( pCandidateExperiments_->length() == 0 ) {
fprintf( pAO, "No possible experiments found\n" );
return;
}
bool bMoreAvailableExperiments;
bool bConsiderMoreExperiments = true;
do {
experiment* pBestExp = pGetBestExperimentCandidate();
if ( pBestExp == NULL ) {
bMoreAvailableExperiments = false;
bConsiderMoreExperiments = false;
}
else {
// found the best experiment
assert( pBestExp->dErrorsFixed_ >= consedParameters::pGetConsedParameters()->dAutoFinishMinNumberOfErrorsFixedByAnExpToUse_ );
// the best experiment is fine--transfer it
pBestExp->nPurpose_ = nFixWeakRegion;
pBestExp->chooseExperiment();
possiblyUpdateContigHighQualitySegment( pBestExp );
pChosenExperiments_->insert( pBestExp );
if ( bUniversalPrimerNotCustomPrimerReads )
pUniversalPrimerExperimentsForLowQualityRegions_->insert( pBestExp );
ConsEd::pGetAssembly()->pAllChosenExperiments_->insert( pBestExp );
adjustExpectedNumberOfErrors( pBestExp );
findWindowOfMostErrorsFixed( pBestExp );
pBestExp->print();
fprintf( pAO, "After this experiment,\ncalculated number of errors in extended contig: %.2f and consensus: %.2f\n\n",
dExpectedNumberOfErrorsInExtendedContig_,
dExpectedNumberOfErrorsInConsensus_ );
if ( !bStillTooManyExpectedErrors() )
bConsiderMoreExperiments = false;
else {
adjustRemainingExperimentCandidates( pBestExp );
// check no bugs
double dErrorsInConsensus;
double dErrorsInExtendedContig;
recalculateCurrentExpectedNumberOfErrorsForThisContig(
dErrorsInConsensus,
dErrorsInExtendedContig );
if (
( ABS( dErrorsInConsensus - dExpectedNumberOfErrorsInConsensus_ ) > .000001 ) ||
( ABS( dErrorsInExtendedContig - dExpectedNumberOfErrorsInExtendedContig_ ) > .000001 )
) {
fprintf( pAO, "Warning: dExpectedNumberOfErrorsInConsensus = %e but calculated is %e and dExpectedNumberOfErrorsInExtendedContig_ = %e but calculated is %e\n",
dExpectedNumberOfErrorsInConsensus_,
dErrorsInConsensus,
dExpectedNumberOfErrorsInExtendedContig_,
dErrorsInExtendedContig );
}
else {
fprintf( pAO, "OK3\n" );
}
// end check no bugs
}
} // if ( pBestExp == NULL ) else
} while( bConsiderMoreExperiments );
// tell user why we stopped
// We could have stopped because:
// 1) the error rate was below the threshold
// 2) there were no more experiments
// 3) there were more experiments, but they weren't worth doing
if ( bStillTooManyExpectedErrors() ) {
fprintf( pAO, "There are still too many expected errors but we stopped because\n" );
if ( bMoreAvailableExperiments )
fprintf( pAO, "all remaining available experiments each solved too \nfew errors\nto be worthwhile doing\n" );
else
fprintf( pAO, "there were no more experiments this program could \nfind to resolve the remaining errors.\n" );
printWorstRemainingRegion();
}
else
fprintf( pAO, "The expected error count for the contig should be acceptable after doing these experiments\n" );
}
void Contig :: printWorstRemainingRegion() {
int nWindowOfWorstErrors = ConsEd::pGetAssembly()->pModelRead_->length() -
ConsEd::pGetAssembly()->nModelReadHighQualitySegmentStart_;
double dConsensusErrors = 0;
double dConsensusErrorsInWindow = 0;
double dMaxErrorsInWindow = -1;
int nUnpaddedMaxLocationLeft;
int nUnpaddedMaxLocationRight;
if ( pUnpaddedErrorProbabilities_->nGetStartIndex() + nWindowOfWorstErrors - 1 >
pUnpaddedErrorProbabilities_->nGetEndIndex() ) {
// this is the case in which the extended consensus is shorter than
// a single read
nUnpaddedMaxLocationLeft = nGetUnpaddedStartIndex();
nUnpaddedMaxLocationRight = nGetUnpaddedEndIndex();
dMaxErrorsInWindow = dExpectedNumberOfErrorsInExtendedContig_;
}
else {
// get the first window's worth
double dErrorsInCurrentWindow = 0.0;
for( int nUnpaddedConsPos = pUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpaddedConsPos <= pUnpaddedErrorProbabilities_->nGetStartIndex() +
nWindowOfWorstErrors - 1;
++nUnpaddedConsPos ) {
dErrorsInCurrentWindow +=
(*pUnpaddedErrorProbabilities_)[ nUnpaddedConsPos ];
}
dMaxErrorsInWindow = dErrorsInCurrentWindow;
nUnpaddedMaxLocationLeft =
pUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpaddedMaxLocationRight =
pUnpaddedErrorProbabilities_->nGetStartIndex() +
nWindowOfWorstErrors - 1;
for( int nUnpaddedConsPosRight =
pUnpaddedErrorProbabilities_->nGetStartIndex() + nWindowOfWorstErrors;
nUnpaddedConsPosRight <= pUnpaddedErrorProbabilities_->nGetEndIndex();
++nUnpaddedConsPosRight ) {
// remove the base just to the left of the leftmost and add the
// base at the right
dErrorsInCurrentWindow =
dErrorsInCurrentWindow
- (*pUnpaddedErrorProbabilities_)[ nUnpaddedConsPosRight - nWindowOfWorstErrors ]
+ (*pUnpaddedErrorProbabilities_)[ nUnpaddedConsPosRight ];
if ( dErrorsInCurrentWindow > dMaxErrorsInWindow ) {
dMaxErrorsInWindow = dErrorsInCurrentWindow;
nUnpaddedMaxLocationLeft =
nUnpaddedConsPosRight - nWindowOfWorstErrors + 1;
nUnpaddedMaxLocationRight = nUnpaddedConsPosRight;
}
}
} // if ( pUnpaddedErrorProbabilities_->nGetStartIndex() ... else
fprintf( pAO, "Worst %d base region that is left unfixed is from %d to %d with %.3f errors\n",
nWindowOfWorstErrors,
nUnpaddedMaxLocationLeft,
nUnpaddedMaxLocationRight,
dMaxErrorsInWindow );
}
void Contig :: createCumulativeErrorArray() {
int nArrayLength =
2* (
consedParameters::pGetConsedParameters()->nAutoFinishNumberOfBasesBetweenContigsAssumed_ ) +
nGetUnpaddedSeqLength();
int nOffset =
-consedParameters::pGetConsedParameters()->nAutoFinishNumberOfBasesBetweenContigsAssumed_
+ Contig :: nGetUnpaddedStartIndex();
pCumUnpaddedErrors_ =
new mbtValVector<double>( nArrayLength, nOffset, 0 );
(*pCumUnpaddedErrors_)[pCumUnpaddedErrors_->nGetStartIndex() ] =
(*pUnpaddedErrorProbabilities_)[ pCumUnpaddedErrors_->nGetStartIndex() ];
for( int nPos = pCumUnpaddedErrors_->nGetStartIndex() + 1;
nPos <= pCumUnpaddedErrors_->nGetEndIndex();
++nPos ) {
(*pCumUnpaddedErrors_)[ nPos ] = (*pUnpaddedErrorProbabilities_)[ nPos ]
+ (*pCumUnpaddedErrors_)[ nPos - 1 ];
}
}
void Contig :: resetCumulativeErrorArray() {
(*pCumUnpaddedErrors_)[pCumUnpaddedErrors_->nGetStartIndex() ] =
(*pUnpaddedErrorProbabilities_)[ pCumUnpaddedErrors_->nGetStartIndex() ];
for( int nPos = pCumUnpaddedErrors_->nGetStartIndex() + 1;
nPos <= pCumUnpaddedErrors_->nGetEndIndex();
++nPos ) {
(*pCumUnpaddedErrors_)[ nPos ] = (*pUnpaddedErrorProbabilities_)[ nPos ]
+ (*pCumUnpaddedErrors_)[ nPos - 1 ];
}
}
void Contig :: setUnpaddedErrorProbabilitiesAndMore() {
setUnpaddedErrorProbabilities();
// now look for tags that indicate not finishing. Whenever you find
// one, set the errors to zero there.
setErrorsToZeroWithinDoNotFinishTags();
handleTagsToNotExtendContig();
lookForCloneEndTags();
setHighQualitySegment();
// special case for people who are using autofinish just to close gaps.
// This is a little inefficient since it first sets all error probabilities
// (above) and then sets most of them to back to zero. It leaves
// the low quality quality errors near the ends of the consensus
// so that autofinish can improve its quality so it can extend
// into the gaps.
if ( !pCP->bAutoFinishCoverLowConsensusQualityRegions_ ) {
if ( nUnpaddedHighQualitySegmentStart_ != -1 &&
nUnpaddedHighQualitySegmentEnd_ != -1 ) {
for( int nUnpadded = nUnpaddedHighQualitySegmentStart_;
nUnpadded <= nUnpaddedHighQualitySegmentEnd_;
++nUnpadded ) {
(*pUnpaddedErrorProbabilities_)[ nUnpadded ] = 0.0;
}
}
}
// save a copy since the other copy will be updated as autofinish
// gets reads (this may change)
pOriginalUnpaddedErrorProbabilities_ = new mbtValVector<double>(
pUnpaddedErrorProbabilities_->length(),
pUnpaddedErrorProbabilities_->nOffset_,
0 );
for( int nPos = pUnpaddedErrorProbabilities_->nGetStartIndex();
nPos <= pUnpaddedErrorProbabilities_->nGetEndIndex();
++nPos ) {
(*pOriginalUnpaddedErrorProbabilities_)[ nPos ] =
(*pUnpaddedErrorProbabilities_)[ nPos ];
}
// set the array of arrays of distributions
// Intially, since there are no reads chosen, at each base, the
// array of distributions will be null.
pDistributionsOfChosenExperimentsArray_ =
new bigArrayOfLittleArraysOfDistributions(
pUnpaddedErrorProbabilities_->length(),
pUnpaddedErrorProbabilities_->nOffset_ );
pSavedDistributionsOfChosenExperimentsArray_ =
new bigArrayOfLittleArraysOfDistributions(
pUnpaddedErrorProbabilities_->length(),
pUnpaddedErrorProbabilities_->nOffset_ );
}
void Contig :: setUnpaddedErrorProbabilities() {
int nLength = nGetUnpaddedSeqLength() +
2 * consedParameters::pGetConsedParameters()-> nAutoFinishNumberOfBasesBetweenContigsAssumed_;
pUnpaddedErrorProbabilities_ = new mbtValVector<double>(
nLength,
-consedParameters::pGetConsedParameters()-> nAutoFinishNumberOfBasesBetweenContigsAssumed_
+ Contig::nGetUnpaddedStartIndex(),
0 );
pSavedUnpaddedErrorProbabilities_
= new mbtValVector<double>( *pUnpaddedErrorProbabilities_ );
// the Contig::nGetUnpaddedStartIndex above makes this array so that position 1 is the same
// as position 1 in the contig in unpadded consensus positions
// We count a gap base as an error *unless* it is off the end of the
// clone end. To figure this out, we need to use the clone end code.
double dErrorOfALeftGapBase;
if ( pCP->bAutoFinishCloseGaps_ ) {
if ( bIsThisEndTheEndOfTheClone( nLeftGap ) )
dErrorOfALeftGapBase = 0.0;
else
dErrorOfALeftGapBase = 1.0;
}
else
dErrorOfALeftGapBase = 0.0;
int nUnpaddedConsPos;
for( nUnpaddedConsPos = pUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpaddedConsPos <= 0;
++nUnpaddedConsPos )
(*pUnpaddedErrorProbabilities_)[ nUnpaddedConsPos ] =
dErrorOfALeftGapBase;
double dErrorOfARightGapBase;
if ( pCP->bAutoFinishCloseGaps_ ) {
if ( bIsThisEndTheEndOfTheClone( nRightGap ) )
dErrorOfARightGapBase = 0.0;
else
dErrorOfARightGapBase = 1.0;
}
else
dErrorOfARightGapBase = 0.0;
// found and fixed bug via desk check
for( nUnpaddedConsPos = nGetUnpaddedEndIndex() + 1;
nUnpaddedConsPos <= pUnpaddedErrorProbabilities_->nGetEndIndex();
++nUnpaddedConsPos )
(*pUnpaddedErrorProbabilities_)[ nUnpaddedConsPos ] =
dErrorOfARightGapBase;
int nUnpaddedPos = Contig::nGetUnpaddedStartIndex();
Ntide nt;
int nQuality;
double dError;
for( int nConsPos = nGetStartConsensusIndex();
nConsPos <= nGetEndConsensusIndex();
++nConsPos ) {
nt = ntGetCons( nConsPos );
if ( !nt.bIsPad() ) {
// changed 3/17/99 since 98 should not be considered a high quality
// dErrorRateFromQuality handles this correctly--see
// setErrorRateFromQuality.cpp
nQuality = nt.qualGetQuality();
(*pUnpaddedErrorProbabilities_)[ nUnpaddedPos ] =
dErrorRateFromQuality[ nQuality ];
++nUnpaddedPos;
}
}
assert( nUnpaddedPos == ( nGetUnpaddedEndIndex() + 1 ) );
}
void Contig :: handleTagsToNotExtendContig() {
determineIfTagsPreventExtendingContig();
if ( bTagsSayDoNotExtendToLeft_ ) {
for( int nUnpaddedConsPos = pUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpaddedConsPos <= Contig::nGetUnpaddedStartIndex() - 1;
++nUnpaddedConsPos )
(*pUnpaddedErrorProbabilities_)[ nUnpaddedConsPos ] = 0.0;
}
if ( bTagsSayDoNotExtendToRight_ ) {
for( int nUnpaddedConsPos = Contig::nGetUnpaddedEndIndex() + 1;
nUnpaddedConsPos <= pUnpaddedErrorProbabilities_->nGetEndIndex();
++nUnpaddedConsPos )
(*pUnpaddedErrorProbabilities_)[ nUnpaddedConsPos ] = 0.0;
}
}
RWCString soDoNotDoPCR( "doNotDoPCR" );
void Contig :: determineIfTagsPreventPCR() {
if ( bAlreadyExaminedContigForDoNotDoPCRTags_ )
return;
bAlreadyExaminedContigForDoNotDoPCRTags_ = true;
fprintf( pAO, "looking for tags indicating do not do PCR in contig %s\n",
soGetName().data() );
// first look in the reads
for( int nRead = 0; nRead < nGetNumberOfFragsInContig(); ++nRead ) {
LocatedFragment* pLocFrag = pLocatedFragmentGet( nRead );
for( int nTag = 0; nTag < pLocFrag->nGetNumberOfTags(); ++nTag ) {
tag* pTag = pLocFrag->pGetTag( nTag );
if ( pTag->soGetTagType() != soDoNotDoPCR )
continue;
int nUnpaddedConsPosRangeStart =
nUnpaddedIndex( pTag->nPaddedConsPosStart_ ) -
pCP->nAutoFinishDoNotExtendContigsIfTagsAreThisCloseToContigEnd_;
int nUnpaddedConsPosRangeEnd =
nUnpaddedIndex( pTag->nPaddedConsPosEnd_ ) +
pCP->nAutoFinishDoNotExtendContigsIfTagsAreThisCloseToContigEnd_;
bool bOverlapsLeftEnd = false;
bool bOverlapsRightEnd = false;
if ( nUnpaddedConsPosRangeStart <= nGetUnpaddedStartIndex() &&
nGetUnpaddedStartIndex() <= nUnpaddedConsPosRangeEnd ) {
bOverlapsLeftEnd = true;
}
if ( nUnpaddedConsPosRangeStart <= nGetUnpaddedEndIndex() &&
nGetUnpaddedEndIndex() <= nUnpaddedConsPosRangeEnd ) {
bOverlapsRightEnd = true;
}
if ( bOverlapsLeftEnd && bOverlapsRightEnd ) {
/* pathological case in which this is a small contig and
we are having trouble figuring out which end the tag
refers to */
int nUnpaddedConsPosStart =
nUnpaddedIndex( pTag->nPaddedConsPosStart_ );
int nUnpaddedConsPosEnd =
nUnpaddedIndex( pTag->nPaddedConsPosEnd_ );
if ( ABS( nUnpaddedConsPosStart - nGetUnpaddedStartIndex() )
<
ABS( nUnpaddedConsPosEnd - nGetUnpaddedEndIndex() ) )
bTagsSayDoNotDoPCRWithLeftEnd_ = true;
else
bTagsSayDoNotDoPCRWithRightEnd_ = true;
}
else {
if ( bOverlapsLeftEnd )
bTagsSayDoNotDoPCRWithLeftEnd_ = true;
else
bTagsSayDoNotDoPCRWithRightEnd_ = true;
}
} // for( int nTag = ...
} // (for( int nRead = ...
// now look at consensus tags
for( int nConsensusTag = 0; nConsensusTag < nGetNumberOfTags();
++nConsensusTag ) {
tag* pTag = pGetTag( nConsensusTag );
if ( pTag->soGetTagType() != soDoNotDoPCR )
continue;
int nUnpaddedConsPosRangeStart =
nUnpaddedIndex( pTag->nPaddedConsPosStart_ )
-
pCP->nAutoFinishDoNotExtendContigsIfTagsAreThisCloseToContigEnd_;
int nUnpaddedConsPosRangeEnd =
nUnpaddedIndex( pTag->nPaddedConsPosEnd_ )
+
pCP->nAutoFinishDoNotExtendContigsIfTagsAreThisCloseToContigEnd_;
bool bOverlapsLeftEnd = false;
bool bOverlapsRightEnd = false;
if ( nUnpaddedConsPosRangeStart <= nGetUnpaddedStartIndex() &&
nGetUnpaddedStartIndex() <= nUnpaddedConsPosRangeEnd )
bOverlapsLeftEnd = true;
if ( nUnpaddedConsPosRangeStart <= nGetUnpaddedEndIndex() &&
nGetUnpaddedEndIndex() <= nUnpaddedConsPosRangeEnd )
bOverlapsRightEnd = true;
if ( bOverlapsLeftEnd && bOverlapsRightEnd ) {
/* pathological case in which this is a small contig and we
are having trouble figuring out which end the tag refers
to */
int nUnpaddedConsPosStart =
nUnpaddedIndex( pTag->nPaddedConsPosStart_ );
int nUnpaddedConsPosEnd =
nUnpaddedIndex( pTag->nPaddedConsPosEnd_ );
if ( ABS( nUnpaddedConsPosStart - nGetUnpaddedStartIndex() )
<
ABS( nUnpaddedConsPosEnd - nGetUnpaddedEndIndex() ) )
bTagsSayDoNotDoPCRWithLeftEnd_ = true;
else
bTagsSayDoNotDoPCRWithRightEnd_ = true;
}
else if ( bOverlapsLeftEnd )
bTagsSayDoNotDoPCRWithLeftEnd_ = true;
else if ( bOverlapsRightEnd )
bTagsSayDoNotDoPCRWithRightEnd_ = true;
}
if ( bTagsSayDoNotDoPCRWithLeftEnd_ )
fprintf( pAO, "tags say do not do pcr with left end\n" );
if ( bTagsSayDoNotDoPCRWithRightEnd_ )
fprintf( pAO, "tags say do not do pcr with right end\n" );
}
void Contig :: determineIfTagsPreventExtendingContig() {
if ( bAlreadyExaminedContigForDoNotExtendTags_ )
return;
bAlreadyExaminedContigForDoNotExtendTags_ = true;
fprintf( pAO, "looking for tags indicating do not extend contig in contig %s\n",
soGetName().data() );
// first look in the reads
for( int nRead = 0; nRead < nGetNumberOfFragsInContig(); ++nRead ) {
LocatedFragment* pLocFrag = pLocatedFragmentGet( nRead );
for( int nTag = 0; nTag < pLocFrag->nGetNumberOfTags(); ++nTag ) {
tag* pTag = pLocFrag->pGetTag( nTag );
if ( pCP->aAutoFinishTagsToNotExtend_.index( pTag->soGetTagType() ) == RW_NPOS ) continue;
int nUnpaddedConsPosStart =
nUnpaddedIndex( pTag->nPaddedConsPosStart_ ) -
pCP->nAutoFinishDoNotExtendContigsIfTagsAreThisCloseToContigEnd_;
int nUnpaddedConsPosEnd =
nUnpaddedIndex( pTag->nPaddedConsPosEnd_ ) +
pCP->nAutoFinishDoNotExtendContigsIfTagsAreThisCloseToContigEnd_;
if ( nUnpaddedConsPosStart <= nGetUnpaddedStartIndex() &&
nGetUnpaddedStartIndex() <= nUnpaddedConsPosEnd ) {
bTagsSayDoNotExtendToLeft_ = true;
fprintf( pAO, "do not extend contig %s to left because of read tag %s on read %s\n",
soGetName().data(),
pTag->soGetTagType().data(),
pLocFrag->soGetName().data() );
}
if ( nUnpaddedConsPosStart <= nGetUnpaddedEndIndex() &&
nGetUnpaddedEndIndex() <= nUnpaddedConsPosEnd ) {
bTagsSayDoNotExtendToRight_ = true;
fprintf( pAO, "do not extend contig %s to right because of read tag %s on read %s\n",
soGetName().data(),
pTag->soGetTagType().data(),
pLocFrag->soGetName().data() );
}
}
}
// now look in consensus tags
for( int nConsensusTag = 0; nConsensusTag < nGetNumberOfTags();
++nConsensusTag ) {
tag* pTag = pGetTag( nConsensusTag );
if ( pCP->aAutoFinishTagsToNotExtend_.index( pTag->soGetTagType() ) ==
RW_NPOS ) continue;
int nUnpaddedConsPosStart = nUnpaddedIndex( pTag->nPaddedConsPosStart_ )
- pCP->nAutoFinishDoNotExtendContigsIfTagsAreThisCloseToContigEnd_;
int nUnpaddedConsPosEnd = nUnpaddedIndex( pTag->nPaddedConsPosEnd_ )
+ pCP->nAutoFinishDoNotExtendContigsIfTagsAreThisCloseToContigEnd_;
if ( nUnpaddedConsPosStart <= nGetUnpaddedStartIndex() &&
nGetUnpaddedStartIndex() <= nUnpaddedConsPosEnd ) {
bTagsSayDoNotExtendToLeft_ = true;
fprintf( pAO, "do not extend contig %s to left because of consensus tag %s\n",
soGetName().data(),
pTag->soGetTagType().data() );
}
if ( nUnpaddedConsPosStart <= nGetUnpaddedEndIndex() &&
nGetUnpaddedEndIndex() <= nUnpaddedConsPosEnd ) {
fprintf( pAO, "do not extend contig %s to right because of consensus tag %s\n",
soGetName().data(),
pTag->soGetTagType().data() );
bTagsSayDoNotExtendToRight_ = true;
}
}
}
void Contig :: lookForCloneEndTags() {
fprintf( pAO, "looking for cloneEnd tags in contig %s\n",
soGetName().data() );
// first look in the reads
for( int nRead = 0; nRead < nGetNumberOfFragsInContig(); ++nRead ) {
LocatedFragment* pLocFrag = pLocatedFragmentGet( nRead );
for( int nTag = 0; nTag < pLocFrag->nGetNumberOfTags(); ++nTag ) {
tag* pTag = pLocFrag->pGetTag( nTag );
if ( pTag->soGetTagType() != "cloneEnd" )
continue;
// found a cloneEnd tag!
dealWithCloneEndTag( pTag );
}
}
// now look in consensus tags
for( int nConsensusTag = 0; nConsensusTag < nGetNumberOfTags();
++nConsensusTag ) {
tag* pTag = pGetTag( nConsensusTag );
if ( pTag->soGetTagType() != "cloneEnd" )
continue;
// found a cloneEnd tag!
dealWithCloneEndTag( pTag );
}
}
void Contig :: dealWithCloneEndTag( tag* pTag ) {
fprintf( pAO, "found clone end tag at %d to %d with insert at %s\n",
nUnpaddedIndex( pTag->nPaddedConsPosStart_ ),
nUnpaddedIndex( pTag->nPaddedConsPosEnd_ ),
pTag->soMiscData_.data() );
int nUnpaddedLeft;
int nUnpaddedRight;
if ( pTag->soMiscData_ == "->" ) {
nUnpaddedLeft = pUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpaddedRight = nUnpaddedIndex( pTag->nPaddedConsPosEnd_ );
bLeftEndOfContigExaminedForEndOfClone_ = true;
bLeftEndOfContigIsEndOfClone_ = true;
}
else if ( pTag->soMiscData_ == "<-" ) {
nUnpaddedLeft = nUnpaddedIndex( pTag->nPaddedConsPosStart_ );
nUnpaddedRight = pUnpaddedErrorProbabilities_->nGetEndIndex();
bRightEndOfContigExaminedForEndOfClone_ = true;
bRightEndOfContigIsEndOfClone_ = true;
}
else {
// something screwed up with this tag
return;
}
// set the error probability to zero within the clone vector
// so it isn't finished.
for( int nUnpadded = nUnpaddedLeft; nUnpadded <= nUnpaddedRight;
++nUnpadded ) {
(*pUnpaddedErrorProbabilities_)[ nUnpadded ] = 0.0;
}
}
void Contig :: setErrorsToZeroWithinDoNotFinishTags() {
// look for do-not-finish tags in both reads and in the consensus
// First look in the reads.
for( int nRead = 0; nRead < nGetNumberOfFragsInContig(); ++nRead ) {
LocatedFragment* pLocFrag = pLocatedFragmentGet( nRead );
for( int nTag = 0; nTag < pLocFrag->nGetNumberOfTags(); ++nTag ) {
tag* pTag = pLocFrag->pGetTag( nTag );
if ( consedParameters::pGetConsedParameters()->aAutoFinishTagsToNotFinish_.index( pTag->soGetTagType() ) == RW_NPOS ) continue;
int nUnpaddedConsPosStart = nUnpaddedIndex( pTag->nPaddedConsPosStart_ );
int nUnpaddedConsPosEnd = nUnpaddedIndex( pTag->nPaddedConsPosEnd_ );
// read tags might be on a read that protrudes far beyond the
// end of the consensus, and thus isn't on the extended contig
int nOverlapLeft;
int nOverlapRight;
if ( bIntersect( nUnpaddedConsPosStart,
nUnpaddedConsPosEnd,
pUnpaddedErrorProbabilities_->nGetStartIndex(),
pUnpaddedErrorProbabilities_->nGetEndIndex(),
nOverlapLeft,
nOverlapRight ) ) {
for( int nUnpadded = nOverlapLeft;
nUnpadded <= nOverlapRight;
++nUnpadded ) {
(*pUnpaddedErrorProbabilities_)[ nUnpadded ] = 0.0;
}
}
} // for( nTag = 0 ...
} // for( int nRead = ...
// now look for consensus tags
for( int nConsensusTag = 0; nConsensusTag < nGetNumberOfTags();
++nConsensusTag ) {
tag* pTag = pGetTag( nConsensusTag );
if ( consedParameters::pGetConsedParameters()->aAutoFinishTagsToNotFinish_.index( pTag->soGetTagType() ) == RW_NPOS ) continue;
int nUnpaddedConsPosStart = nUnpaddedIndex( pTag->nPaddedConsPosStart_ );
int nUnpaddedConsPosEnd = nUnpaddedIndex( pTag->nPaddedConsPosEnd_ );
for( int nUnpadded = nUnpaddedConsPosStart;
nUnpadded <= nUnpaddedConsPosEnd;
++nUnpadded ) {
(*pUnpaddedErrorProbabilities_)[ nUnpadded ] = 0.0;
}
} // for( int nConsensusTag = 0 ...
}
void Contig :: setUnpaddedReadsThatCouldImproveQualities() {
int nArrayLength =
2* (
consedParameters::pGetConsedParameters()->nAutoFinishNumberOfBasesBetweenContigsAssumed_ ) +
nGetUnpaddedSeqLength();
int nOffset =
-consedParameters::pGetConsedParameters()->nAutoFinishNumberOfBasesBetweenContigsAssumed_
+ nGetUnpaddedStartIndex();
pUnpaddedReadTypesCoveringEachBase_ = new mbtValVector<unsigned char>(
nArrayLength,
nOffset,
'\0' );
pOriginalUnpaddedReadTypesCoveringEachBase_ = new mbtValVector<unsigned char>(
nArrayLength,
nOffset );
pSavedUnpaddedReadTypesCoveringEachBase_
= new mbtValVector<unsigned char>( nArrayLength, nOffset );
unsigned char ucZero = 0;
for( int nPos = pUnpaddedReadTypesCoveringEachBase_->nGetStartIndex();
nPos <= pUnpaddedReadTypesCoveringEachBase_->nGetEndIndex();
++nPos ) {
(*pUnpaddedReadTypesCoveringEachBase_)[ nPos ] = ucZero;
}
for( int nRead = 0; nRead < nGetNumberOfFragsInContig();
++nRead ) {
LocatedFragment* pLocFrag = pLocatedFragmentGet( nRead );
// such a read doesn't belong where it is--do not
// consider it as covering the region
if ( pLocFrag->bWholeReadIsUnaligned_ ) continue;
int nUnpaddedLeft;
int nUnpaddedRight;
// avoid subscript errors
if ( bIntersect(
nUnpaddedIndex( pLocFrag->nAlignClipStart_ ),
nUnpaddedIndex( pLocFrag->nAlignClipEnd_ ),
pUnpaddedReadTypesCoveringEachBase_->nGetStartIndex(),
pUnpaddedReadTypesCoveringEachBase_->nGetEndIndex(),
nUnpaddedLeft,
nUnpaddedRight ) ) {
unsigned char ucStrandAndChemistry = pLocFrag->ucGetReadType();
for( int nUnpaddedPos = nUnpaddedLeft;
nUnpaddedPos <= nUnpaddedRight; ++nUnpaddedPos ) {
(*pUnpaddedReadTypesCoveringEachBase_)[ nUnpaddedPos ] |=
ucStrandAndChemistry;
}
}
}
// copy to pOriginalUnpaddedReadTypes_
for( int nUnpadded = pUnpaddedReadTypesCoveringEachBase_->nGetStartIndex();
nUnpadded <= pUnpaddedReadTypesCoveringEachBase_->nGetEndIndex();
++nUnpadded ) {
(*pOriginalUnpaddedReadTypesCoveringEachBase_)[ nUnpadded ] =
(*pUnpaddedReadTypesCoveringEachBase_)[ nUnpadded ];
}
}
void Contig :: createUniversalPrimerExperimentCandidatesList() {
if ( pCP->bAutoFinishAllowResequencingReads_ ) {
fprintf( pAO, "creating list of possible experiments with universal primers..." );
fflush( pAO );
createExperimentsWithResequencingUniversalPrimers();
fprintf( pAO, "done\n" );
fflush( pAO );
}
if ( pCP->bAutoFinishAllowDeNovoUniversalPrimerSubcloneReads_ ) {
fprintf( pAO, "creating list of possible experiments of reverses where there is only a forward or forwards where there is only a reverse..." );
fflush( pAO );
createExperimentsWithFirstForwardsOrFirstReverses();
fprintf( pAO, "done\n\n" );
fflush( pAO );
}
}
void Contig :: createCustomPrimerExperimentCandidatesList() {
if ( pCP->bAutoFinishAllowCustomPrimerSubcloneReads_ ||
pCP->bAutoFinishAllowWholeCloneReads_ ) {
fprintf( pAO, "creating list of possible experiments with custom primers..." );
fflush( pAO );
createExperimentsWithCustomPrimers();
fprintf( pAO, "done\n\n" );
fflush( pAO );
}
}
void Contig :: createExperimentsWithCustomPrimers() {
// now create experiments with custom primers
pCP->setParametersForSequencingPrimersNotPCRPrimers( true );
// create an experiment if it will solve more than
// nAutoFinishMinimumErrorThresholdOfRead
for( int nReadType = 0; nReadType < 2; ++nReadType ) {
readTypes ucRead;
// for now, just consider terminator reads
if (nReadType == 0 )
ucRead = TOP_STRAND_TERMINATOR_READS;
else if ( nReadType == 1)
ucRead = BOTTOM_STRAND_TERMINATOR_READS;
else
assert( false );
bool bTopStrandNotBottomStrand = bIsTopStrandRead( ucRead );
int nFirstExperimentHighQualityLeftEnd;
int nLastExperimentHighQualityLeftEnd;
// g=gap C=consensus base
// ggggggCCCCCCCCCCCCCCCCCCCCCCCCCCCgggggg
// top strand reads from leftmost to rightmost
// ppppHHHHHHllll ppppppHHHHHHHllllll
// bottom strand reads
// lllllHHHHpppp llllHHHHHpppp
// That is, we try to make reads wherever there is
// room on the consensus for the primer.
// Primers can be of varying lengths. Shall we save enough
// room for the longest?
int nFirstExperimentFirstBaseOfReadUnpadded;
int nLastExperimentFirstBaseOfReadUnpadded;
if ( bTopStrandNotBottomStrand ) {
// we must be able to fit a primer onto the consensus.
// ggggggCCCCCCCCCCCCCCCCCCCCCCCCCCCgggggg
// ppppHHHHHHllll ppppppHHHHHHHllllll
nFirstExperimentFirstBaseOfReadUnpadded = nGetUnpaddedStartIndex() +
pCP->nPrimersMinimumLengthOfAPrimerToUse_;
nLastExperimentFirstBaseOfReadUnpadded = nGetUnpaddedEndIndex() -
pCP->nPrimersMinimumLengthOfAPrimerToUse_ + 1;
}
else {
// ggggggCCCCCCCCCCCCCCCCCCCCCCCCCCCgggggg
// lllllHHHHpppp llllHHHHHpppp
nFirstExperimentFirstBaseOfReadUnpadded = nGetUnpaddedStartIndex() - 1;
nLastExperimentFirstBaseOfReadUnpadded = nGetUnpaddedEndIndex() -
pCP->nPrimersMinimumLengthOfAPrimerToUse_;
}
double dMinQuality =
consedParameters::pGetConsedParameters()->dAutoFinishMinNumberOfErrorsFixedByAnExpToUse_;
// all unpadded positions
// nFirstBaseOfReadUnpadded is the most 5' base of the read, after the primer
// Thus it is the left end of the read for top strand reads,
// and the right end of the read for bottom strand reads.
for( int nFirstBaseOfReadUnpadded = nFirstExperimentFirstBaseOfReadUnpadded;
nFirstBaseOfReadUnpadded <= nLastExperimentFirstBaseOfReadUnpadded;
++nFirstBaseOfReadUnpadded ) {
// see if this experiment will solve more than the minimum # of errors
// I will do that by asking a simpler question: "is there more than
// the minimum number of errors in the entire region covered by
// the potential read?" If not, skip this experiment.
double dMaxErrors = dGetInitialMaxErrorsFixed2(
nFirstBaseOfReadUnpadded,
bTopStrandNotBottomStrand );
if ( dMaxErrors < pCP->dAutoFinishMinNumberOfErrorsFixedByAnExp_ )
continue;
// However, the quality might be too
// low to create a primer. If a primer can't be created,
// there is no point in considering this experiment further.
// So check that next:
if ( !bCustomPrimerQualityTest( nFirstBaseOfReadUnpadded,
bTopStrandNotBottomStrand ) ) continue;
if ( pCP->bAutoFinishAllowWholeCloneReads_ ) {
maybeCreateExperimentWithCustomPrimer(
ucRead,
bTopStrandNotBottomStrand,
nFirstBaseOfReadUnpadded,
true, // bCloneTemplateNotSubcloneTemplate
false ); // bForCoveringSingleSubcloneRegions
}
if ( pCP->bAutoFinishAllowCustomPrimerSubcloneReads_ ) {
maybeCreateExperimentWithCustomPrimer(
ucRead,
bTopStrandNotBottomStrand,
nFirstBaseOfReadUnpadded,
false, // bCloneTemplateNotSubcloneTemplate
false ); // bForCoveringSingleSubcloneRegions
}
} // for( int nFirstBaseOfReadUnpadded = ...
} // for (int nReadType...
// no longer sorting this list--just pick out top one each time.
// Much more efficient.
}
bool Contig :: bCustomPrimerQualityTest( const int nFirstBaseOfReadUnpadded,
const bool bTopStrandNotBottomStrand ) {
// nFirstBaseOfReadUnpadded is in 1based unpadded coordinates, but primer
// arrays in in 0 based coordinates so convert
int nZeroBasedUnpaddedLeftMost;
int nZeroBasedUnpaddedRightMost;
if ( bTopStrandNotBottomStrand ) {
// subtract 1 to convert to 0 based coordinates.
// then subtract 1 again to move from read to primer
nZeroBasedUnpaddedRightMost = nFirstBaseOfReadUnpadded - 2;
nZeroBasedUnpaddedLeftMost = nZeroBasedUnpaddedRightMost -
consedParameters::pGetConsedParameters()->nPrimersMinimumLengthOfAPrimerToUse_
+ 1;
}
else {
// substract 1 to convert to 0-based coordinates
// then add 1 to move from read to primer
nZeroBasedUnpaddedLeftMost = nFirstBaseOfReadUnpadded;
nZeroBasedUnpaddedRightMost = nZeroBasedUnpaddedLeftMost +
consedParameters::pGetConsedParameters()->nPrimersMinimumLengthOfAPrimerToUse_
- 1;
}
int n1BasedLeft = nZeroBasedUnpaddedLeftMost + 1;
int n1BasedRight = nZeroBasedUnpaddedRightMost + 1;
// if the primer is off of the consensus, then can't make a custom oligo
// Note that this only checks that the minimum length primer is on the
// consensus. It is possible that a longer length primer will be
// off the consensus. We will need to check for that when we create
// the primers.
if ( n1BasedLeft < nGetUnpaddedStartIndex() ) {
return( false );
}
if ( n1BasedRight > nGetUnpaddedEndIndex() ) {
return( false );
}
unsigned char ucMinAcceptableQuality =
consedParameters::pGetConsedParameters()->nPrimersMinQuality_;
for( int nPrimerPos = nZeroBasedUnpaddedLeftMost;
nPrimerPos <= nZeroBasedUnpaddedRightMost;
++nPrimerPos ) {
if ( pUnpaddedContig_->pUnpaddedQualityArray_[ nPrimerPos ] <
ucMinAcceptableQuality )
return( false );
}
return( true );
} // bCustomPrimerQualityTest
void Contig :: maybeCreateExperimentWithCustomPrimer(
const readTypes ucRead,
const bool bTopStrandNotBottomStrand,
const int nFirstBaseOfReadUnpadded,
// 1 -based coordinates
const bool bCloneTemplateNotSubcloneTemplate,
const bool bForCoveringSingleSubcloneRegions
) {
int nUnpaddedPosOf3PrimeEnd;
if ( bTopStrandNotBottomStrand ) {
// subtract 1 to move from read to primer
// (both are in 1-based coordinates)
nUnpaddedPosOf3PrimeEnd = nFirstBaseOfReadUnpadded - 1;
}
else {
// add 1 to move from read to primer
// (both are in 1-based coordinates)
nUnpaddedPosOf3PrimeEnd = nFirstBaseOfReadUnpadded + 1;
}
primerType* pPrimerArray;
int nNumberOfPrimers;
bool bSomeAcceptablePrimers = true;
createPrimerCandidatesWith3PrimeEnd( this,
nUnpaddedPosOf3PrimeEnd,
bTopStrandNotBottomStrand,
pPrimerArray,
nNumberOfPrimers,
bSomeAcceptablePrimers );
if (!bSomeAcceptablePrimers ) return;
// only check the shortest--if it fails, they all will fail
primerType* pShortestPrimer =
pFindShortestAcceptablePrimer( pPrimerArray, nNumberOfPrimers );
// the following must be true or else bPrimerHasMatchesElsewhere won't
// do any checking
if ( pShortestPrimer->nUnpaddedLength_ !=
consedParameters::pGetConsedParameters()->nPrimersMinimumLengthOfAPrimerToUse_
) {
cout << "primer strange length = " << pShortestPrimer->nUnpaddedLength_
<< endl;
cout.flush();
assert( false );
}
bool bFoundPrimer = false;
for( int nPrimer = 0; nPrimer < nNumberOfPrimers; ++nPrimer ) {
primerType* pPrimer = pPrimerArray + nPrimer;
if ( ConsEd::pGetAssembly()->bFoundAutoFinishDebugPrimer(
pPrimer ) ) {
bFoundPrimer = true;
break;
}
}
if ( bFoundPrimer ) {
cerr << "found primer before bPrimerHasMatchesElsewhere, shortest version with same 3' end is ";
for( int n = 0; n < pShortestPrimer->nUnpaddedLength_; ++n ) {
cerr << pShortestPrimer->szPrimer_[ n ];
}
cerr << endl;
}
if ( bPrimerHasMatchesElsewhere( pShortestPrimer,
bTopStrandNotBottomStrand,
bCloneTemplateNotSubcloneTemplate
) )
return;
// only check the shortest--if it fails, they all will fail
if ( bPrimerHasMatchesAgainstSequencesInFile( pShortestPrimer,
bCloneTemplateNotSubcloneTemplate) )
return;
whichPrimersHaveAcceptableMeltingTemp( pPrimerArray, nNumberOfPrimers,
bSomeAcceptablePrimers );
if ( !bSomeAcceptablePrimers ) return;
whichPrimersHaveAcceptableSelfMatch( pPrimerArray, nNumberOfPrimers,
bSomeAcceptablePrimers );
if ( !bSomeAcceptablePrimers ) return;
int nNumberOfAcceptablePrimers;
checkPrimersForMononucleotideRepeats( pPrimerArray, nNumberOfPrimers,
nNumberOfAcceptablePrimers );
if ( nNumberOfAcceptablePrimers == 0 ) return;
checkPrimersForACGT( pPrimerArray, nNumberOfPrimers,
nNumberOfAcceptablePrimers );
if ( nNumberOfAcceptablePrimers == 0 ) return;
// if reached here, there is at least 1 acceptable primer for this
// position. There may be more, in which case I think that we
// should just take the shortest.
primerType* pPrimer = pFindShortestAcceptablePrimer( pPrimerArray,
nNumberOfPrimers );
// if this is using a subclone template, we need to be sure
// that there is some acceptable template
if ( !bCloneTemplateNotSubcloneTemplate ) {
// isn't pPrimer->pContig_ always this? (DG, Dec 1998)
if ( ! pPrimer->pContig_->bFindTemplatesForPrimer(
pPrimer,
ucRead
) ) {
return;
}
}
if ( bFoundPrimer ) {
cerr << "found primer after all checks except errors fixed ";
for( int n = 0; n < pShortestPrimer->nUnpaddedLength_; ++n ) {
cerr << pShortestPrimer->szPrimer_[ n ];
}
cerr << endl;
}
// before we just checked roughly if the read would fix enough errors.
// Now let's check precisely
// this needs to take into account the template lengths for subclone
// reads. For whole clone reads, we will assume that the read will
// be the whole length of the model read. In the future, we may
// downgrade this quality.
// find the read in unpadded consensus coordinates
double dErrorsFixed;
if ( bForCoveringSingleSubcloneRegions ) {
dErrorsFixed = -2.0; // no purpose served by calculating this.
// since it isn't used.
}
else {
if ( bCloneTemplateNotSubcloneTemplate ) {
dErrorsFixed = dGetErrorsFixedForWholeCloneRead(
bTopStrandNotBottomStrand,
nFirstBaseOfReadUnpadded,
false, // bChooseExperiment
false, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
}
else {
// subclone template(s)
// This calculation will consider all of the different templates
dErrorsFixed = dGetErrorsFixedForCustomPrimerSubcloneRead(
pPrimer,
bTopStrandNotBottomStrand,
nFirstBaseOfReadUnpadded,
false, // bChooseExperiment
false, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
}
if (
dErrorsFixed < pCP->dAutoFinishMinNumberOfErrorsFixedByAnExp_ )
return;
} // if ( bForCoveringSingleSubcloneRegions ) else
if ( bFoundPrimer ) {
cerr << "found primer before createExperimentWithCustomPrimer ";
for( int n = 0; n < pShortestPrimer->nUnpaddedLength_; ++n ) {
cerr << pShortestPrimer->szPrimer_[ n ];
}
cerr << " dErrorsFixed = " << dErrorsFixed << endl;
cerr << endl;
}
createExperimentWithCustomPrimer(
pPrimer,
ucRead,
bTopStrandNotBottomStrand,
nFirstBaseOfReadUnpadded,
bCloneTemplateNotSubcloneTemplate,
bForCoveringSingleSubcloneRegions,
dErrorsFixed );
}
void Contig :: createExperimentWithCustomPrimer(
primerType* pPrimer,
const readTypes ucRead,
const bool bTopStrandNotBottomStrand,
const int nFirstBaseOfReadUnpadded,
const bool bCloneTemplateNotSubcloneTemplate,
const bool bForCoveringSingleSubcloneRegions,
const double dErrorsFixed ) {
// we will get rid of the primer's memory so save primer info
primerType* pSavedPrimer = (primerType*) malloc( sizeof( primerType ) );
memcpy( pSavedPrimer, pPrimer, sizeof( primerType ) );
experiment* pExp = new experiment();
pExp->dErrorsFixed_ = dErrorsFixed;
pExp->ucStrandAndChemistry_ = ucRead;
// pExp->bCustomPrimerNotUniversalPrimer_ = true;
pExp->nReadType_ = nWalk;
pExp->bCloneTemplateNotSubcloneTemplate_ =
bCloneTemplateNotSubcloneTemplate;
pExp->fCost_ = ( bCloneTemplateNotSubcloneTemplate ) ?
consedParameters::pGetConsedParameters()->dAutoFinishCostOfCustomPrimerCloneReaction_ :
consedParameters::pGetConsedParameters()->dAutoFinishCostOfCustomPrimerSubcloneReaction_;
pExp->pCustomPrimer_ = pSavedPrimer;
pExp->pContig_ = this;
pExp->pSubcloneTemplate_ = NULL;
pExp->nFirstBaseOfReadUnpadded_ = nFirstBaseOfReadUnpadded;
// Note: these positions do *not* take into account the extents of
// the subclone templates
if ( bTopStrandNotBottomStrand ) {
pExp->nUnpaddedLeft_ = nFirstBaseOfReadUnpadded;
pExp->nUnpaddedRight_ = nFirstBaseOfReadUnpadded +
ConsEd::pGetAssembly()->pModelRead_->length() - 1;
}
else {
pExp->nUnpaddedRight_ = nFirstBaseOfReadUnpadded;
pExp->nUnpaddedLeft_ = nFirstBaseOfReadUnpadded -
ConsEd::pGetAssembly()->pModelRead_->length() + 1;
}
pExp->bTopStrandNotBottomStrand_ = bTopStrandNotBottomStrand;
pExp->nState_ = nAvailableExperiment;
// note in the below calculations that the nAutoFinishPotential...
// parameters are 0-based
if ( !bForCoveringSingleSubcloneRegions )
setVariousErrorData( pExp );
else {
orderTemplatesAndSetSingleSubcloneValues( pExp );
if ( pExp->nSingleSubcloneBasesFixed_ <
consedParameters::pGetConsedParameters()->nAutoFinishMinNumberOfSingleSubcloneBasesFixedByAnExp_ ) {
// this experiment doesn't fix enough single subclone bases--reject it
delete pExp;
return;
}
// in the case of custom primer subclone experiments for
// single subclone regions, the errors fixed has not yet been set
// (I actually don't know if it matters, because this value is not
// printed out--let's see if we can get away without even calculating
// it.
// if ( pExp->bIsCustomPrimerSubcloneRead() ) {
// pExp->dErrorsFixed_ = dGetErrorsFixedForCustomPrimerSubcloneRead(
// pPrimer,
// bTopStrandNotBottomStrand,
// nFirstBaseOfReadUnpadded,
// false, // bChooseExperiment
// false, // bJustConsensusNotGaps
// false ); // bFindWindowOfMostErrorsFixed
// }
}
pCandidateExperiments_->insert( pExp );
}
void Contig :: setVariousErrorData( experiment* pExp ) {
calculateErrorsFixedJustInConsensus( pExp );
if ( ( nGetUnpaddedStartIndex() <= pExp->nUnpaddedLeft_ )
&&
( pExp->nUnpaddedRight_ <= nGetUnpaddedEndIndex() ) ) {
pExp->nExtendsIntoWhichGap_ = nNoGap;
}
else {
if ( pExp->nUnpaddedLeft_ < nGetUnpaddedStartIndex() )
pExp->nExtendsIntoWhichGap_ = nLeftGap;
else
pExp->nExtendsIntoWhichGap_ = nRightGap;
}
pExp->dErrorsFixedPerDollar_ =
pExp->dErrorsFixed_ / pExp->fCost_;
// let's just drop the whole business (DG, 10/98)
// The .00000001 business is because, due to rounding error, the
// errors fixed in the consensus can be larger
// if ( ! (pExp->dErrorsFixedJustInConsensus_ <=
// pExp->dErrorsFixed_ + 0.0000000001 ) ) {
// printf( "problem with the following experiment: " );
// printf( "diff = %e\n",
// pExp->dErrorsFixedJustInConsensus_ -
// pExp->dErrorsFixed_ );
// printf( " pExp->dErrorsFixedJustInConsensus_ = %f \npExp->dErrorsFixed_ = %f\n",
// pExp->dErrorsFixedJustInConsensus_,
// pExp->dErrorsFixed_ );
// pExp->print();
// assert( false );
// }
}
bool Contig :: bStillTooManyExpectedErrors() {
if ( dExpectedNumberOfErrorsInExtendedContig_ > dTargetNumberOfErrors_ )
return true;
else
return false;
}
void Contig :: calculateErrorsFixedJustInConsensus( experiment* pExp ) {
// now calculate the dErrorsFixedJustInConsensus_
// If the high quality region and low quality regions are on the
// consensus entirely, then
// this is the same as the dErrorsFixed_
// Otherwise, we have to calculate it explicitly.
if (
( nGetUnpaddedStartIndex() <= pExp->nUnpaddedLeft_ )
&&
( pExp->nUnpaddedRight_ <= nGetUnpaddedEndIndex() )
)
pExp->dErrorsFixedJustInConsensus_ = pExp->dErrorsFixed_;
else {
// case in which the read extends out over the end of the consensus
// and thus only some of the errors fixed improve the consensus
// accuracy
if ( pExp->bIsWholeCloneRead() ) {
pExp->dErrorsFixedJustInConsensus_ =
dGetErrorsFixedForWholeCloneRead(
pExp->bTopStrandNotBottomStrand_,
pExp->nFirstBaseOfReadUnpadded_,
false, // bChooseExperiment
true, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
}
else if ( pExp->bIsCustomPrimerSubcloneRead() ) {
pExp->dErrorsFixedJustInConsensus_ =
dGetErrorsFixedForCustomPrimerSubcloneRead(
pExp->pCustomPrimer_,
pExp->bTopStrandNotBottomStrand_,
pExp->nFirstBaseOfReadUnpadded_,
false, // bChooseExperiment
true, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
}
else if ( pExp->bIsUniversalPrimerRead() ) {
pExp->dErrorsFixedJustInConsensus_ =
dGetErrorsFixedForSubcloneRead(
&(pExp->pSubcloneTemplate_),
1, // nNumberOfSubclones
pExp->bTopStrandNotBottomStrand_,
pExp->nFirstBaseOfReadUnpadded_,
true, // bExcludeVectorBasesAtBeginningOfRead
false, // bChooseExperiment
true, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
}
else
assert( false );
} // if ( ( nGetUnpaddedStartIndex() <= ... else
pExp->dErrorsFixedJustInConsensusPerDollar_ =
pExp->dErrorsFixedJustInConsensus_ / pExp->fCost_;
}
void Contig :: adjustExpectedNumberOfErrors( experiment* pExp ) {
dExpectedNumberOfErrorsInExtendedContig_ -= pExp->dErrorsFixed_;
dExpectedNumberOfErrorsInConsensus_ -= pExp->dErrorsFixedJustInConsensus_;
}
experiment* Contig :: pGetBestExperimentCandidateButMightHaveBeenTriedBefore() {
// this has been written so that it only returns experiments
// that fix the min required number of errors in the consensus
experiment* pBestSoFar = NULL;
double dErrorsFixedPerDollarBestSoFar = 0.0;
for( int nCan = 0; nCan < pCandidateExperiments_->length(); ++nCan ) {
experiment* pExp = (*pCandidateExperiments_)[ nCan ];
if ( pExp->nState_ == nChosenExperiment ) continue;
if ( pExp->nState_ == nRejectedExperiment ) continue;
// commented out because rounding error can make this not true
// assert( pExp->dErrorsFixedJustInConsensus_ <=
// pExp->dErrorsFixed_ + 0.0000000001 );
// I don't believe we should have any requirement about fixing a minimum
// number of errors in the consensus. DG, Jan 2000. But we should
// have a minimum requirement about fixing total errors.
if ( pExp->dErrorsFixed_ <
consedParameters::pGetConsedParameters()->dAutoFinishMinNumberOfErrorsFixedByAnExpToUse_ )
continue;
if ( pExp->dErrorsFixedPerDollar_ > dErrorsFixedPerDollarBestSoFar ) {
// check that, if the read is a subclone custom primer read,
// it isn't too close to another subclone custom primer read
// on the same strand that has already been chosen
if ( pCP->bAutoFinishDoNotAllowSubcloneCustomPrimerReadsCloseTogether_ ) {
if ( pExp->nReadType_ == nWalk &&
!pExp->bCloneTemplateNotSubcloneTemplate_ &&
bThereIsAnotherSubcloneCustomPrimerReadTooCloseOnTheSameStrand( pExp ) ) {
// this experiment can never be considered again in this
// round of autofinish, so reject it
pExp->nState_ = nRejectedExperiment;
continue;
}
}
if ( pCP->bAutoFinishDoNotAllowWholeCloneCustomPrimerReadsCloseTogether_ &&
pExp->nReadType_ == nWalk &&
pExp->bCloneTemplateNotSubcloneTemplate_ &&
bThereIsAnotherWholeCloneCustomPrimerReadTooCloseOnTheSameStrand( pExp ) ) {
pExp->nState_ = nRejectedExperiment;
continue;
}
// If the experiment is partly off the extended high quality
// segment of the
// consensus, then check that it overlaps be enough
if (
( pExp->nUnpaddedLeft_ < nUpdatedUnpaddedHighQualitySegmentStart_ )
||
( nUnpaddedHighQualitySegmentEnd_ < pExp->nUnpaddedRight_ ) ) {
if ( nGetAmountOfOverlap( pExp->nUnpaddedLeft_,
pExp->nUnpaddedRight_,
nUnpaddedHighQualitySegmentStart_,
nUnpaddedHighQualitySegmentEnd_ ) <
pCP->nAutoFinishMinBaseOverlapBetweenAReadAndHighQualitySegmentOfConsensus_ ) {
fprintf( pAO, "\nrejecting experiment from %d to %d because it doesn't overlap the high quality segment %d to %d by at least %d bases\n",
pExp->nUnpaddedLeft_,
pExp->nUnpaddedRight_,
nUnpaddedHighQualitySegmentStart_,
nUnpaddedHighQualitySegmentEnd_,
pCP->nAutoFinishMinBaseOverlapBetweenAReadAndHighQualitySegmentOfConsensus_ );
continue;
}
}
pBestSoFar = pExp;
dErrorsFixedPerDollarBestSoFar = pExp->dErrorsFixedPerDollar_;
}
}
// note: pBestSoFar might be null if all experiments are chosen
return( pBestSoFar );
}
experiment* Contig :: pGetBestExperimentCandidate() {
experiment* pBestExperiment = NULL;
bool bTryAgain = false;
do {
pBestExperiment = pGetBestExperimentCandidateButMightHaveBeenTriedBefore();
if ( pBestExperiment ) {
if ( pBestExperiment->bHasThisExperimentOrOneLikeItBeenTriedBefore() ) {
pBestExperiment->nState_ = nRejectedExperiment;
bTryAgain = true;
}
else
bTryAgain = false;
}
else
bTryAgain = false;
} while( bTryAgain );
return( pBestExperiment );
}
void Contig :: adjustRemainingExperimentCandidates(
experiment* pExperimentJustChosen ) {
// first, adjust the ErrorsFixed arrays
adjustConsensusDistributionsAndArrays( pExperimentJustChosen );
// don't bother adjusting the 1st experiment, since that is
// the one we just chose
if ( pCandidateExperiments_->length() < 2 )
return;
for( int nExperiment = 0;
nExperiment < pCandidateExperiments_->length();
++nExperiment ) {
experiment* pExperimentToBeAdjusted = (*pCandidateExperiments_)[ nExperiment ];
if (pExperimentToBeAdjusted == pExperimentJustChosen ) continue;
// note that once an experiment falls below the minimum number of
// errors fixed, it is rejected and won't be adjusted or considered
// again
if (pExperimentToBeAdjusted->nState_ != nAvailableExperiment )
continue;
if ( pExperimentJustChosen->bOverlaps( pExperimentToBeAdjusted ) ) {
adjustErrorsFixed( pExperimentToBeAdjusted );
}
}
}
// After having adjusted the pUnpaddedErrorProbabilities_ array,
// recalculate the errors fixed by each experiment. Note that it is
// only necessary to do this to the experiments that overlap the one
// that was just chosen, since any other experiment will fix the same
// number of errors that it did before.
void Contig :: adjustErrorsFixed( experiment* pExperimentToBeAdjusted ) {
if ( pExperimentToBeAdjusted->bIsCustomPrimerSubcloneRead() ) {
// if ( pExperimentToBeAdjusted->nUnpaddedLeft_ == -177 ) {
// char szPrimer[ 200 ];
// for( int n = 0; n < pExperimentToBeAdjusted->pCustomPrimer_->nUnpaddedLength_; ++n )
// szPrimer[n] = pExperimentToBeAdjusted->pCustomPrimer_->szPrimer_[n];
// szPrimer[ pExperimentToBeAdjusted->pCustomPrimer_->nUnpaddedLength_ ] = '\0';
// fprintf( pAO, "found experiment at -177\n" );
// if ( strcmp( szPrimer, "cctccacacctcccc" ) == 0 ) {
// fprintf( pAO, "found experiment we think should not have been chosen:" );
// pExperimentToBeAdjusted->print();
// bFoundDebug = true;
// }
// }
pExperimentToBeAdjusted->dErrorsFixed_ =
dGetErrorsFixedForCustomPrimerSubcloneRead(
pExperimentToBeAdjusted->pCustomPrimer_,
pExperimentToBeAdjusted->bTopStrandNotBottomStrand_,
pExperimentToBeAdjusted->nFirstBaseOfReadUnpadded_,
false, // bChooseExperiment
false, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
if ( ConsEd::pGetAssembly()->bFoundAutoFinishDebugPrimer( pExperimentToBeAdjusted->pCustomPrimer_ ) ) {
fprintf( pAO, "\nfound primer after adjustErrorsFixed--printing:\n" );
pExperimentToBeAdjusted->print();
}
}
else if ( pExperimentToBeAdjusted->bIsWholeCloneRead() ) {
// whole clone reads
pExperimentToBeAdjusted->dErrorsFixed_ =
dGetErrorsFixedForWholeCloneRead(
pExperimentToBeAdjusted->bTopStrandNotBottomStrand_,
pExperimentToBeAdjusted->nFirstBaseOfReadUnpadded_,
false, // bChooseExperiment
false, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
if ( ConsEd::pGetAssembly()->bFoundAutoFinishDebugPrimer( pExperimentToBeAdjusted->pCustomPrimer_ ) ) {
fprintf( pAO, "\nfound debug read after adjustErrorsFixed--printing:\n" );
pExperimentToBeAdjusted->print();
}
}
else if ( pExperimentToBeAdjusted->bIsUniversalPrimerRead() ) {
// universal primer reads
bool bJustCountConsensusErrorsFixed = false;
if ( !pExperimentToBeAdjusted->bDeNovoUniversalPrimerRead_ &&
!pCP->bAutoFinishAllowResequencingReadsToExtendContigs_ )
bJustCountConsensusErrorsFixed = true;
pExperimentToBeAdjusted->dErrorsFixed_ =
dGetErrorsFixedForSubcloneRead(
&(pExperimentToBeAdjusted->pSubcloneTemplate_),
1, // nNumberOfSubclones
pExperimentToBeAdjusted->bTopStrandNotBottomStrand_,
pExperimentToBeAdjusted->nFirstBaseOfReadUnpadded_,
true, // bExcludeVectorBasesAtBeginningOfRead
false, // bChooseExperiment
bJustCountConsensusErrorsFixed, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
if ( ConsEd::pGetAssembly()->bFoundAutoFinishDebugRead( pExperimentToBeAdjusted ) ) {
fprintf( pAO, "\nfound debug read after adjustErrorsFixed--printing:\n" );
pExperimentToBeAdjusted->print();
}
}
else
assert( false );
calculateErrorsFixedJustInConsensus( pExperimentToBeAdjusted );
pExperimentToBeAdjusted->dErrorsFixedPerDollar_ =
pExperimentToBeAdjusted->dErrorsFixed_ / pExperimentToBeAdjusted->fCost_;
// note: I need to check where this routine is used. If it is only used
// by parts of the code that won't pick an experiment below the
// error rate threshold, there is no point in recalculating the errors
// here if the error rate is already below the threshold--in that
// case this should be near the top of the routine. (DG 11/99)
// This routine is only used for available experiments that overlap
// an experiment just picked. Thus the errors fixed do need to be
// recalculated.
if ( pExperimentToBeAdjusted->dErrorsFixedJustInConsensus_ < consedParameters::pGetConsedParameters()->dAutoFinishMinNumberOfErrorsFixedByAnExpToUse_ ) {
pExperimentToBeAdjusted->nState_ = nRejectedExperiment;
}
}
// When an experiment is chosen, adjust the pUnpaddedErrorProbabilities_
// to reflect the fact that some of the errors have been fixed in
// the region covered by the chosen experiment.
void Contig :: adjustConsensusDistributionsAndArrays(
experiment* pExperimentJustChosen ) {
// if ( pExperimentJustChosen->nUnpaddedLeft_ == -177 ) {
// char szPrimer[ 200 ];
// for( int n = 0; n < pExperimentJustChosen->pCustomPrimer_->nUnpaddedLength_; ++n )
// szPrimer[n] = pExperimentJustChosen->pCustomPrimer_->szPrimer_[n];
// szPrimer[ pExperimentJustChosen->pCustomPrimer_->nUnpaddedLength_ ] = '\0';
// fprintf( pAO, "found experiment at -177\n" );
// if ( strcmp( szPrimer, "cctccacacctcccc" ) == 0 ) {
// fprintf( pAO, "found experiment we think should not have been chosen:" );
// pExperimentJustChosen->print();
// bFoundDebug = true;
// }
// }
// 3 cases: custom primer subclone, whole clone, univ primer
double dErrorsFixed;
if ( pExperimentJustChosen->bIsUniversalPrimerRead() ) {
bool bJustCountConsensusErrorsFixed = false;
if ( !pExperimentJustChosen->bDeNovoUniversalPrimerRead_ &&
!pCP->bAutoFinishAllowResequencingReadsToExtendContigs_ )
bJustCountConsensusErrorsFixed = true;
dErrorsFixed = dGetErrorsFixedForSubcloneRead(
&(pExperimentJustChosen->pSubcloneTemplate_),
1, // nNumberOfSubclones
pExperimentJustChosen->bTopStrandNotBottomStrand_,
pExperimentJustChosen->nFirstBaseOfReadUnpadded_,
true, // bExcludeVectorBasesAtBeginningOfRead
true, // bChooseExperiment
bJustCountConsensusErrorsFixed, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
if ( ConsEd::pGetAssembly()->bFoundAutoFinishDebugRead( pExperimentJustChosen ) ) {
fprintf( pAO, "\nfound debug read in adjustConsensusDistributionsAndArrays, fixed %g--printing:\n",
dErrorsFixed );
pExperimentJustChosen->print();
}
}
else if ( pExperimentJustChosen->bIsCustomPrimerSubcloneRead() ) {
dErrorsFixed = dGetErrorsFixedForCustomPrimerSubcloneRead(
pExperimentJustChosen->pCustomPrimer_,
pExperimentJustChosen->bTopStrandNotBottomStrand_,
pExperimentJustChosen->nFirstBaseOfReadUnpadded_,
true, // bChooseExperiment
false, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
}
else if ( pExperimentJustChosen->bIsWholeCloneRead() ) {
dErrorsFixed = dGetErrorsFixedForWholeCloneRead(
pExperimentJustChosen->bTopStrandNotBottomStrand_,
pExperimentJustChosen->nFirstBaseOfReadUnpadded_,
true, // bChooseExperiment
false, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
}
else
assert( false ); // shouldn't be any other kind of experiment
if (
( ABS( dErrorsFixed - pExperimentJustChosen->dErrorsFixed_ ) >
.000001 ) &&
( pExperimentJustChosen->nPurpose_ != nCoverSingleSubcloneBases ) ) {
fprintf( pAO, "debugging: adjustConsensusDistributionsAndArrays, dErrorsFixed = %e but pExperimentJustChosen->dErrorsFixed_ = %e for experiment from %d to %d id start: %d\n",
dErrorsFixed,
pExperimentJustChosen->dErrorsFixed_,
pExperimentJustChosen->nUnpaddedLeft_,
pExperimentJustChosen->nUnpaddedRight_,
pExperimentJustChosen->aExpID_[0] );
cerr << "alert--see \"debugging\" in .out file--there is a bug" << endl;
}
}
void Contig :: printNoExperiments() {
fprintf( pAO, "CONTIG: %s does not need any experiments because it is already below the error threshold ",
(char*) soGetName().data() );
fprintf( pAO, "ACCEPTABLE ERRORS: %.1f CURRENT ERRORS EXTENDED CONTIG: %.1f CONSENSUS: %.1f\n ",
dTargetNumberOfErrors_,
dExpectedNumberOfErrorsInExtendedContig_,
dExpectedNumberOfErrorsInConsensus_ );
}
void Contig :: setTargetNumberOfErrorsForThisContig() {
// this includes gap errors
int nBasesToCareAbout = 0;
for( int nUnpadded = pOriginalUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpadded <= pOriginalUnpaddedErrorProbabilities_->nGetEndIndex();
++nUnpadded ) {
if ( (*pOriginalUnpaddedErrorProbabilities_)[ nUnpadded ] > 0.0 )
++nBasesToCareAbout;
}
dTargetNumberOfErrors_ =
( (double) consedParameters::pGetConsedParameters()->nAutoFinishMaxAcceptableErrorsPerMegabase_ * (double) nBasesToCareAbout ) / (double) 1000000.0;
}
void Contig :: setExpectedNumberOfErrorsForThisContig() {
// this includes errors both in the consensus and in the gaps
dExpectedNumberOfErrorsInExtendedContig_ = 0.0;
int nUnpadded;
for( nUnpadded = pOriginalUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpadded <= pOriginalUnpaddedErrorProbabilities_->nGetEndIndex();
++nUnpadded )
dExpectedNumberOfErrorsInExtendedContig_ +=
(*pOriginalUnpaddedErrorProbabilities_)[ nUnpadded ];
dExpectedNumberOfErrorsInConsensus_ = 0.0;
for( nUnpadded = nGetUnpaddedStartIndex();
nUnpadded <= nGetUnpaddedEndIndex();
++nUnpadded )
dExpectedNumberOfErrorsInConsensus_ +=
(*pOriginalUnpaddedErrorProbabilities_)[ nUnpadded ];
}
void Contig :: recalculateCurrentExpectedNumberOfErrorsForThisContig(
double& dErrorsInConsensus,
double& dErrorsInExtendedContig ) {
dErrorsInExtendedContig = 0.0;
int nUnpadded;
for( nUnpadded = pUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpadded <= pUnpaddedErrorProbabilities_->nGetEndIndex();
++nUnpadded ) {
dErrorsInExtendedContig +=
(*pUnpaddedErrorProbabilities_)[ nUnpadded ];
}
dErrorsInConsensus = 0.0;
for( nUnpadded = nGetUnpaddedStartIndex();
nUnpadded <= nGetUnpaddedEndIndex();
++nUnpadded ) {
dErrorsInConsensus += (*pUnpaddedErrorProbabilities_)[ nUnpadded ];
}
}
double Contig :: dGetExpectedNumberOfErrorsPerMegabaseForConsensus() {
// this includes only errors in the consensus
int nUnpaddedBases = 0;
double dExpectedNumberOfErrors = 0.0;
double dError;
double dQuality;
for( int nPos = nGetStartConsensusIndex();
nPos <= nGetEndConsensusIndex();
++nPos ) {
// tried ntGetCons( nPos ).bIsPad() but that still called the
// operator=
Ntide nt = ntGetCons( nPos );
if ( !nt.bIsPad() ) {
Quality qual = nt.qualGetQuality();
// do not include edited bases in the quality
// calculation
if (qual == ucQualityLowEdited ) continue;
if (qual == ucQualityHighEdited ) continue;
dError = dErrorRateFromQuality[ qual ];
dExpectedNumberOfErrors += dError;
++nUnpaddedBases;
}
}
if ( nUnpaddedBases > 0 )
dExpectedNumberOfErrors =
(dExpectedNumberOfErrors * 1000000.0)/nUnpaddedBases;
else
dExpectedNumberOfErrors = -1;
return( dExpectedNumberOfErrors );
}
bool Contig :: bFindTemplatesForPrimer( primerType* pPrimer,
const readTypes ucRead ) {
// Then we use the template's
// error rates (and the redundancy) to see which template will be
// best and will cover. (I am not going to use Pat's
// parameter for length of not--there might be situations in which
// we would only have a very short template, but it would be good
// enough for the problem at hand.)
// We are only interested in templates that cover this area:
// (primer) (high quality)
// ---> llllllHHHHHHHHHHHHHHllllllllllllllllllllllllllllll
// a d b where a = nReadConsPosLeft and b = nReadConsPosRight
// ----- (minimum acceptable extents of template--one more
// base than primer)
// try templates like this:
// ----------- (not ok)
// ------------------------------ (ok)
// -------- (not ok)
// ----------------- ok
// ------------------------------
// ------------------------------
// ------------------------------
// ------------------------------
// all of the ones below are past the primer, and all further
// ones will be also, so we can't stop looking at templates
// ------------------------------
// ------------------------------
// ------------------------------
//
// or for bottom strand templates:
//
// llllllllllllHHHHHHHHHllll <---
// ----- (minimum acceptable extents
// of template--one more base
// than primer)
// try templates like this:
// -------------------------------- (not ok)
// --------------------------------------- ok
// ------------------------------ (not ok)
// -------------------------------------- ok
// --------------------------------------
// --------------------------------------
// --------------------------------------
// --------------------------------------
// --------------------------------------
// .
// .
// .
// --------------------------------------
// all of the ones below are no good because they don't overlap
// the high quality or low quality part of the read at all. Since
// the templates are sorted by left end, we can stop looking at templates
//
// --------------------------------------
// --------------------------------------
// --------------------------------------
// --------------------------------------
bool bTopStrandNotBottomStrand = bIsTopStrandRead( ucRead );
int nReadUnpaddedLeft;
int nReadUnpaddedRight;
if ( bTopStrandNotBottomStrand ) {
// note that even though the read won't start until
// pPrimer->nUnpaddedEnd_ + 1, the template must start at
// pPrimer->nUnpaddedStart_ or to the left of it. Otherwise
// the primer won't bind to the template.
nReadUnpaddedLeft = pPrimer->nUnpaddedStart_;
nReadUnpaddedRight = pPrimer->nUnpaddedStart_ +
pCP->nPrimersWhenChoosingATemplateMinPotentialReadLength_;
}
else {
// for reverse primers, nUnpaddedEnd_ is still
// the rightmost base of the primer
nReadUnpaddedRight = pPrimer->nUnpaddedEnd_;
nReadUnpaddedLeft = pPrimer->nUnpaddedEnd_ -
pCP->nPrimersWhenChoosingATemplateMinPotentialReadLength_;
}
// since the templates are sorted by left-end, let's be generous and
// look back further than is possible for a template to start and
// still cover this read (if the user is honest with
// the maximum template size)
subcloneTTemplate* pSubTT = new subcloneTTemplate();
pSubTT->nUnpaddedStart_ = nReadUnpaddedRight -
consedParameters::pGetConsedParameters()->nPrimersMaxInsertSizeOfASubclone_;
int nTemplate = aSubcloneTemplates_.nFindIndexOfMatchOrSuccessor( pSubTT );
delete pSubTT; // this was just used to find the first one
if ( nTemplate == RW_NPOS ) {
pPrimer->bAcceptable_ = false;
pPrimer->nWhyIsPrimerNotAcceptable_ = BAD_PRIMER_HAS_NO_TEMPLATE;
return( false ); // no templates cover the read
}
bool bTemplatesStillMightCover = true;
for( ; nTemplate < aSubcloneTemplates_.entries() && bTemplatesStillMightCover; ++nTemplate ) {
pSubTT = aSubcloneTemplates_[ nTemplate ];
if ( pSubTT->nUnpaddedStart_ > nReadUnpaddedLeft )
bTemplatesStillMightCover = false;
else {
if ( !pSubTT->bOKToUseThisTemplate_ ) continue;
// check if the end of the template is to the right of the
// read we want to make
if ( pSubTT->nUnpaddedEnd_ < nReadUnpaddedRight ) continue;
// if reached here, the template covers the desired read
// check that this template can be used for a read in the
// desired orientation
if ( bTopStrandNotBottomStrand ) {
if ( !pSubTT->bOKToUseForTopStrandReads_ ) {
continue;
}
}
else {
if ( !pSubTT->bOKToUseForBottomStrandReads_ ) {
continue;
}
}
pSubTT->tryToAddTemplateToPrimer( pPrimer );
}
}
int nNumberOfTemplates = 0;
for( nTemplate = 0; nTemplate < nNUMBER_OF_TEMPLATES; ++nTemplate ) {
if ( pPrimer->pSubcloneTemplate_[ nTemplate ] )
++nNumberOfTemplates;
else
break;
}
// check to see that at least one template was found
if ( nNumberOfTemplates == 0 ) {
pPrimer->bAcceptable_ = false;
pPrimer->nWhyIsPrimerNotAcceptable_ = BAD_PRIMER_HAS_NO_TEMPLATE;
return( false );
}
if ( nNumberOfTemplates < pCP->nPrimersMinNumberOfTemplatesForPrimers_ ) {
pPrimer->bAcceptable_ = false;
pPrimer->nWhyIsPrimerNotAcceptable_ = BAD_PRIMER_NOT_ENOUGH_TEMPLATES;
return( false );
}
return( true );
}
void Contig :: createExperimentsWithResequencingUniversalPrimers() {
for( int nTemplate = 0; nTemplate < aSubcloneTemplates_.entries();
++nTemplate ) {
subcloneTTemplate* pSubcloneTemplate = aSubcloneTemplates_[ nTemplate ];
if ( !pSubcloneTemplate->bOKToUseThisTemplate_ ) continue;
LocatedFragment* pLocFragOfExistingTopStrandUniversalPrimer = NULL;
LocatedFragment* pLocFragOfExistingBottomStrandUniversalPrimer = NULL;
for( int nRead = 0; nRead < pSubcloneTemplate->aReads_.entries();
++nRead ) {
LocatedFragment* pLocFrag = pSubcloneTemplate->aReads_[ nRead ];
// Below has the effect of creating an experiment for each
// universal primer dye primer read. This is ok, but
// inefficient in 2 ways:
// 1) there may already by a dye
// terminator read and thus the new read won't fix any errors
// and
// 2) there may be more than one dye primer read.
//
// There is a better way--just check to see if there is any
// dye primer read and no dye terminator read. This would
// not be helpful, though, in the redundant phase. For now,
// I'm not changing it since efficiency is not a problem
// currently (8/99 DG)
// DG 11/99: Now that the initial and redundant phases are done
// at once, we can do the following: just check if there is any
// top strand read from this template. If so, do another. Check
// if there is any bottom strand read from this template. If so,
// do another. Don't worry about whether the existing reads
// are dye primer or dye terminator reads. That will be taken
// care of by the error probability arrays--reads of redundant
// chemistry/strand will fix fewer errors.
// fixed error 9/02 Jody Schwartz, LBL,
// in which the existing read was a pcr end read
// Autofinish tried to do a resequence from it and got a
// segmentation fault.
if ( !pLocFrag->bIsUniversalPrimerRead() ) continue;
// this can happen since some templates have reads in
// different contigs
if ( pLocFrag->pGetContig() != this ) continue;
if ( pLocFrag->bComp() ) {
if ( !pLocFragOfExistingBottomStrandUniversalPrimer )
pLocFragOfExistingBottomStrandUniversalPrimer = pLocFrag;
else {
// which read is better quality?
if ( pLocFrag->nGetHighQualitySegmentLengthSortOfUnpadded() >
pLocFragOfExistingBottomStrandUniversalPrimer->nGetHighQualitySegmentLengthSortOfUnpadded() )
pLocFragOfExistingBottomStrandUniversalPrimer = pLocFrag;
}
}
else {
if ( !pLocFragOfExistingTopStrandUniversalPrimer )
pLocFragOfExistingTopStrandUniversalPrimer = pLocFrag;
else {
if ( pLocFrag->nGetHighQualitySegmentLengthSortOfUnpadded() >
pLocFragOfExistingTopStrandUniversalPrimer->nGetHighQualitySegmentLengthSortOfUnpadded() )
pLocFragOfExistingTopStrandUniversalPrimer = pLocFrag;
}
}
} // for( int nRead = 0; ...
if ( pLocFragOfExistingTopStrandUniversalPrimer )
createExperimentWithResequencingUniversalPrimer(
pLocFragOfExistingTopStrandUniversalPrimer,
pSubcloneTemplate );
if ( pLocFragOfExistingBottomStrandUniversalPrimer )
createExperimentWithResequencingUniversalPrimer(
pLocFragOfExistingBottomStrandUniversalPrimer,
pSubcloneTemplate );
}
}
void Contig :: createExperimentWithResequencingUniversalPrimer(
LocatedFragment* pExistingLocFrag,
subcloneTTemplate* pSubcloneTemplate ) {
bool bTopStrandNotBottomStrand = !pExistingLocFrag->bComp();
int nFirstBaseOfReadUnpadded;
// these used to be off by 1 so that you would see 79 399 in the
// autoFinishExpTag's (DG 7/29/99)
// DG 11/99: The location of the beginning of the new read is best
// determined by an existing universal primer read--*not* by the location
// of the vector/insert junction since universal primer reads start
// *way* before (80 or so bases) the vector/insert junction. I changed
// the above so we check which read is highest quality and use that one
// if there is more than one
if ( bTopStrandNotBottomStrand ) {
nFirstBaseOfReadUnpadded = pExistingLocFrag->nGetAlignStartUnpadded();
}
else {
nFirstBaseOfReadUnpadded = pExistingLocFrag->nGetAlignEndUnpadded();
}
readTypes ucStrandAndChemistry =
( bTopStrandNotBottomStrand ?
TOP_STRAND_TERMINATOR_READS :
BOTTOM_STRAND_TERMINATOR_READS );
// if we know where the end of the template is, we should use
// that information
// if we have already run into vector, use that information
// to shorten the length of the possible read
int nUnpaddedLeft;
int nUnpaddedRight;
if ( bTopStrandNotBottomStrand ) {
nUnpaddedLeft = nFirstBaseOfReadUnpadded;
nUnpaddedRight = nFirstBaseOfReadUnpadded +
ConsEd::pGetAssembly()->pModelRead_->length() - 1;
}
else {
nUnpaddedRight = nFirstBaseOfReadUnpadded;
nUnpaddedLeft = nFirstBaseOfReadUnpadded -
ConsEd::pGetAssembly()->pModelRead_->length() + 1;
}
// the read starts at the end of the X's rather than
// at the start of the existing read
if ( pSubcloneTemplate->bUnpaddedStartKnownPrecisely_ ) {
nUnpaddedLeft =
MAX( nUnpaddedLeft,
pSubcloneTemplate->nUnpaddedStart_ );
}
if ( pSubcloneTemplate->bUnpaddedEndKnownPrecisely_ ) {
nUnpaddedRight =
MIN( nUnpaddedRight,
pSubcloneTemplate->nUnpaddedEnd_ );
}
if ( nUnpaddedLeft > nUnpaddedRight ) return;
if ( pCP->bAutoFinishAllowResequencingReadsOnlyForRunsAndStops_ ) {
bool bCrossesARun;
bool bCrossesAStop;
if ( !bDoesRegionCrossARunOrStop( nUnpaddedLeft,
nUnpaddedRight,
bCrossesARun,
bCrossesAStop ) )
return;
}
// since dGetErrorsFixed is expensive, let's first see if we can
// quickly eliminate this read cheaply
double dMaxErrorsFixed =
dGetInitialMaxErrorsFixed( nUnpaddedLeft,
nUnpaddedRight );
if ( dMaxErrorsFixed < pCP->dAutoFinishMinNumberOfErrorsFixedByAnExpToUse_ )
return;
// now, if the read has made the cut, then check it more rigorously
bool bJustCountConsensusErrorsFixed =
pCP->bAutoFinishAllowResequencingReadsToExtendContigs_ ?
false : true;
double dErrorsFixed = dGetErrorsFixedForSubcloneRead(
&pSubcloneTemplate,
1, // only 1 subclone
bTopStrandNotBottomStrand,
nFirstBaseOfReadUnpadded,
true, // exclude vector bases
false, // bChooseExperiment
bJustCountConsensusErrorsFixed,
// bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
// don't suggest an experiment if the errors fixed
// is not sufficiently high
if ( dErrorsFixed <
consedParameters::pGetConsedParameters()->dAutoFinishMinNumberOfErrorsFixedByAnExpToUse_ )
return;
experiment* pExp = new experiment();
pExp->dErrorsFixed_ = dErrorsFixed;
pExp->ucStrandAndChemistry_ = ucStrandAndChemistry;
// pExp->bCustomPrimerNotUniversalPrimer_ = false;
pExp->nReadType_ = pExistingLocFrag->nReadType_;
pExp->bCloneTemplateNotSubcloneTemplate_ = false;
pExp->fCost_ =
consedParameters::pGetConsedParameters()->dAutoFinishCostOfResequencingUniversalPrimerSubcloneReaction_;
pExp->pCustomPrimer_ = NULL;
pExp->pContig_ = this;
pExp->pSubcloneTemplate_ = pSubcloneTemplate;
pExp->nFirstBaseOfReadUnpadded_ = nFirstBaseOfReadUnpadded;
pExp->nUnpaddedLeft_ = nUnpaddedLeft;
pExp->nUnpaddedRight_ = nUnpaddedRight;
pExp->bTopStrandNotBottomStrand_ = bTopStrandNotBottomStrand;
pExp->nState_ = nAvailableExperiment;
assert( pExp->nUnpaddedLeft_ <= pExp->nUnpaddedRight_ );
setVariousErrorData( pExp );
pCandidateExperiments_->insert( pExp );
if ( ConsEd::pGetAssembly()->bFoundAutoFinishDebugRead( pExp ) ) {
fprintf( pAO, "\nfound debug read in createExperimentWithResequencingUniversalPrimer--printing:\n" );
pExp->print();
}
}
typedef struct {
LocatedFragment* pLocFrag_;
int nVectorInsertJunction_; // mark on not-X base
} vectorReadType;
static int nSortByXNotXJunction( const vectorReadType* pRead1,
const vectorReadType* pRead2 ) {
if ( pRead1->nVectorInsertJunction_ < pRead2->nVectorInsertJunction_ )
return( -1 );
else if ( pRead1->nVectorInsertJunction_ > pRead2->nVectorInsertJunction_ )
return( 1 );
else
return( 0 );
}
bool Contig :: bIsThisEndTheEndOfTheClone( const int nWhichEnd ) {
if ( nWhichEnd == nLeftGap ) {
if ( !bLeftEndOfContigExaminedForEndOfClone_ )
determineWhetherEndOfClone( nWhichEnd );
return( bLeftEndOfContigIsEndOfClone_ );
}
else {
if ( !bRightEndOfContigExaminedForEndOfClone_ )
determineWhetherEndOfClone( nWhichEnd );
return( bRightEndOfContigIsEndOfClone_ );
}
}
void Contig :: determineWhetherEndOfClone( const int nWhichEnd ) {
bool bEndOfClone;
if ( pCP->bAutoFinishCDNANotGenomic_ )
bEndOfClone = bIsThisEndTheEndOfTheCloneCDNA( nWhichEnd );
else
bEndOfClone = bIsThisEndTheEndOfTheCloneGenomic( nWhichEnd );
if ( nWhichEnd == nLeftGap ) {
bLeftEndOfContigExaminedForEndOfClone_ = true;
bLeftEndOfContigIsEndOfClone_ = bEndOfClone;
}
else {
bRightEndOfContigExaminedForEndOfClone_ = true;
bRightEndOfContigIsEndOfClone_ = bEndOfClone;
}
}
bool Contig :: bIsThisEndTheEndOfTheCloneCDNA( const int nWhichEnd ) {
const int nAmountIntoConsensus = 50;
const int nAmountOffConsensus = 50;
int nWindowToLookLeft;
int nWindowToLookRight;
bool bLookForTopStrandNotBottomStrandRead;
if ( nWhichEnd == nLeftGap ) {
nWindowToLookLeft = nGetStartConsensusIndex() - nAmountOffConsensus;
nWindowToLookRight = nGetStartConsensusIndex() + nAmountIntoConsensus;
bLookForTopStrandNotBottomStrandRead = true;
}
else {
nWindowToLookLeft = nGetEndConsensusIndex() - nAmountIntoConsensus;
nWindowToLookRight = nGetEndConsensusIndex() + nAmountOffConsensus;
bLookForTopStrandNotBottomStrandRead = false;
}
ContigView* pContigView = pGetContigView( nWindowToLookLeft,
nWindowToLookRight );
bool bThisIsCloneEnd = false;
for( int nRead = 0; nRead < pContigView->nGetNumberOfFragments();
++nRead ) {
LocatedFragment* pLocFrag = pContigView->pLocatedFragmentGet( nRead );
if ( bLookForTopStrandNotBottomStrandRead == pLocFrag->bComp() ) continue;
if ( pLocFrag->nReadType_ == nUniversalForward ||
pLocFrag->nReadType_ == nUniversalReverse ) {
bThisIsCloneEnd = true;
break;
}
}
return( bThisIsCloneEnd );
}
bool Contig :: bIsThisEndTheEndOfTheCloneGenomic( const int nWhichEnd ) {
const int nAmountIntoConsensus = 50;
const int nAmountOffConsensus = 10;
int nWindowToLookLeft;
int nWindowToLookRight;
if ( nWhichEnd == nLeftGap ) {
nWindowToLookLeft = nGetStartConsensusIndex() - nAmountOffConsensus;
nWindowToLookRight = nGetStartConsensusIndex() + nAmountIntoConsensus;
}
else {
nWindowToLookLeft = nGetEndConsensusIndex() - nAmountIntoConsensus;
nWindowToLookRight = nGetEndConsensusIndex() + nAmountOffConsensus;
}
// get all reads in a big window (60 bases) around the end of the consensus
ContigView* pContigView = pGetContigView( nWindowToLookLeft,
nWindowToLookRight );
if ( pContigView->nGetNumberOfFragments() < 2 )
return( false );
vectorReadType* pVectorReadArray = ( vectorReadType*) calloc( pContigView->nGetNumberOfFragments() + 1,
sizeof( vectorReadType ) );
int nNumberOfXNotXReads = 0;
const int nLookingForX = 1;
const int nLookingForNotX = 2;
int nRead;
for( nRead = 0; nRead < pContigView->nGetNumberOfFragments();
++nRead ) {
LocatedFragment* pFrag = pContigView->pLocatedFragmentGet( nRead );
if ( pFrag->bIsWholeReadLowQuality() ) continue;
int nHighQualityStart = pFrag->nGetHighQualityStart();
int nHighQualityEnd = pFrag->nGetHighQualityEnd();
if ( !bIntervalsIntersect( nHighQualityStart,
nHighQualityEnd,
nWindowToLookLeft,
nWindowToLookRight ) )
continue;
// if reached here, we have a read whose high quality segment
// intersects the window around the end of the consensus
// widen this window a little to possibly include some low quality
// bases as well, but not larger than the window around the
// end of the consensus
int nConsPosStart;
int nConsPosEnd;
assert( bIntersect( nWindowToLookLeft,
nWindowToLookRight,
pFrag->nGetAlignStart(),
pFrag->nGetAlignEnd(),
nConsPosStart,
nConsPosEnd ) );
// find first/last x
bool bFoundX = false;
bool bFoundNotX = false;
int nXNotXJunction; // mark on not-X base
if ( nWhichEnd == nLeftGap ) {
int nState = nLookingForX;
// want to find where XXXXBBBBB junction is
bool bStillLooking = true;
for( int nConsPos = nConsPosStart; nConsPos <= nConsPosEnd && bStillLooking; ++nConsPos ) {
char c = pFrag->ntGetFragFromConsPos( nConsPos ).cGetBase();
if ( c == '*' ) continue;
if ( nState == nLookingForX ) {
if ( c == 'x' ) {
bFoundX = true;
nState = nLookingForNotX;
}
}
else if ( nState == nLookingForNotX ) {
if ( c != 'x' ) {
bFoundNotX = true;
nXNotXJunction = nConsPos;
bStillLooking = false;
}
}
} // for( int nConsPos
} // if ( nWhichEnd == nLeftGap ) {
else {
int nState = nLookingForNotX;
// want to find where (notX)XXXXX junction is
bool bStillLooking = true;
for( int nConsPos = nConsPosStart; nConsPos <= nConsPosEnd &&
bStillLooking; ++nConsPos ) {
char c = pFrag->ntGetFragFromConsPos( nConsPos ).cGetBase();
if ( c == '*' ) continue;
if ( nState == nLookingForNotX ) {
if ( c != 'x' ) {
bFoundNotX = true;
nState = nLookingForX;
}
}
else if ( nState == nLookingForX ) {
if ( c == 'x' ) {
bFoundX = true;
nXNotXJunction = nConsPos - 1;
bStillLooking = false;
}
}
} // for (int nConsPos ..
}
// for this read to be informative, there must be an X/notX junction
if ( bFoundX && bFoundNotX ) {
// add this read and its junction to the list
pVectorReadArray[ nNumberOfXNotXReads ].pLocFrag_ = pFrag;
pVectorReadArray[ nNumberOfXNotXReads ].nVectorInsertJunction_
= nXNotXJunction;
++nNumberOfXNotXReads;
}
} // for( int nRead...
qsort( pVectorReadArray,
nNumberOfXNotXReads,
sizeof( vectorReadType ),
( ( int(*) ( const void*, const void* ) ) nSortByXNotXJunction ) );
// look for a pair of reads in which one of the reads that both have
// junction close to each other and one of them is over 100 bases from
// the beginning
const int nLittleWindow = 6;
bool bFoundPairOfReads = false;
int nVectorInsertJunction = 0;
for( nRead = 0; nRead < nNumberOfXNotXReads - 1; ++nRead ) {
int nLocDiff =
pVectorReadArray[ nRead + 1 ].nVectorInsertJunction_
-
pVectorReadArray[ nRead ].nVectorInsertJunction_;
assert( nLocDiff >= 0 );
if ( nLocDiff <= nLittleWindow ) {
// found possible vector/insert junction
// Let's see if one of these reads has the junction
// not in the 1st hundred bases
// If the gap is the left end, then the read must be top strand
// (pointing in) and have the xxxxxAAAAAA be over base position 100.
// This is because bottom strand reads could have
// xxxxxxxAAAAAAA over base position 100 due to a short insert
// size and the sequence running into sequencing vector.
// If the gap is on the right end, then the read must be bottom
// (pointing in) and have the AAAAAAxxxxx after base 100.
// This is because top strand (pointing out) reads could have
// AAAAAAxxxxxx due to short insert size. I actually saw this.
// Phil had me remove the above condition and instead have
// a condition that there is no read that has all high quality
// non-X in a window around the junction.
// I want to check that the 2 reads are not from the same template.
// If they are, then the fact that 2 reads have XXXX's starting
// in the same place could just be due to sequencing (e.g., M13)
// vector. Phil agreed 3/99
if ( pVectorReadArray[ nRead ].pLocFrag_->soGetTemplateName() ==
pVectorReadArray[ nRead + 1 ].pLocFrag_->soGetTemplateName() )
continue;
if (
( pVectorReadArray[ nRead ].pLocFrag_->nOrientedUnpaddedFragPosFromConsPos( pVectorReadArray[ nRead ].nVectorInsertJunction_ ) > 100 )
||
( pVectorReadArray[ nRead + 1].pLocFrag_->nOrientedUnpaddedFragPosFromConsPos( pVectorReadArray[ nRead + 1 ].nVectorInsertJunction_ ) > 100 ) ) {
// I think we've found a vector/insert junction
bFoundPairOfReads = true;
nVectorInsertJunction = pVectorReadArray[ nRead ].nVectorInsertJunction_;
}
} // if ( nLocDiff <= 6 )
} // for( int nRead ...
// now we have a putative vector/insert junction
free( pVectorReadArray );
if ( !bFoundPairOfReads )
return( false );
// now do this additional check--make sure that all other reads
// are either:
// 1) low quality
// 2) have some X's pointing out of the contig if we look far enough
// That is, if a read is high quality at the junction, then it must
// not be all non-X's in a window.
// That window used to be only 6 bases. This failed in a case like this:
// BBBBBBBBBBBBBBBBBBB(contig) (consensus)
// XXXXXXXXXXXXXXXXXXBBBBBBBBBBBBBBBBBBBBBBBBBBB (read 1)
// xxxxxxxxxxxxxxxxxxxbbbbbbbbbbbbbbbbbb
// xxxxxxxxxxxxxxxxxxxbbbbbbbbbbbbbbbbbb
// read 1 was the high quality read.
// It's vector wasn't masked all the way to the vector/insert junction.
// So now I look 15 bases to see if the read is all high quality out of
// the contig. If I have no X's in that 15 base high quality window,
// then I assume not a clone end.
const int nWindowForLookingForHighQuality = 15;
nWindowToLookLeft = ( nWhichEnd == nLeftGap ) ?
( nVectorInsertJunction - nWindowForLookingForHighQuality ) :
( nVectorInsertJunction ) ;
nWindowToLookRight = ( nWhichEnd == nLeftGap ) ?
( nVectorInsertJunction ) :
( nVectorInsertJunction + nWindowForLookingForHighQuality );
for( nRead = 0; nRead < pContigView->nGetNumberOfFragments();
++nRead ) {
LocatedFragment* pLocFrag = pContigView->pLocatedFragmentGet( nRead );
if ( pLocFrag->bIsWholeReadLowQuality() ) continue;
int nHighQualityStart = pLocFrag->nGetHighQualityStart();
int nHighQualityEnd = pLocFrag->nGetHighQualityEnd();
int nIntersect1;
int nIntersect2;
if ( !bIntersect( nHighQualityStart,
nHighQualityEnd,
nWindowToLookLeft,
nWindowToLookRight,
nIntersect1,
nIntersect2
) )
continue;
if ( ( nIntersect1 != nWindowToLookLeft ) ||
( nIntersect2 != nWindowToLookRight ) )
continue;
// if reached here, we have a read whose high quality segment
// contains the window just outside the putative vector/insert
// junction.
// The condition that will show that this is *not* the
// cloning site is if all these bases are high quality not X's
bool bFoundX = false;
for( int nConsPos = nIntersect1; nConsPos <= nIntersect2 && !bFoundX; ++nConsPos )
if ( pLocFrag->ntGetFragFromConsPos( nConsPos ).cGetBase() == 'x' )
bFoundX = true;
if ( !bFoundX ) {
// oh, oh. There is a read with all non-X's outside the putative
// insert. So we must not have found the vector/insert junction.
return( false );
}
}// for( int nRead...
// so if passed all these tests, then it must be a cloning site
return( true );
}
void Contig :: createExperimentsWithFirstForwardsOrFirstReverses() {
double dMinQuality =
consedParameters::pGetConsedParameters()->dAutoFinishMinNumberOfErrorsFixedByAnExpToUse_;
for( int nTemplate = 0; nTemplate < aSubcloneTemplates_.entries();
++nTemplate ) {
subcloneTTemplate* pSub = aSubcloneTemplates_[ nTemplate ];
bool bTopNotBottomStrand;
int nReadType;
if ( !pSub->bOKToCallDeNovoUniversalPrimerRead( bTopNotBottomStrand,
nReadType,
this ) )
continue;
if ( pSub->bChosenPriorToMakingListOfExperimentCandidates_ )
continue;
int nFirstBaseOfReadUnpadded;
int nUnpaddedLeft;
int nUnpaddedRight;
if ( bTopNotBottomStrand ) {
// case like this:
// xxxx<--------xxxx( existing read)
// -------->xxxx( reverse to be called)
// or this:
// --------> <--------xxxx(existing read)
// <-------insert size-------->
nFirstBaseOfReadUnpadded = pSub->nUnpaddedStart_ -
pCP->nAutoFinishNumberOfVectorBasesAtBeginningOfAUniveralPrimerRead_;
nUnpaddedLeft = pSub->nUnpaddedStart_;
nUnpaddedRight = nFirstBaseOfReadUnpadded +
ConsEd::pGetAssembly()->pModelRead_->length() - 1;
}
else {
// ----------> (existing univ primer read )
// <-------------------- (mate to suggest)
// ^ ( nConsPosStart_ )
// (first base of read)^
// case like this:
// xxxx----------->xxxxx (existing read)
// <------- (candidate reverse)
// or case like this:
// xxxx-----------> (existing read)
// <--------- (candidate reverse)
// <-------insert size------>
nFirstBaseOfReadUnpadded = pSub->nUnpaddedEnd_ +
pCP->nAutoFinishNumberOfVectorBasesAtBeginningOfAUniveralPrimerRead_;
nUnpaddedLeft = nFirstBaseOfReadUnpadded -
ConsEd::pGetAssembly()->pModelRead_->length() + 1;
nUnpaddedRight = pSub->nUnpaddedEnd_;
}
// case in which very short model read and long parameter
// pCP->nAutoFinishNumberOfVectorBasesAtBeginningOfAUniveralPrimerRead_;
if ( nUnpaddedLeft > nUnpaddedRight ) continue;
// is it possible to have a low quality region on the extended contig
// but not a high quality region? Yes:
// (contig)
// ------------------------------
// ------> (existing read)
// <llllllHHHHHHHHH---
// (reverse)
double dErrorsFixed = dGetInitialMaxErrorsFixed2(
nFirstBaseOfReadUnpadded,
bTopNotBottomStrand );
// see if this experiment will solve more than the minimum # of errors
if ( dErrorsFixed
< dMinQuality ) {
continue;
}
// if passes that test, then find errors exactly
dErrorsFixed = dGetErrorsFixedForSubcloneRead(
&pSub,
1, // only 1 subclone
bTopNotBottomStrand,
nFirstBaseOfReadUnpadded,
true, // exclude vector bases
false, // bChooseExperiment
false, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
if ( dErrorsFixed < dMinQuality ) {
continue;
}
// We have a possible experiment
experiment* pExp = new experiment();
pExp->dErrorsFixed_ = dErrorsFixed;
readTypes ucChemistryAndStrand = ( bTopNotBottomStrand ?
TOP_STRAND_TERMINATOR_READS :
BOTTOM_STRAND_TERMINATOR_READS );
pExp->ucStrandAndChemistry_ = ucChemistryAndStrand;
pExp->nReadType_ = nReadType;
pExp->bCloneTemplateNotSubcloneTemplate_ = false;
pExp->fCost_ =
pCP->dAutoFinishCostOfDeNovoUniversalPrimerSubcloneReaction_;
pExp->pCustomPrimer_ = 0;
// must figure out the errors on just the consensus
pExp->pContig_ = this;
pExp->pSubcloneTemplate_ = pSub;
pExp->nFirstBaseOfReadUnpadded_ = nFirstBaseOfReadUnpadded;
pExp->nUnpaddedLeft_ = nUnpaddedLeft;
pExp->nUnpaddedRight_ = nUnpaddedRight;
// case in which the model read is shorter than the
// parameter nAutoFinishNumberOfVectorBasesAtBeginningOfAUniveralPrimerRead_
if ( pExp->nUnpaddedLeft_ > pExp->nUnpaddedRight_ ) continue;
pExp->bTopStrandNotBottomStrand_ = bTopNotBottomStrand;
pExp->nState_ = nAvailableExperiment;
pExp->bDeNovoUniversalPrimerRead_ = true;
setVariousErrorData( pExp );
pCandidateExperiments_->insert( pExp );
} // for( int nTemplate = 0; nTemplate < aSubcloneTemplates_.entries();
}
// Note: this can only be used before pUnpaddedErrorProbabilities_
// starts being modified
double Contig :: dGetOriginalConsensusErrors( const int nUnpaddedLeft,
const int nUnpaddedRight ) {
int nUnpaddedLeftOnConsensus;
int nUnpaddedRightOnConsensus;
if ( !bIntersect( nUnpaddedLeft,
nUnpaddedRight,
nGetUnpaddedStartIndex(),
nGetUnpaddedEndIndex(),
nUnpaddedLeftOnConsensus,
nUnpaddedRightOnConsensus )) {
// the region is completely off the consensus
return( 0.0 );
}
double dOriginalErrors = 0.0;
for( int nUnpadded = nUnpaddedLeftOnConsensus;
nUnpadded <= nUnpaddedRightOnConsensus;
++nUnpadded ) {
dOriginalErrors += (*pOriginalUnpaddedErrorProbabilities_)[ nUnpadded ];
}
return( dOriginalErrors );
}
double Contig :: dGetInitialMaxErrorsFixed2( const int nFirstBaseOfReadUnpadded,
const bool bTopStrandNotBottomStrand ) {
int nUnpaddedLeft;
int nUnpaddedRight;
if ( bTopStrandNotBottomStrand ) {
nUnpaddedLeft = nFirstBaseOfReadUnpadded;
nUnpaddedRight = nUnpaddedLeft +
ConsEd::pGetAssembly()->pModelRead_->length() - 1;
}
else {
nUnpaddedRight = nFirstBaseOfReadUnpadded;
nUnpaddedLeft = nUnpaddedRight -
ConsEd::pGetAssembly()->pModelRead_->length() + 1;
}
return( dGetInitialMaxErrorsFixed( nUnpaddedLeft, nUnpaddedRight ) );
}
double Contig :: dGetInitialMaxErrorsFixed( const int nUnpaddedLeft,
const int nUnpaddedRight ) {
int nRightEndOfRead = nUnpaddedRight;
int nJustToLeftOfLeftEndOfRead = nUnpaddedLeft - 1;
// I think this method may not take into account the error probability
// of the leftmost base if the read is hanging off the left end
// (DG 1/00). Yes, it does, since the pCumUnpaddedErrors starts
// one base to the left of the leftmost extended consensus error probability
// array. (DG 1/00).
int nIntersectLeft;
int nIntersectRight;
if ( !bIntersect( nJustToLeftOfLeftEndOfRead,
nRightEndOfRead,
pCumUnpaddedErrors_->nGetStartIndex(),
pCumUnpaddedErrors_->nGetEndIndex(),
nIntersectLeft,
nIntersectRight ) )
return( 0.0 );
return( (*pCumUnpaddedErrors_)[ nIntersectRight ] -
(*pCumUnpaddedErrors_)[ nIntersectLeft ] );
}
void Contig :: dumpContigErrorProbabilities( const int nFileID ) {
FileName filDump = "error_probabilities." + RWCString( (long) nFileID );
fprintf( pAO, "%s\n", filDump.data() );
FILE* pDump = fopen( filDump.data(), "w" );
if ( !pDump ) {
PANIC_OST( ost ) << "In dumpContigErrorProbablities: cannot open " <<
filDump << ends;
throw SysRequestFailed( ost.str() );
}
for( int nUnpadded = 2851;
nUnpadded <= 2860;
++nUnpadded ) {
double dErrorProbability = (*pUnpaddedErrorProbabilities_)[ nUnpadded ];
double dQuality = -10.0 * log10( dErrorProbability );
fprintf( pDump, "%d %.2f %.6e\n", nUnpadded, dQuality,
dErrorProbability );
}
fclose( pDump );
}
bool Contig :: bThereIsAnotherSubcloneCustomPrimerReadTooCloseOnTheSameStrand(
experiment* pExperimentCandidate ) {
for( int nExp = 0; nExp < pChosenExperiments_->length(); ++nExp ) {
experiment* pExp = (*pChosenExperiments_)[ nExp ];
if ( ABS( pExp->nFirstBaseOfReadUnpadded_ - pExperimentCandidate->nFirstBaseOfReadUnpadded_ )
< pCP->nAutoFinishDoNotAllowSubcloneCustomPrimerReadsCloserThanThisManyBases_ ) {
if (
( pExp->bTopStrandNotBottomStrand_ == pExperimentCandidate->bTopStrandNotBottomStrand_ ) &&
( pExp->nReadType_ == nWalk ) &&
!pExp->bCloneTemplateNotSubcloneTemplate_
) {
return( true );
}
}
}
return( false );
}
bool Contig :: bThereIsAnotherWholeCloneCustomPrimerReadTooCloseOnTheSameStrand(
experiment* pExperimentCandidate ) {
for( int nExp = 0; nExp < pChosenExperiments_->length(); ++nExp ) {
experiment* pExp = (*pChosenExperiments_)[ nExp ];
if ( ABS( pExp->nFirstBaseOfReadUnpadded_ -
pExperimentCandidate->nFirstBaseOfReadUnpadded_ )
<
pCP->nAutoFinishDoNotAllowWholeCloneCustomPrimerReadsCloserThanThisManyBases_ ) {
if (
( pExp->bTopStrandNotBottomStrand_ == pExperimentCandidate->bTopStrandNotBottomStrand_ ) &&
( pExp->nReadType_ == nWalk ) &&
pExp->bCloneTemplateNotSubcloneTemplate_
) {
return( true );
}
}
}
return( false );
}
void Contig :: possiblyUpdateContigHighQualitySegment( experiment* pChosenExp ) {
if ( pChosenExp->nUnpaddedLeft_ < nUpdatedUnpaddedHighQualitySegmentStart_ )
nUpdatedUnpaddedHighQualitySegmentStart_ =
pChosenExp->nUnpaddedLeft_;
if ( pChosenExp->nUnpaddedRight_ > nUpdatedUnpaddedHighQualitySegmentEnd_ )
nUpdatedUnpaddedHighQualitySegmentEnd_ =
pChosenExp->nUnpaddedRight_;
}
class endpointsOfSubcloneOfRead {
public:
endpointsOfSubcloneOfRead( const int nUnpaddedLeftEndpoint,
const int nUnpaddedRightEndpoint,
subcloneTTemplate* pSub ) :
nUnpaddedLeftEndpoint_( nUnpaddedLeftEndpoint ),
nUnpaddedRightEndpoint_( nUnpaddedRightEndpoint ),
pSub_( pSub ) {}
// default ctor to keep RWTValOrderedVector happy
endpointsOfSubcloneOfRead() {}
// also needed to keep RWTValOrderedVector happy on hp-ux
bool operator==( const endpointsOfSubcloneOfRead& end ) {
if ( nUnpaddedLeftEndpoint_ == end.nUnpaddedLeftEndpoint_ &&
nUnpaddedRightEndpoint_ == end.nUnpaddedRightEndpoint_ &&
pSub_ == end.pSub_ )
return( true );
else
return( false );
}
public:
int nUnpaddedLeftEndpoint_;
int nUnpaddedRightEndpoint_;
subcloneTTemplate* pSub_;
};
static int nCompareLeftEndpoints( endpointsOfSubcloneOfRead* pEndPoints1,
endpointsOfSubcloneOfRead* pEndPoints2 ) {
if ( pEndPoints1->nUnpaddedLeftEndpoint_ < pEndPoints2->nUnpaddedLeftEndpoint_ )
return( -1 );
else if ( pEndPoints1->nUnpaddedLeftEndpoint_ == pEndPoints2->nUnpaddedLeftEndpoint_ )
return( 0 );
else
return( 1 );
}
static int nCompareRightEndpoints( endpointsOfSubcloneOfRead* pEndPoints1,
endpointsOfSubcloneOfRead* pEndPoints2 ) {
if ( pEndPoints1->nUnpaddedRightEndpoint_ < pEndPoints2->nUnpaddedRightEndpoint_ )
return( -1 );
else if ( pEndPoints1->nUnpaddedRightEndpoint_ == pEndPoints2->nUnpaddedRightEndpoint_ )
return( 0 );
else
return( 1 );
}
double Contig :: dGetErrorsFixedForCustomPrimerSubcloneRead(
primerType* pPrimer,
const bool bTopStrandNotBottomStrand,
const int nFirstBaseOfReadUnpadded,
const bool bChooseExperiment,
const bool bJustConsensusNotGaps,
const bool bFindWindowOfMostErrorsFixed ) {
int nNumberOfTemplatesToUse = 0;
for( int nTemplate = 0;
nTemplate < pCP->nAutoFinishHowManyTemplatesYouIntendToUseForCustomPrimerSubcloneReactions_;
++nTemplate ) {
if ( pPrimer->pSubcloneTemplate_[ nTemplate ] )
++nNumberOfTemplatesToUse;
}
return( dGetErrorsFixedForSubcloneRead(
pPrimer->pSubcloneTemplate_,
nNumberOfTemplatesToUse,
bTopStrandNotBottomStrand,
nFirstBaseOfReadUnpadded,
false, // bExcludeVectorBaseAtBeginningOfRead
bChooseExperiment, // bChooseExperiment
bJustConsensusNotGaps,
bFindWindowOfMostErrorsFixed ) );
}
static int nWindowOfMostErrorsFixedUnpaddedRight;
static double dErrorsFixedInWindowOfMostErrorsFixed;
static double dErrorsFixedInCurrentWindow;
static double dErrorsFixed[10];
static int nPointerToLastAddedBase;
static int n10thBaseUnpadded;
static float fProductDistributionSameStrand[ nMAX_QUALITY_LEVEL2 + 1 ];
static float fProductDistributionOtherStrand[ nMAX_QUALITY_LEVEL2 + 1 ];
bool bProductDistributionSameStrandAllOnes = true;
bool bProductDistributionOtherStrandAllOnes = true;
double Contig :: dGetErrorsFixedForSubcloneRead(
subcloneTTemplate** ppSubcloneTemplateArray,
const int nNumberOfTemplates,
const bool bTopStrandNotBottomStrand,
const int nFirstBaseOfReadUnpadded,
const bool bExcludeVectorBasesAtBeginningOfRead,
const bool bChooseExperiment,
const bool bJustConsensusNotGaps,
const bool bFindWindowOfMostErrorsFixed ) {
// fixing bug in which, if the entire read is off the extended consensus,
// such as when finding all reads that might close a gap, the
// window of most errors fixed was garbage
if ( bFindWindowOfMostErrorsFixed ) {
// these will be changed if we are able to calculate them
nWindowOfMostErrorsFixedUnpaddedRight = 0;
dErrorsFixedInWindowOfMostErrorsFixed = 0.0;
}
// nLeftmostBaseOfRead and nRightmostBaseOfRead and the
// extents of the read(s) created with this primer if the read
// were not to be limited by the primers.
int nLeftmostBaseOfRead;
int nRightmostBaseOfRead;
if ( bTopStrandNotBottomStrand ) {
if ( bExcludeVectorBasesAtBeginningOfRead )
nLeftmostBaseOfRead = nFirstBaseOfReadUnpadded +
pCP->nAutoFinishNumberOfVectorBasesAtBeginningOfAUniveralPrimerRead_;
else
nLeftmostBaseOfRead = nFirstBaseOfReadUnpadded;
nRightmostBaseOfRead = nFirstBaseOfReadUnpadded +
ConsEd::pGetAssembly()->pModelRead_->length()
- 1;
}
else {
nLeftmostBaseOfRead = nFirstBaseOfReadUnpadded -
ConsEd::pGetAssembly()->pModelRead_->length()
+ 1;
if ( bExcludeVectorBasesAtBeginningOfRead )
nRightmostBaseOfRead = nFirstBaseOfReadUnpadded -
pCP->nAutoFinishNumberOfVectorBasesAtBeginningOfAUniveralPrimerRead_;
else
nRightmostBaseOfRead = nFirstBaseOfReadUnpadded;
}
// this would happen if this project had short reads--shorter than
// the beginning vector on universal primer reads
if ( ! ( nLeftmostBaseOfRead <= nRightmostBaseOfRead ) )
return( 0.0 );
// intersect this with the extended consensus--we are only interested
// in the portion of the read on the extended consensus since no credit
// is given for any part of the read that sticks out beyond this.
// (usually this will be all of it)
int nOverlapReadLeft;
int nOverlapReadRight;
if ( ! bIntersect( pUnpaddedErrorProbabilities_->nGetStartIndex(),
pUnpaddedErrorProbabilities_->nGetEndIndex(),
nLeftmostBaseOfRead,
nRightmostBaseOfRead,
nOverlapReadLeft,
nOverlapReadRight ) )
return( 0.0 );
if ( bJustConsensusNotGaps ) {
if ( ! bIntersect( nOverlapReadLeft,
nOverlapReadRight,
nGetUnpaddedStartIndex(),
nGetUnpaddedEndIndex(),
nOverlapReadLeft,
nOverlapReadRight ) ) {
return( 0.0 );
}
}
RWTValOrderedVector<endpointsOfSubcloneOfRead> aEndpointsOfTransitionsBetweenTemplateCoverage(
(size_t) nNumberOfTemplates );
for( int nTemplate = 0;
nTemplate < nNumberOfTemplates;
++nTemplate ) {
subcloneTTemplate* pSub = ppSubcloneTemplateArray[ nTemplate ];
int nTemplateReadLeft;
int nTemplateReadRight;
if ( bIntersect( nOverlapReadLeft,
nOverlapReadRight,
pSub->nUnpaddedStart_,
pSub->nUnpaddedEnd_,
nTemplateReadLeft,
nTemplateReadRight ) ) {
aEndpointsOfTransitionsBetweenTemplateCoverage.insert(
endpointsOfSubcloneOfRead( nTemplateReadLeft,
nTemplateReadRight,
pSub )
);
}
} // for( int nTemplate = 0;
// sort this array
void* pArray =
(void*) aEndpointsOfTransitionsBetweenTemplateCoverage.pArray_;
size_t nNumberOfElements =
aEndpointsOfTransitionsBetweenTemplateCoverage.nCurrentLength_;
size_t nSizeOfAnElement = sizeof( endpointsOfSubcloneOfRead );
if ( bTopStrandNotBottomStrand )
qsort( pArray, nNumberOfElements, nSizeOfAnElement,
( (int(*) ( const void*, const void*) ) nCompareRightEndpoints ) );
else
qsort( pArray, nNumberOfElements, nSizeOfAnElement,
( (int(*) ( const void*, const void*) ) nCompareLeftEndpoints ) );
// now for each base, see how many errors will be fixed
// Note that each read here will be the same strand and chemistry
// now for each *unpadded* consensus base, we will compute how many
// errors are fixed. A read gets no credit for fixing anything
// where there is a pad in the consensus.
readTypes ucNewReadsStrandAndChemistry =
( bTopStrandNotBottomStrand ?
TOP_STRAND_TERMINATOR_READS :
BOTTOM_STRAND_TERMINATOR_READS );
readTypes ucOtherStrandAndChemistry =
( !bTopStrandNotBottomStrand ?
TOP_STRAND_TERMINATOR_READS :
BOTTOM_STRAND_TERMINATOR_READS );
if ( bFindWindowOfMostErrorsFixed ) {
dErrorsFixedInWindowOfMostErrorsFixed = 0.0;
dErrorsFixedInCurrentWindow = 0.0;
// find 2 critical points:
// ---------------------------------------
// [ ]
// ^ 10th base
// ^ last base
// before the 10th base, we just accumulate errors fixed
// at the 10th base, we set the first best window
// after the 10th base, when we add a base, we take one away
n10thBaseUnpadded =
aEndpointsOfTransitionsBetweenTemplateCoverage[ 0 ].nUnpaddedLeftEndpoint_ + 9;
nPointerToLastAddedBase = -1; // prepare to increment to 0
}
double dTotalErrorsFixed = 0.0;
for( int nTemplateTransitions = 0;
nTemplateTransitions < aEndpointsOfTransitionsBetweenTemplateCoverage.length();
++nTemplateTransitions ) {
int nUnpaddedLeft;
int nUnpaddedRight;
int nNumberOfTemplatesInThisInterval;
// top strand reads have templates like this:
// ------ (template 1)
// ---------- (template 2)
// ------------ (template 3)
// <----><--><>
// A B C
// So when the nNumberOfTemplatesInThisInterval is 2 (region B), we use
// template 2 and 3
// bottom strand reads have templates like this:
// -------------------- template 1
// ---------------- template 2
// ----------- template 3
// <--><---><--------->
// A B C
// so when nNumberOfTemplatesInThisInterval is 2 (region B), we use templates
// 1 and 2
// nUnpaddedLeft and nUnpaddedRight (below) are the boundaries
// of the regions A, B, C in turn
if ( bTopStrandNotBottomStrand ) {
nNumberOfTemplatesInThisInterval =
aEndpointsOfTransitionsBetweenTemplateCoverage.length()
- nTemplateTransitions;
if ( nTemplateTransitions == 0 )
nUnpaddedLeft =
aEndpointsOfTransitionsBetweenTemplateCoverage[ 0 ].nUnpaddedLeftEndpoint_;
else
nUnpaddedLeft =
aEndpointsOfTransitionsBetweenTemplateCoverage[ nTemplateTransitions - 1 ].nUnpaddedRightEndpoint_
+ 1;
nUnpaddedRight =
aEndpointsOfTransitionsBetweenTemplateCoverage[ nTemplateTransitions ].nUnpaddedRightEndpoint_;
}
else {
// bottom strand
nNumberOfTemplatesInThisInterval = nTemplateTransitions + 1;
nUnpaddedLeft = aEndpointsOfTransitionsBetweenTemplateCoverage[ nTemplateTransitions ].nUnpaddedLeftEndpoint_;
if ( nTemplateTransitions ==
( aEndpointsOfTransitionsBetweenTemplateCoverage.length() - 1 ) )
nUnpaddedRight =
aEndpointsOfTransitionsBetweenTemplateCoverage[ nTemplateTransitions ].nUnpaddedRightEndpoint_;
else
nUnpaddedRight =
aEndpointsOfTransitionsBetweenTemplateCoverage[ nTemplateTransitions + 1 ].nUnpaddedLeftEndpoint_ - 1;
}
MBTValOrderedOffsetVector<unsigned char>* pModelReadWithThisNumberOfTemplates =
(*(ConsEd::pGetAssembly()->pArrayOfModelReadsWithSeveralTemplates_))[ nNumberOfTemplatesInThisInterval ];
// set up for walking along read
int nIncrement;
int nSeqPosition;
if ( bTopStrandNotBottomStrand ) {
nSeqPosition = nUnpaddedLeft - nFirstBaseOfReadUnpadded + 1;
nIncrement = 1;
// --------------> read
// ----------- extended consensus
// ^ nOverlapLeft (start here and nSeqPosition then increases
//
}
else {
// <----------------- read
// --------------- extended consensus
// ^ nFirstBaseOfReadUnpadded
// ^ nOverlapRight
// ^ nOverlapLeft (start here and nSeqPosition then decreases
nSeqPosition = nFirstBaseOfReadUnpadded - nUnpaddedLeft + 1;
nIncrement = -1;
}
// walk along read
double dErrorsFixedThisBase2;
for( int nUnpadded = nUnpaddedLeft; nUnpadded <= nUnpaddedRight;
++nUnpadded, nSeqPosition += nIncrement ) {
unsigned char ucMedianMaxQualityOfNewReads;
RWTPtrOrderedVector<distribution>* pArrayOfDistributions
= (*pDistributionsOfChosenExperimentsArray_)[ nUnpadded ];
if ( !pArrayOfDistributions ) {
// case of no experiments chosen at this consensus location
// This is no longer true since there may be a minilibrary
// called which will make pUnpaddedErrorProbabilities zero
// where it covers the consensus.
// assert( (*pUnpaddedErrorProbabilities_)[ nUnpadded ] ==
// (*pOriginalUnpaddedErrorProbabilities_)[ nUnpadded ] );
ucMedianMaxQualityOfNewReads =
(*pModelReadWithThisNumberOfTemplates)[ nSeqPosition ];
dErrorsFixedThisBase2 =
dGetErrorsFixedAtBase(
(*pUnpaddedErrorProbabilities_)[ nUnpadded ],
dErrorRateFromQuality[ ucMedianMaxQualityOfNewReads ],
ucNewReadsStrandAndChemistry,
(*pUnpaddedReadTypesCoveringEachBase_)[ nUnpadded ] );
dTotalErrorsFixed += dErrorsFixedThisBase2;
// bChooseExperiment && nUnpadded >= -28 && nUnpadded <= -26
// if ( bFoundDebug && bChooseExperiment ) {
// fprintf( pAO, "no exp yet, nUnpadded = %d and total errors fixed = %g median quality = %d seq pos = %d\n",
// nUnpadded, dErrorsFixedAtBase2, (int) ucMedianMaxQualityOfNewReads,
// nSeqPosition );
// }
if ( bChooseExperiment ) {
// update the consensus error probabilities
(*pUnpaddedErrorProbabilities_)[ nUnpadded ] -=
dErrorsFixedThisBase2;
(*pUnpaddedReadTypesCoveringEachBase_)[ nUnpadded ] |=
ucNewReadsStrandAndChemistry;
// must put an array of distributions on the consensus
RWTPtrOrderedVector<distribution>* pArrayOfDistributions =
new RWTPtrOrderedVector<distribution>( (size_t) 3 );
(*pDistributionsOfChosenExperimentsArray_)[ nUnpadded ] =
pArrayOfDistributions;
// now make a distribution for each template add it to
// the array of distributions
int nTemplateStart;
int nTemplateEnd;
if ( bTopStrandNotBottomStrand ) {
nTemplateStart =
aEndpointsOfTransitionsBetweenTemplateCoverage.length() -
nNumberOfTemplatesInThisInterval;
nTemplateEnd =
aEndpointsOfTransitionsBetweenTemplateCoverage.length() - 1;
}
else {
nTemplateStart = 0;
nTemplateEnd =
nNumberOfTemplatesInThisInterval - 1;
}
float* pModelReadDistribution =
(*ConsEd::pGetAssembly()->pModelReadDistributions_)[ nSeqPosition ];
for( int nTemplate = nTemplateStart; nTemplate <= nTemplateEnd;
++nTemplate ) {
distribution* pDis = new distribution();
// fixed bug 5/8/2000
pDis->pSub_ = aEndpointsOfTransitionsBetweenTemplateCoverage[ nTemplate ].pSub_;
pDis->bTopStrandNotBottomStrand_ = bTopStrandNotBottomStrand;
for( int nQuality = 0; nQuality <= nMAX_QUALITY_LEVEL2;
++nQuality ) {
pDis->fCumulativeQuality_[ nQuality ] =
pModelReadDistribution[ nQuality ];
}
pArrayOfDistributions->insert( pDis );
}
}
}
else {
// all hell and fury
float* pModelReadDistribution =
(*ConsEd::pGetAssembly()->pModelReadDistributions_)[ nSeqPosition ];
// top strand reads have templates like this:
// ------ (template 1)
// ---------- (template 2)
// ------------ (template 3)
// <----><--><>
// A B C
// So when the nNumberOfTemplatesInThisInterval is 2 (region B), we use
// template 2 and 3
// bottom strand reads have templates like this:
// -------------------- template 1
// ---------------- template 2
// ----------- template 3
// <--><---><--------->
// A B C
// so when nNumberOfTemplatesInThisInterval is 2 (region B), we use templates
// 1 and 2
int nTemplateStart;
int nTemplateEnd;
if ( bTopStrandNotBottomStrand ) {
nTemplateStart =
aEndpointsOfTransitionsBetweenTemplateCoverage.length() -
nNumberOfTemplatesInThisInterval;
nTemplateEnd =
aEndpointsOfTransitionsBetweenTemplateCoverage.length() - 1;
}
else {
nTemplateStart = 0;
nTemplateEnd =
nNumberOfTemplatesInThisInterval - 1;
}
// make a little array of bool that indicates which of the
// consensus templates we have already used
RWTValVector<bool> aChosenTemplatesCombined(
pArrayOfDistributions->length(),
false
);
bProductDistributionSameStrandAllOnes = true;
bProductDistributionOtherStrandAllOnes = true;
for( int nTemplate = nTemplateStart;
nTemplate <= nTemplateEnd; ++nTemplate ) {
// found bug 5/8/2000
subcloneTTemplate* pSubOfConsideredRead =
aEndpointsOfTransitionsBetweenTemplateCoverage[ nTemplate ].pSub_;
// check if this template is found in one of the chosen
// reads at this base position
distribution* pDistributionWithSameTemplate = NULL;
int nDistributionWithSameTemplate;
for( int nDis = 0; nDis < pArrayOfDistributions->length();
++nDis ) {
if ( pSubOfConsideredRead ==
(*pArrayOfDistributions)[ nDis ]->pSub_ ) {
if ( bTopStrandNotBottomStrand ==
(*pArrayOfDistributions)[ nDis ]->bTopStrandNotBottomStrand_ ) {
pDistributionWithSameTemplate =
(*pArrayOfDistributions)[ nDis ];
nDistributionWithSameTemplate = nDis;
break;
}
}
}
if (!pDistributionWithSameTemplate ) {
// easy case--just multiply this distribution
if ( bProductDistributionSameStrandAllOnes ) {
assignDistributions( fProductDistributionSameStrand,
pModelReadDistribution );
bProductDistributionSameStrandAllOnes = false;
}
else
multiplyDistributions( fProductDistributionSameStrand,
pModelReadDistribution );
if ( bChooseExperiment ) {
// must add the distribution for this new template
distribution* pNewDistribution = new distribution();
pNewDistribution->pSub_ = pSubOfConsideredRead;
pNewDistribution->bTopStrandNotBottomStrand_ =
bTopStrandNotBottomStrand;
for( int nQuality = 0; nQuality <= nMAX_QUALITY_LEVEL2;
++nQuality )
pNewDistribution->fCumulativeQuality_[ nQuality ] =
pModelReadDistribution[ nQuality ];
pArrayOfDistributions->insert( pNewDistribution );
}
}
else {
// case in which there is a chosen read with the same
// template. In this case, find the max of the two
// distributions, and then multiply it.
aChosenTemplatesCombined[ nDistributionWithSameTemplate ]
= true;
if ( bProductDistributionSameStrandAllOnes ) {
makeDistributionAllOnes( fProductDistributionSameStrand );
bProductDistributionSameStrandAllOnes = false;
}
multiplyDistributionWithMinimum(
fProductDistributionSameStrand,
pModelReadDistribution,
pDistributionWithSameTemplate->fCumulativeQuality_ );
if ( bChooseExperiment ) {
for( int nQuality = 0; nQuality <= nMAX_QUALITY_LEVEL2;
++nQuality ) {
pDistributionWithSameTemplate->fCumulativeQuality_[ nQuality ] =
MIN( pModelReadDistribution[ nQuality ],
pDistributionWithSameTemplate->fCumulativeQuality_[ nQuality ] );
}
} // if ( bChooseExperiment ) {...
}
} // for( int nTemplate = nTemplateStart;
// We have figured out the product distribution of all
// templates in the considered read. We need now to include
// any distributions on chosen reads that have have templates
// distinct from any on the considered read.
for( int nChosenTemplate = 0;
nChosenTemplate < aChosenTemplatesCombined.length();
++nChosenTemplate ) {
if ( !aChosenTemplatesCombined[ nChosenTemplate ] ) {
float* pCumulativeQualityArray =
(*pArrayOfDistributions)[ nChosenTemplate ]->fCumulativeQuality_;
if ( (*pArrayOfDistributions)[ nChosenTemplate ]->bTopStrandNotBottomStrand_
== bTopStrandNotBottomStrand ) {
if ( bProductDistributionSameStrandAllOnes ) {
assignDistributions( fProductDistributionSameStrand,
pCumulativeQualityArray );
bProductDistributionSameStrandAllOnes = false;
}
else
multiplyDistributions2( fProductDistributionSameStrand,
pCumulativeQualityArray );
}
else {
if ( bProductDistributionOtherStrandAllOnes ) {
assignDistributions( fProductDistributionOtherStrand,
pCumulativeQualityArray );
bProductDistributionOtherStrandAllOnes = false;
}
else
multiplyDistributions3( fProductDistributionOtherStrand,
pCumulativeQualityArray );
}
}
} // for( int nChosenTemplate = 0; ...
// have figured out the product distribution at this
// base position for the considered read and the chosen reads
// of the same strand and of the other strand.
// figure out the median of the distributions of each strand
// and then combine these medians with the original
// consensus quality to get the new consensus quality.
unsigned char ucMedianQualitySameStrand;
if ( bProductDistributionSameStrandAllOnes )
ucMedianQualitySameStrand = 0;
else {
for( int nQuality = 0; nQuality <= nMAX_QUALITY_LEVEL2;
++nQuality ) {
if ( fProductDistributionSameStrand[ nQuality ] >= 0.5 ) {
ucMedianQualitySameStrand = nQuality;
break;
}
}
}
// there should be some way of speeding this up when there
// is *no* reads on the opposite strand. I don't know how
// often this will be the case--perhaps more than half the time?
// Sped up.
unsigned char ucMedianQualityOtherStrand;
if ( bProductDistributionOtherStrandAllOnes )
ucMedianQualityOtherStrand = 0;
else {
for( int nQuality = 0; nQuality <= nMAX_QUALITY_LEVEL2;
++nQuality ) {
if ( fProductDistributionOtherStrand[ nQuality ] >= 0.5 ) {
ucMedianQualityOtherStrand = nQuality;
break;
}
}
}
double dErrorsFixedByOtherStrand =
dGetErrorsFixedAtBase(
(*pOriginalUnpaddedErrorProbabilities_)[ nUnpadded ],
dErrorRateFromQuality[ ucMedianQualityOtherStrand ],
ucOtherStrandAndChemistry,
(*pOriginalUnpaddedReadTypesCoveringEachBase_)[ nUnpadded ] );
double dErrorsFixedBySameStrand =
dGetErrorsFixedAtBase(
(*pOriginalUnpaddedErrorProbabilities_)[ nUnpadded ] - dErrorsFixedByOtherStrand,
dErrorRateFromQuality[ ucMedianQualitySameStrand ],
ucNewReadsStrandAndChemistry,
(*pOriginalUnpaddedReadTypesCoveringEachBase_)[ nUnpadded ] );
double dErrorsFixedThisBase = dErrorsFixedBySameStrand +
dErrorsFixedByOtherStrand;
// this represents the total errors fixed at this base by all
// reads chosen plus the read currently considered.
// But we are only interested in the *additional* errors
// fixed by this read, not counting the errors fixed by
// previous reads.
double dAdditionalErrorsFixed = dErrorsFixedThisBase2 =
dErrorsFixedThisBase +
(*pUnpaddedErrorProbabilities_)[ nUnpadded ]
-
(*pOriginalUnpaddedErrorProbabilities_)[ nUnpadded ];
// note that (*pOriginalErrorProbabilities_)[ nUnpadded ] -
// (*pUnpaddedErrorProbabilities_)[ nUnpadded ]
// is the errors fixed by already chosen experiments
dTotalErrorsFixed += dAdditionalErrorsFixed;
// if ( bChooseExperiment && nUnpadded >= -28 && nUnpadded <= -26 ) {
// if ( bFoundDebug && bChooseExperiment ) {
// fprintf( pAO, "choosing additional experiment, nUnpadded = %d and total errors fixed = %g additional errors fixed = %g seq pos = %d dTotalErrors fixed = %g\n",
// nUnpadded, dErrorsFixedThisBase,
// dAdditionalErrorsFixed,
// nSeqPosition,
// dTotalErrorsFixed );
// if ( nUnpadded == 180 ) {
// if ( !bProductDistributionSameStrandAllOnes ) {
// fprintf( pAO, "product dist same strand: " );
// dumpDistribution( fProductDistributionSameStrand );
// }
// else
// fprintf( pAO, "product dist same strand all ones\n" );
// if ( !bProductDistributionOtherStrandAllOnes ) {
// fprintf( pAO, "product dist other strand: " );
// dumpDistribution( fProductDistributionOtherStrand );
// }
// else
// fprintf( pAO, "product dist other strand all ones\n" );
// fprintf( pAO, "model read distribution: " );
// dumpDistribution( pModelReadDistribution );
// }
// }
if ( bChooseExperiment ) {
(*pUnpaddedErrorProbabilities_)[ nUnpadded ] -=
dAdditionalErrorsFixed;
(*pUnpaddedReadTypesCoveringEachBase_)[ nUnpadded ] |=
ucNewReadsStrandAndChemistry;
}
} // if ( !pArrayOfDistributions ) else
if ( bFindWindowOfMostErrorsFixed ) {
if ( nUnpadded <= n10thBaseUnpadded ) {
++nPointerToLastAddedBase;
dErrorsFixed[ nPointerToLastAddedBase ] = dErrorsFixedThisBase2;
dErrorsFixedInCurrentWindow +=
dErrorsFixedThisBase2;
if ( nUnpadded == n10thBaseUnpadded ) {
dErrorsFixedInWindowOfMostErrorsFixed =
dErrorsFixedInCurrentWindow;
nWindowOfMostErrorsFixedUnpaddedRight = nUnpadded;
}
}
else {
++nPointerToLastAddedBase;
if ( nPointerToLastAddedBase > 9 )
nPointerToLastAddedBase = 0;
dErrorsFixedInCurrentWindow
-= dErrorsFixed[ nPointerToLastAddedBase ];
dErrorsFixedInCurrentWindow +=
dErrorsFixedThisBase2;
dErrorsFixed[ nPointerToLastAddedBase ] = dErrorsFixedThisBase2;
if ( dErrorsFixedInCurrentWindow >
dErrorsFixedInWindowOfMostErrorsFixed ) {
dErrorsFixedInWindowOfMostErrorsFixed =
dErrorsFixedInCurrentWindow;
nWindowOfMostErrorsFixedUnpaddedRight = nUnpadded;
}
}
} // if ( bFindWindowOfMostErrorsFixed ) {
} // int nUnpadded =...
} // for( int nTemplateTransitions = 0;
return( dTotalErrorsFixed );
} // dGetErrorsFixedForCustomPrimerRead
double Contig :: dGetErrorsFixedForWholeCloneRead(
const bool bTopStrandNotBottomStrand,
const int nFirstBaseOfReadUnpadded,
const bool bChooseExperiment,
const bool bJustConsensusNotGaps,
const bool bFindWindowOfMostErrorsFixed ) {
return( dGetErrorsFixedForSubcloneRead(
&pFakeSubcloneForWholeCloneReads_,
1,
bTopStrandNotBottomStrand,
nFirstBaseOfReadUnpadded,
false, // bExcludeVectorBasesAtBeginningOfRead,
bChooseExperiment,
bJustConsensusNotGaps,
bFindWindowOfMostErrorsFixed ) );
}
void Contig :: createFakeSubcloneTTemplatesForWholeCloneReads() {
pFakeSubcloneForWholeCloneReads_ = new subcloneTTemplate();
int nUnpaddedMin = pUnpaddedErrorProbabilities_->nGetStartIndex();
int nUnpaddedMax = pUnpaddedErrorProbabilities_->nGetEndIndex();
pFakeSubcloneForWholeCloneReads_->nUnpaddedStart_ =
nUnpaddedMin;
pFakeSubcloneForWholeCloneReads_->nUnpaddedEnd_ =
nUnpaddedMax;
pFakeSubcloneForWholeCloneReads_->bOKToUseThisTemplate_ = true;
// pFakeSubcloneForWholeCloneReads_->bTemplateIsDoubleStranded_ = true;
pFakeSubcloneForWholeCloneReads_->pContig_ = this;
pFakeSubcloneForWholeCloneReads_->bUnpaddedStartKnownPrecisely_ = true;
pFakeSubcloneForWholeCloneReads_->bUnpaddedEndKnownPrecisely_ = true;
pFakeSubcloneForWholeCloneReads_->bOKToUseForTopStrandReads_ = true;
pFakeSubcloneForWholeCloneReads_->bOKToUseForBottomStrandReads_ = true;
}
void Contig :: findWindowOfMostErrorsFixed( experiment* pExp ) {
double dErrorsFixed;
if ( pExp->bIsUniversalPrimerRead() ) {
bool bJustCountConsensusErrorsFixed = false;
if ( !pExp->bDeNovoUniversalPrimerRead_ &&
!pCP->bAutoFinishAllowResequencingReadsToExtendContigs_ )
bJustCountConsensusErrorsFixed = true;
dErrorsFixed = dGetErrorsFixedForSubcloneRead(
&(pExp->pSubcloneTemplate_),
1, // nNumberOfSubclones
pExp->bTopStrandNotBottomStrand_,
pExp->nFirstBaseOfReadUnpadded_,
true, // bExcludeVectorBasesAtBeginningOfRead
false, // bChooseExperiment
bJustCountConsensusErrorsFixed, // bJustConsensusNotGaps
true ); // bFindWindowOfMostErrorsFixed
}
else if ( pExp->bIsCustomPrimerSubcloneRead() ) {
dErrorsFixed = dGetErrorsFixedForCustomPrimerSubcloneRead(
pExp->pCustomPrimer_,
pExp->bTopStrandNotBottomStrand_,
pExp->nFirstBaseOfReadUnpadded_,
false, // bChooseExperiment
false, // bJustConsensusNotGaps
true ); // bFindWindowOfMostErrorsFixed
}
else if ( pExp->bIsWholeCloneRead() ) {
dErrorsFixed = dGetErrorsFixedForWholeCloneRead(
pExp->bTopStrandNotBottomStrand_,
pExp->nFirstBaseOfReadUnpadded_,
false, // bChooseExperiment
false, // bJustConsensusNotGaps
true ); // bFindWindowOfMostErrorsFixed
}
else
assert( false ); // shouldn't be any other kind of experiment
// the result is stored in these static variables
pExp->nWindowOfMostErrorsFixedUnpaddedRight_ =
nWindowOfMostErrorsFixedUnpaddedRight;
pExp->dErrorsFixedInWindow_ =
dErrorsFixedInWindowOfMostErrorsFixed;
}
void Contig :: clearExperimentCandidatesList() {
for( int nExp = 0; nExp < pCandidateExperiments_->length(); ++nExp ) {
experiment* pExp = (*pCandidateExperiments_)[ nExp ];
if ( pExp->nState_ != nChosenExperiment )
delete pExp;
}
pCandidateExperiments_->clear();
}
void Contig :: saveStateOfConsensus() {
*pSavedUnpaddedErrorProbabilities_
= *pUnpaddedErrorProbabilities_;
*pSavedUnpaddedReadTypesCoveringEachBase_
= *pUnpaddedReadTypesCoveringEachBase_;
*pSavedDistributionsOfChosenExperimentsArray_ =
*pDistributionsOfChosenExperimentsArray_;
}
void Contig :: restoreStateOfConsensus() {
*pUnpaddedReadTypesCoveringEachBase_ =
*pSavedUnpaddedReadTypesCoveringEachBase_;
// restore old distribution array
*pDistributionsOfChosenExperimentsArray_ =
*pSavedDistributionsOfChosenExperimentsArray_;
// now restore (partly) the pUnpaddedErrorProbabilities_ array
double dRedundancy = pCP->dAutoFinishRedundancy_;
if ( dRedundancy < 1.0 )
dRedundancy = 1.0;
if ( dRedundancy > 2.0 )
dRedundancy = 2.0;
for( int nUnpadded = pUnpaddedErrorProbabilities_->nGetStartIndex();
nUnpadded <= pUnpaddedErrorProbabilities_->nGetEndIndex();
++nUnpadded ) {
if ( (*pUnpaddedErrorProbabilities_)[ nUnpadded ] == 0.0 ||
(*pSavedUnpaddedErrorProbabilities_)[ nUnpadded ] == 0.0 )
(*pUnpaddedErrorProbabilities_)[ nUnpadded ] = 0.0;
else {
// create affine sum of error probabilities before and
// after the first pass of universal primer reads
// Use the affine sum of logs instead of error probabilities
// since this will reflect more closely the number of additional
// reads to be called.
(*pUnpaddedErrorProbabilities_)[ nUnpadded ] =
exp(
(2.0 - dRedundancy) *
log( (*pUnpaddedErrorProbabilities_)[ nUnpadded ] )
+
( dRedundancy - 1.0 ) *
log( (*pSavedUnpaddedErrorProbabilities_)[ nUnpadded ] )
);
}
// if ( dRedundancy == 2.0 ) {
// if ( ABS( (*pSavedUnpaddedErrorProbabilities_)[ nUnpadded ] -
// (*pUnpaddedErrorProbabilities_)[ nUnpadded ] )
// > .001 ) {
// fprintf( pAO, "too big of difference\n" );
// fprintf( pAO, "%.9g ",
// (*pSavedUnpaddedErrorProbabilities_)[ nUnpadded ] -
// (*pUnpaddedErrorProbabilities_)[ nUnpadded ] );
// }
// }
// else if ( dRedundancy == 1.0 ) {
// fprintf( pAO, "%.9g ",
// (*pUnpaddedErrorProbabilities_)[ nUnpadded ] - d );
// }
}
}
void Contig :: recalculateConsensusUsingAllChosenUniversalPrimerExperiments() {
for( int nExp = 0;
nExp < pUniversalPrimerExperimentsForLowQualityRegions_->length();
++nExp ) {
experiment* pExp =
(*pUniversalPrimerExperimentsForLowQualityRegions_)[ nExp ];
double d =
dGetErrorsFixedForSubcloneRead(
&(pExp->pSubcloneTemplate_),
1, // nNumberOfSubclones
pExp->bTopStrandNotBottomStrand_,
pExp->nFirstBaseOfReadUnpadded_,
true, // bExcludeVectorBasesAtBeginningOfRead
true, // bChooseExperiment
false, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
}
}
void Contig :: recalculateErrorsFixedOfExperimentCandidates() {
for( int nCan = 0; nCan < pCandidateExperiments_->length(); ++nCan ) {
experiment* pExp = (*pCandidateExperiments_)[ nCan ];
if ( pExp->nState_ == nChosenExperiment ) continue;
// make all experiments fair game again
if ( pExp->nState_ == nRejectedExperiment )
pExp->nState_ = nAvailableExperiment;
bool bJustCountConsensusErrorsFixed = false;
// optionally do not allow resequencing reads to close gaps
if ( pExp->bIsUniversalPrimerRead() &&
!pExp->bDeNovoUniversalPrimerRead_ &&
!pCP->bAutoFinishAllowResequencingReadsToExtendContigs_ )
bJustCountConsensusErrorsFixed = true;
pExp->dErrorsFixed_ =
dGetErrorsFixedForSubcloneRead(
&(pExp->pSubcloneTemplate_),
1, // nNumberOfSubclones
pExp->bTopStrandNotBottomStrand_,
pExp->nFirstBaseOfReadUnpadded_,
true, // bExcludeVectorBasesAtBeginningOfRead
false, // bChooseExperiment
bJustCountConsensusErrorsFixed, // bJustConsensusNotGaps
false ); // bFindWindowOfMostErrorsFixed
setVariousErrorData( pExp );
}
}
void Contig :: chooseExperimentsToCoverLowQualityRegionsAndGapsFor9_66() {
fprintf( pAO, "now attempting to solve weak regions in consensus and extend the consensus...\n\n" );
// Possible problem: if the consensus is high quality, but you want
// autofinish to extend the consensus, it might not do so.
// (DG, Feb 2000)
if ( dExpectedNumberOfErrorsInExtendedContig_ <= dTargetNumberOfErrors_ ) {
fprintf( pAO, "Estimated number of errors in extended contig %.2f is already less than the maximum %.2f for this contig\n",
dExpectedNumberOfErrorsInExtendedContig_,
dTargetNumberOfErrors_ );
return;
}
chooseExperimentsToCoverLowQualityRegionsAndGapsFor9_66_2();
if ( pCP->bAutoFinishCloseGaps_ )
warnUserIfCouldNotExtendContig();
}
void Contig :: chooseExperimentsToCoverLowQualityRegionsAndGapsFor9_66_2() {
if ( pCandidateExperiments_->length() == 0 ) {
fprintf( pAO, "No possible experiments found\n" );
return;
}
bool bMoreAvailableExperiments;
bool bConsiderMoreExperiments = true;
do {
experiment* pBestExp = pGetBestExperimentCandidate();
if ( pBestExp == NULL ) {
bMoreAvailableExperiments = false;
bConsiderMoreExperiments = false;
}
else {
// found the best experiment
assert( pBestExp->dErrorsFixed_ >= consedParameters::pGetConsedParameters()->dAutoFinishMinNumberOfErrorsFixedByAnExpToUse_ );
// the best experiment is fine--transfer it
pBestExp->nPurpose_ = nFixWeakRegion;
pBestExp->chooseExperiment();
possiblyUpdateContigHighQualitySegment( pBestExp );
pChosenExperiments_->insert( pBestExp );
ConsEd::pGetAssembly()->pAllChosenExperiments_->insert( pBestExp );
adjustExpectedNumberOfErrors( pBestExp );
findWindowOfMostErrorsFixed( pBestExp );
pBestExp->print();
fprintf( pAO, "After this experiment,\ncalculated number of errors in extended contig: %.2f and consensus: %.2f\n\n",
dExpectedNumberOfErrorsInExtendedContig_,
dExpectedNumberOfErrorsInConsensus_ );
if ( !bStillTooManyExpectedErrors() )
bConsiderMoreExperiments = false;
else {
adjustRemainingExperimentCandidates( pBestExp );
// check no bugs
double dErrorsInConsensus;
double dErrorsInExtendedContig;
recalculateCurrentExpectedNumberOfErrorsForThisContig(
dErrorsInConsensus,
dErrorsInExtendedContig );
if (
( ABS( dErrorsInConsensus - dExpectedNumberOfErrorsInConsensus_ ) > .000001 ) ||
( ABS( dErrorsInExtendedContig - dExpectedNumberOfErrorsInExtendedContig_ ) > .000001 )
) {
fprintf( pAO, "Warning: dExpectedNumberOfErrorsInConsensus = %e but calculated is %e and dExpectedNumberOfErrorsInExtendedContig_ = %e but calculated is %e\n",
dExpectedNumberOfErrorsInConsensus_,
dErrorsInConsensus,
dExpectedNumberOfErrorsInExtendedContig_,
dErrorsInExtendedContig );
}
else {
fprintf( pAO, "OK3\n" );
}
// end check no bugs
}
} // if ( pBestExp == NULL ) else
} while( bConsiderMoreExperiments );
// tell user why we stopped
// We could have stopped because:
// 1) the error rate was below the threshold
// 2) there were no more experiments
// 3) there were more experiments, but they weren't worth doing
if ( bStillTooManyExpectedErrors() ) {
fprintf( pAO, "There are still too many expected errors but we stopped because\n" );
if ( bMoreAvailableExperiments )
fprintf( pAO, "all remaining available experiments each solved too \nfew errors\nto be worthwhile doing\n" );
else
fprintf( pAO, "there were no more experiments this program could \nfind to resolve the remaining errors.\n" );
printWorstRemainingRegion();
}
else
fprintf( pAO, "The expected error count for the contig should be acceptable after doing these experiments\n" );
}
void Contig :: createExperimentCandidatesListFor9_66() {
fprintf( pAO, "creating list of possible experiments with universal primers..." );
fflush( pAO );
createExperimentsWithResequencingUniversalPrimers();
fprintf( pAO, "done\n" );
fflush( pAO );
if ( pCP->bAutoFinishAllowCustomPrimerSubcloneReads_ ||
pCP->bAutoFinishAllowWholeCloneReads_ ) {
fprintf( pAO, "creating list of possible experiments with custom primers..." );
fflush( pAO );
createExperimentsWithCustomPrimers();
fprintf( pAO, "done\n\n" );
fflush( pAO );
}
fprintf( pAO, "creating list of possible experiments of reverses where there is only a forward or forwards where there is only a reverse..." );
fflush( pAO );
createExperimentsWithFirstForwardsOrFirstReverses();
fprintf( pAO, "done\n\n" );
fflush( pAO );
}
bool mbtValVectorOfBool::bMasksSet_ = false;
unsigned char
mbtValVectorOfBool::ucSetBitMasks_[nNumberOfBitsInMbtValVectorOfBoolElement];
unsigned char
mbtValVectorOfBool::ucClearBitMasks_[nNumberOfBitsInMbtValVectorOfBoolElement];
void Contig :: determineLocationsOfRunsAndStops() {
int nLength = nGetUnpaddedSeqLength();
int nOffset = nGetUnpaddedStartIndex();
pUnpaddedLocationsOfRuns_ = new mbtValVectorOfBool( nLength,
nOffset,
"Contig::pUnpaddedLocationsOfRuns_" );
pUnpaddedLocationsOfStops_ = new mbtValVectorOfBool( nLength,
nOffset,
"Contig::pUnpaddedLocationsOfStops_" );
gotoList* pStopsGotoList = new gotoList();
addStopsCausingLCQs( pStopsGotoList,
false ); // bSetPaddedPositionsArray
const int nWindowSizeOfStops = 20;
int nGotoItem;
for( nGotoItem = 0;
nGotoItem < pStopsGotoList->nGetNumGotoItems();
++nGotoItem ) {
gotoItem* pGotoItem = pStopsGotoList->pGetGotoItem( nGotoItem );
int nUnpaddedLeft;
int nUnpaddedRight;
if ( pGotoItem->soInitialDescription_ == " stop left of here" ) {
nUnpaddedRight = pGotoItem->nUnpaddedGotoItemEnd_;
nUnpaddedLeft =
pGotoItem->nUnpaddedGotoItemEnd_ - nWindowSizeOfStops + 1;
}
else if ( pGotoItem->soInitialDescription_ == " stop right of here" ) {
nUnpaddedLeft = pGotoItem->nUnpaddedGotoItemStart_;
nUnpaddedRight = nUnpaddedLeft + nWindowSizeOfStops - 1;
}
// don't be off the consensus
if ( bIntersect( nUnpaddedLeft,
nUnpaddedRight,
nGetUnpaddedStartIndex(), // consensus
nGetUnpaddedEndIndex(),
nUnpaddedLeft,
nUnpaddedRight ) ) {
for( int nUnpadded = nUnpaddedLeft; nUnpadded <= nUnpaddedRight;
++nUnpadded ) {
pUnpaddedLocationsOfStops_->setValue( nUnpadded, 1 );
}
}
} // for( int nGotoItem = 0;
delete pStopsGotoList;
gotoList* pRunsGotoList = new gotoList();
addRunsAndMicrosatellitesCausingLCQs( pRunsGotoList,
false ); // bSetPaddedPositionsArray
for( nGotoItem = 0;
nGotoItem < pRunsGotoList->nGetNumGotoItems();
++nGotoItem ) {
gotoItem* pGotoItem = pRunsGotoList->pGetGotoItem( nGotoItem );
int nUnpaddedLeft = pGotoItem->nUnpaddedGotoItemStart_;
int nUnpaddedRight = pGotoItem->nUnpaddedGotoItemEnd_;
// don't be off the consensus
if ( bIntersect( nUnpaddedLeft,
nUnpaddedRight,
nGetUnpaddedStartIndex(), // consensus
nGetUnpaddedEndIndex(),
nUnpaddedLeft,
nUnpaddedRight ) ) {
for( int nUnpadded = nUnpaddedLeft; nUnpadded <= nUnpaddedRight;
++nUnpadded ) {
pUnpaddedLocationsOfRuns_->setValue( nUnpadded, 1 );
}
}
} // for( int nGotoItem = 0;
delete pRunsGotoList;
}
bool Contig :: bDoesRegionCrossARunOrStop( const int nUnpaddedLeft,
const int nUnpaddedRight,
bool& bCrossesRun,
bool& bCrossesStop ) {
bCrossesRun = false;
bCrossesStop = false;
int nUnpaddedLeftOnConsensus;
int nUnpaddedRightOnConsensus;
if ( bIntersect( nUnpaddedLeft,
nUnpaddedRight,
nGetUnpaddedStartIndex(),
nGetUnpaddedEndIndex(),
nUnpaddedLeftOnConsensus,
nUnpaddedRightOnConsensus ) ) {
for( int nUnpadded = nUnpaddedLeftOnConsensus;
nUnpadded <= nUnpaddedRightOnConsensus; ++nUnpadded ) {
if ( (*pUnpaddedLocationsOfStops_)[ nUnpadded ] )
bCrossesStop = true;
if ( (*pUnpaddedLocationsOfRuns_)[ nUnpadded ] )
bCrossesRun = true;
}
}
if ( bCrossesRun || bCrossesStop ) {
return( true );
}
else
return( false );
}
void Contig :: determineWhetherThisExperimentNeedsSpecialChemistry(
experiment* pExp ) {
bDoesRegionCrossARunOrStop( pExp->nUnpaddedLeft_,
pExp->nUnpaddedRight_,
pExp->bCrossesRun_,
pExp->bCrossesStop_ );
}