diff --git a/docs/content/tools/genomecov.rst b/docs/content/tools/genomecov.rst index e0c73a6e..4e4e8d7d 100755 --- a/docs/content/tools/genomecov.rst +++ b/docs/content/tools/genomecov.rst @@ -45,6 +45,7 @@ Usage and option summary **-ibam** | BAM file as input for coverage. Each BAM alignment in A added to the total coverage for the genome. | Use "stdin" or simply "-" if passing it with a UNIX pipe: For example: | ``samtools view -b | genomeCoverageBed -ibam stdin -g hg18.genome`` +**-ibams** | The input list file each line is filename in BAM format, Report each file chrom by chrom and final all together. **-d** Report the depth at each genome position with 1-based coordinates. **-dz** Report the depth at each genome position with 0-based coordinates. **-bg** Report depth in BedGraph format. For details, see: http://genome.ucsc.edu/goldenPath/help/bedgraph.html diff --git a/src/genomeCoverageBed/genomeCoverageBed.cpp b/src/genomeCoverageBed/genomeCoverageBed.cpp index 621c461d..be1f6a03 100644 --- a/src/genomeCoverageBed/genomeCoverageBed.cpp +++ b/src/genomeCoverageBed/genomeCoverageBed.cpp @@ -13,7 +13,7 @@ Licenced under the GNU General Public License 2.0 license. #include "genomeCoverageBed.h" -BedGenomeCoverage::BedGenomeCoverage(string bedFile, string genomeFile, +BedGenomeCoverage::BedGenomeCoverage(vector bedFiles, string genomeFile, bool eachBase, bool startSites, bool bedGraph, bool bedGraphAll, int max, float scale, @@ -23,7 +23,7 @@ BedGenomeCoverage::BedGenomeCoverage(string bedFile, string genomeFile, bool eachBaseZeroBased, bool add_gb_track_line, string gb_track_line_opts) { - _bedFile = bedFile; + _bedFiles = bedFiles; _genomeFile = genomeFile; _eachBase = eachBase; _eachBaseZeroBased = eachBaseZeroBased; @@ -42,20 +42,23 @@ BedGenomeCoverage::BedGenomeCoverage(string bedFile, string genomeFile, _gb_track_line_opts = gb_track_line_opts; _currChromName = ""; _currChromSize = 0 ; + _currBedFile = 0 ; + _lastBedFile = -1 ; + - if (_bamInput == false) { _genome = new GenomeFile(genomeFile); } - + PrintTrackDefinitionLine(); if (_bamInput == false) { + string bedFile = _bedFiles.front(); _bed = new BedFile(bedFile); CoverageBed(); } else { - CoverageBam(_bedFile); + CoverageBam(_bedFiles); } } @@ -98,7 +101,7 @@ void BedGenomeCoverage::StartNewChrom(const string& newChrom) { // empty the previous chromosome and reserve new std::vector().swap(_currChromCoverage); - if (_visitedChromosomes.find(newChrom) != _visitedChromosomes.end()) { + if (_visitedChromosomes.find(newChrom) != _visitedChromosomes.end() && _currChromName != "") { cerr << "Input error: Chromosome " << _currChromName << " found in non-sequential lines. This suggests that the input file is not sorted correctly." << endl; @@ -143,6 +146,7 @@ void BedGenomeCoverage::AddBlockedCoverage(const vector &bedBlocks) { } + void BedGenomeCoverage::CoverageBed() { BED a; @@ -194,87 +198,154 @@ void BedGenomeCoverage::CoverageBed() { void BedGenomeCoverage::PrintFinalCoverage() { + if (_currChromName.length() > 0) { + ReportChromCoverage(_currChromCoverage, _currChromSize, + _currChromName, _currChromDepthHist); + } - - // process the results of the last chromosome. - ReportChromCoverage(_currChromCoverage, _currChromSize, - _currChromName, _currChromDepthHist); if (_eachBase == false && _bedGraph == false && _bedGraphAll == false) { ReportGenomeCoverage(_currChromDepthHist); } } -void BedGenomeCoverage::CoverageBam(string bamFile) { +void BedGenomeCoverage::CoverageBam(vector &_bedFiles) { ResetChromCoverage(); - // open the BAM file - BamReader reader; - if (!reader.Open(bamFile)) { - cerr << "Failed to open BAM file " << bamFile << endl; - exit(1); + std::vector::iterator bedFileIt = _bedFiles.begin(); + std::vector::iterator bedFileFend = _bedFiles.end(); +// int nfile = _bedFiles.size(); +// vector bedpos(nfile); // no seek tell in BamReader, later. need add function in Bamtools api and store readed position for speed up + + for(bedFileIt = _bedFiles.begin(); bedFileIt != bedFileFend; bedFileIt++){ + string bamFile = *bedFileIt; + + // open the BAM file + BamReader reader; + if (!reader.Open(bamFile)) { + cerr << "Failed to open BAM file " << bamFile << endl; + exit(1); + } + + // get header & reference information + string header = reader.GetHeaderText(); + RefVector refs = reader.GetReferenceData(); + + for(unsigned int ichrom=0;ichrom < refs.size(); ichrom++){ + _tovisitChromosomes.insert(refs.at(ichrom).RefName); + } + reader.Close(); } - // get header & reference information - string header = reader.GetHeaderText(); - RefVector refs = reader.GetReferenceData(); - - // load the BAM header references into a BEDTools "genome file" - _genome = new GenomeFile(refs); - // convert each aligned BAM entry to BED - // and compute coverage on B - BamAlignment bam; - while (reader.GetNextAlignment(bam)) { - // skip if the read is unaligned - if (bam.IsMapped() == false) - continue; - - // skip if we care about strands and the strand isn't what - // the user wanted - if ( (_filterByStrand == true) && - ((_requestedStrand == "-") != bam.IsReverseStrand()) ) - continue; - - // extract the chrom, start and end from the BAM alignment - string chrom(refs.at(bam.RefID).RefName); - CHRPOS start = bam.Position; - CHRPOS end = bam.GetEndPosition(false, false) - 1; - - // are we on a new chromosome? - if ( chrom != _currChromName ) - StartNewChrom(chrom); - - // add coverage accordingly. - if (!_only_5p_end && !_only_3p_end) { - bedVector bedBlocks; - // we always want to split blocks when a D CIGAR op is found. - // if the user invokes -split, we want to also split on N ops. - if (_obeySplits) { // "D" true, "N" true - GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true, true); + // iterative through chroms(refernce) + set::iterator iterChrom; + for(iterChrom = _tovisitChromosomes.begin(); iterChrom != _tovisitChromosomes.end(); iterChrom++){ + string _chromtoread = *iterChrom; + + for(bedFileIt = _bedFiles.begin(); bedFileIt != bedFileFend; bedFileIt++){ + _currBedFile ++; + string bamFile = *bedFileIt; + + // open the BAM file + BamReader reader; + if (!reader.Open(bamFile)) { + cerr << "Failed to open BAM file " << bamFile << endl; + exit(1); } - else { // "D" true, "N" false - GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true, false); + + // get header & reference information + string header = reader.GetHeaderText(); + RefVector refs = reader.GetReferenceData(); + + // load the BAM header references into a BEDTools "genome file" + _genome = new GenomeFile(refs); + + // skip if no chrom(reference) find in bam file + bool nochrominfile=0; + for(unsigned int irefid=0;irefid < refs.size(); irefid++){ + string tmpchrom(refs.at(irefid).RefName); + if(_chromtoread == tmpchrom){ + nochrominfile=1; + } + } + if(!nochrominfile){ + if(reader.IsOpen()){ + reader.Close(); + } + continue; + } + + // are we on a new chromosome? will reset Coverage + if ( _chromtoread.c_str() != _currChromName) + StartNewChrom(_chromtoread.c_str()); + + // convert each aligned BAM entry to BED + // and compute coverage on B + BamAlignment bam; + bool readedblock = 0; + while (reader.GetNextAlignment(bam)) { + // skip if the read is unaligned + if (bam.IsMapped() == false) + continue; + + // skip if we care about strands and the strand isn't what + // the user wanted + if ( (_filterByStrand == true) && + ((_requestedStrand == "-") != bam.IsReverseStrand()) ) + continue; + + // extract the chrom, start and end from the BAM alignment + string chrom(refs.at(bam.RefID).RefName); + + // skip rest of file if new chrom read over + if(chrom != _chromtoread && readedblock){ + reader.Close(); + break; + } + CHRPOS start = bam.Position; + CHRPOS end = bam.GetEndPosition(false, false) - 1; + + // add coverage accordingly. + if (!_only_5p_end && !_only_3p_end) { + bedVector bedBlocks; + // we always want to split blocks when a D CIGAR op is found. + // if the user invokes -split, we want to also split on N ops. + if (_obeySplits) { // "D" true, "N" true + GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true, true); + } + else { // "D" true, "N" false + GetBamBlocks(bam, refs.at(bam.RefID).RefName, bedBlocks, true, false); + } + AddBlockedCoverage(bedBlocks); + } + else if (_only_5p_end) { + int pos = ( !bam.IsReverseStrand() ) ? start : end; + AddCoverage(pos,pos); + } + else if (_only_3p_end) { + int pos = ( bam.IsReverseStrand() ) ? start : end; + AddCoverage(pos,pos); + } + // flag at least readed one for this chrom + if (!readedblock) { + readedblock = 1; + } + } + + if(reader.IsOpen()){ + // close the BAM + reader.Close(); } - AddBlockedCoverage(bedBlocks); - } - else if (_only_5p_end) { - int pos = ( !bam.IsReverseStrand() ) ? start : end; - AddCoverage(pos,pos); - } - else if (_only_3p_end) { - int pos = ( bam.IsReverseStrand() ) ? start : end; - AddCoverage(pos,pos); } } - // close the BAM - reader.Close(); PrintFinalCoverage(); } void BedGenomeCoverage::ReportChromCoverage(const vector &chromCov, const int &chromSize, const string &chrom, chromHistMap &chromDepthHist) { + chromHistMap tmpchromDepthHist; // use tmpchromDepthHist so not accumulate when report everyfile. if (_eachBase) { int depth = 0; // initialize the depth int offset = (_eachBaseZeroBased)?0:1; @@ -303,15 +374,17 @@ void BedGenomeCoverage::ReportChromCoverage(const vector &chromCov, const // maximum bin requested, then readjust the depth to be the max if (depth >= _max) { chromDepthHist[chrom][_max]++; + tmpchromDepthHist[chrom][_max]++; } else { chromDepthHist[chrom][depth]++; + tmpchromDepthHist[chrom][depth]++; } depth = depth - chromCov[pos].ends; } // report the histogram for each chromosome - histMap::const_iterator depthIt = chromDepthHist[chrom].begin(); - histMap::const_iterator depthEnd = chromDepthHist[chrom].end(); + histMap::const_iterator depthIt = tmpchromDepthHist[chrom].begin(); + histMap::const_iterator depthEnd = tmpchromDepthHist[chrom].end(); for (; depthIt != depthEnd; ++depthIt) { int depth = depthIt->first; unsigned int numBasesAtDepth = depthIt->second; diff --git a/src/genomeCoverageBed/genomeCoverageBed.h b/src/genomeCoverageBed/genomeCoverageBed.h index f2a877fd..f4fc110a 100644 --- a/src/genomeCoverageBed/genomeCoverageBed.h +++ b/src/genomeCoverageBed/genomeCoverageBed.h @@ -41,7 +41,7 @@ class BedGenomeCoverage { public: // constructor - BedGenomeCoverage(string bedFile, string genomeFile, + BedGenomeCoverage(vector bedFiles, string genomeFile, bool eachBase, bool startSites, bool bedGraph, bool bedGraphAll, int max, float scale, @@ -57,7 +57,7 @@ class BedGenomeCoverage { private: // data (parms) - string _bedFile; + vector _bedFiles; string _genomeFile; bool _bamInput; bool _eachBase; @@ -82,14 +82,18 @@ class BedGenomeCoverage { chromDepthMap _chromCov; string _currChromName ; vector _currChromCoverage; + vector _allChromCoverage; chromHistMap _currChromDepthHist; int _currChromSize ; + unsigned int _currBedFile ; + int _lastBedFile ; set _visitedChromosomes; + set _tovisitChromosomes; // methods void CoverageBed(); - void CoverageBam(string bamFile); + void CoverageBam(vector &_bamFiles); void LoadBamHeaderIntoGenomeFile(const string &bamFile); void ReportChromCoverage(const vector &, const int &chromSize, const string &chrom, chromHistMap&); void ReportGenomeCoverage(chromHistMap &chromDepthHist); diff --git a/src/genomeCoverageBed/genomeCoverageMain.cpp b/src/genomeCoverageBed/genomeCoverageMain.cpp index ca221ae6..47869769 100644 --- a/src/genomeCoverageBed/genomeCoverageMain.cpp +++ b/src/genomeCoverageBed/genomeCoverageMain.cpp @@ -30,7 +30,9 @@ int genomecoverage_main(int argc, char* argv[]) { bool showHelp = false; // input files + vector bedFiles; string bedFile; + string bedFilelist; string genomeFile; int max = INT_MAX; float scale = 1.0; @@ -75,6 +77,7 @@ int genomecoverage_main(int argc, char* argv[]) { if ((i+1) < argc) { haveBed = true; bedFile = argv[i + 1]; + bedFiles.push_back(bedFile); i++; } } @@ -83,6 +86,23 @@ int genomecoverage_main(int argc, char* argv[]) { haveBed = true; bamInput = true; bedFile = argv[i + 1]; + bedFiles.push_back(bedFile); + i++; + } + } + else if(PARAMETER_CHECK("-ibams", 6, parameterLength)) { + if ((i+1) < argc) { + haveBed = true; + bamInput = true; + bedFilelist = argv[i + 1]; + ifstream listFile; + listFile.open(bedFilelist.c_str()); + string line; + if (listFile.is_open()) { + while (getline(listFile, line)) { + bedFiles.push_back(line); + } + } i++; } } @@ -197,7 +217,7 @@ int genomecoverage_main(int argc, char* argv[]) { } if (!showHelp) { - BedGenomeCoverage *bc = new BedGenomeCoverage(bedFile, genomeFile, eachBase, + BedGenomeCoverage *bc = new BedGenomeCoverage(bedFiles, genomeFile, eachBase, startSites, bedGraph, bedGraphAll, max, scale, bamInput, obeySplits, filterByStrand, requestedStrand, @@ -222,7 +242,9 @@ void genomecoverage_help(void) { cerr << "Options: " << endl; + cerr << "\t-is\t\t" << "The input list file each line is . Report each file chrom by chrom and final all together." << endl; cerr << "\t-ibam\t\t" << "The input file is in BAM format." << endl; + cerr << "\t-ibams\t\t" << "The input list file each line is filename in BAM format. Report each file chrom by chrom and final all together." << endl; cerr << "\t\t\tNote: BAM _must_ be sorted by position" << endl << endl; cerr << "\t-d\t\t" << "Report the depth at each genome position (with one-based coordinates)." << endl; diff --git a/test/genomecov/test-genomecov.sh b/test/genomecov/test-genomecov.sh index 2c00c24d..5ed2c87d 100644 --- a/test/genomecov/test-genomecov.sh +++ b/test/genomecov/test-genomecov.sh @@ -144,4 +144,4 @@ check obs exp rm obs exp -rm *.bam \ No newline at end of file +#rm *.bam diff --git a/test/genomecov/three_blocks.sam b/test/genomecov/three_blocks.sam index 9bf1a627..0758b34a 100644 --- a/test/genomecov/three_blocks.sam +++ b/test/genomecov/three_blocks.sam @@ -1,3 +1,8 @@ @HD VN:1.0 GO:none SO:coordinate @SQ SN:chr1 LN:1000 +@SQ SN:chr10 LN:1000 +@SQ SN:chr19 LN:1000 three_blocks 16 chr1 1 40 10M10N10M10N10M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50 +three_blocks 16 chr1 1 40 30M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50 +three_blocks 16 chr10 1 40 10M10N10M10N10M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50 +three_blocks 16 chr19 1 40 10M10N10M10N10M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50 diff --git a/test/genomecov/three_blocks_match_1bp.bed b/test/genomecov/three_blocks_match_1bp.bed index de8b10d2..eed2f93a 100644 --- a/test/genomecov/three_blocks_match_1bp.bed +++ b/test/genomecov/three_blocks_match_1bp.bed @@ -1 +1,2 @@ chr1 10 60 three_blocks_nomatch 0 + 0 0 0 3 11,10,10, 0,20,40, +chr19 1 60 three_blocks_nomatch 0 + 0 0 0 3 11,10,10, 0,20,40, diff --git a/test/genomecov/two_blocks.sam b/test/genomecov/two_blocks.sam index a28928ab..62c4d18c 100644 --- a/test/genomecov/two_blocks.sam +++ b/test/genomecov/two_blocks.sam @@ -1,3 +1,7 @@ @HD VN:1.0 GO:none SO:coordinate @SQ SN:chr1 LN:1000 -two_blocks 16 chr1 1 40 15M10N15M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50 \ No newline at end of file +@SQ SN:chr2 LN:1000 +two_blocks 16 chr1 1 40 15M10N15M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50 +two_blocks 16 chr1 1 40 30M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50 +two_blocks 16 chr2 1 40 30M * 0 0 GAAGGCCACCGCCGCGGTTATTTTCCTTCA CCCDDB?=FJIIJIGFJIJHIJJJJJJJJI MD:Z:50 +