7
$\begingroup$

I am constructing a bit of software which pipes the outputs of the bam file via samtools view into a script for parsing. My goal is to (somehow) make this process more efficient, and faster than samtools view.

I am only using 3-4 fields in the bam. Naturally, if I only wanted those fields, I could do something like samtools view file1.bam | cut -f ##.

However, I'm trying to be faster than this approach. More specifically:

(1) Given the bam file format, is it theoretically possible to output several requested SAM fields faster than samtools view|cut -f?

(2) Is there software available to achieve (1), or would one need to hack samtools?

Improve this question
$\endgroup$
5
  • 1
    $\begingroup$ Which fields do you need exactly? Using htslib is relatively easy if you know C, but if you can get by with something like sambamba then that'd be vastly easier still. $\endgroup$
    – Devon Ryan
    Aug 22, 2017 at 21:41
  • $\begingroup$ @DevonRyan Outputting the bam could be faster with sambambda? These are "mandatory" fields, so they're present in each bam "row" $\endgroup$
    – ShanZhengYang
    Aug 22, 2017 at 22:22
  • 1
    $\begingroup$ samtools view, especially in the way it is used here, produces a text representation of a BAM file. Using a library (or tool that taps into a library) allows faster access at the binary level by avoiding full text conversion of each alignment prior to filtering. $\endgroup$
    – gringer ♦
    Aug 22, 2017 at 22:39
  • $\begingroup$ @gringer I'm not following this exactly. Could you provide more details/an explanation? $\endgroup$
    – ShanZhengYang
    Aug 22, 2017 at 22:52
  • $\begingroup$ Sambamba is generally faster than samtools, or at least that used to be the case. Samtools is a reference implementation, it doesn't always have to be the fastest. $\endgroup$
    – Devon Ryan
    Aug 23, 2017 at 4:53

3 Answers 3

Reset to default
6
$\begingroup$

The BAM file format is not a text-based format. It has a specific binary structure, specified in reasonable detail in the SAM file format specification. Whenever this information is displayed on a screen as text, it needs to be converted from the binary format to a text format, which takes a bit of time and processing power.

As this question suggests, if only a specific field from alignments is needed (or fields), then it will probably be better to extract just those fields and do any necessary conversion only on those fields. While this can be done by writing a BAM parser from scratch, many developers have already written software libraries to process BAM files in this way.

Devon Ryan has suggested htslib, which is a C library written by a group that includes the people who wrote the SAM/BAM/CRAM file format specifications. There's also pysam, which is a python wrapper around htslib.

The particular tool that is used will depend on your familiarity with programming and the specific thing that you want to do. If you just want to "output a BAM file", then cat is the quickest:

cat file1.bam

... but you probably don't want to do that, because it seems like you want to process a text-based representation with a script. Because you haven't specified which fields you are interested in, it's not possible to suggest the best thing to use. In the end, I expect that an "efficient" solution to your problem would involve htslib in some form. This is not really hacking samtools, it's using the backend of samtools to process BAM data.

However, it's useful to also think about the time cost of coding. What application do you have which means that the text-based processing method is not sufficient? It takes time to write code, and a lot more time to debug that code to make sure it's doing the right thing in all situations. If this is for a one-off thing, then samtools view output fed into your script may be the quickest solution.

Improve this answer
$\endgroup$
2
  • $\begingroup$ I'm a bit familiar with pysam. It seems like this would be very inefficient, as I would need to input the bam into a python data structure, e.g. ` pysam.AlignmentFile()`. As for the script, yes, I currently am parsing the text-format of the outputted bam file. I have done nothing with compression/binary yet. $\endgroup$
    – ShanZhengYang
    Aug 23, 2017 at 4:02
  • 1
    $\begingroup$ Whether or not it is inefficient depends on the way the library is implemented. It's possible for data structures to contain pointers to files (or places within files) and only retrieve information from the files when requested. There's no need to load the entire file into memory to explore just one aspect of it. $\endgroup$
    – gringer ♦
    Aug 23, 2017 at 21:49
6
$\begingroup$

I modified your original question: as you are extracting 4 fields, you are not outputting BAM. The answer to the modified question is: yes, you can write a C program with htslib (or with bamtools, bioD, bioGo or rust-bio). Formatting an entire SAM is fairly expensive. You can see this by comparing samtools view aln.bam > /dev/null and samtools view -u aln.bam > /dev/null. With a C program, you can select fields to output. This will give you a noticeable performance boost, depending on the fields you extract.

In general, if you really care about performance, don't use script. The time you save from optimizing samtools view|cut -f will be quickly neutralized by the inefficiency of most scripting languages.

Improve this answer
$\endgroup$
2
  • $\begingroup$ > With a C program, you can select fields to output. Can you detail a bit more how this would be done? The BAM format is in binary, so one would have to convert to text first and then extract the fields desired? $\endgroup$
    – ShanZhengYang
    Aug 29, 2017 at 10:05
  • 2
    $\begingroup$ @ShanZhengYang The library parses the file record for you and places it in memory. From there, you print (or otherwise process) from the in-memory data only your desired fields in whatever format you want. There is no text parsing required. $\endgroup$
    – mattm
    Aug 29, 2017 at 15:36
3
$\begingroup$

Another library you can use for this purpose is the htsjdk, which is written in java. Using htsjdk with java is analogous to using htslib with C; the BAM format is already handled by the library and you can manipulate fields in your own code. The same basic analysis applies to java as C; you do not need to convert the BAM file to text and then parse it. Generally, a program implemented in java will be slower than C but still significantly faster than scripting. You do get language features of java: portability, array bounds checking, garbage collection, etc.

Improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.