本文整理汇总了C++中BamMultiReader类的典型用法代码示例。如果您正苦于以下问题:C++ BamMultiReader类的具体用法?C++ BamMultiReader怎么用?C++ BamMultiReader使用的例子?那么, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了BamMultiReader类的11个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: if
bool ConvertTool::ConvertToolPrivate::Run(void) {
// ------------------------------------
// initialize conversion input/output
// set to default input if none provided
if ( !m_settings->HasInput )
m_settings->InputFiles.push_back(Options::StandardIn());
// open input files
BamMultiReader reader;
if ( !reader.Open(m_settings->InputFiles) ) {
cerr << "bamtools convert ERROR: could not open input BAM file(s)... Aborting." << endl;
return false;
}
// if input is not stdin & a region is provided, look for index files
if ( m_settings->HasInput && m_settings->HasRegion ) {
if ( !reader.LocateIndexes() ) {
cerr << "bamtools convert ERROR: could not locate index file(s)... Aborting." << endl;
return false;
}
}
// retrieve reference data
m_references = reader.GetReferenceData();
// set region if specified
BamRegion region;
if ( m_settings->HasRegion ) {
if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
if ( reader.HasIndexes() ) {
if ( !reader.SetRegion(region) ) {
cerr << "bamtools convert ERROR: set region failed. Check that REGION describes a valid range" << endl;
reader.Close();
return false;
}
}
} else {
cerr << "bamtools convert ERROR: could not parse REGION: " << m_settings->Region << endl;
cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
<< endl;
reader.Close();
return false;
}
}
// if output file given
ofstream outFile;
if ( m_settings->HasOutput ) {
// open output file stream
outFile.open(m_settings->OutputFilename.c_str());
if ( !outFile ) {
cerr << "bamtools convert ERROR: could not open " << m_settings->OutputFilename
<< " for output" << endl;
return false;
}
// set m_out to file's streambuf
m_out.rdbuf(outFile.rdbuf());
}
// -------------------------------------
// do conversion based on format
bool convertedOk = true;
// pileup is special case
// conversion not done per alignment, like the other formats
if ( m_settings->Format == FORMAT_PILEUP )
convertedOk = RunPileupConversion(&reader);
// all other formats
else {
bool formatError = false;
// set function pointer to proper conversion method
void (BamTools::ConvertTool::ConvertToolPrivate::*pFunction)(const BamAlignment&) = 0;
if ( m_settings->Format == FORMAT_BED ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBed;
else if ( m_settings->Format == FORMAT_FASTA ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFasta;
else if ( m_settings->Format == FORMAT_FASTQ ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFastq;
else if ( m_settings->Format == FORMAT_JSON ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintJson;
else if ( m_settings->Format == FORMAT_SAM ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintSam;
else if ( m_settings->Format == FORMAT_YAML ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintYaml;
else {
cerr << "bamtools convert ERROR: unrecognized format: " << m_settings->Format << endl;
cerr << "Please see documentation for list of supported formats " << endl;
formatError = true;
convertedOk = false;
}
// if format selected ok
if ( !formatError ) {
// if SAM format & not omitting header, print SAM header first
if ( (m_settings->Format == FORMAT_SAM) && !m_settings->IsOmittingSamHeader )
//.........这里部分代码省略.........
示例2: Run
int CountTool::Run(int argc, char* argv[]) {
// parse command line arguments
Options::Parse(argc, argv, 1);
// if no '-in' args supplied, default to stdin
if ( !m_settings->HasInput )
m_settings->InputFiles.push_back(Options::StandardIn());
// open reader without index
BamMultiReader reader;
if (!reader.Open(m_settings->InputFiles, false, true)) {
cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
return 1;
}
// alignment counter
BamAlignment al;
int alignmentCount(0);
// if no region specified, count entire file
if ( !m_settings->HasRegion ) {
while ( reader.GetNextAlignmentCore(al) )
++alignmentCount;
}
// otherwise attempt to use region as constraint
else {
// if region string parses OK
BamRegion region;
if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
// attempt to re-open reader with index files
reader.Close();
bool openedOK = reader.Open(m_settings->InputFiles, true, true );
// if error
if ( !openedOK ) {
cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
return 1;
}
// if index data available, we can use SetRegion
if ( reader.IsIndexLoaded() ) {
// attempt to use SetRegion(), if failed report error
if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
cerr << "ERROR: Region requested, but could not set BamReader region to REGION: " << m_settings->Region << " Aborting." << endl;
reader.Close();
return 1;
}
// everything checks out, just iterate through specified region, counting alignments
while ( reader.GetNextAlignmentCore(al) )
++alignmentCount;
}
// no index data available, we have to iterate through until we
// find overlapping alignments
else {
while( reader.GetNextAlignmentCore(al) ) {
if ( (al.RefID >= region.LeftRefID) && ( (al.Position + al.Length) >= region.LeftPosition ) &&
(al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
{
++alignmentCount;
}
}
}
}
// error parsing REGION string
else {
cerr << "ERROR: Could not parse REGION - " << m_settings->Region << endl;
cerr << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
reader.Close();
return 1;
}
}
// print results
cout << alignmentCount << endl;
// clean & exit
reader.Close();
return 0;
}
示例3: if
bool ConvertTool::ConvertToolPrivate::Run(void) {
// ------------------------------------
// initialize conversion input/output
// set to default input if none provided
if ( !m_settings->HasInput )
m_settings->InputFiles.push_back(Options::StandardIn());
// open input files
BamMultiReader reader;
if ( !m_settings->HasInput ) { // don't attempt to open index for stdin
if ( !reader.Open(m_settings->InputFiles, false) ) {
cerr << "Could not open input files" << endl;
return false;
}
} else {
if ( !reader.Open(m_settings->InputFiles, true) ) {
if ( !reader.Open(m_settings->InputFiles, false) ) {
cerr << "Could not open input files" << endl;
return false;
} else {
cerr << "Opened reader without index file, jumping is disabled." << endl;
}
}
}
m_references = reader.GetReferenceData();
// set region if specified
BamRegion region;
if ( m_settings->HasRegion ) {
if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
if ( !reader.SetRegion(region) ) {
cerr << "Could not set BamReader region to REGION: " << m_settings->Region << endl;
return false;
}
} else {
cerr << "Could not parse REGION: " << m_settings->Region << endl;
return false;
}
}
// if output file given
ofstream outFile;
if ( m_settings->HasOutput ) {
// open output file stream
outFile.open(m_settings->OutputFilename.c_str());
if ( !outFile ) {
cerr << "Could not open " << m_settings->OutputFilename << " for output" << endl;
return false;
}
// set m_out to file's streambuf
m_out.rdbuf(outFile.rdbuf());
}
// -------------------------------------
// do conversion based on format
bool convertedOk = true;
// pileup is special case
// conversion not done per alignment, like the other formats
if ( m_settings->Format == FORMAT_PILEUP )
convertedOk = RunPileupConversion(&reader);
// all other formats
else {
bool formatError = false;
// set function pointer to proper conversion method
void (BamTools::ConvertTool::ConvertToolPrivate::*pFunction)(const BamAlignment&) = 0;
if ( m_settings->Format == FORMAT_BED ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBed;
else if ( m_settings->Format == FORMAT_BEDGRAPH ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintBedGraph;
else if ( m_settings->Format == FORMAT_FASTA ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFasta;
else if ( m_settings->Format == FORMAT_FASTQ ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintFastq;
else if ( m_settings->Format == FORMAT_JSON ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintJson;
else if ( m_settings->Format == FORMAT_SAM ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintSam;
else if ( m_settings->Format == FORMAT_WIGGLE ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintWiggle;
else if ( m_settings->Format == FORMAT_YAML ) pFunction = &BamTools::ConvertTool::ConvertToolPrivate::PrintYaml;
else {
cerr << "Unrecognized format: " << m_settings->Format << endl;
cerr << "Please see help|README (?) for details on supported formats " << endl;
formatError = true;
convertedOk = false;
}
// if format selected ok
if ( !formatError ) {
// if SAM format & not omitting header, print SAM header first
if ( (m_settings->Format == FORMAT_SAM) && !m_settings->IsOmittingSamHeader )
m_out << reader.GetHeaderText();
// iterate through file, doing conversion
BamAlignment a;
while ( reader.GetNextAlignment(a) )
(this->*pFunction)(a);
//.........这里部分代码省略.........
示例4: main
int main ( int argc, char *argv[] ) {
struct parameters *param = 0;
param = interface(param, argc, argv);
//bam input and generate index if not yet
//-------------------------------------------------------------------------------------------------------+
// BAM input (file or filenames?) |
//-------------------------------------------------------------------------------------------------------+
char *fof = param->mapping_f;
FILE *IN=NULL;
char linefof[5000];
int filecount=0;
vector <string> fnames;
if (strchr(fof,' ')!=NULL) {
char *ptr;
ptr=strtok(fof," ");
while (ptr!=NULL) {
fnames.push_back(ptr);
filecount++;
ptr=strtok(NULL," ");
}
} else {
IN=fopen(fof,"rt");
if (IN!=NULL) {
long linecount=0;
while (fgets(linefof,5000-1,IN)!=NULL) {
linecount++;
if (linefof[0]!='#' && linefof[0]!='\n') {
char *ptr=strchr(linefof,'\n');
if (ptr!=NULL && ptr[0]=='\n') {
ptr[0]='\0';
}
FILE *dummy=NULL;
dummy=fopen(linefof,"rt");
if (dummy!=NULL) { // seems to be a file of filenames...
fclose(dummy);
fnames.push_back(linefof);
filecount++;
} else if (filecount==0 || linecount>=1000-1) { // seems to be a single file
fnames.push_back(fof);
filecount++;
break;
}
}
}
fclose(IN);
}
} //file or file name decided and stored in vector "fnames"
cerr << "the input mapping files are:" << endl;
vector <string>::iterator fit = fnames.begin();
for(; fit != fnames.end(); fit++) {
cerr << *fit << endl;
}
//-------------------------------------------------------------------------------------------------------+
// end of file or filenames |
//-------------------------------------------------------------------------------------------------------+
// open the BAM file(s)
BamMultiReader reader;
reader.Open(fnames);
// get header & reference information
string header = reader.GetHeaderText();
RefVector refs = reader.GetReferenceData();
// attempt to open BamWriter
BamWriter writer;
string outputBam = param->writer;
if ( outputBam != "" ) {
if ( !writer.Open(param->writer, header, refs) ) {
cerr << "Could not open output BAM file" << endl;
exit(0);
}
}
BamAlignment bam;
while (reader.GetNextAlignment(bam)) { //change RG
string rg = "RG";
string rgType = "Z";
string rgValue = "1";
bam.EditTag(rg,rgType,rgValue);
writer.SaveAlignment(bam);
} // read a bam
return 0;
} //main
示例5: filelist
bool FilterTool::FilterToolPrivate::Run(void) {
// set to default input if none provided
if ( !m_settings->HasInput && !m_settings->HasInputFilelist )
m_settings->InputFiles.push_back(Options::StandardIn());
// add files in the filelist to the input file list
if ( m_settings->HasInputFilelist ) {
ifstream filelist(m_settings->InputFilelist.c_str(), ios::in);
if ( !filelist.is_open() ) {
cerr << "bamtools filter ERROR: could not open input BAM file list... Aborting." << endl;
return false;
}
string line;
while ( getline(filelist, line) )
m_settings->InputFiles.push_back(line);
}
// initialize defined properties & user-specified filters
// quit if failed
if ( !SetupFilters() )
return false;
// open reader without index
BamMultiReader reader;
if ( !reader.Open(m_settings->InputFiles) ) {
cerr << "bamtools filter ERROR: could not open input files for reading." << endl;
return false;
}
// retrieve reader header & reference data
const string headerText = reader.GetHeaderText();
filterToolReferences = reader.GetReferenceData();
// determine compression mode for BamWriter
bool writeUncompressed = ( m_settings->OutputFilename == Options::StandardOut() &&
!m_settings->IsForceCompression );
BamWriter::CompressionMode compressionMode = BamWriter::Compressed;
if ( writeUncompressed ) compressionMode = BamWriter::Uncompressed;
// open BamWriter
BamWriter writer;
writer.SetCompressionMode(compressionMode);
if ( !writer.Open(m_settings->OutputFilename, headerText, filterToolReferences) ) {
cerr << "bamtools filter ERROR: could not open " << m_settings->OutputFilename << " for writing." << endl;
reader.Close();
return false;
}
// if no region specified, filter entire file
BamAlignment al;
if ( !m_settings->HasRegion ) {
while ( reader.GetNextAlignment(al) ) {
if ( CheckAlignment(al) )
writer.SaveAlignment(al);
}
}
// otherwise attempt to use region as constraint
else {
// if region string parses OK
BamRegion region;
if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
// attempt to find index files
reader.LocateIndexes();
// if index data available for all BAM files, we can use SetRegion
if ( reader.HasIndexes() ) {
// attempt to use SetRegion(), if failed report error
if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
cerr << "bamtools filter ERROR: set region failed. Check that REGION describes a valid range" << endl;
reader.Close();
return false;
}
// everything checks out, just iterate through specified region, filtering alignments
while ( reader.GetNextAlignment(al) )
if ( CheckAlignment(al) )
writer.SaveAlignment(al);
}
// no index data available, we have to iterate through until we
// find overlapping alignments
else {
while ( reader.GetNextAlignment(al) ) {
if ( (al.RefID >= region.LeftRefID) && ((al.Position + al.Length) >= region.LeftPosition) &&
(al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
{
if ( CheckAlignment(al) )
writer.SaveAlignment(al);
}
}
}
}
//.........这里部分代码省略.........
示例6: CropBam
int CropBamTool::CropBam()
{
// open bam files
BamMultiReader bamReader;
bamReader.Open(bamFiles);
// the dictionary of chromosomes
RefVector genome = bamReader.GetReferenceData();
// get the scanning window
vector<tuple<int,int,int>> windows;
int numWindows = GenericRegionTools::toScanWindow(genome, regionStrings, windows);
unordered_set<string> readpool;
// temporary struct for sequence object
typedef struct {
string name;
int head_soft_clip;
int tail_soft_clip;
string seq;
string qual;
}cropbam_seq_t;
// temporary struct for unique seqs
map<string,list<cropbam_seq_t>> uniqueSeqPool;
// lambda expression for output
auto Output = [this](cropbam_seq_t &a){
if (this->outFormat=="fasta"){
cout << ">" << a.name << "\t"
<< "head_soft_clip=" << a.head_soft_clip << "\t"
<< "tail_soft_clip=" << a.tail_soft_clip << "\t"
<< endl
<< a.seq << endl;
}
if (this->outFormat=="fastq"){
cout << "@" << a.name << "\t"
<< "head_soft_clip=" << a.head_soft_clip << "\t"
<< "tail_soft_clip=" << a.tail_soft_clip << "\t"
<< endl
<< a.seq << endl;
cout << "+" << endl
<< a.qual << endl;
}
};
// loop over windows
omp_set_dynamic(0);
omp_set_num_threads(numThreads);
#pragma omp parallel for shared(genome)
for (int i=0; i<numWindows; i++)
{
clock_t tStart = clock();
bamReader.Open(bamFiles);
int wId = get<0>(windows[i]);
int wLp = get<1>(windows[i]);
int wRp = get<2>(windows[i]);
if (verbose>=1) Verbose("process the window " + genome[wId].RefName + ":" + to_string(wLp+1) + "-" + to_string(wRp));
// rewind the bam reader
bamReader.Rewind();
// set the region
bamReader.SetRegion(wId, wLp, wId, wRp);
int numReads = 0;
// retrieve the alignment
BamAlignment aln;
while (bamReader.GetNextAlignment(aln))
{
// skip the alignment if it doesn't overlap the window
if (aln.Position>=wRp || aln.GetEndPosition()<=wLp)
continue;
// skip the invalid alignment
if (!isValidAlignment(aln, readLenThres, mapQualThres, alnFlagMarker))
continue;
// skip the alignment harboring too many mismatches
if (!GenericBamAlignmentTools::validReadIdentity(aln, 1-alnIdenThres))
continue;
stringstream keyss;
keyss << GenericBamAlignmentTools::getBamAlignmentName(aln) << "-"
<< wId << "-" << wLp << "-" << wRp;
string key = keyss.str();
auto ptr = readpool.find(key);
if (ptr!=readpool.end())
continue;
readpool.emplace(key);
// get the partial read
string readSegment, readQualSegment, genomeSegment;
GenericBamAlignmentTools::getLocalAlignment(aln, wLp, wRp-wLp, readSegment, readQualSegment, genomeSegment);
// add soft clip
int hsc=0;
//.........这里部分代码省略.........
示例7: bbctools_create
int bbctools_create( BbcUtils::OptParser &optParser ) {
const vector<string> cmdArgs = optParser.getArgs();
// remove .bbc extension from bbc file root, if present
string bbcfileRoot = optParser.getOptValue( "bbc" );
int i = bbcfileRoot.size() - 4;
if( i > 0 && bbcfileRoot.substr(i,4) == ".bbc" ) {
bbcfileRoot = bbcfileRoot.substr(0,i);
}
bool f_bci = optParser.getOptBoolean("index");
bool f_cbc = optParser.getOptBoolean("coarse");
string targetRegions = optParser.getOptValue("regions");
string annotationFields = optParser.getOptValue( "annotationFields");
vector<string> auxRegionSplit = BbcUtils::mapToPairList(annotationFields);
string sumstatsFile = optParser.getOptValue("sumStats");
string covstatsFile = optParser.getOptValue("covStats");
string readOrigFile = optParser.getOptValue("readOrigin");
string readType = optParser.getOptValue("readType");
string covDepths = optParser.getOptValue("covDepths","-");
double minPcCov = optParser.getOptNumber("minPcCov");
int32_t primerLength = optParser.getOptInteger( "primerLength", (readType == "AmpliSeq" ? 30 : 0) );
int32_t maxE2eEndGap = optParser.getOptInteger( "e2eGap", (readType == "AmpliSeq" ? 2 : 0) );
bool autoCreateBamIndex = optParser.getOptBoolean("autoCreateBamIndex");
bool samdepth = optParser.getOptBoolean("samdepth");
int32_t filterQuality = optParser.getOptInteger("minMAPQ");
int32_t minAlignLength = optParser.getOptInteger("minAlignLength");
bool filterDuplicates = optParser.getOptBoolean("noDups");
bool filterUnique = optParser.getOptBoolean("unique");
uint32_t skipFlag = filterDuplicates ? 0x704 : 0x304;
uint16_t minMapQuality = filterUnique ? 1 : filterQuality;
bool onlyOnTargetReads = optParser.getOptBoolean("onTargetReads");
bool onlyOnTargetBases = optParser.getOptBoolean("onTargetBases");
// possible future options
bool invertOnTarget = false;
// check basic valid argument values and combinations
int numOuts = !bbcfileRoot.empty() + !covstatsFile.empty() + !sumstatsFile.empty() + !readOrigFile.empty();
int numPipes = (bbcfileRoot == "-") + (covstatsFile == "-") + (sumstatsFile == "-") + (readOrigFile == "-");
if( numOuts == 0 && !f_bci && !f_cbc ) {
bbcfileRoot = "-"; // default if no other output specified
} else if( numPipes > 1 ) {
cerr << "Error: bbctools create: Only one file output (--covStats, --sumStats, --readOrigin or --bbc) may be piped to STDOUT." << endl;
return -1;
} else if( samdepth && numOuts ) {
cerr << "Error: bbctools create: --samdepth (-s) option may only be used without other output options." << endl;
return -1;
}
// check if single argument is a BBC file and leave open for reading if so
BbcView bbcView;
bool haveBbcFile = cmdArgs.size() == 1 && bbcView.Open( cmdArgs[0], true );
bbcView.SelectPrintStream( samdepth ? "SAMDEPTH" : "BBCVIEW" );
// check distinction between default and explicit no target regions - only for BBC input
bool explicitNoTargetRegions = false;
if( targetRegions == "-" ) {
explicitNoTargetRegions = haveBbcFile;
targetRegions = "";
}
if( targetRegions.empty() ) {
if( onlyOnTargetBases && explicitNoTargetRegions && !invertOnTarget ) {
cerr << "Warning: bbctools create --onTargetBases (-b) option with --regions '-' produces no coverage." << endl;
} else if( onlyOnTargetReads ) {
cerr << "Error: bbctools create --onTargetReads (-r) option requires a --regions file." << endl;
return -1;
}
}
// check for legal BBC create options
if( f_bci || f_cbc ) {
if( (bbcfileRoot.empty() || bbcfileRoot == "-") && !haveBbcFile ) {
string opt = f_bci ? "--index (-i)" : "--coarse (-c)";
cerr << "Error: bbctools create "+opt+" option requires the --bbc (-B) option or a BBC source file." << endl;
return -1;
}
}
BamMultiReader bamReader;
if( haveBbcFile ) {
// warn for options that do not work with BBC input
if( filterQuality > 0 || filterDuplicates || filterUnique || minAlignLength ) {
cerr << "Warning: SAM flag, alignment length and MAPQ filters ignored for BBC source file." << endl;
}
if( samdepth ) {
cerr << "Error: --samdepth option is not supported for BBC source files." << endl;
return -1;
}
if( !readOrigFile.empty() ) {
cerr << "Error: --readOrigin option is not supported for BBC source files." << endl;
return -1;
}
} else {
// check / open for multiple BAM file inputs
if ( !bamReader.Open(cmdArgs) ) {
if( cmdArgs.size() == 1 ) cerr << "ERROR: Could not read input BAM file:";
else cerr << "ERROR: Could not read all input BAM files:";
// get and clean up bamtools error msg
string errMsg = bamReader.GetErrorString();
//.........这里部分代码省略.........
示例8: file
bool MergeTool::MergeToolPrivate::Run(void) {
// set to default input if none provided
if ( !m_settings->HasInputBamFilename )
m_settings->InputFiles.push_back(Options::StandardIn());
// opens the BAM files (by default without checking for indexes)
BamMultiReader reader;
if ( !reader.Open(m_settings->InputFiles) ) {
cerr << "bamtools merge ERROR: could not open input BAM file(s)... Aborting." << endl;
return false;
}
// retrieve header & reference dictionary info
std::string mergedHeader = reader.GetHeaderText();
RefVector references = reader.GetReferenceData();
// determine compression mode for BamWriter
bool writeUncompressed = ( m_settings->OutputFilename == Options::StandardOut() &&
!m_settings->IsForceCompression );
BamWriter::CompressionMode compressionMode = BamWriter::Compressed;
if ( writeUncompressed ) compressionMode = BamWriter::Uncompressed;
// open BamWriter
BamWriter writer;
writer.SetCompressionMode(compressionMode);
if ( !writer.Open(m_settings->OutputFilename, mergedHeader, references) ) {
cerr << "bamtools merge ERROR: could not open "
<< m_settings->OutputFilename << " for writing." << endl;
reader.Close();
return false;
}
// if no region specified, store entire contents of file(s)
if ( !m_settings->HasRegion ) {
BamAlignment al;
while ( reader.GetNextAlignmentCore(al) )
writer.SaveAlignment(al);
}
// otherwise attempt to use region as constraint
else {
// if region string parses OK
BamRegion region;
if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
// attempt to find index files
reader.LocateIndexes();
// if index data available for all BAM files, we can use SetRegion
if ( reader.HasIndexes() ) {
// attempt to use SetRegion(), if failed report error
if ( !reader.SetRegion(region.LeftRefID,
region.LeftPosition,
region.RightRefID,
region.RightPosition) )
{
cerr << "bamtools merge ERROR: set region failed. Check that REGION describes a valid range"
<< endl;
reader.Close();
return false;
}
// everything checks out, just iterate through specified region, storing alignments
BamAlignment al;
while ( reader.GetNextAlignmentCore(al) )
writer.SaveAlignment(al);
}
// no index data available, we have to iterate through until we
// find overlapping alignments
else {
BamAlignment al;
while ( reader.GetNextAlignmentCore(al) ) {
if ( (al.RefID >= region.LeftRefID) && ( (al.Position + al.Length) >= region.LeftPosition ) &&
(al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
{
writer.SaveAlignment(al);
}
}
}
}
// error parsing REGION string
else {
cerr << "bamtools merge ERROR: could not parse REGION - " << m_settings->Region << endl;
cerr << "Check that REGION is in valid format (see documentation) and that the coordinates are valid"
<< endl;
reader.Close();
writer.Close();
return false;
}
}
// clean & exit
reader.Close();
writer.Close();
return true;
//.........这里部分代码省略.........
示例9: main
int main ( int argc, char *argv[] ) {
struct parameters *param = 0;
param = interface(param, argc, argv);
//region file input (the region file should be sorted as the same way as the bam file)
ifstream region_f;
region_f.open(param->region_f, ios_base::in); // the region file is opened
//bam input and generate index if not yet
//-------------------------------------------------------------------------------------------------------+
// BAM input (file or filenames?) |
//-------------------------------------------------------------------------------------------------------+
char *fof = param->mapping_f;
FILE *IN=NULL;
char linefof[5000];
int filecount=0;
vector <string> fnames;
if (strchr(fof,' ')!=NULL) {
char *ptr;
ptr=strtok(fof," ");
while (ptr!=NULL) {
fnames.push_back(ptr);
filecount++;
ptr=strtok(NULL," ");
}
} else {
IN=fopen(fof,"rt");
if (IN!=NULL) {
long linecount=0;
while (fgets(linefof,5000-1,IN)!=NULL) {
linecount++;
if (linefof[0]!='#' && linefof[0]!='\n') {
char *ptr=strchr(linefof,'\n');
if (ptr!=NULL && ptr[0]=='\n') {
ptr[0]='\0';
}
FILE *dummy=NULL;
dummy=fopen(linefof,"rt");
if (dummy!=NULL) { // seems to be a file of filenames...
fclose(dummy);
fnames.push_back(linefof);
filecount++;
} else if (filecount==0 || linecount>=1000-1) { // seems to be a single file
fnames.push_back(fof);
filecount++;
break;
}
}
}
fclose(IN);
}
} //file or file name decided and stored in vector "fnames"
cerr << "the input mapping files are:" << endl;
vector <string>::iterator fit = fnames.begin();
for(; fit != fnames.end(); fit++) {
cerr << *fit << endl;
}
//-------------------------------------------------------------------------------------------------------+
// end of file or filenames |
//-------------------------------------------------------------------------------------------------------+
// open the BAM file(s)
BamMultiReader reader;
reader.Open(fnames);
// get header & reference information
string header = reader.GetHeaderText();
RefVector refs = reader.GetReferenceData();
if ( ! reader.LocateIndexes() ) // opens any existing index files that match our BAM files
reader.CreateIndexes(); // creates index files for BAM files that still lack one
// locus bias
struct lb empty_profile = {0,0,0,0};
vector <struct lb> locus_b(1000, empty_profile);
// output locus bias file
string locus_bias_set = param->lbias;
ofstream locus_bias;
if ( locus_bias_set != "" ) {
locus_bias.open(param->lbias);
if ( !locus_bias ) {
cerr << "can not open locus_bias file.\n";
exit(0);
}
}
//should decide which chromosome
string line;
string old_chr = "SRP";
string type = param->type;
//whether do some position-level pile-up stuff
bool posc = false;
ofstream posc_f;
ofstream chrmap_f;
//.........这里部分代码省略.........
示例10: ogeNameThread
int FileReader::runInternal()
{
ogeNameThread("am_FileReader");
if(!format_specified)
format = deduceFileFormat();
if(format == FORMAT_BAM)
{
BamMultiReader reader;
if(!reader.Open(filenames)) {
cerr << "Error opening BAM files." << endl;
reader.Close();
return -1;
}
header = reader.GetHeader();
references = reader.GetReferenceData();
open = true;
BamAlignment * al;
while(true)
{
if(load_string_data)
al = reader.GetNextAlignment();
else
al = reader.GetNextAlignmentCore();
if(!al)
break;
putOutputAlignment(al);
}
reader.Close();
} else if(format == FORMAT_SAM) {
vector<SamReader> readers;
SamHeader first_header;
// before doing any reading, open the files to
// verify they are the right format, etc.
for(int i = 0; i < filenames.size(); i++) {
SamReader reader;
if(!reader.Open(filenames[i])) {
cerr << "Error opening SAM file: " << filenames[i] << endl;
return -1;
}
if(filenames.size() > 1 && i == 0)
first_header = header;
// TODO: We can probably find a better way to deal with multiple SAM file headers,
// but for now we should disallow different headers to avoid issues.
if(i > 0 && header.ToString() != first_header.ToString())
cerr << "Warning! SAM input files have different headers." << endl;
reader.Close();
}
for(int i = 0; i < filenames.size(); i++) {
SamReader reader;
if(!reader.Open(filenames[i])) {
cerr << "Error opening SAM file: " << filenames[i] << endl;
return -1;
}
header = reader.GetHeader();
references = reader.GetReferenceData();
open = true;
if(filenames.size() > 1 && i == 0)
first_header = header;
BamAlignment * al = NULL;
while(true)
{
al = reader.GetNextAlignment();
if(NULL == al)
break;
putOutputAlignment(al);
}
reader.Close();
}
} else {
cerr << "FileReader couldn't detect file format. Aborting." << endl;
exit(-1);
return -1;
}
return 0;
}
示例11: Run
int MergeTool::Run(int argc, char* argv[]) {
// parse command line arguments
Options::Parse(argc, argv, 1);
// set to default input if none provided
if ( !m_settings->HasInputBamFilename )
m_settings->InputFiles.push_back(Options::StandardIn());
// opens the BAM files (by default without checking for indexes)
BamMultiReader reader;
if ( !reader.Open(m_settings->InputFiles, false, true) ) {
cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
return 1;
}
// retrieve header & reference dictionary info
std::string mergedHeader = reader.GetHeaderText();
RefVector references = reader.GetReferenceData();
// open writer
BamWriter writer;
bool writeUncompressed = ( m_settings->OutputFilename == Options::StandardOut() && !m_settings->IsForceCompression );
if ( !writer.Open(m_settings->OutputFilename, mergedHeader, references, writeUncompressed) ) {
cerr << "ERROR: Could not open BAM file " << m_settings->OutputFilename << " for writing... Aborting." << endl;
reader.Close();
return 1;
}
// if no region specified, store entire contents of file(s)
if ( !m_settings->HasRegion ) {
BamAlignment al;
while ( reader.GetNextAlignmentCore(al) )
writer.SaveAlignment(al);
}
// otherwise attempt to use region as constraint
else {
// if region string parses OK
BamRegion region;
if ( Utilities::ParseRegionString(m_settings->Region, reader, region) ) {
// attempt to re-open reader with index files
reader.Close();
bool openedOK = reader.Open(m_settings->InputFiles, true, true );
// if error
if ( !openedOK ) {
cerr << "ERROR: Could not open input BAM file(s)... Aborting." << endl;
return 1;
}
// if index data available, we can use SetRegion
if ( reader.IsIndexLoaded() ) {
// attempt to use SetRegion(), if failed report error
if ( !reader.SetRegion(region.LeftRefID, region.LeftPosition, region.RightRefID, region.RightPosition) ) {
cerr << "ERROR: Region requested, but could not set BamReader region to REGION: " << m_settings->Region << " Aborting." << endl;
reader.Close();
return 1;
}
// everything checks out, just iterate through specified region, storing alignments
BamAlignment al;
while ( reader.GetNextAlignmentCore(al) )
writer.SaveAlignment(al);
}
// no index data available, we have to iterate through until we
// find overlapping alignments
else {
BamAlignment al;
while ( reader.GetNextAlignmentCore(al) ) {
if ( (al.RefID >= region.LeftRefID) && ( (al.Position + al.Length) >= region.LeftPosition ) &&
(al.RefID <= region.RightRefID) && ( al.Position <= region.RightPosition) )
{
writer.SaveAlignment(al);
}
}
}
}
// error parsing REGION string
else {
cerr << "ERROR: Could not parse REGION - " << m_settings->Region << endl;
cerr << "Be sure REGION is in valid format (see README) and that coordinates are valid for selected references" << endl;
reader.Close();
writer.Close();
return 1;
}
}
// clean & exit
reader.Close();
writer.Close();
return 0;
}