本文整理汇总了C++中SamFile类的典型用法代码示例。如果您正苦于以下问题：C++ SamFile类的具体用法？C++ SamFile怎么用？C++ SamFile使用的例子？那么, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了SamFile类的15个代码示例，这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞，您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: testReadBam
void testReadBam()
{
SamFile inSam;
assert(inSam.OpenForRead("testFiles/testBam.bam"));
// Call generic test which since the sam and bam are identical, should
// contain the same results.
testRead(inSam);
inSam.Close();
testFlagRead("testFiles/testBam.bam");
}示例2: CalcClusters
void GenomeRegionSeqStats::CalcClusters(String &bamFile, int minMapQuality)
{
SamFile sam;
SamRecord samRecord;
SamFileHeader samHeader;
if(!sam.OpenForRead(bamFile.c_str()))
error("Open BAM file %s failed!\n", bamFile.c_str());
if(!sam.ReadHeader(samHeader)) {
error("Read BAM file header %s failed!\n", bamFile.c_str());
}
if(depth.size()==0) depth.resize(referencegenome.sequenceLength());
String contigLabel;
uint32_t start;
uint32_t gstart;
Reset();
while(sam.ReadRecord(samHeader, samRecord))
{
nReads++;
if(samRecord.getFlag() & SamFlag::UNMAPPED) nUnMapped++;
if(samRecord.getMapQuality() < minMapQuality) continue;
CigarRoller cigar(samRecord.getCigar());
int nonClipSequence = 0;
if(cigar.size()!=0 && cigar[0].operation==Cigar::softClip)
nonClipSequence = cigar[0].count;
contigLabel = samRecord.getReferenceName();
start = nonClipSequence + samRecord.get0BasedPosition(); // start is 0-based
gstart = referencegenome.getGenomePosition(contigLabel.c_str(), start);
if(IsInRegions(contigLabel, start, start+samRecord.getReadLength())) continue;
for(uint32_t i=gstart; i<gstart+samRecord.getReadLength(); i++)
if(depth[i]<MAXDP)
depth[i]++;
nMappedOutTargets++;
}
}示例3: assert
void modify::modifyTags()
{
assert(samIn.OpenForRead(myFilename.c_str()));
// Read the sam header.
assert(samIn.ReadHeader(samHeader));
SamFile samOut;
SamFile bamOut;
std::string inputType = myFilename.substr(myFilename.find_last_of('.'));
std::string outFileBase = "results/updateTagFrom";
if(inputType == ".bam")
{
outFileBase += "Bam";
}
else
{
outFileBase += "Sam";
}
std::string outFile = outFileBase + ".sam";
assert(samOut.OpenForWrite(outFile.c_str()));
outFile = outFileBase + ".bam";
assert(bamOut.OpenForWrite(outFile.c_str()));
assert(samOut.WriteHeader(samHeader));
assert(bamOut.WriteHeader(samHeader));
int count = 0;
// Read the records.
while(samIn.ReadRecord(samHeader, samRecord))
{
if(count == 0)
{
assert(samRecord.rmTag("MD", 'Z'));
}
else if(count == 2)
{
assert(samRecord.rmTags("XT:A;MD:Z;AB:c;NM:i"));
}
else if(count == 4)
{
assert(samRecord.rmTags("MD:Z,AB:c,NM:i"));
}
assert(bamOut.WriteRecord(samHeader, samRecord));
assert(samOut.WriteRecord(samHeader, samRecord));
++count;
}
}示例4: CalcRegionStats
void GenomeRegionSeqStats::CalcRegionStats(String &bamFile)
{
SamFile sam;
SamRecord samRecord;
SamFileHeader samHeader;
if(!sam.OpenForRead(bamFile.c_str()))
error("Open BAM file %s failed!\n", bamFile.c_str());
if(!sam.ReadHeader(samHeader)) {
error("Read BAM file header %s failed!\n", bamFile.c_str());
}
String contigLabel;
int start, end;
Reset();
while(sam.ReadRecord(samHeader, samRecord))
{
nReads++;
if(samRecord.getFlag() & SamFlag::UNMAPPED) nUnMapped++;
if(contigFinishedCnt>=contigs.size()) continue;
CigarRoller cigar(samRecord.getCigar());
int nonClipSequence = 0;
if(cigar.size()!=0 && cigar[0].operation==Cigar::softClip)
nonClipSequence = cigar[0].count;
contigLabel = samRecord.getReferenceName();
start = nonClipSequence + samRecord.get0BasedPosition(); // start is 0-based
end = start + samRecord.getReadLength() - 1;
if(UpdateRegionStats(contigLabel, start, end)) nMapped2Targets++;
}
CalcRegionReadCountInGCBins();
CalcGroupReadCountInGCBins();
std::cout << "Total reads : " << nReads << std::endl;
}示例5: readCoordRecord
SamStatus::Status ClipOverlap::readCoordRecord(SamFile& samIn,
SamRecord** recordPtr,
MateMapByCoord& mateMap,
SamCoordOutput* outputBufferPtr)
{
// Null pointer, so get a new pointer.
if(*recordPtr == NULL)
{
*recordPtr = myPool.getRecord();
if(*recordPtr == NULL)
{
// Failed to allocate a new record.
// Try to free up records from the mate map
if(!forceRecordFlush(mateMap, outputBufferPtr))
{
std::cerr << "Failed to flush the output buffer.\n";
return(SamStatus::FAIL_IO);
}
// Try to get a new record, one should have been cleared.
*recordPtr = myPool.getRecord();
if(*recordPtr == NULL)
{
std::cerr << "Failed to allocate any records.\n";
return(SamStatus::FAIL_MEM);
}
}
}
// RecordPtr is set.
if(!samIn.ReadRecord(mySamHeader, **recordPtr))
{
// Nothing to process, so return.
return(samIn.GetStatus());
}
return(SamStatus::SUCCESS);
}示例6: MEI
/*
if a discordant read is mapped to MEI (overlap with MEI coord)
add centor of ( anchor_end + 3*avr_ins_var )
skip unmap & clip
check 3 types at the same time
*/
void Sites::AddDiscoverSiteFromSingleBam( SingleSampleInfo & si )
{
/*for(int i=0;i<NMEI; i++) {
std::cout << "m=" << i << ": ";
for(map<string, map<int, bool> >::iterator t=meiCoord[m].begin(); t!=meiCoord[m].end(); t++)
std::cout << t->first << "->" << t->second.size() << " ";
std::cout << std::endl;
}*/
avr_read_length = si.avr_read_length;
avr_ins_size = si.avr_ins_size;
min_read_length = avr_read_length / 3;
current_depth = si.depth;
// total_depth += current_depth;
resetDepthAddFlag();
SamFile bam;
SamFileHeader bam_header;
OpenBamAndBai( bam,bam_header, si.bam_name );
for( int i=0; i<pchr_list->size(); i++ ) {
string chr = (*pchr_list)[i];
// if ( !single_chr.empty() && chr.compare(single_chr)!=0 )
// continue;
if ( siteList.find(chr) == siteList.end() )
siteList[chr].clear();
// map<string, map<int, SingleSite> >::iterator pchr = siteList[m].find(chr);
// map<string, map<int, bool> >::iterator coord_chr_ptr = meiCoord[m].find(chr);
// if (coord_chr_ptr == meiCoord[m].end())
// continue;
bool section_status;
if (range_start<0) { // no range
section_status = bam.SetReadSection( chr.c_str() );
if (!section_status) {
string str = "Cannot set read section at chr " + chr;
morphWarning( str );
}
}
else { // set range
section_status = bam.SetReadSection( chr.c_str(), range_start, range_end );
if (!section_status) {
string str = "Cannot set read section at chr " + chr + " " + std::to_string(range_start) + "-" + std::to_string(range_end);
morphWarning( str );
}
}
// DO ADDING
// if (siteList[chr].empty())
// p_reach_last = 1;
// else {
// p_reach_last = 0;
pnearest = siteList[chr].begin();
// }
SingleSite new_site; // temporary cluster. will be added to map later.
new_site.depth = current_depth;
bool start_new = 1; // check if need to add new_site to map and start new new_site
SamRecord rec;
int between = 0; // count #reads after new_site.end. If end changed, add it to rcount and reset to zero
while( bam.ReadRecord( bam_header, rec ) ) {
if (!start_new) {
if (rec.get1BasedPosition() >= new_site.end)
between++;
else
new_site.rcount++;
}
if (rec.getFlag() & 0x2)
continue;
if ( OneEndUnmap( rec.getFlag() ) )
continue;
if ( IsSupplementary(rec.getFlag()) )
continue;
if ( rec.getReadLength() < min_read_length )
continue;
if ( rec.getMapQuality() < MIN_QUALITY )
continue;
if (chr.compare(rec.getMateReferenceName())==0 && rec.getInsertSize() < abs(avr_ins_size*2))
continue;
bool is_mei = 0;
vector<bool> is_in_coord;
is_in_coord.resize(3, 0);
for(int m=0; m<NMEI; m++) {
map<string, map<int, bool> >::iterator coord_chr_ptr = meiCoord[m].find(rec.getMateReferenceName());
if (coord_chr_ptr == meiCoord[m].end())
is_in_coord[m] = 0;
else
is_in_coord[m] = isWithinCoord( rec.get1BasedMatePosition(), coord_chr_ptr->second ); // within MEI coord
if (is_in_coord[m])
is_mei = 1;
}
if (!is_mei)
continue;
if (start_new) {
setNewCluster( is_in_coord, new_site,rec);
start_new = 0;
//.........这里部分代码省略.........
示例7: setReadCountInSection
void setReadCountInSection( vector<int> & raw_counts, string & chr, int center, SamFile & samIn, SamFileHeader & samHeader, vector<RefSeq*> & REF_SEQ )
{
int st = center - WIN/2;
int ed = center + WIN/2;
bool section_status = samIn.SetReadSection( chr.c_str(), st, ed );
if (!section_status) {
std::cerr << "Warning: Unable to set read section: " << chr << ": " << st << "-" << ed << ". Set section depth = 0!" << std::endl;
return;
}
// proper reads
map<string, vector< DiscPair > > disc_recs; // record where the disc come from
ProperDeck pDeck( REF_SEQ );
SamRecord sam_rec;
while( samIn.ReadRecord(samHeader, sam_rec) ) {
bool pass_qc = PassQC( sam_rec );
if ( !pass_qc )
continue;
if ( sam_rec.getFlag() & 0x2 ) {
if ( sam_rec.getInsertSize() > 0 )
pDeck.Add( sam_rec );
else {
RetrievedIndex rv = pDeck.RetrieveIndex( sam_rec );
int index = getRetrievedIndexToRawCounts( rv );
if (index >= 0) { // for read partially in window, only add clip
if ( sam_rec.get1BasedPosition() < st || sam_rec.get1BasedAlignmentEnd() > ed ) {
if ( index >= 2 )
raw_counts[ index ]++;
}
else
raw_counts[ index ]++;
}
}
}
else { // disc: rec info and wait to reset section later
if ( sam_rec.getReadLength() < 30 )
continue;
string mate_chr = sam_rec.getMateReferenceName();
// check if this one is valid as anchor
DiscPair new_pair( 1, 0, sam_rec, REF_SEQ );
disc_recs[mate_chr].push_back( new_pair );
}
}
// disc reads
for( map<string, vector< DiscPair > >::iterator chr_it = disc_recs.begin(); chr_it != disc_recs.end(); chr_it++ ) {
for( vector< DiscPair >::iterator dp_it = chr_it->second.begin(); dp_it != chr_it->second.end(); dp_it++ ) {
bool section_status = samIn.SetReadSection( chr_it->first.c_str(), dp_it->GetSecondAlignPosition(), dp_it->GetSecondAlignPosition() + WIN );
if (!section_status) {
std::cerr << "ERROR: Unable to set read section: " << chr << ": " << st << "-" << ed << std::endl;
exit(1);
}
SamRecord sam_rec;
while( samIn.ReadRecord(samHeader, sam_rec) ) {
bool pass_qc = PassQC( sam_rec );
if ( !pass_qc )
continue;
if ( !DiscSamPass(sam_rec) )
continue;
int position = sam_rec.get1BasedPosition();
if ( position > dp_it->GetFirstAlignPosition() )
break;
if ( sam_rec.getFlag() & 0x2 )
continue;
bool same_pair = dp_it->IsSamePair( sam_rec );
if ( !same_pair )
continue;
// now add ro raw stats: always use first as anchor
dp_it->AddSecondToPair( sam_rec );
Loci loci = dp_it->GetFirstLoci();
if ( dp_it->FirstValid() ) {
int index = getLociToRawCounts( loci );
raw_counts[ index ]++;
// clear & break
break;
}
}
}
}
}示例8: BEGIN_LONG_PARAMETERS
int Convert::execute(int argc, char **argv)
{
// Extract command line arguments.
String inFile = "";
String outFile = "";
String refFile = "";
bool lshift = false;
bool noeof = false;
bool params = false;
bool useBases = false;
bool useEquals = false;
bool useOrigSeq = false;
bool recover = false;
ParameterList inputParameters;
BEGIN_LONG_PARAMETERS(longParameterList)
LONG_STRINGPARAMETER("in", &inFile)
LONG_STRINGPARAMETER("out", &outFile)
LONG_STRINGPARAMETER("refFile", &refFile)
LONG_PARAMETER("lshift", &lshift)
LONG_PARAMETER("noeof", &noeof)
LONG_PARAMETER("recover", &recover)
LONG_PARAMETER("params", ¶ms)
LONG_PARAMETER_GROUP("SequenceConversion")
EXCLUSIVE_PARAMETER("useBases", &useBases)
EXCLUSIVE_PARAMETER("useEquals", &useEquals)
EXCLUSIVE_PARAMETER("useOrigSeq", &useOrigSeq)
LONG_PHONEHOME(VERSION)
END_LONG_PARAMETERS();
inputParameters.Add(new LongParameters ("Input Parameters",
longParameterList));
// parameters start at index 2 rather than 1.
inputParameters.Read(argc, argv, 2);
// If no eof block is required for a bgzf file, set the bgzf file type to
// not look for it.
if(noeof)
{
// Set that the eof block is not required.
BgzfFileType::setRequireEofBlock(false);
}
// Check to see if the in file was specified, if not, report an error.
if(inFile == "")
{
printUsage(std::cerr);
inputParameters.Status();
// In file was not specified but it is mandatory.
std::cerr << "--in is a mandatory argument, "
<< "but was not specified" << std::endl;
return(-1);
}
if(outFile == "")
{
printUsage(std::cerr);
inputParameters.Status();
// In file was not specified but it is mandatory.
std::cerr << "--out is a mandatory argument, "
<< "but was not specified" << std::endl;
return(-1);
}
// Check to see if the ref file was specified.
// Open the reference.
GenomeSequence* refPtr = NULL;
if(refFile != "")
{
refPtr = new GenomeSequence(refFile);
}
SamRecord::SequenceTranslation translation;
if((useBases) && (refPtr != NULL))
{
translation = SamRecord::BASES;
}
else if((useEquals) && (refPtr != NULL))
{
translation = SamRecord::EQUAL;
}
else
{
useOrigSeq = true;
translation = SamRecord::NONE;
}
if(params)
{
inputParameters.Status();
}
// Open the input file for reading.
SamFile samIn;
if(recover) samIn.setAttemptRecovery(true);
samIn.OpenForRead(inFile);
//.........这里部分代码省略.........
示例9: main
//.........这里部分代码省略.........
else {
gpLogger->write_log("\t--RG [%s]",vsRGHeaders[0].c_str());
}
if ( vsPGHeaders.empty() ) {
gpLogger->write_log("\t--PG []");
}
else {
for(uint32_t i=0; i < vsPGHeaders.size(); ++i) {
gpLogger->write_log("\t--PG [%s]",vsPGHeaders[i].c_str());
}
}
if ( (vsHDHeaders.empty() ) && ( vsRGHeaders.empty() ) && ( vsPGHeaders.empty() ) && ( !bClear ) && ( sFasta.empty() ) ) {
gpLogger->warning("No option is in effect for modifying BAM files. The input and output files will be identical");
}
if ( ( vsHDHeaders.size() > 1 ) || ( vsRGHeaders.size() > 1 ) ) {
gpLogger->error("HD and RG headers cannot be multiple");
}
FastaFile fastaFile;
if ( ! sFasta.empty() ) {
if ( fastaFile.open(sFasta.c_str()) ) {
gpLogger->write_log("Reading the reference file %s",sFasta.c_str());
fastaFile.readThru();
fastaFile.close();
gpLogger->write_log("Finished reading the reference file %s",sFasta.c_str());
}
else {
gpLogger->error("Failed to open reference file %s",sFasta.c_str());
}
}
SamFile samIn;
SamFile samOut;
if ( ! samIn.OpenForRead(sInFile.c_str()) ) {
gpLogger->error("Cannot open BAM file %s for reading - %s",sInFile.c_str(), SamStatus::getStatusString(samIn.GetStatus()) );
}
if ( ! samOut.OpenForWrite(sOutFile.c_str()) ) {
gpLogger->error("Cannot open BAM file %s for writing - %s",sOutFile.c_str(), SamStatus::getStatusString(samOut.GetStatus()) );
}
SamFileHeader samHeader;
SamHeaderRecord* pSamHeaderRecord;
samIn.ReadHeader(samHeader);
// check the sanity of SQ file
// make sure the SN and LN matches, with the same order
if ( bCheckSQ ) {
unsigned int numSQ = 0;
while( (pSamHeaderRecord = samHeader.getNextHeaderRecord()) != NULL ) {
if ( pSamHeaderRecord->getType() == SamHeaderRecord::SQ ) {
++numSQ;
}
}
if ( numSQ != fastaFile.vsSequenceNames.size() ) {
gpLogger->error("# of @SQ tags are different from the original BAM and the reference file");
}
// iterator over all @SQ objects
for(unsigned int i=0; i < numSQ; ++i) {
pSamHeaderRecord = samHeader.getSQ(fastaFile.vsSequenceNames[i].c_str());
if ( fastaFile.vsSequenceNames[i].compare(pSamHeaderRecord->getTagValue("SN")) != 0 ) {
gpLogger->error("SequenceName is not identical between fasta and input BAM file");示例10: usage
// main function
int TrimBam::execute(int argc, char ** argv)
{
SamFile samIn;
SamFile samOut;
int numTrimBaseL = 0;
int numTrimBaseR = 0;
bool noeof = false;
bool ignoreStrand = false;
bool noPhoneHome = false;
std::string inName = "";
std::string outName = "";
if ( argc < 5 ) {
usage();
std::cerr << "ERROR: Incorrect number of parameters specified\n";
return(-1);
}
inName = argv[2];
outName = argv[3];
static struct option getopt_long_options[] = {
// Input options
{ "left", required_argument, NULL, 'L'},
{ "right", required_argument, NULL, 'R'},
{ "ignoreStrand", no_argument, NULL, 'i'},
{ "noeof", no_argument, NULL, 'n'},
{ "noPhoneHome", no_argument, NULL, 'p'},
{ "nophonehome", no_argument, NULL, 'P'},
{ "phoneHomeThinning", required_argument, NULL, 't'},
{ "phonehomethinning", required_argument, NULL, 'T'},
{ NULL, 0, NULL, 0 },
};
int argIndex = 4;
if(argv[argIndex][0] != '-')
{
// This is the number of bases to trim off both sides
// so convert to a number.
numTrimBaseL = atoi(argv[argIndex]);
numTrimBaseR = numTrimBaseL;
++argIndex;
}
int c = 0;
int n_option_index = 0;
// Process any additional parameters
while ( ( c = getopt_long(argc, argv,
"L:R:in", getopt_long_options, &n_option_index) )
!= -1 )
{
switch(c)
{
case 'L':
numTrimBaseL = atoi(optarg);
break;
case 'R':
numTrimBaseR = atoi(optarg);
break;
case 'i':
ignoreStrand = true;
break;
case 'n':
noeof = true;
break;
case 'p':
case 'P':
noPhoneHome = true;
break;
case 't':
case 'T':
PhoneHome::allThinning = atoi(optarg);
break;
default:
fprintf(stderr,"ERROR: Unrecognized option %s\n",
getopt_long_options[n_option_index].name);
return(-1);
}
}
if(!noPhoneHome)
{
PhoneHome::checkVersion(getProgramName(), VERSION);
}
if(noeof)
{
// Set that the eof block is not required.
BgzfFileType::setRequireEofBlock(false);
}
if ( ! samIn.OpenForRead(inName.c_str()) ) {
fprintf(stderr, "***Problem opening %s\n",inName.c_str());
return(-1);
}
if(!samOut.OpenForWrite(outName.c_str())) {
fprintf(stderr, "%s\n", samOut.GetStatusMessage());
return(samOut.GetStatus());
}
//.........这里部分代码省略.........
示例11: processFile
int GapInfo::processFile(const char* inputFileName, const char* outputFileName,
const char* refFile, bool detailed,
bool checkFirst, bool checkStrand)
{
// Open the file for reading.
SamFile samIn;
samIn.OpenForRead(inputFileName);
// Read the sam header.
SamFileHeader samHeader;
samIn.ReadHeader(samHeader);
SamRecord samRecord;
GenomeSequence* refPtr = NULL;
if(strcmp(refFile, "") != 0)
{
refPtr = new GenomeSequence(refFile);
}
IFILE outFile = ifopen(outputFileName, "w");
// Map for summary.
std::map<int, int> gapInfoMap;
// Keep reading records until ReadRecord returns false.
while(samIn.ReadRecord(samHeader, samRecord))
{
uint16_t samFlags = samRecord.getFlag();
if((!SamFlag::isMapped(samFlags)) ||
(!SamFlag::isMateMapped(samFlags)) ||
(!SamFlag::isPaired(samFlags)) ||
(samFlags & SamFlag::SECONDARY_ALIGNMENT) ||
(SamFlag::isDuplicate(samFlags)) ||
(SamFlag::isQCFailure(samFlags)))
{
// unmapped, mate unmapped, not paired,
// not the primary alignment,
// duplicate, fails vendor quality check
continue;
}
// No gap info if the chromosome names are different or
// are unknown.
int32_t refID = samRecord.getReferenceID();
if((refID != samRecord.getMateReferenceID()) || (refID == -1))
{
continue;
}
int32_t readStart = samRecord.get0BasedPosition();
int32_t mateStart = samRecord.get0BasedMatePosition();
// If the mate starts first, then the pair was processed by
// the mate.
if(mateStart < readStart)
{
continue;
}
if((mateStart == readStart) && (SamFlag::isReverse(samFlags)))
{
// read and mate start at the same position, so
// only process the forward strand.
continue;
}
// Process this read pair.
int32_t readEnd = samRecord.get0BasedAlignmentEnd();
int32_t gapSize = mateStart - readEnd - 1;
if(detailed)
{
// Output the gap info.
ifprintf(outFile, "%s\t%d\t%d",
samRecord.getReferenceName(), readEnd+1, gapSize);
// Check if it is not the first or if it is not the forward strand.
if(checkFirst && !SamFlag::isFirstFragment(samFlags))
{
ifprintf(outFile, "\tNotFirst");
}
if(checkStrand && SamFlag::isReverse(samFlags))
{
ifprintf(outFile, "\tReverse");
}
ifprintf(outFile, "\n");
}
else
{
// Summary.
// Skip reads that are not the forward strand.
if(SamFlag::isReverse(samFlags))
{
// continue
continue;
}
//.........这里部分代码省略.........
示例12: BEGIN_LONG_PARAMETERS
//.........这里部分代码省略.........
if(myRefID != UNSET_REF && myRefName.Length() != 0)
{
std::cerr << "Can't specify both refID and refName" << std::endl;
inputParameters.Status();
return(-1);
}
if(myRefID != UNSET_REF && bed.Length() != 0)
{
std::cerr << "Can't specify both refID and bed" << std::endl;
inputParameters.Status();
return(-1);
}
if(myRefName.Length() != 0 && bed.Length() != 0)
{
std::cerr << "Can't specify both refName and bed" << std::endl;
inputParameters.Status();
return(-1);
}
if(!bed.IsEmpty())
{
myBedFile = ifopen(bed, "r");
}
if(params)
{
inputParameters.Status();
}
// Open the file for reading.
mySamIn.OpenForRead(inFile);
// Open the output file for writing.
SamFile samOut;
samOut.OpenForWrite(outFile);
// Open the bam index file for reading if a region was specified.
if((myRefName.Length() != 0) || (myRefID != UNSET_REF) || (myBedFile != NULL))
{
mySamIn.ReadBamIndex(indexFile);
}
// Read & write the sam header.
mySamIn.ReadHeader(mySamHeader);
samOut.WriteHeader(mySamHeader);
// Read the sam records.
SamRecord samRecord;
// Track the status.
int numSectionRecords = 0;
// Set returnStatus to success. It will be changed
// to the failure reason if any of the writes fail.
SamStatus::Status returnStatus = SamStatus::SUCCESS;
while(getNextSection())
{
// Keep reading records until they aren't anymore.
while(mySamIn.ReadRecord(mySamHeader, samRecord))
{
if(!readName.IsEmpty())
{
// Check for readname.
if(strcmp(samRecord.getReadName(), readName.c_str()) != 0)
{
// not a matching read name, so continue to the next record.示例13: testIndex
void testIndex(BamIndex& bamIndex)
{
assert(bamIndex.getNumMappedReads(1) == 2);
assert(bamIndex.getNumUnMappedReads(1) == 0);
assert(bamIndex.getNumMappedReads(0) == 4);
assert(bamIndex.getNumUnMappedReads(0) == 1);
assert(bamIndex.getNumMappedReads(23) == -1);
assert(bamIndex.getNumUnMappedReads(23) == -1);
assert(bamIndex.getNumMappedReads(-1) == 0);
assert(bamIndex.getNumUnMappedReads(-1) == 2);
assert(bamIndex.getNumMappedReads(-2) == -1);
assert(bamIndex.getNumUnMappedReads(-2) == -1);
assert(bamIndex.getNumMappedReads(22) == 0);
assert(bamIndex.getNumUnMappedReads(22) == 0);
// Get the chunks for reference id 1.
Chunk testChunk;
SortedChunkList chunkList;
assert(bamIndex.getChunksForRegion(1, -1, -1, chunkList) == true);
assert(!chunkList.empty());
testChunk = chunkList.pop();
assert(chunkList.empty());
assert(testChunk.chunk_beg == 0x4e7);
assert(testChunk.chunk_end == 0x599);
// Get the chunks for reference id 0.
assert(bamIndex.getChunksForRegion(0, -1, -1, chunkList) == true);
assert(!chunkList.empty());
testChunk = chunkList.pop();
assert(chunkList.empty());
assert(testChunk.chunk_beg == 0x360);
assert(testChunk.chunk_end == 0x4e7);
// Get the chunks for reference id 2.
assert(bamIndex.getChunksForRegion(2, -1, -1, chunkList) == true);
assert(!chunkList.empty());
testChunk = chunkList.pop();
assert(chunkList.empty());
assert(testChunk.chunk_beg == 0x599);
assert(testChunk.chunk_end == 0x5ea);
// Get the chunks for reference id 3.
// There isn't one for this ref id, but still successfully read the file,
// so it should return true, but the list should be empty.
assert(bamIndex.getChunksForRegion(3, -1, -1, chunkList) == true);
assert(chunkList.empty());
// Test reading an indexed bam file.
SamFile inFile;
assert(inFile.OpenForRead("testFiles/sortedBam.bam"));
inFile.setSortedValidation(SamFile::COORDINATE);
assert(inFile.ReadBamIndex("testFiles/sortedBam.bam.bai"));
SamFileHeader samHeader;
assert(inFile.ReadHeader(samHeader));
SamRecord samRecord;
// Test getting num mapped/unmapped reads.
assert(inFile.getNumMappedReadsFromIndex(1) == 2);
assert(inFile.getNumUnMappedReadsFromIndex(1) == 0);
assert(inFile.getNumMappedReadsFromIndex(0) == 4);
assert(inFile.getNumUnMappedReadsFromIndex(0) == 1);
assert(inFile.getNumMappedReadsFromIndex(23) == -1);
assert(inFile.getNumUnMappedReadsFromIndex(23) == -1);
assert(inFile.getNumMappedReadsFromIndex(-1) == 0);
assert(inFile.getNumUnMappedReadsFromIndex(-1) == 2);
assert(inFile.getNumMappedReadsFromIndex(-2) == -1);
assert(inFile.getNumUnMappedReadsFromIndex(-2) == -1);
assert(inFile.getNumMappedReadsFromIndex(22) == 0);
assert(inFile.getNumUnMappedReadsFromIndex(22) == 0);
assert(inFile.getNumMappedReadsFromIndex("2", samHeader) == 2);
assert(inFile.getNumUnMappedReadsFromIndex("2", samHeader) == 0);
assert(inFile.getNumMappedReadsFromIndex("1", samHeader) == 4);
assert(inFile.getNumUnMappedReadsFromIndex("1", samHeader) == 1);
assert(inFile.getNumMappedReadsFromIndex("22", samHeader) == 0);
assert(inFile.getNumUnMappedReadsFromIndex("22", samHeader) == 0);
assert(inFile.getNumMappedReadsFromIndex("", samHeader) == 0);
assert(inFile.getNumUnMappedReadsFromIndex("*", samHeader) == 2);
assert(inFile.getNumMappedReadsFromIndex("unknown", samHeader) == -1);
assert(inFile.getNumUnMappedReadsFromIndex("unknown", samHeader) == -1);
assert(inFile.getNumMappedReadsFromIndex("X", samHeader) == 0);
assert(inFile.getNumUnMappedReadsFromIndex("X", samHeader) == 0);
// Section -1 = Ref *: 2 records (8 & 10 from testSam.sam that is reflected
// in the validation.
assert(inFile.SetReadSection(-1));
assert(inFile.ReadRecord(samHeader, samRecord));
validateRead8(samRecord);
assert(inFile.ReadRecord(samHeader, samRecord));
validateRead10(samRecord);
assert(inFile.ReadRecord(samHeader, samRecord) == false);
// Section 2 = Ref 3: 1 records (9 from testSam.sam that is reflected
// in the validation.
assert(inFile.SetReadSection(2));
assert(inFile.ReadRecord(samHeader, samRecord));
validateRead9(samRecord);
assert(inFile.ReadRecord(samHeader, samRecord) == false);
//.........这里部分代码省略.........
示例14: testRead
void testRead(SamFile &inSam)
{
// Read the SAM Header.
SamFileHeader samHeader;
assert(inSam.ReadHeader(samHeader));
validateHeader(samHeader);
testCopyHeader(samHeader);
testModHeader(samHeader);
SamRecord samRecord;
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead1(samRecord);
// Set a new quality and get the buffer.
samRecord.setQuality("ABCDE");
validateRead1ModQuality(samRecord);
// void* buffer = samRecord.getRecordBuffer();
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead2(samRecord);
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead3(samRecord);
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead4(samRecord);
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead5(samRecord);
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead6(samRecord);
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead7(samRecord);
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead8(samRecord);
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead9(samRecord);
assert(inSam.ReadRecord(samHeader, samRecord) == true);
validateRead10(samRecord);
}示例15: handleSortedByReadName
SamStatus::Status ClipOverlap::handleSortedByReadName(SamFile& samIn,
SamFile* samOutPtr)
{
// Set returnStatus to success. It will be changed
// to the failure reason if any of the writes fail.
SamStatus::Status returnStatus = SamStatus::SUCCESS;
// Read the sam records.
SamRecord* prevSamRecord = NULL;
SamRecord* samRecord = new SamRecord;
SamRecord* tmpRecord = new SamRecord;
if((samRecord == NULL) || (tmpRecord == NULL))
{
std::cerr << "Failed to allocate a SamRecord, so exit.\n";
return(SamStatus::FAIL_MEM);
}
// Keep reading records until ReadRecord returns false.
while(samIn.ReadRecord(mySamHeader, *samRecord))
{
int16_t flag = samRecord->getFlag();
if((flag & myIntExcludeFlags) != 0)
{
// This read should not be checked for overlaps.
// Check if there is a previous SamRecord.
if(prevSamRecord != NULL)
{
// There is a previous record.
// If it has a different read name, write it.
if(strcmp(samRecord->getReadName(),
prevSamRecord->getReadName()) != 0)
{
// Different read name, so write the previous record.
if((samOutPtr != NULL) && !myOverlapsOnly)
{
if(!samOutPtr->WriteRecord(mySamHeader, *prevSamRecord))
{
// Failed to write a record.
fprintf(stderr, "%s\n", samOutPtr->GetStatusMessage());
returnStatus = samOutPtr->GetStatus();
}
}
// Clear the previous record info.
tmpRecord = prevSamRecord;
prevSamRecord = NULL;
}
// If it has the same read name, leave it in case there is another read with that name
}
// This record is not being checked for overlaps, so just write it and continue
if((samOutPtr != NULL) && !myOverlapsOnly)
{
if(!samOutPtr->WriteRecord(mySamHeader, *samRecord))
{
// Failed to write a record.
fprintf(stderr, "%s\n", samOutPtr->GetStatusMessage());
returnStatus = samOutPtr->GetStatus();
}
}
continue;
}
if(prevSamRecord == NULL)
{
// Nothing to compare this record to, so set this
// record to the previous, and the next record.
prevSamRecord = samRecord;
samRecord = tmpRecord;
tmpRecord = NULL;
continue;
}
// Check if the read name matches the previous read name.
if(strcmp(samRecord->getReadName(),
prevSamRecord->getReadName()) == 0)
{
bool overlap = false;
// Same Read Name, so check clipping.
OverlapHandler::OverlapInfo prevClipInfo =
myOverlapHandler->getOverlapInfo(*prevSamRecord);
OverlapHandler::OverlapInfo curClipInfo =
myOverlapHandler->getOverlapInfo(*samRecord);
// If either indicate a complete clipping, clip both.
if((prevClipInfo == OverlapHandler::NO_OVERLAP_WRONG_ORIENT) ||
(curClipInfo == OverlapHandler::NO_OVERLAP_WRONG_ORIENT))
{
overlap = true;
myOverlapHandler->handleNoOverlapWrongOrientation(*prevSamRecord);
// Don't update stats since this is the 2nd in the pair
myOverlapHandler->handleNoOverlapWrongOrientation(*samRecord,
false);
}
else if((prevClipInfo == OverlapHandler::OVERLAP) ||
(prevClipInfo == OverlapHandler::SAME_START))
{
// The previous read starts at or before the current one.
overlap = true;
myOverlapHandler->handleOverlapPair(*prevSamRecord,
*samRecord);
//.........这里部分代码省略.........