/*****************************************************************************
# 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 "autoReport.h"
#include "contig.h"
#include "locatedFragment.h"
#include "rwcstring.h"
#include "assembly.h"
#include "consed.h"
#include "consedParameters.h"
#include "terminateIfNoPhdDir.h"
#include "soAddCommas.h"
#include "filename.h"
#include "rwctokenizer.h"
#include "fileDefines.h"
#include "openLogFiles.h"
#include "printAutoFinishMiscInfo.h"
#include "soLine.h"
#include "soGetErrno.h"
#include "listOfLibraries.h"
#include "lowconsqual.h"
#include "highQualityDiscrepancyGotoList.h"
#include "singleSubcloneRegionsGotoList.h"
#include <stdio.h>
#include "primateSpecies.h"
#include "fwdRevPair2.h"
#include "fwdRevPair.h"
autoReport :: autoReport( const FileName& filAceFileToOpen ) :
filAceFileToOpen_( filAceFileToOpen ) {}
void autoReport :: doIt() {
pCP->bOKToUseGui_ = false;
ConsEd* pConsed = new ConsEd();
pConsed->pAutoReport_ = this;
openLogFiles( filAceFileToOpen_ );
// this uses the output file to report the error
// so it must have openAutoFinishOutputFile before it
terminateIfNoPhdDir();
printAutoFinishMiscInfo( filAceFileToOpen_, "" );
// I want to read the list of libraries before reading the ace file
// to make debugging by the user easier
pCP->pListOfLibraries_ = new listOfLibraries();
pCP->pListOfLibraries_->aLibraries_.soName_ =
"pCP->pListOfLibraries_->aLibraries_ in autoEdit.cpp";
pCP->pListOfLibraries_->parseLibraryFile();
parseSpecies();
// tried moving this to above to before openAutoFinishOutputFiles
// so that, if
// the user specifies a bogus ace file, such as -ace -ace,
// we don't create several zero-length files that have filenames
// that the shell can't handle. However, if I do that, then any
// error encountered by openAssemblyFile will try to write to
// pAO and give a segmentation fault
// this does "new Assembly"
ConsEd::pGetConsEd()->openAssemblyFile( filAceFileToOpen_ );
Assembly* pAssembly = ConsEd::pGetAssembly();
pAssembly->setAllPaddedPositionsArrays();
checkPrimateSpeciesAndConstants( this );
if ( pCP->bAutoReportPrintSpeciesAlignment_ ) {
dumpSpeciesAlignment();
}
if ( pCP->bAutoReportPrintReadAlignment_ ) {
printReadAlignment();
}
if ( pCP->bAutoReportPrintDiscrepantRegions_ ) {
printDiscrepantRegions();
}
if ( pCP->bAutoReportPrintBasesInDiscrepantRegions_ ) {
printBasesInDiscrepantRegions();
}
if ( pCP->bAutoReportPrintMinimumQualityHistogram_ ) {
printMinimumQualityHistogram();
}
if ( pCP->bAutoReportPrintNumberOfIsolatedPads_ ) {
printNumberOfIsolatedPads();
}
if ( pCP->bAutoReportPrintNumberOfIsolatedPadsForEachSpecies_ ) {
printNumberOfIsolatedPadsForEachSpecies();
}
if ( pCP->bAutoReportPrintAgreeDisagreeBetweenPairsOfSpecies_ ) {
printAgreeDisagreeBetweenPairsOfSpecies();
}
if ( pCP->bAutoReportPrintAgreeDisagreeBetweenPairsOfSpecies2_ ) {
printAgreeDisagreeBetweenPairsOfSpecies2();
}
if ( pCP->bAutoReportPrintScaffolds_ ) {
printScaffolds();
}
if ( pCP->bAutoReportCalculateErrorProbabilitiesByComparingPTroPPan_ ) {
calculateErrorProbabilitiesByComparingPTroPPan();
}
if ( pCP->bAutoReportPrintIfReadsAreCorrectlyAligned_ ) {
printIfReadsAreCorrectlyAligned();
}
if ( pCP->bAutoReportPrintLengthsOfUnalignedHighQualitySegmentsOfReads_ ) {
printLengthsOfUnalignedHighQualitySegments();
}
if ( pCP->bAutoReportPrintLengthsOfAlignedSegmentsOfReads_ ) {
printLengthsOfAlignedSegments();
}
if ( pCP->bAutoReportCompareTopAndBottomStrands_ ) {
compareTopAndBottomStrands();
}
if ( pCP->bAutoReportCompareTopAndBottomStrands2_ ) {
compareTopAndBottomStrands2();
}
if ( pCP->bAutoReportCompareTopAndBottomStrands3_ ) {
compareTopAndBottomStrands3();
}
if ( pCP->bAutoReportCompareTopAndBottomStrands4_ ) {
compareTopAndBottomStrands4();
}
if ( pCP->bAutoReportSingleSignalInfo_ ) {
singleSignalInfo();
}
if ( pCP->bAutoReportSingleSignalInfo2_ ) {
singleSignalInfo2();
}
if ( pCP->bAutoReportCompareTopAndBottomStrandsWithHuman_ ) {
compareTopAndBottomStrandsWithHuman();
}
if ( pCP->bAutoReportCountColumnsForGroupsOfSpecies_ ) {
countColumnsForGroupsOfSpecies();
}
if ( pCP->bAutoReportCompareHQSWithLQS_ ) {
compareHQSWithLQS();
}
if ( pCP->bAutoReportLowQualityBasesInHQS_ ) {
lowQualityBasesInHQS();
}
if ( pCP->bAutoReportSingleSignalOrQuality_ ) {
singleSignalOrQuality();
}
if ( pCP->bAutoReportDiscrepancyRateInFlankedRegions_ ) {
discrepancyRateInFlankedRegions();
}
if ( pCP->bAutoReportDiscrepancyRateInFlankedRegions2_ ) {
discrepancyRateInFlankedRegions2();
}
if ( pCP->bAutoReportDiscrepancyRateInFlankedRegions4_ ) {
discrepancyRateInFlankedRegions4();
}
if ( pCP->bAutoReportDiscrepancyRateInFlankedRegions5_ ) {
discrepancyRateInFlankedRegions5();
}
if ( pCP->bAutoReportGoodReadsBug_ ) {
countReadsDueToBugInGoodReads();
}
if ( pCP->bAutoReportCompareTopAndBottomStrandsNoHuman_ ) {
compareTopAndBottomStrandsNoHuman();
}
if ( pCP->bAutoReportHighQualitySegmentData_ ) {
highQualitySegmentData();
}
if ( pCP->bAutoReportPrintFlankedColumns_ ) {
printFlankedColumns();
}
if ( pCP->bAutoReportPrintFlankedColumns2_ ) {
printFlankedColumns2();
}
if ( pCP->bAutoReportPrintFlankedColumns3_ ) {
printFlankedColumns3();
}
if ( pCP->bAutoReportPrintFlankedColumns4_ ) {
printFlankedColumns4();
}
if ( pCP->bAutoReportCountAcceptableColumnsWithNoneOnLeft_ ) {
countAcceptableColumnsWithNoneOnLeft();
}
if ( pCP->bAutoReportCountAllMutations_ ) {
countAllMutations();
}
if ( pCP->bAutoReportCountAllMutationsML_ ) {
countAllMutationsML();
}
if ( pCP->bAutoReportPrintHighQualityDiscrepancies_ ) {
printHighQualityDiscrepancies();
}
if ( pCP->bAutoReportPrintLowConsensusQualityRegions_ ) {
printLowConsensusQualityRegions();
}
if ( pCP->bAutoReportPrintSingleSubcloneRegions_ ) {
printSingleSubcloneRegions();
}
if ( pCP->bAutoReportPrintMutationsWithContext_ ) {
printMutationsWithContext();
}
if ( pCP->bAutoReportPrintCrudeChimpHumanMutations_ ) {
printCrudeChimpHumanMutations();
}
if ( pCP->bAutoReportPrintToCompareToReich_ ) {
printToCompareToReich();
}
if ( pCP->bAutoReportPrintLinkingForwardReversePairs_ ) {
printLinkingForwardReversePairs();
}
if ( pCP->bAutoReportPrintFilteredInconsistentForwardReversePairs_ ) {
printFilteredInconsistentForwardReversePairs();
}
if ( pCP->bAutoReportPrintHighlyDiscrepantRegions_ ) {
printHighlyDiscrepantRegions();
}
if ( pCP->bAutoReportPrintReadNamesInRegion_ ) {
ConsEd::pGetAssembly()->printReadNamesInRegion();
}
if ( pCP->bAutoReportPrintAssemblySummary_ ) {
ConsEd::pGetAssembly()->printAssemblySummary();
}
fclose( pAO );
cerr << "see " << pCP->filAutoFinishFullOutput_ << endl;
}
void autoReport :: dumpSpeciesAlignment() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
nContig++) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
dumpSpeciesAlignmentForOneContig( pContig );
}
}
void autoReport :: dumpSpeciesAlignmentForOneContig( Contig* pContig ) {
// create arrays for each species with length equal to this
// contig length (unpadded )
aArrayOfPtrsToArraysOfBases_.clearAndDestroy();
aArrayOfPtrsToArraysOfQualities_.clearAndDestroy();
aArrayOfPtrsToArraysOfHighQualitySegmentFlags_.clearAndDestroy();
aArrayOfPtrsToArraysOfLocatedFragments_.clearAndDestroy();
aArrayOfPtrsToArraysOfReadPositionsInDirectionOfSequencing_.clearAndDestroy();
int nUnpaddedLength = pContig->nGetUnpaddedSeqLength();
int nSpecies;
for( nSpecies = 0; nSpecies < aArrayOfSpecies_.length(); ++nSpecies ) {
RWCString soSpecies = aArrayOfSpecies_[ nSpecies ];
RWTValVectorOfBases* pArrayOfBases =
new RWTValVector<char>( nUnpaddedLength, '0' );
RWTValVectorOfQualities* pArrayOfQualities =
new RWTValVector<signed char>( nUnpaddedLength, -1 );
RWCString soArrayName = soSpecies + " bit array of high quality segment";
RWTValVectorOfHighQualitySegmentFlags* pArrayOfHighQualitySegmentFlags =
new RWTValVector<bool>( nUnpaddedLength, 0 );
pArrayOfHighQualitySegmentFlags->soName_ = soArrayName;
RWTPtrVectorOfLocatedFragment* pArrayOfReads =
new RWTPtrVector<LocatedFragment>( nUnpaddedLength, NULL );
RWTValVectorOfUnpaddedReadPositionsInDirectionOfSequencing*
pArrayOfUnpaddedReadPositionsInDirectionOfSequencing =
new RWTValVector<int>( nUnpaddedLength, -1 );
pArrayOfUnpaddedReadPositionsInDirectionOfSequencing->soName_ =
soSpecies +
" RWTValVectorOfUnpaddedReadPositionsInDirectionOfSequencing";
aArrayOfPtrsToArraysOfBases_.insert( pArrayOfBases );
aArrayOfPtrsToArraysOfQualities_.insert( pArrayOfQualities );
aArrayOfPtrsToArraysOfHighQualitySegmentFlags_.insert(
pArrayOfHighQualitySegmentFlags );
aArrayOfPtrsToArraysOfLocatedFragments_.insert( pArrayOfReads );
aArrayOfPtrsToArraysOfReadPositionsInDirectionOfSequencing_.insert(
pArrayOfUnpaddedReadPositionsInDirectionOfSequencing );
// find the reads for these species
for( int nRead = 0; nRead < pContig->nGetNumberOfFragsInContig();
++nRead ) {
LocatedFragment* pLocFrag = pContig->pLocatedFragmentGet( nRead );
if ( !pLocFrag->soGetName().bContains( soSpecies ) ) continue;
// if reached here, the read is of the species we want
if ( pLocFrag->bIsWholeReadUnaligned() ) continue;
// figure out the unpadded read position in direction of
// sequencing
int nUnpaddedReadPosInDirectionOfSequencing = -666;
bool bFirstTime = true;
for( int nConsPos = pLocFrag->nGetAlignClipStart();
nConsPos <= pLocFrag->nGetAlignClipEnd();
++nConsPos ) {
Ntide nt = pLocFrag->ntGetFragFromConsPos( nConsPos );
// keeping track of read position
if ( bFirstTime ) {
bFirstTime = false;
nUnpaddedReadPosInDirectionOfSequencing =
pLocFrag->nOrientedUnpaddedFragPosFromConsPos(
pLocFrag->nGetAlignClipStart() );
}
else {
if ( !nt.bIsPad() ) {
if ( pLocFrag->bComp() ) {
--nUnpaddedReadPosInDirectionOfSequencing;
}
else {
++nUnpaddedReadPosInDirectionOfSequencing;
}
}
}
// if the human sequence has a pad in this position, ignore
// this position. Since the consensus is the same as the
// human sequence, use it.
if ( pContig->ntGetCons( nConsPos ).bIsPad() ) continue;
int n0UnpaddedPos = pContig->nUnpaddedIndex( nConsPos ) - 1;
bool bReplaceBase = false;
if ( (*pArrayOfReads)[ n0UnpaddedPos ] == NULL ) {
// this is the first read of this species
bReplaceBase = true;
}
else {
// case in which this is the 2nd read of this species.
LocatedFragment* pLocFragB = (*pArrayOfReads)[ n0UnpaddedPos ];
if ( nt.cGetBase() == (*pArrayOfBases)[ n0UnpaddedPos ] ) {
// update various arrays
// Note that this means that the
// RWTPtrVectorOfLocatedFragment array will not necessarily
// contain a read of this quality and high quality segment
// since one read may be in the high quality segment
// and the other higher quality. We will report
// the best of both.
// update quality array
if ( nt.qualGetQualityWithout98or99() >
(*pArrayOfQualities)[ n0UnpaddedPos ] ) {
(*pArrayOfQualities)[ n0UnpaddedPos ] =
nt.qualGetQualityWithout98or99();
}
// update the high quality segment array
if ( pLocFrag->bIsInHighQualitySegmentOfRead( nConsPos )
&&
( (*pArrayOfHighQualitySegmentFlags)[ n0UnpaddedPos ]
== false ) ) {
(*pArrayOfHighQualitySegmentFlags)[ n0UnpaddedPos ] =
true;
}
}
else {
// case in which the 2 reads disagree
// choose the one that is in the high quality segment
bool bReadAInHighQualitySegment =
pLocFrag->bIsInHighQualitySegmentOfRead( nConsPos );
bool bReadBInHighQualitySegment =
pLocFragB->bIsInHighQualitySegmentOfRead( nConsPos );
if ( bReadAInHighQualitySegment &&
!bReadBInHighQualitySegment ) {
bReplaceBase = true;
}
else if ( !bReadAInHighQualitySegment &&
bReadBInHighQualitySegment ) {
continue; // keep using the existing read
}
else {
assert( bReadAInHighQualitySegment ==
bReadBInHighQualitySegment );
// both reads are in high quality segment or
// both reads are in low quality segment
// choose the read based on the average quality
// in a window
const int nWindowSize = 11;
float fQualityReadA =
pLocFrag->fGetAverageQualityInWindow( nConsPos,
nWindowSize );
float fQualityReadB =
pLocFragB->fGetAverageQualityInWindow( nConsPos,
nWindowSize );
if ( fQualityReadA > fQualityReadB ) {
bReplaceBase = true;
}
}
} // // if ( nt.cGetBase() != ... else
} // if ( (*pArrayOfReads)[ n0UnpaddedPos ] == NULL ) else
if ( bReplaceBase ) {
(*pArrayOfBases)[ n0UnpaddedPos ] = nt.cGetBase();
(*pArrayOfQualities)[ n0UnpaddedPos ] = nt.qualGetQualityWithout98or99();
(*pArrayOfHighQualitySegmentFlags)[ n0UnpaddedPos ] =
pLocFrag->bIsInHighQualitySegmentOfRead( nConsPos );
(*pArrayOfReads)[ n0UnpaddedPos ] = pLocFrag;
(*pArrayOfUnpaddedReadPositionsInDirectionOfSequencing)[ n0UnpaddedPos ] =
nUnpaddedReadPosInDirectionOfSequencing;
}
} // for( int nConsPos = pLocFrag->nGetAlignClipStart();
} // for( int nRead = 0; nRead < pContig->nGetNumberOfFragsInContig();
} // for( int nSpecies = 0; nSpecies < aArrayOfSpecies.length(); ++nSpecies ) {
printThisContig( pContig );
}
void autoReport :: parseSpecies() {
aArrayOfSpecies_.clear();
RWCTokenizer tok( pCP->soAutoReportSpecies_ );
RWCString soSpecies;
while( !( soSpecies = tok() ).isNull() ) {
aArrayOfSpecies_.insert( soSpecies );
}
}
void autoReport :: printThisContig( Contig* pContig ) {
pContig->setStartNumberingUnpaddedConsensusAtUserDefinedIfNecessary();
// print headings
fprintf( pAO, "\n\nAutoReportPrintSpeciesAlignment {\n" );
fprintf( pAO, "CONTIG: %s\n\n",
pContig->soGetName().data() );
for( int nSpecies = 0; nSpecies < aArrayOfSpecies_.length(); ++nSpecies ) {
fprintf( pAO, " %10s", aArrayOfSpecies_[ nSpecies ].data() );
}
fprintf( pAO, "\n" );
RWCString soLineToOutput( (size_t) 1000 );
for( int n0Unpadded = 0; n0Unpadded < pContig->nGetUnpaddedSeqLength();
++n0Unpadded ) {
soLineToOutput = pCP->soAutoReportPrefix_;
soLineToOutput += " ";
soLineToOutput +=
soAddCommas( n0Unpadded + pContig->nStartNumberingUnpaddedConsensusAtUserDefined_ );
soLineToOutput.padOnLeft( 12 );
int nSpecies;
for( nSpecies = 0; nSpecies < aArrayOfSpecies_.length();
++nSpecies ) {
RWTValVectorOfBases* pArrayOfBases =
aArrayOfPtrsToArraysOfBases_[ nSpecies ];
soLineToOutput += " ";
if ( (*pArrayOfBases)[ n0Unpadded ] == '*' )
soLineToOutput += "-";
else
soLineToOutput += (*pArrayOfBases)[ n0Unpadded ];
}
for( nSpecies = 0; nSpecies < aArrayOfSpecies_.length();
++nSpecies ) {
RWTValVectorOfQualities* pArrayOfQualities =
aArrayOfPtrsToArraysOfQualities_[ nSpecies ];
soLineToOutput += " ";
RWCString soTemp( (long) (*pArrayOfQualities)[ n0Unpadded ] );
soTemp.padOnLeft( 2 );
soLineToOutput += soTemp;
}
for( nSpecies = 0; nSpecies < aArrayOfSpecies_.length();
++nSpecies ) {
RWTValVectorOfHighQualitySegmentFlags*
pArrayOfHighQualitySegmentFlags =
aArrayOfPtrsToArraysOfHighQualitySegmentFlags_[ nSpecies ];
soLineToOutput += " ";
soLineToOutput +=
( (*pArrayOfHighQualitySegmentFlags)[ n0Unpadded ] ? "1" : "0" );
}
if ( pCP->bAutoReportPrintReadPositions_ ) {
for( nSpecies = 0; nSpecies < aArrayOfSpecies_.length();
++nSpecies ) {
RWTValVectorOfUnpaddedReadPositionsInDirectionOfSequencing*
pArrayOfUnpaddedReadPositionsInDirectionOfSequencing =
aArrayOfPtrsToArraysOfReadPositionsInDirectionOfSequencing_[
nSpecies ];
soLineToOutput += " ";
int nReadPositionInDirectionOfSequencing =
(*pArrayOfUnpaddedReadPositionsInDirectionOfSequencing)[ n0Unpadded ];
soLineToOutput += RWCString( (long) nReadPositionInDirectionOfSequencing );
}
}
if ( pCP->bAutoReportPrintChosenReadName_ ) {
for( nSpecies = 0; nSpecies < aArrayOfSpecies_.length();
++nSpecies ) {
RWTPtrVectorOfLocatedFragment* pArrayOfReads =
aArrayOfPtrsToArraysOfLocatedFragments_[ nSpecies ];
LocatedFragment* pLocFrag = (*pArrayOfReads)[ n0Unpadded];
RWCString soLastPartOfReadName;
if ( pLocFrag ) {
soLastPartOfReadName =
pLocFrag->soGetName().soGetLastPart(
pCP->nAutoReportNumbersOfCharactersOfChosenReadNameToBePrinted_ );
}
else {
soLastPartOfReadName = "0";
}
soLineToOutput += " ";
soLineToOutput += soLastPartOfReadName;
}
}
fprintf( pAO, "%s\n", soLineToOutput.data() );
}
fprintf( pAO, "\n} AutoReportPrintSpeciesAlignment\n\n" );
}
void autoReport :: printReadAlignment() {
readReadFile();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
nContig++) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
printReadAlignmentForOneContig( pContig );
}
}
void autoReport :: printReadAlignmentForOneContig( Contig* pContig ) {
// print headings
fprintf( pAO, "\n\nAutoReportPrintReadAlignment {\n" );
fprintf( pAO, "CONTIG: %s\n\n",
pContig->soGetName().data() );
fprintf( pAO, "Columns {\n" );
fprintf( pAO, "prefix (typically chromosome)\n" );
fprintf( pAO, "user-defined position\n" );
fprintf( pAO, "bases (0 if unaligned)\n" );
fprintf( pAO, "qualities (-1 if unaligned)\n" );
fprintf( pAO, "high quality segment flags\n" );
fprintf( pAO, "read position in direction of sequencing (-1 if unaligned)\n" );
fprintf( pAO, "} Columns\n" );
// should print the list of reads
RWTPtrOrderedVector<LocatedFragment> aArrayOfReadsToPrintInThisContig;
findListOfReadsInThisContig( pContig, &aArrayOfReadsToPrintInThisContig );
fprintf( pAO, "Reads {\n" );
int nRead;
for( nRead = 0; nRead < aArrayOfReadsToPrintInThisContig.length();
++nRead ) {
LocatedFragment* pLocFrag = aArrayOfReadsToPrintInThisContig[ nRead ];
fprintf( pAO, "%s\n", pLocFrag->soGetName().data() );
}
fprintf( pAO, "} Reads\n" );
// no point in printing out all of the positions before the first aligned
// read and after the last aligned read
int nFirstAlignedConsPos;
int nLastAlignedConsPos;
bool bFirst = true;
for( nRead = 0; nRead < aArrayOfReadsToPrintInThisContig.length(); ++nRead ) {
LocatedFragment* pLocFrag = aArrayOfReadsToPrintInThisContig[ nRead ];
if ( bFirst ) {
bFirst = false;
nFirstAlignedConsPos = pLocFrag->nGetAlignClipStart();
nLastAlignedConsPos = pLocFrag->nGetAlignClipEnd();
}
else {
nFirstAlignedConsPos = MIN( nFirstAlignedConsPos,
pLocFrag->nGetAlignClipStart() );
nLastAlignedConsPos = MAX( nLastAlignedConsPos,
pLocFrag->nGetAlignClipEnd() );
}
}
RWTValVector<int> aLastReadPosition( aArrayOfReadsToPrintInThisContig.length(), -666 );
aLastReadPosition.soName_ = "autoReport::aLastReadPositions";
RWTValVector<bool> aLastReadPositionSet( aArrayOfReadsToPrintInThisContig.length(), false );
aLastReadPositionSet.soName_ = "autoReport::aLastReadPositionSet";
for( int nConsPos = nFirstAlignedConsPos;
nConsPos <= nLastAlignedConsPos; ++nConsPos ) {
// problem here with the read positions. Need to think of efficient
// way to print them. Perhaps we can keep an array of the current
// read positions for each read. Then as move through each
// consensus position, could go through each read and recover its "state"
// (past read position) and get its new position based on whether
// it has a pad or not and its orientation.
soLine = pCP->soAutoReportPrefix_;
soLine += " ";
soLine +=
soAddCommas(
pContig->nUnpaddedIndexStartingAtUserDefinedPosition( nConsPos ) );
int nRead;
// bases
for( nRead = 0; nRead < aArrayOfReadsToPrintInThisContig.length(); ++nRead ) {
LocatedFragment* pLocFrag = aArrayOfReadsToPrintInThisContig[ nRead ];
if ( pLocFrag->bIsInAlignedPartOfRead( nConsPos ) ) {
Ntide nt = pLocFrag->ntGetFragFromConsPos( nConsPos );
soLine += " ";
soLine += nt.cGetBase();
}
else {
soLine += " 0";
}
}
// qualities
for( nRead = 0; nRead < aArrayOfReadsToPrintInThisContig.length(); ++nRead ) {
LocatedFragment* pLocFrag = aArrayOfReadsToPrintInThisContig[ nRead ];
if ( pLocFrag->bIsInAlignedPartOfRead( nConsPos ) ) {
Ntide nt = pLocFrag->ntGetFragFromConsPos( nConsPos );
soLine += " ";
RWCString soTemp =
RWCString( (long) nt.qualGetQualityWithout98or99() );
soTemp.padOnLeft( 2 );
soLine += soTemp;
}
else {
soLine += " -1";
}
}
// high quality segment
for( nRead = 0; nRead < aArrayOfReadsToPrintInThisContig.length(); ++nRead ) {
LocatedFragment* pLocFrag = aArrayOfReadsToPrintInThisContig[ nRead ];
// can't decide whether this should be aligned as well
if ( pLocFrag->bIsInHighQualitySegmentOfRead( nConsPos )
&&
pLocFrag->bIsInAlignedPartOfRead( nConsPos ) ) {
soLine += " 1";
}
else {
soLine += " 0";
}
}
for( nRead = 0; nRead < aArrayOfReadsToPrintInThisContig.length(); ++nRead ) {
LocatedFragment* pLocFrag = aArrayOfReadsToPrintInThisContig[ nRead ];
if ( pLocFrag->bIsInAlignedPartOfRead( nConsPos ) ) {
Ntide nt = pLocFrag->ntGetFragFromConsPos( nConsPos );
if ( aLastReadPositionSet[ nRead ] ) {
if ( !nt.bIsPad() ) {
if ( pLocFrag->bComp() ) {
--aLastReadPosition[ nRead ];
}
else {
++aLastReadPosition[ nRead ];
}
}
}
else {
aLastReadPositionSet[ nRead ] = true;
aLastReadPosition[ nRead ] =
pLocFrag->nOrientedUnpaddedFragPosFromConsPos( nConsPos );
}
// if this is a pad, then shouldn't we print something
// else other than the most recently printed position?
// I guess that any program reading this file would
// know it was a pad by looking at the read base, but still...
soLine += " ";
soLine += RWCString( (long) aLastReadPosition[ nRead ] );
}
else {
// the read is unaligned at this position
soLine += " -1";
}
} // for( nRead =
fprintf( pAO, "%s\n", soLine.data() );
} // for( int nConsPos = nFirstAlignedConsPos;
fprintf( pAO, "} AutoReportPrintReadAlignment\n\n" );
}
void autoReport :: readReadFile() {
aArrayOfReads_.clear();
FILE* pReadFile = fopen( pCP->filAutoReportPrintTheseReads_.data(),
"r" );
if ( pReadFile == NULL ) {
RWCString soErrorMessage = "could not open consed.autoReportPrintTheseReads: ";
soErrorMessage += pCP->filAutoReportPrintTheseReads_;
soErrorMessage += soGetErrno();
THROW_ERROR( soErrorMessage );
}
while( fgets( soLine.data(), nMaxLineSize, pReadFile ) != NULL ) {
soLine.nCurrentLength_ = strlen( soLine.data() );
soLine.stripAllWhitespaceExceptInternal();
aArrayOfReads_.insert( soLine );
}
fclose( pReadFile );
if ( aArrayOfReads_.length() == 0 ) {
RWCString soErrorMessage = "there are no reads in consed.autoReportPrintTheseReads: ";
soErrorMessage += pCP->filAutoReportPrintTheseReads_;
THROW_ERROR( soErrorMessage );
}
// might want these printed in a particular order
// aArrayOfReads.removeDuplicates();
}
void autoReport :: findListOfReadsInThisContig(
Contig* pContig,
RWTPtrOrderedVector<LocatedFragment>*
pArrayOfReadsInThisContigToPrint ) {
pArrayOfReadsInThisContigToPrint->clear();
for( int nRead = 0; nRead < pContig->nGetNumberOfFragsInContig(); ++nRead ) {
LocatedFragment* pLocFrag = pContig->pLocatedFragmentGet( nRead );
for( int nRead2 = 0; nRead2 < aArrayOfReads_.length(); ++nRead2 ) {
if ( pLocFrag->soGetName() == aArrayOfReads_[ nRead2 ] ) {
pArrayOfReadsInThisContigToPrint->insert( pLocFrag );
}
}
}
}
void autoReport :: printDiscrepantRegions() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
nContig++) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printDiscrepantRegions( );
}
}
void autoReport :: printMinimumQualityHistogram() {
RWTValVector<int> aMinimumQualityHistogram( 101, // capacity
0 ); // initial value
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->addToMinimumQualityHistogram( aMinimumQualityHistogram );
}
fprintf( pAO, "MINIMUM_QUALITY_HISTOGRAM {\n" );
for( int nQuality = 0; nQuality <= 99; ++nQuality ) {
fprintf( pAO, "%d %d\n",
nQuality, aMinimumQualityHistogram[ nQuality ] );
}
fprintf( pAO, "} MINIMUM_QUALITY_HISTOGRAM\n" );
}
void autoReport :: printNumberOfIsolatedPads() {
mbtValVector<int> aIsolatedPrimatePads( 101, 0, 0 );
// nCapacity, nOffset, initialValue
mbtValVector<int> aIsolatedHumanPads( 101, 0, 0 );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->addNumberOfIsolatedPads( aIsolatedPrimatePads,
aIsolatedHumanPads );
}
fprintf( pAO, "ISOLATED_PRIMATE_PADS {\n" );
int nQuality;
for( nQuality = 0; nQuality <= 99; ++nQuality ) {
fprintf( pAO, "%d %d\n", nQuality, aIsolatedPrimatePads[ nQuality ] );
}
fprintf( pAO, "} ISOLATED_PRIMATE_PADS\n" );
fprintf( pAO, "ISOLATED_HUMAN_PADS {\n" );
for( nQuality = 0; nQuality <= 99; ++nQuality ) {
fprintf( pAO, "%d %d\n", nQuality, aIsolatedHumanPads[ nQuality ] );
}
fprintf( pAO, "} ISOLATED_HUMAN_PADS\n" );
}
void autoReport :: printNumberOfIsolatedPadsForEachSpecies() {
RWTValOrderedVector<RWCString> aSpecies;
RWTValOrderedVector<int> aNumberOfIsolatedPadsInPrimateButNotInHuman;
RWTValOrderedVector<int> aNumberOfIsolatedPadsInHumanButNotInPrimate;
RWCString soSpecies;
RWCTokenizer tok( pCP->soAutoReportSpecies_ );
while( !( soSpecies = tok() ).isNull() ) {
aSpecies.insert( soSpecies );
}
aNumberOfIsolatedPadsInPrimateButNotInHuman.clear();
aNumberOfIsolatedPadsInHumanButNotInPrimate.clear();
int nSpecies;
for( nSpecies = 0; nSpecies <= aSpecies.length();
++nSpecies ) {
aNumberOfIsolatedPadsInPrimateButNotInHuman.insert( 0 );
aNumberOfIsolatedPadsInHumanButNotInPrimate.insert( 0 );
}
int nPossibleSites = 0;
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->countPadsForEachSpecies(
aSpecies,
aNumberOfIsolatedPadsInPrimateButNotInHuman,
aNumberOfIsolatedPadsInHumanButNotInPrimate,
nPossibleSites );
}
fprintf( pAO, "POSSIBLE_SITES: %d\n", nPossibleSites );
fprintf( pAO, "ISOLATED_PADS_FOR_EACH_SPECIES {\n" );
for( nSpecies = 0; nSpecies < aSpecies.length(); ++nSpecies ) {
fprintf( pAO, "%s: %d %d %d\n",
aSpecies[ nSpecies ].data(),
aNumberOfIsolatedPadsInPrimateButNotInHuman[ nSpecies ],
aNumberOfIsolatedPadsInHumanButNotInPrimate[ nSpecies ],
aNumberOfIsolatedPadsInPrimateButNotInHuman[ nSpecies ] -
aNumberOfIsolatedPadsInHumanButNotInPrimate[ nSpecies ] );
}
fprintf( pAO, "} ISOLATE_PADS_FOR_EACH_SPECIES\n" );
}
void autoReport :: printBasesInDiscrepantRegions() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printBasesInDiscrepantRegions( this );
}
}
void autoReport :: printAgreeDisagreeBetweenPairsOfSpecies() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printAgreeDisagreeBetweenPairsOfSpecies( this );
}
}
void autoReport :: printAgreeDisagreeBetweenPairsOfSpecies2() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printAgreeDisagreeBetweenPairsOfSpecies2( this );
}
}
void autoReport :: calculateErrorProbabilitiesByComparingPTroPPan() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->calculateErrorProbabilitiesByComparingPTroPPan( this );
}
}
void autoReport :: printScaffolds() {
RWCString soContigMap = ConsEd::pGetAssembly()->soGetContigMap();
soContigMap.translate( ",", "\n" );
fprintf( pAO, "%s\n", soContigMap.data() );
}
void autoReport :: printIfReadsAreCorrectlyAligned() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->areReadsCorrectlyAligned( this );
}
}
void autoReport :: printLengthsOfUnalignedHighQualitySegments() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printLengthsOfUnalignedHighQualitySegments();
}
}
void autoReport :: printLengthsOfAlignedSegments() {
fprintf( pAO, "AlignedSegmentsOfReads{\n" );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printLengthsOfAlignedSegments();
}
fprintf( pAO, "}AlignedSegmentsOfReads\n" );
}
void autoReport :: compareTopAndBottomStrands() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->compareTopAndBottomStrands( this );
}
}
void autoReport :: singleSignalInfo() {
fprintf( pAO, "singleSignalInfo {\n" );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->singleSignalInfo( this );
}
fprintf( pAO, "} singleSignalInfo\n" );
}
void autoReport :: singleSignalInfo2() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->singleSignalInfo2();
}
}
void autoReport :: compareTopAndBottomStrandsWithHuman() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->compareTopAndBottomStrandsWithHuman( this );
}
}
void autoReport :: compareTopAndBottomStrands2() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->compareTopAndBottomStrands2( this );
}
}
void autoReport :: compareTopAndBottomStrands3() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->compareTopAndBottomStrands3( this );
}
}
void autoReport :: compareTopAndBottomStrands4() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->compareTopAndBottomStrands4( this );
}
}
void autoReport :: countColumnsForGroupsOfSpecies() {
parseSpecies();
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->countColumnsForGroupsOfSpecies( this );
}
}
void autoReport :: compareHQSWithLQS() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->compareHQSWithLQS( this );
}
}
void autoReport :: lowQualityBasesInHQS() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->lowQualityBasesInHQS( this );
}
}
void autoReport :: singleSignalOrQuality() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->singleSignalOrQuality( this );
}
}
void autoReport :: discrepancyRateInFlankedRegions() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->discrepancyRateInFlankedRegions( this );
}
}
void autoReport :: discrepancyRateInFlankedRegions2() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->discrepancyRateInFlankedRegions2( this );
}
}
void autoReport :: discrepancyRateInFlankedRegions4() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->discrepancyRateInFlankedRegions4( this );
}
}
void autoReport :: discrepancyRateInFlankedRegions5() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->discrepancyRateInFlankedRegions5( this );
}
}
void autoReport :: countReadsDueToBugInGoodReads() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->countReadsDueToBugInGoodReads();
}
}
void autoReport :: compareTopAndBottomStrandsNoHuman() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->compareTopAndBottomStrandsNoHuman( this );
}
}
void autoReport :: highQualitySegmentData() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->highQualitySegmentData();
}
}
void autoReport :: printFlankedColumns() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printFlankedColumns( this );
}
}
void autoReport :: printFlankedColumns2() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printFlankedColumns2( this );
}
}
void autoReport :: printFlankedColumns3() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printFlankedColumns3( this );
}
}
void autoReport :: printFlankedColumns4() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printFlankedColumns4( this );
}
}
void autoReport :: countAcceptableColumnsWithNoneOnLeft() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->countAcceptableColumnsWithNoneOnLeft( this );
}
}
void autoReport :: countAllMutations() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->countAllMutations( this );
}
}
void autoReport :: countAllMutationsML() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->countAllMutationsML( this );
}
}
void autoReport :: printHighQualityDiscrepancies() {
fprintf( pAO, "highQualityDiscrepancies {\n" );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
highQualityDiscrepancyGotoList myGotoList(
pContig,
pCP->bAutoReportHighQualityDiscrepanciesExcludeCompressionOrG_dropoutTags_,
pCP->bAutoReportHighQualityDiscrepanciesExcludeMostPads_ );
for( int nGotoItem = 0; nGotoItem < myGotoList.nGetNumGotoItems();
++nGotoItem ) {
gotoItem* pGotoItem = myGotoList.pGetGotoItem( nGotoItem );
fprintf( pAO, "%s\n",
pGotoItem->soFullDescription_.data() );
}
}
fprintf( pAO, "} highQualityDiscrepancies\n" );
}
void autoReport :: printLowConsensusQualityRegions() {
fprintf( pAO, "lowConsensusQualityRegions {\n" );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
LowConsQualGotoList myGotoList( pContig );
for( int nGotoItem = 0; nGotoItem < myGotoList.nGetNumGotoItems();
++nGotoItem ) {
gotoItem* pGotoItem = myGotoList.pGetGotoItem( nGotoItem );
fprintf( pAO, "%s\n",
pGotoItem->soFullDescription_.data() );
}
}
fprintf( pAO, "} lowConsensusQualityRegions\n" );
}
void autoReport :: printSingleSubcloneRegions() {
fprintf( pAO, "singleSubcloneRegions {\n" );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
singleSubcloneRegionsGotoList myGotoList( pContig );
for( int nGotoItem = 0; nGotoItem < myGotoList.nGetNumGotoItems();
++nGotoItem ) {
gotoItem* pGotoItem = myGotoList.pGetGotoItem( nGotoItem );
fprintf( pAO, "%s\n",
pGotoItem->soFullDescription_.data() );
}
}
fprintf( pAO, "} singleSubcloneRegions\n" );
}
void autoReport :: printLinkingForwardReversePairs() {
fprintf( pAO, "linkingForwardReversePairs {\n" );
ConsEd::pGetAssembly()->figureOutContigOrderAndOrientation();
fprintf( pAO, "} linkingForwardReversePairs\n" );
}
void autoReport :: printFilteredInconsistentForwardReversePairs() {
Assembly* pAssembly = ConsEd::pGetAssembly();
pAssembly->setContigTemplateArrays();
pAssembly->filterInconsistentFwdRevPairs();
fprintf( pAO, "filteredInconsistentFwdRevPairs {\n" );
fprintf( pAO, "read1 read2 contig1 contig2 (-> or <- for read1) (-> or <- for read2) (why inconsistent)\n" );
for( int nPair = 0;
nPair < pAssembly->pArrayOfConfirmedInconsistentFwdRevPairs_->length();
++nPair ) {
fwdRevPair* pFwdRevPair =
(*(pAssembly->pArrayOfConfirmedInconsistentFwdRevPairs_))[ nPair ];
// how shall we print these?
// read1 read2 contig1 contig2 -> <- problem
RWCString soArrow1;
if ( pFwdRevPair->pLocFrag1_->bComp() ) {
soArrow1 = "<-";
}
else {
soArrow1 = "->";
}
RWCString soArrow2;
if ( pFwdRevPair->pLocFrag2_->bComp() ) {
soArrow2 = "<-";
}
else {
soArrow2 = "->";
}
// convert pFwdRevPair->cProblem_ to text
fwdRevPair2 dummy( pFwdRevPair->pLocFrag1_,
pFwdRevPair->pLocFrag2_,
NULL, // pGctToUse
' ', // cType
pFwdRevPair->cProblem_ );
RWCString soProblem = dummy.soGetProblem();
if ( soProblem.bContains( "->" ) ||
soProblem.bContains( "<-" ) ) {
soProblem = "wrong orientation";
}
fprintf( pAO, "%s %s %s %s %s %s %s\n",
pFwdRevPair->pLocFrag1_->soGetName().data(),
pFwdRevPair->pLocFrag2_->soGetName().data(),
pFwdRevPair->pContigOfRead1_->soGetName().data(),
pFwdRevPair->pContigOfRead2_->soGetName().data(),
soArrow1.data(),
soArrow2.data(),
soProblem.data() );
}
fprintf( pAO, "} filteredInconsistentFwdRevPairs\n" );
}
void autoReport :: printMutationsWithContext() {
fprintf( pAO, "printMutationsWithContext {\n" );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printMutationsWithContext( this );
}
fprintf( pAO, "} printMutationsWithContext\n" );
}
void autoReport :: printCrudeChimpHumanMutations() {
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printCrudeChimpHumanMutations( this );
}
}
void autoReport :: printToCompareToReich() {
fprintf( pAO, "printToCompareToReich {\n" );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
pContig->printToCompareToReich( this );
}
fprintf( pAO, "} printToCompareToReich\n" );
}
void autoReport :: printHighlyDiscrepantRegions() {
fprintf( pAO, "printHighlyDiscrepantRegions {\n" );
RWCString soTitle;
RWCString soFirstLine;
RWCString soSecondLine;
soTitle = "Highly Discrepant Positions";
soFirstLine = "min # of discrepant reads: " +
RWCString( (long) pCP->nNavigateByHighlyDiscrepantPositionsMinDiscrepantReads_ ) +
" min quality: " + RWCString( (long) pCP->nNavigateByHighlyDiscrepantPositionsIgnoreBasesBelowThisQuality_ ) +
" \"r\": base of reference seq";
soSecondLine = "max depth of coverage: " +
RWCString( (long) pCP->nNavigateByHighlyDiscrepantPositionsMaxDepthOfCoverage_ );
if ( pCP->bNavigateByHighlyDiscrepantPositionsJustListIndels_ ) {
soSecondLine += " just indels";
}
soSecondLine += " and ignoring reference seq";
RWCString soThirdLine = " A C G T * pos contig";
fprintf( pAO, "%s\n", soTitle.data() );
fprintf( pAO, "%s\n", soFirstLine.data() );
fprintf( pAO, "%s\n", soSecondLine.data() );
fprintf( pAO, "%s\n", soThirdLine.data() );
for( int nContig = 0; nContig < ConsEd::pGetAssembly()->nNumContigs();
++nContig ) {
Contig* pContig = ConsEd::pGetAssembly()->pGetContig( nContig );
gotoList* pGotoList = new gotoList();
pContig->navigateByHighlyDiscrepantPositions(
pCP->nNavigateByHighlyDiscrepantPositionsMinDiscrepantReads_,
pCP->nNavigateByHighlyDiscrepantPositionsMaxDepthOfCoverage_,
pCP->nNavigateByHighlyDiscrepantPositionsIgnoreBasesBelowThisQuality_,
pCP->bNavigateByHighlyDiscrepantPositionsJustListIndels_,
pCP->bNavigateByHighlyDiscrepantPositionsIgnoreOtherReadsStartingAtSameLocation_,
pCP->bNavigateByHighlyDiscrepantPositionsIgnoreIfListedBasesInConsensus_,
pCP->soNavigateByHighlyDiscrepantPositionsIgnoreIfTheseBasesInConsensus_,
true, // bGotoListNotArrays
pGotoList );
// print out
for( int nGotoItem = 0; nGotoItem < pGotoList->nGetNumGotoItems();
++nGotoItem ) {
gotoItem* pGotoItem = pGotoList->pGetGotoItem( nGotoItem );
fprintf( pAO, "%s\n",
pGotoItem->soFullDescription_.data() );
}
}
fprintf( pAO, "} printHighlyDiscrepantRegions\n" );
}