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

// get next alignment among all files without parsing character data from alignments
bool BamMultiReader::GetNextAlignmentCore(BamAlignment& nextAlignment) {

    // bail out if we are at EOF in all files, means no more alignments to process
    if (!HasOpenReaders())
        return false;

    // when all alignments have stepped into a new target sequence, update our
    // current reference sequence id
    UpdateReferenceID();

    // our lowest alignment and reader will be at the front of our alignment index
    BamAlignment* alignment = alignments.begin()->second.second;
    BamReader* reader = alignments.begin()->second.first;

    // now that we have the lowest alignment in the set, save it by copy to our argument
    nextAlignment = BamAlignment(*alignment);
    //memcpy(&nextAlignment, alignment, sizeof(BamAlignment));

    // remove this alignment index entry from our alignment index
    alignments.erase(alignments.begin());

    // and add another entry if we can get another alignment from the reader
    if (reader->GetNextAlignmentCore(*alignment)) {
        alignments.insert(make_pair(make_pair(alignment->RefID, alignment->Position), 
                                    make_pair(reader, alignment)));
    } else { // do nothing
        //cerr << "reached end of file " << lowestReader->GetFilename() << endl;
    }

    return true;

}

// 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;
}

// 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;
}

void StructuralVariations::findTranslocationsOnTheFly(string bamFileName, bool outtie, float meanCoverage, string outputFileHeader, map<string,int> SV_options) {
	size_t start = time(NULL);
	//open the bam file
	BamReader bamFile;
	bamFile.Open(bamFileName);
	//Information from the header is needed to initialize the data structure
	SamHeader head = bamFile.GetHeader();
	// now create Translocation on the fly
	Window *window;

	window = new Window(bamFileName,outtie,meanCoverage,outputFileHeader,SV_options);
	window->initTrans(head);
	//expands a vector so that it is large enough to hold reads from each contig in separate elements
	window->eventReads.resize(SV_options["contigsNumber"]);
	window->eventSplitReads.resize(SV_options["contigsNumber"]);

	window-> binnedCoverage.resize(SV_options["contigsNumber"]);
	window-> linksFromWin.resize(SV_options["contigsNumber"]);
	
	window -> numberOfEvents = 0;

	string line;
	string coverageFile=outputFileHeader+".tab";
	ifstream inputFile( coverageFile.c_str() );
	int line_number=0;
	while (std::getline( inputFile, line )){
		if(line_number > 0){
			vector<string> splitline;
    		std::stringstream ss(line);
    		std::string item;
    		while (std::getline(ss, item, '\t')) {
        		splitline.push_back(item);
    		}
			window -> binnedCoverage[window -> contig2position[splitline[0]]].push_back(atof(splitline[3].c_str()));
		}
		line_number += 1;
	}
	inputFile.close();


	//Initialize bam entity
	BamAlignment currentRead;
	//now start to iterate over the bam file
	int counter = 0;
	while ( bamFile.GetNextAlignmentCore(currentRead) ) {
	  if(currentRead.IsMapped()) {
	    window->insertRead(currentRead);
	  }
	}
	for(int i=0;i< window-> eventReads.size();i++){
	  if(window -> eventReads[i].size() >= window -> minimumPairs){
	    window->computeVariations(i);
	  }
	  window->eventReads[i]=queue<BamAlignment>();
	  window->eventSplitReads[i] = vector<BamAlignment>();
	}
	  
	window->interChrVariationsVCF.close();
	window->intraChrVariationsVCF.close();
	printf ("variant calling time consumption= %lds\n", time(NULL) - start);
}

// generates mutiple sorted temp BAM files from single unsorted BAM file
bool SortTool::SortToolPrivate::GenerateSortedRuns(void) {
    
    // open input BAM file
    BamReader reader;
    if ( !reader.Open(m_settings->InputBamFilename) ) {
        cerr << "bamtools sort ERROR: could not open " << m_settings->InputBamFilename
             << " for reading... Aborting." << endl;
        return false;
    }
    
    // get basic data that will be shared by all temp/output files 
    SamHeader header = reader.GetHeader();
    header.SortOrder = ( m_settings->IsSortingByName
                       ? Constants::SAM_HD_SORTORDER_QUERYNAME
                       : Constants::SAM_HD_SORTORDER_COORDINATE );
    m_headerText = header.ToString();
    m_references = reader.GetReferenceData();
    
    // set up alignments buffer
    BamAlignment al;
    vector<BamAlignment> buffer;
    buffer.reserve( (size_t)(m_settings->MaxBufferCount*1.1) );
    bool bufferFull = false;
    
    // if sorting by name, we need to generate full char data
    // so can't use GetNextAlignmentCore()
    if ( m_settings->IsSortingByName ) {

        // iterate through file
        while ( reader.GetNextAlignment(al)) {

            // check buffer's usage
            bufferFull = ( buffer.size() >= m_settings->MaxBufferCount );

            // store alignments until buffer is "full"
            if ( !bufferFull )
                buffer.push_back(al);

            // if buffer is "full"
            else {

                // push any unmapped reads into buffer,
                // don't want to split these into a separate temp file
                if ( !al.IsMapped() )
                    buffer.push_back(al);

                // "al" is mapped, so create a sorted temp file with current buffer contents
                // then push "al" into fresh buffer
                else {
                    CreateSortedTempFile(buffer);
                    buffer.push_back(al);
                }
            }
        }
    }

    // sorting by position, can take advantage of GNACore() speedup
    else {

        // iterate through file
        while ( reader.GetNextAlignmentCore(al) ) {

            // check buffer's usage
            bufferFull = ( buffer.size() >= m_settings->MaxBufferCount );

            // store alignments until buffer is "full"
            if ( !bufferFull )
                buffer.push_back(al);

            // if buffer is "full"
            else {

                // push any unmapped reads into buffer,
                // don't want to split these into a separate temp file
                if ( !al.IsMapped() )
                    buffer.push_back(al);

                // "al" is mapped, so create a sorted temp file with current buffer contents
                // then push "al" into fresh buffer
                else {
                    CreateSortedTempFile(buffer);
                    buffer.push_back(al);
                }
            }
        }
    }

    // handle any leftover buffer contents
    if ( !buffer.empty() )
        CreateSortedTempFile(buffer);
    
    // close reader & return success
    reader.Close();
    return true;
}

int main(int argc, char** argv) {
    FILE *fsp = 0, *fep = 0;

    int flag;       // 0
    char chrom[255];    // 1
    int pos;        // 2
    int quality;    // 3
    char *cigar;    // 4

    size_t len = 0;
    
    char chromcopy[255] = { 0 };

    char fn[255];

    int i;
    int forward, reverse;

    bool direction; // 0
    bool leftflag;

    const char *prefix = "";

    int cnt = 0;
    BamReader reader;

    //if (argc < 2 || argc > 4) usage();

    for (i=1; i<argc; i++) {
        if ((strcmp(argv[i], "-p") == 0) && (i != argc-1)) {
            prefix = argv[++i];
        } else {
            if (!reader.IsOpen()) {
                try {
                    reader.Open(argv[i]);
                } catch (exception& e) {
                    cout << e.what();
                    throw;
                }
            } else {
                usage();
            }
        }
    }
    if (!reader.IsOpen()) usage();

    const SamHeader header = reader.GetHeader();
    const RefVector references = reader.GetReferenceData();
    
    BamAlignment al;
    while(reader.GetNextAlignmentCore(al)) {
        flag = al.AlignmentFlag;
        if ((flag & 256) || (flag & 4) || (flag & 512) || (flag & 1024)) continue;

        direction = ( (flag & 16) == 16);
        strcpy(chrom, references[al.RefID].RefName.c_str());
        pos = al.Position+1;
        quality = al.MapQuality;
        if (quality == 0) continue;

        if (strcmp(chromcopy, chrom) != 0) {
            if (fsp) fclose(fsp);
            if (fep) fclose(fep);
            fn[0] = 0;
            if (prefix != "") {
                strcat(fn, prefix);
            }
            strcat(fn, chrom);
            strcat(fn, "_forward.txt");
            fsp = fopen(fn, "w");
            fn[0] = 0;
            if (prefix != "") {
                strcat(fn, prefix);
            }
            strcat(fn, chrom);
            strcat(fn, "_reverse.txt");
            fep = fopen(fn, "w");
            strcpy(chromcopy, chrom);
        }

        vector<CigarOp> cigars = al.CigarData;
        forward = pos;
        reverse = pos;
        leftflag = true;
        for (i = 0; i < cigars.size(); i++) {
            if (cigars[i].Type == 'S' || cigars[i].Type == 'H') {
                if (leftflag == false) break;
                continue;
            } else {
                if (cigars[i].Type != 'I') reverse += cigars[i].Length;
                leftflag = false;
            }
        }
        if (direction == 0)
            fprintf(fsp, "%d\n", forward);
        else
            fprintf(fep, "%d\n", reverse-1);

        cnt++;
        if (cnt % 500000 == 0)
//.........这里部分代码省略.........