/*****************************************************************************
# 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 "searchForSingleSubcloneRegions.h"
#include "guiMultiContigNavigator.h"
#include "contig.h"
#include "locatedFragment.h"
#include "consedParameters.h"
#include "rwtvalvector.h"
#include "rwtptrvector.h"
#include <stdio.h>
void searchForSingleSubcloneRegions( Widget widGuiContigWin,
Contig* pContig,
ContigWin* pContigWin ) {
RWCString soColumnHeadings1 = "Contig Start End";
RWCString soColumnHeadings2 = " Pos Pos";
guiMultiContigNavigator* pGuiMultiContigNav =
new guiMultiContigNavigator(
"Regions Covered by Only 1 Subclone (or None at All)",
soColumnHeadings1,
soColumnHeadings2,
100,
"", // extra button
NULL, // extra button
NULL, // extra button
widGuiContigWin,
NULL, // pGotoList
pContigWin
);
int nPaddedLength = pContig->nGetDisplayableContigLength();
RWTPtrVector<LocatedFragment> aPaddedContigOfReadPtrs( nPaddedLength, NULL );
RWTValVector<bool> aPaddedContigOfFlags( nPaddedLength, false );
for( int nRead = 0; nRead < pContig->nGetNumberOfFragsInContig();
++nRead ) {
// get pointer to this located frag from contig
LocatedFragment* pLocFrag = pContig->pLocatedFragmentGet( nRead );
// Phil said 980611 that:
// > David: We *do* consider it covered by the unaligned part of a read if the
// > user has turn off 'dim unaligned ends of read'. To be consistent, we
// > should then say that in this case it is covered by a totally unaligned
// > read.
// Phil: No, this case is very different, since the read doesn't belong there at
// all. (In the case where part of the read is aligned, the unaligned part is
// still "covering" that region.)
if ( pLocFrag->bIsWholeReadUnaligned() )
continue;
RWCString soSubcloneName = pLocFrag->soGetTemplateName();
int nConsPosStart = pLocFrag->nGetAlignStart();
int nConsPosEnd = pLocFrag->nGetAlignEnd();
if ( !consedParameters::pGetConsedParameters()->bSearchFunctionsUseLowQualityEndsOfReads_ ) {
if (pLocFrag->bIsWholeReadLowQuality() )
continue; // this read is no help
nConsPosStart = pLocFrag->nGetHighQualityStart();
nConsPosEnd = pLocFrag->nGetHighQualityEnd();
}
if ( !consedParameters::pGetConsedParameters()->bSearchFunctionsUseUnalignedEndsOfReads_ ) {
if ( nConsPosStart < pLocFrag->nGetAlignClipStart() )
nConsPosStart = pLocFrag->nGetAlignClipStart();
if ( pLocFrag->nGetAlignClipEnd() < nConsPosEnd )
nConsPosEnd = pLocFrag->nGetAlignClipEnd();
}
for( int nConsPos = nConsPosStart; nConsPos <= nConsPosEnd; ++nConsPos ) {
if ( ! aPaddedContigOfReadPtrs[ nConsPos - pContig->nGetFirstDisplayableContigPos() ] ) {
aPaddedContigOfReadPtrs[ nConsPos - pContig->nGetFirstDisplayableContigPos() ] =
pLocFrag;
}
else {
// case in which there is already a subclone in that position
// Let's see if it is the same subclone
// If this base is already double covered by subclones, the
// additional read doesn't matter
if (aPaddedContigOfFlags[ nConsPos - pContig->nGetFirstDisplayableContigPos() ] )
continue;
LocatedFragment* pLocFragAlreadyFound =
aPaddedContigOfReadPtrs[ nConsPos - pContig->nGetFirstDisplayableContigPos() ];
if ( pLocFragAlreadyFound->soGetTemplateName() !=
soSubcloneName ) {
// so this base is covered by more than one subclone
aPaddedContigOfFlags[ nConsPos - pContig->nGetFirstDisplayableContigPos() ] = true;
}
} // if ( ! aPaddedContig...
} // for( int nConsPos
} // for( int nRead = 0; ...
// now find the single stranded segments
int nStart;
int nEnd;
bool bWithinSingleSubcloneRegion = false;
for( int nConsPos = pContig->nGetStartConsensusIndex();
nConsPos <= pContig->nGetEndConsensusIndex();
++nConsPos ) {
Ntide nt = pContig->ntGetCons( nConsPos );
// ignore pads
if ( nt.bIsPad() ) continue;
bool bDouble = aPaddedContigOfFlags[ nConsPos - pContig->nGetFirstDisplayableContigPos() ];
if ( !bDouble ) {
// in single stranded region
if (bWithinSingleSubcloneRegion ) {
// nEnd keeps growing with each base added
nEnd = nConsPos;
}
else {
// start a single stranded region
bWithinSingleSubcloneRegion = true;
nStart = nConsPos;
nEnd = nConsPos;
}
}
else {
if ( bWithinSingleSubcloneRegion ) {
// just finished a single stranded region
char szLength[200];
sprintf( szLength, "%5d bp", nEnd - nStart + 1 );
RWCString soDescription( szLength );
gotoItem* pGotoItem = new gotoItem(
pContig,
NULL, // consensus location--not read
nStart,
nEnd,
pContig->nUnpaddedIndex( nStart ),
pContig->nUnpaddedIndex( nEnd ),
soDescription );
bWithinSingleSubcloneRegion = false;
pGuiMultiContigNav->appendToList( pGotoItem );
}
} // if ( aPaddedContig ... else
} // for( int nConsPos
if ( bWithinSingleSubcloneRegion ) {
// just finished a single subclone region
char szLength[200];
sprintf( szLength, "%5d bp", nEnd - nStart + 1 );
RWCString soDescription( szLength );
gotoItem* pGotoItem = new gotoItem(
pContig,
NULL, // consensus not read
nStart,
nEnd,
pContig->nUnpaddedIndex( nStart ),
pContig->nUnpaddedIndex( nEnd ),
soDescription );
pGuiMultiContigNav->appendToList( pGotoItem );
}
pGuiMultiContigNav->displayList( true ); // sort by left end positions
}