MATLAB: Is BioMap slow for indexing BAM files in R2016b

bamtoolsbioBioinformatics Toolboxcomputational

I am using "BioMap" to index very large "BAM" files. However, for files with several thousand references, it is taking longer than a day to finish. Why is this happening?

Best Answer

"BAM" files (and "SAM" files, their ASCII encoded cousin) can be very large. As a consequence they must be indexed and analyzed on disk. After indexing, there are two additional files that the "BioMap" class uses to access the data efficiently from disk, which have the extensions "BAI" and "LINEARINDEX". However, if both of these files are not present, "BioMap" must construct them. This is very computationally intensive for files containing many reference sequences, as their is a large amount of I/O that must occur in addition to the indexing analysis.

If the references in the "BAM" file are ordered, and the "BAI" file is present, it may be possible to achieve faster results using the "bamread" function.

Related Solutions

MATLAB: How to convert a BAM file to SAM and GFF file to GTF for use in featurecount

The featurecount function currently only supports SAM and GTF files. The only workaround is to convert your files to these formats.

Using a BioMap object is currently the most straightforward way to convert from SAM to BAM using the Bioinformatics Toolbox.

b = BioMap('test.bam');
write(b,'test.sam','format','sam')

A direct conversion from GFF3 to GTF may not be possible without some simplifying assumptions, and there is currently no way to do this using the Bioinformatics Toolbox. Since GTF is identical to GFF2, you may be able to give the filename of a GFF2 file directly to featurecount.

If you have a third party conversion tool, you can call this tool from within MATLAB by using the ! character, as in this example.

if exist('test.sam') ~= 2
    !thirdpartyconversiontool test.bam test.sam
end
if exist('test.gtf') ~= 2
    !thirdpartyconversiontool test.gff test.gtf
end

MATLAB: How to create/find a bam index file for using the bamread function

You can use a BioMap object to load the information from the BAM file without the BAI file. See instructions here: http://www.mathworks.com/help/bioinfo/ug/storing-and-managing-short-read-sequence-data-in-objects.html

For example if you want to load the data from chromosome 1, you do this:

info = baminfo(bamFilename,'ScanDictionary',true);
bm = BioMap(bamFilename, 'SelectReference', 'ChrI');

Instead of 'ChrI' you can put any of the references that were used for mapping. These are stored here:

info.ScannedDictionary

The BioMap line will also generate a BAI file which you can use in bamread but that usually takes a lot more memory and it is slower.

Best Answer

Related Solutions

MATLAB: How to convert a BAM file to SAM and GFF file to GTF for use in featurecount

MATLAB: How to create/find a bam index file for using the bamread function

Related Question