pezmaster31/bamtools — Issue #38: Slow performance when using FASTA reference genome.
Author: colinhercus
Created: 2012-03-08 03:32:40+00:00
State: open

Hi,

I was running bamtools pileup with a fasta input and it was really slow. I tracked this down to fasta module which retrieved the full reference sequence each time one base was required. I mod'd this so it would calculate the offset of the base req'd and then seek and getc to retrieve just one byte. It improved the performance of pileup dramatically when using Human genome.

Kind Regards, Colin

```
diff --git a/src/utils/bamtools_fasta.cpp b/src/utils/bamtools_fasta.cpp
index d3ad080..f3d27fe 100644
--- a/src/utils/bamtools_fasta.cpp
+++ b/src/utils/bamtools_fasta.cpp
@@ -251,21 +251,25 @@ bool Fasta::FastaPrivate::GetBase(const int& refId, const int& position, char& b
return false;
}

- // seek to beginning of sequence data
- if ( fseeko(Stream, referenceData.Offset, SEEK_SET) != 0 ) {
- cerr << "FASTA error : could not sek in file" << endl;
- return false;
- }
-
- // retrieve sequence
- string sequence = "";
- if ( !GetNextSequence(sequence) ) {
- cerr << "FASTA error : could not retrieve base from FASTA file" << endl;
+ int64_t lines = position / referenceData.LineLength;
+ int64_t lineOffset = position % referenceData.LineLength;
+ int64_t SeekTo = referenceData.Offset + lines * referenceData.ByteLength + lineOffset;
+
+ // seek to required base
+ if ( fseeko(Stream, SeekTo, SEEK_SET) != 0 ) {
+ cerr << "FASTA error : could not seek in file" << endl;
return false;
}
-
- // set base & return success
- base = sequence.at(position);
+ base = getc(Stream);
return true;
}