Skip to main content

Breaking BAM

It's fun to mess with our bioinformatics tools and laugh at ourselves. The BAM format carries a query name string. In an idle moment I wondered, how long a string can I put here before bad things happen to the BAM?

Generally this string carries information related to the device that produced the read, lane number, run number and so and so forth. It's completely specification free and everyone encodes information here in their own way, so you would think it's an arbitrary length string. Almost.

Try out the short Python snippet below. It generates a dummy fasta file (test.fa) and a dummy BAM file (test.bam) with one 'aligned' read. The only funky thing is that you can make the qname field of the read as long as you want (by varying the lone input parameter to the script).

import sys
import pysam

  qnamelen = int(sys.argv[1])
  qnamelen = 255

with open('test.fa', 'w') as f:

bam_hdr = {'HD': {'VN': '1.4'},
           'SQ': [{'LN': 100, 'SN': 'test'}]}
bf = pysam.Samfile('test2.bam', 'wb', header=bam_hdr)  # Write binary BAM with header
ar = pysam.AlignedRead()
ar.qname = 'Q'*qnamelen  # Gives Tablet headaches when str len > 254
ar.seq = 'A'*50
ar.qual = '~'*50
ar.mapq = 100
ar.pos = 0
ar.cigarstring = '50M'
pysam.sort('test2.bam', 'test')

Now open the pair of fasta and BAM files in samtools tview, Tablet and IGV (and indeed in anything else). Laugh loudly until your office mates banish you.

Only samtools tview - the least fancy one - loads the alignment without batting an eyelid.

The specs allow an arbitrary length null terminated string.

(Oh, as an aside, try this. Create a .fasta file with this simple corruption: Remove the newline from the sequence id in the fasta file (so that the sequence looks something like 'testAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA'. Load the sequence in IGV. Now fix the sequence and try and convince IGV to load the new corrected sequence. Hint: it won't go away until you use igvtools to reindex the fasta file. running bwa index won't help)


Popular posts from this blog

Flowing text in inkscape (Poster making)

You can flow text into arbitrary shapes in inkscape. (From a hint here).

You simply create a text box, type your text into it, create a frame with some drawing tool, select both the text box and the frame (click and shift) and then go to text->flow into frame.


The omnipresent anonymous asked:
Trying to enter sentence so that text forms the number three...any ideas?
The solution:
Type '3' using the text toolConvert to path using object->pathSize as necessaryRemove fillUngroupType in actual text in new text boxSelect the text and the '3' pathFlow the text

Pandas panel = collection of tables/data frames aligned by index and column

Pandas panel provides a nice way to collect related data frames together while maintaining correspondence between the index and column values:

import pandas as pd, pylab #Full dimensions of a slice of our panel index = ['1','2','3','4'] #major_index columns = ['a','b','c'] #minor_index df = pd.DataFrame(pylab.randn(4,3),columns=columns,index=index) #A full slice of the panel df2 = pd.DataFrame(pylab.randn(3,2),columns=['a','c'],index=['1','3','4']) #A partial slice df3 = pd.DataFrame(pylab.randn(2,2),columns=['a','b'],index=['2','4']) #Another partial slice df4 = pd.DataFrame(pylab.randn(2,2),columns=['d','e'],index=['5','6']) #Partial slice with a new column and index pn = pd.Panel({'A': df}) pn['B'] = df2 pn['C'] = df3 pn['D'] = df4 for key in pn.items: print pn[key] -> output …

Drawing circles using matplotlib

Use the pylab.Circle command

import pylab #Imports matplotlib and a host of other useful modules cir1 = pylab.Circle((0,0), radius=0.75, fc='y') #Creates a patch that looks like a circle (fc= face color) cir2 = pylab.Circle((.5,.5), radius=0.25, alpha =.2, fc='b') #Repeat (alpha=.2 means make it very translucent) ax = pylab.axes(aspect=1) #Creates empty axes (aspect=1 means scale things so that circles look like circles) ax.add_patch(cir1) #Grab the current axes, add the patch to it ax.add_patch(cir2) #Repeat