C++ BamReader::Open方法代码示例

// opens BAM files
bool BamMultiReader::Open(const vector<string> filenames, bool openIndexes, bool coreMode, bool useDefaultIndex) {

    // for filename in filenames
    fileNames = filenames; // save filenames in our multireader
    for (vector<string>::const_iterator it = filenames.begin(); it != filenames.end(); ++it) {
        string filename = *it;
        BamReader* reader = new BamReader;

        bool openedOK = true;
        if (openIndexes) {
            if (useDefaultIndex)
                openedOK = reader->Open(filename, filename + ".bai");
            else
                openedOK = reader->Open(filename, filename + ".bti");
        } else {
            openedOK = reader->Open(filename); // for merging, jumping is disallowed
        }

        // if file opened ok, check that it can be read
        if ( openedOK ) {

            bool fileOK = true;
            BamAlignment* alignment = new BamAlignment;
            if (coreMode) {
                fileOK &= reader->GetNextAlignmentCore(*alignment);
            } else {
                fileOK &= reader->GetNextAlignment(*alignment);
            }

            if (fileOK) {
                readers.push_back(make_pair(reader, alignment)); // store pointers to our readers for cleanup
                alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position),
                                            make_pair(reader, alignment)));
            } else {
                cerr << "WARNING: could not read first alignment in " << filename << ", ignoring file" << endl;
                // if only file available & could not be read, return failure
                if ( filenames.size() == 1 ) return false;
            }

        }

        // TODO; any more error handling on openedOK ??
        else
            return false;
    }

    // files opened ok, at least one alignment could be read,
    // now need to check that all files use same reference data
    ValidateReaders();
    return true;
}

void BamToFastq::SingleFastq() {
    // open the 1st fastq file for writing
    ofstream fq(_fastq1.c_str(), ios::out);
    if ( !fq ) {
        cerr << "Error: The first fastq file (" << _fastq1 << ") could not be opened.  Exiting!" << endl;
        exit (1);
    }
    // open the BAM file
    BamReader reader;
    reader.Open(_bamFile);
    BamAlignment bam;
    while (reader.GetNextAlignment(bam)) {
        // extract the sequence and qualities for the BAM "query"
        string seq  = bam.QueryBases;
        string qual = bam.Qualities;
        if (bam.IsReverseStrand() == true) {
            reverseComplement(seq);
            reverseSequence(qual);
        }
        fq << "@" << bam.Name << endl;
        fq << seq << endl;
        fq << "+" << endl;
        fq << qual << endl;
    }
}

bool RevertTool::RevertToolPrivate::Run(void) {
  
    // opens the BAM file without checking for indexes
    BamReader reader;
    if ( !reader.Open(m_settings->InputFilename) ) {
        cerr << "Could not open input BAM file... quitting." << endl;
        return false;
    }

    // get BAM file metadata
    const string& headerText = reader.GetHeaderText();
    const RefVector& references = reader.GetReferenceData();
    
    // open writer
    BamWriter writer;
    bool writeUncompressed = ( m_settings->OutputFilename == Options::StandardOut() && !m_settings->IsForceCompression );
    if ( !writer.Open(m_settings->OutputFilename, headerText, references, writeUncompressed) ) {
        cerr << "Could not open " << m_settings->OutputFilename << " for writing." << endl;
        return false;
    }

    // plow through file, reverting alignments
    BamAlignment al;
    while ( reader.GetNextAlignment(al) ) {
        RevertAlignment(al);
        writer.SaveAlignment(al);
    }
    
    // clean and exit
    reader.Close();
    writer.Close();
    return true; 
}

void BedIntersectPE::IntersectBamPE(string bamFile) {

    // load the "B" bed file into a map so
    // that we can easily compare "A" to it for overlaps
    _bedB->loadBedFileIntoMap();

    // open the BAM file
    BamReader reader;
    BamWriter writer;
    reader.Open(bamFile);

    // get header & reference information
    string bamHeader = reader.GetHeaderText();
    RefVector refs   = reader.GetReferenceData();

    // open a BAM output to stdout if we are writing BAM
    if (_bamOutput == true) {
        // set compression mode
        BamWriter::CompressionMode compressionMode = BamWriter::Compressed;
        if ( _isUncompressedBam ) compressionMode = BamWriter::Uncompressed;
        writer.SetCompressionMode(compressionMode);
        // open our BAM writer
        writer.Open("stdout", bamHeader, refs);
    }

    // track the previous and current sequence
    // names so that we can identify blocks of
    // alignments for a given read ID.
    string prevName, currName;
    prevName = currName = "";

    vector<BamAlignment> alignments;        // vector of BAM alignments for a given ID in a BAM file.
    alignments.reserve(100);

    _bedA->bedType = 10;                    // it's a full BEDPE given it's BAM

    // rip through the BAM file and convert each mapped entry to BEDPE
    BamAlignment bam1, bam2;
    while (reader.GetNextAlignment(bam1)) {
        // the alignment must be paired
        if (bam1.IsPaired() == true) {
            // grab the second alignment for the pair.
            reader.GetNextAlignment(bam2);

            // require that the alignments are from the same query
            if (bam1.Name == bam2.Name) {
                ProcessBamBlock(bam1, bam2, refs, writer);
            }
            else {
                cerr << "*****ERROR: -bedpe requires BAM to be sorted or grouped by query name. " << endl;
                exit(1);
            }
        }
    }
    // close up
    reader.Close();
    if (_bamOutput == true) {
        writer.Close();
    }
}

int main (int argc, char *argv[]) {
    
     string bamfiletopen = string(argv[1]);
     BamReader reader;

     if ( !reader.Open(bamfiletopen) ) {
    	cerr << "Could not open input BAM files." << endl;
    	return 1;
     }

    BamAlignment al;
    while ( reader.GetNextAlignment(al) ) {
	string reconstructedReference = reconstructRef(&al);
	cout<<al.QueryBases<<endl;
	cout<<reconstructedReference<<endl;

	
	pair< string, vector<int> >  reconP = reconstructRefWithPos(&al);
	for(unsigned int i=0;i<reconP.first.size();i++){
	    cout<<reconP.first[i]<<"\t"<<reconP.second[i]<<endl;
	}

    }
    reader.Close();

    return 0;
}

/*
Description:
	Load all the bam into memory at one time if no parameters set, otherwise load the needed part of the bam.
	Save the parsed info into vector. 
*/
bool BamParse::parseAlignment(int chrom1, int chrom1_begin, int chrom2, int chrom2_end)
{
	BamReader reader;
    if ( !reader.Open(filename) ) {
        cerr << "Bamtools ERROR: could not open input BAM file: " << filename << endl;
        return false;
    }
		
	//check whether need to set a region.
	if(chrom1>-1 && chrom1_begin>-1 && chrom2>-1 && chrom2_end>-1)
	{
		this->loadIndex(reader);
		BamRegion br(chrom1,chrom1_begin,chrom2,chrom2_end);
		bool is_set=reader.SetRegion(br);
		if(is_set==false)
		{
			return false;//cannot set the region.
		}
	}

	//process input data
    BamAlignment al;   
	while ( reader.GetNextAlignment(al) )
	{
		if(al.Position<0) continue;

		BamAlignmentRecord* bar=new BamAlignmentRecord();
		setAlignmentRecord(al,bar);
		bam_aln_records.push_back(bar);
	}

	reader.Close();
	return true;
}

// opens BAM files
bool BamMultiReaderPrivate::Open(const vector<string>& filenames) {

    m_errorString.clear();

    // put all current readers back at beginning (refreshes alignment cache)
    if ( !Rewind() ) {
        const string currentError = m_errorString;
        const string message = string("unable to rewind existing readers: \n\t") + currentError;
        SetErrorString("BamMultiReader::Open", message);
        return false;
    }

    // iterate over filenames
    bool errorsEncountered = false;
    vector<string>::const_iterator filenameIter = filenames.begin();
    vector<string>::const_iterator filenameEnd  = filenames.end();
    for ( ; filenameIter != filenameEnd; ++filenameIter ) {
        const string& filename = (*filenameIter);
        if ( filename.empty() ) continue;

        // attempt to open BamReader
        BamReader* reader = new BamReader;
        const bool readerOpened = reader->Open(filename);

        // if opened OK, store it
        if ( readerOpened )
            m_readers.push_back( MergeItem(reader, new BamAlignment) );

        // otherwise store error & clean up invalid reader
        else {
            m_errorString.append(1, '\t');
            m_errorString += string("unable to open file: ") + filename;
            m_errorString.append(1, '\n');
            errorsEncountered = true;

            delete reader;
            reader = 0;
        }
    }

    // check for errors while opening
    if ( errorsEncountered ) {
        const string currentError = m_errorString;
        const string message = string("unable to open all files: \t\n") + currentError;
        SetErrorString("BamMultiReader::Open", message);
        return false;
    }

    // check for BAM file consistency
    if ( !ValidateReaders() ) {
        const string currentError = m_errorString;
        const string message = string("unable to open inconsistent files: \t\n") + currentError;
        SetErrorString("BamMultiReader::Open", message);
        return false;
    }

    // update alignment cache
    return UpdateAlignmentCache();
}

/**
 * Main work method.  Reads the BAM file once and collects sorted information about
 * the 5' ends of both ends of each read (or just one end in the case of pairs).
 * Then makes a pass through those determining duplicates before re-reading the
 * input file and writing it out with duplication flags set correctly.
 */
int MarkDuplicates::runInternal() {
    
    ogeNameThread("am_MarkDuplicates");

    if(verbose)
        cerr << "Reading input file and constructing read end information." << endl;
    buildSortedReadEndLists();
    
    generateDuplicateIndexes();
    
    if(verbose)
        cerr << "Marking " << numDuplicateIndices << " records as duplicates." << endl;
    
    BamReader in;

    in.Open(getBufferFileName());

    // Now copy over the file while marking all the necessary indexes as duplicates
    long recordInFileIndex = 0;
    
    long written = 0;
    while (true) {
        BamAlignment * prec = in.GetNextAlignment();
        if(!prec)
            break;
        
        if (prec->IsPrimaryAlignment()) {
            if (duplicateIndexes.count(recordInFileIndex) == 1)
                prec->SetIsDuplicate(true);
            else
                prec->SetIsDuplicate(false);
        }
        recordInFileIndex++;
        
        if (removeDuplicates && prec->IsDuplicate()) {
            // do nothing
        }
        else {
            putOutputAlignment(prec);
            if (verbose && read_count && ++written % 100000 == 0) {
                cerr << "\rWritten " << written << " records (" << written * 100 / read_count <<"%)." << std::flush;
            }
        }
    }

    if (verbose && read_count)
        cerr << "\rWritten " << written << " records (" << written * 100 / read_count <<"%)." << endl;

    in.Close();

    remove(getBufferFileName().c_str());
    
    return 0;
}

void BedCoverage::CollectCoverageBam(string bamFile) {

    // load the "B" bed file into a map so
    // that we can easily compare "A" to it for overlaps
    _bedB->loadBedCovFileIntoMap();

    // open the BAM file
    BamReader reader;
    reader.Open(bamFile);

    // get header & reference information
    string header = reader.GetHeaderText();
    RefVector refs = reader.GetReferenceData();

    // convert each aligned BAM entry to BED
    // and compute coverage on B
    BamAlignment bam;
    while (reader.GetNextAlignment(bam)) {
        if (bam.IsMapped()) {
            // treat the BAM alignment as a single "block"
            if (_obeySplits == false) {
                // construct a new BED entry from the current BAM alignment.
                BED a;
                a.chrom  = refs.at(bam.RefID).RefName;
                a.start  = bam.Position;
                a.end    = bam.GetEndPosition(false, false);
                a.strand = "+";
                if (bam.IsReverseStrand()) a.strand = "-";

                _bedB->countHits(a, _sameStrand, _diffStrand, _countsOnly);
            }
            // split the BAM alignment into discrete blocks and
            // look for overlaps only within each block.
            else {
                // vec to store the discrete BED "blocks" from a
                bedVector bedBlocks;
                // since we are counting coverage, we do want to split blocks when a
                // deletion (D) CIGAR op is encountered (hence the true for the last parm)
                GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, false, true);
                // use countSplitHits to avoid over-counting each split chunk
                // as distinct read coverage.
                _bedB->countSplitHits(bedBlocks, _sameStrand, _diffStrand, _countsOnly);
            }
        }
    }
    // report the coverage (summary or histogram) for BED B.
    if (_countsOnly == true)
        ReportCounts();
    else 
        ReportCoverage();
    // close the BAM file
    reader.Close();
}

void parser::fetchLines(vector<BamAlignment>& result, uint32_t n,
        const std::string& file) {
    BamReader bam;
    BamAlignment read;
    Guarded<FileNotGood> g(!(bam.Open(file)), file.c_str());
    const RefVector refvec = bam.GetReferenceData();
    while (bam.GetNextAlignment(read) && n) {
        result.push_back(read);
//        cout << "read " << n << "\t" << read << "\n";
        n--;
    }
}

bool CoverageTool::CoverageToolPrivate::Run(void) {  
  
    // if output filename given
    ofstream outFile;
    if ( m_settings->HasOutputFile ) {
      
        // open output file stream
        outFile.open(m_settings->OutputFilename.c_str());
        if ( !outFile ) {
            cerr << "bamtools coverage ERROR: could not open " << m_settings->OutputFilename
                 << " for output" << endl;
            return false; 
        }
        
        // set m_out to file's streambuf
        m_out.rdbuf(outFile.rdbuf()); 
    } 
    
    //open our BAM reader
    BamReader reader;
    if ( !reader.Open(m_settings->InputBamFilename) ) {
        cerr << "bamtools coverage ERROR: could not open input BAM file: " << m_settings->InputBamFilename << endl;
        return false;
    }

    // retrieve references
    m_references = reader.GetReferenceData();
    
    // set up our output 'visitor'
    CoverageVisitor* cv = new CoverageVisitor(m_references, &m_out);
    
    // set up pileup engine with 'visitor'
    PileupEngine pileup;
    pileup.AddVisitor(cv);
    
    // process input data
    BamAlignment al;    
    while ( reader.GetNextAlignment(al) ) 
        pileup.AddAlignment(al);
    
    // clean up 
    reader.Close();
    if ( m_settings->HasOutputFile )
        outFile.close();
    delete cv;
    cv = 0;
    
    // return success
    return true;
}

int main (int argc, char *argv[]) {

     if( (argc== 1) ||
    	(argc== 2 && string(argv[1]) == "-h") ||
    	(argc== 2 && string(argv[1]) == "-help") ||
    	(argc== 2 && string(argv[1]) == "--help") ){
	 cout<<"Usage:setAsUnpaired [in bam] [outbam]"<<endl<<"this program takes flags all paired sequences as singles"<<endl;
    	return 1;
    }

     string bamfiletopen = string(argv[1]);
     string bamFileOUT   = string(argv[2]);

     BamReader reader;
     BamWriter writer;

     if ( !reader.Open(bamfiletopen) ) {
    	cerr << "Could not open input BAM files." << endl;
    	return 1;
     }
    const SamHeader header = reader.GetHeader();
    const RefVector references = reader.GetReferenceData();
    if ( !writer.Open(bamFileOUT,header,references) ) {
    	cerr << "Could not open output BAM file "<<bamFileOUT << endl;
    	return 1;
    }

    BamAlignment al;
 
    while ( reader.GetNextAlignment(al) ) {
	if(al.IsMapped()){
	    cerr << "Cannot yet handle mapped reads " << endl;
	    return 1;
	}

	
	al.SetIsPaired (false);
	
	writer.SaveAlignment(al);    

    } //while al

    reader.Close();
    writer.Close();

    return 0;
}

int main (int argc, char *argv[]) {

     if( (argc== 1) ||
    	(argc== 2 && string(argv[1]) == "-h") ||
    	(argc== 2 && string(argv[1]) == "-help") ||
    	(argc== 2 && string(argv[1]) == "--help") ){
	 cout<<"Usage:editDist [in bam]"<<endl<<"this program returns the NM field of all aligned reads"<<endl;
	 return 1;
     }

     string bamfiletopen = string(argv[1]);
     // cout<<bamfiletopen<<endl;
     BamReader reader;
     // cout<<"ok"<<endl;
     if ( !reader.Open(bamfiletopen) ) {
	 cerr << "Could not open input BAM files." << endl;
	 return 1;
     }

     BamAlignment al;
     // cout<<"ok"<<endl;
     while ( reader.GetNextAlignment(al) ) {
	 // cout<<al.Name<<endl;
	 if(!al.IsMapped())
	     continue;

	 if(al.HasTag("NM") ){
	     int editDist;
	     if(al.GetTag("NM",editDist) ){
		 cout<<editDist<<endl;
	     }else{
		 cerr<<"Cannot retrieve NM field for "<<al.Name<<endl;
		 return 1;
	     }
	 }else{
	     cerr<<"Warning: read "<<al.Name<<" is aligned but has no NM field"<<endl;
	 }

		    

     } //while al

     reader.Close();

     return 0;
}

int IndexTool::Run(int argc, char* argv[]) {
  
    // parse command line arguments
    Options::Parse(argc, argv, 1);
    
    // open our BAM reader
    BamReader reader;
    reader.Open(m_settings->InputBamFilename);
    
    // create index for BAM file
    bool useDefaultIndex = !m_settings->IsUsingBamtoolsIndex;
    reader.CreateIndex(useDefaultIndex);
    
    // clean & exit
    reader.Close();
    return 0;
}

// opens BAM files
bool BamMultiReader::Open(const vector<string>& filenames, bool openIndexes, bool coreMode, bool preferStandardIndex) {
    
    // for filename in filenames
    fileNames = filenames; // save filenames in our multireader
    for (vector<string>::const_iterator it = filenames.begin(); it != filenames.end(); ++it) {

        const string filename = *it;
        BamReader* reader = new BamReader;

        bool openedOK = true;
        openedOK = reader->Open(filename, "", openIndexes, preferStandardIndex);
        
        // if file opened ok, check that it can be read
        if ( openedOK ) {
           
            bool fileOK = true;
            BamAlignment* alignment = new BamAlignment;
            fileOK &= ( coreMode ? reader->GetNextAlignmentCore(*alignment) : reader->GetNextAlignment(*alignment) );
            
            if (fileOK) {
                readers.push_back(make_pair(reader, alignment)); // store pointers to our readers for cleanup
                alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position),
                                            make_pair(reader, alignment)));
            } else {
                cerr << "WARNING: could not read first alignment in " << filename << ", ignoring file" << endl;
                // if only file available & could not be read, return failure
                if ( filenames.size() == 1 ) return false;
            }
        } 
       
        // TODO; any further error handling when openedOK is false ??
        else 
            return false;
    }

    // files opened ok, at least one alignment could be read,
    // now need to check that all files use same reference data
    ValidateReaders();
    return true;
}