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

Python: Multiprocessing: passing multiple arguments to a function

Write a wrapper function to unpack the arguments before calling the real function. Lambda won't work, for some strange un-Pythonic reason.

import multiprocessing as mp def myfun(a,b): print a + b def mf_wrap(args): return myfun(*args) p = mp.Pool(4) fl = [(a,b) for a in range(3) for b in range(2)] #mf_wrap = lambda args: myfun(*args) -> this sucker, though more pythonic and compact, won't work, fl)

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

Running a task in a separate thread in a Tkinter app.

Use Queues to communicate between main thread and sub-threadUse wm_protocol/protocol to handle quit eventUse Event to pass a message to sub-threadimport Tkinter as tki, threading, Queue, time def thread(q, stop_event): """q is a Queue object, stop_event is an Event. stop_event from """ while(not stop_event.is_set()): if q.empty(): q.put(time.strftime('%H:%M:%S')) class App(object): def __init__(self): self.root = tki.Tk() = tki.Text(self.root, undo=True, width=10, height=1)'left') self.queue = Queue.Queue(maxsize=1) self.poll_thread_stop_event = threading.Event() self.poll_thread = threading.Thread(target=thread, name='Thread', args=(self.queue,self.poll_thread_stop_event)) self.poll_thread.start() self.poll_interval = 250 self.poll() self.root.wm_protocol("WM_DELETE…