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

void getRefAltCounts(BamReader& reader, const BamRegion &r, int& nReadsRef, int& nReadsAlt, vcfRecord record)
{
	// Return the number of reads aligning reference and alternate allele
	if (!reader.IsOpen())
	{
		cout << "Can't access BAM file" << endl;
		return;
	}
	if (!reader.HasIndex())
	{
		cout << "Can't access BAM file index" << endl;
		return;
	}
	
	
	BamAlignment al;
	assert(reader.SetRegion(r));
		
	while ( reader.GetNextAlignment(al) )
	{
		// int baseAlignment = al.Position;
		
		bool found = false;
		int index = 0;
		int indexPosition = al.Position;
		for ( size_t i = 0; i < al.CigarData.size(); i++ )
		{
			switch ( al.CigarData.at(i).Type )
			{
				case 'H' : break; // Hard clip
				case 'S' : break; // Soft clip
				case 'M' : // Match or mismatch
				case 'D' :
				case 'N' :
				case 'P' : 
					if ( indexPosition + al.CigarData.at(i).Length > (record.pos-1) )
					{
						index += (record.pos-1) - indexPosition;
						found = true;
					}
					else
					{
						index += al.CigarData.at(i).Length;
						indexPosition += al.CigarData.at(i).Length;
					}
					break;
				case 'I' :
					index += al.CigarData.at(i).Length;
					break;
			}
			if ( found )
				break;
		}

		if ( !found || index >= al.AlignedBases.size())
		{
			cerr << "Problem with finding variant in read:\n";
			for ( unsigned int i = 0; i < al.AlignedBases.size(); i++ )
				cout << al.AlignedBases.at(i) << " ";
			cout << endl;
			for ( unsigned int i = 0; i < al.CigarData.size(); i++ )
				cout << al.CigarData.at(i).Type << al.CigarData.at(i).Length << " ";
			cout << endl;
			cout << "Read position: " << al.RefID+1 << " " << al.Position+1 << endl;
			cout << "Variant position: " << record.chr << " " << record.pos << endl;
			exit(1);
		}
		
		if ( al.AlignedBases.at(index) == record.ref.at(0) )
			nReadsRef++;
		else if ( al.AlignedBases.at(index) == record.alt.at(0) )
			nReadsAlt++;
	}
}

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)
//.........这里部分代码省略.........