Automated Plotting Of Seq. Reads In A Bam File
6
14
Entering edit mode
10.3 years ago
sa9 ▴ 840

Given a sorted BAM file, I need to plot the sequencing reads that cover every locus in a list of ~60-80 loci (with flanking regions).

Currently, I'm doing this manually using IGV (or any other BAM GUI viewers) by taking screen prints but this is not practical.

What would you recommend to generate such plots in a non-GUI mode? This is going to be part of a pipeline written in Python. I know IGV has some command line tools (called IGVtools) but they are for indexing and counting reads but nothing for viewing BAM file outside the GUI.

bam viewer visualization • 9.0k views
0
Entering edit mode

I think this is a great question, I asked in similar terms a while ago and didn't get a definitive answer: Quick Visual Inspection Of Mapping Images For A List Of Regions

7
Entering edit mode
10.3 years ago
Ryan Dale 4.9k

From Python directly, you could use this simple IGV wrapper from @brentp. Copied from the example usage of that script:

#example usage:

>>> igv = IGV()
>>> igv.genome('hg19')
'OK'

'OK'
>>> igv.go('chr1:45,600-45,800')
'OK'

#save as svg, png, or jpg
>>> igv.save('/tmp/r/region.svg')
'OK'
>>> igv.save('/tmp/r/region.png')
'OK'

# go to a gene name.
>>> igv.go('muc5b')
'OK'
>>> igv.sort()
'OK'
>>> igv.save('muc5b.png')
'OK'

# get a list of commands that will work as an IGV batch script.
>>> print "\n".join(igv.commands)
snapshotDirectory /tmp/igv
genome hg19
goto chr1:45,600-45,800
snapshotDirectory /tmp/r
snapshot region.svg
snapshot region.png
goto muc5b
sort base
snapshot muc5b.png

# Note, there will be some delay as the browser has to load the annotations
# at each step.

1
Entering edit mode

With Perl (or any language) perhaps the simplest approach is to programmatically construct a text file to be used as a script for IGV (http://www.broadinstitute.org/software/igv/batch). The Python stuff above just makes this more convenient from within Python code.

0
Entering edit mode

Simple and effective. Many thanks daler.

0
Entering edit mode

Glad it helps, but I just pointed it out -- all the credit goes to Brent Pedersen

0
Entering edit mode

Sorry, I'm not familiar at all with Phyton, just wondering how to read a file with say, 100 different positions at the genome and then input them into a loop to get the 100 different plots, ... If you know any Perl option to do this would be nice to know it, Many thanks in advance,

0
Entering edit mode

Sorry, I'm not familiar at all with Phyton, just wondering how to read a file with say, 100 different positions at the genome and then input them into a loop to get the 100 different plots, ... Also, If you know any Perl option to do this would be nice to know it, Many thanks in advanc

5
Entering edit mode
8.9 years ago
Rm 8.1k

Download the IGV : You can automate it without manually opening the IGV: If starting file is bed file, use "bedToIgv" from bedtools to generate the igv_batch.file for snapshots:

For example in the Linux: using indexed bam file

bash igv.sh input.bam -g hg19 -b igv_batch_script.txt


#sample igv_batch_script.txt ##

new
genome hg19
snapshotDirectory /path/to/igv_screenshots
goto chr19:59418052-59418053
sort base
collapse
snapshot chr19_59418052-59418053.png
goto chr1:15680529-15680532 chr3:5680521-5680533
sort base
collapse
snapshot chr1_15680529-15680532.chr3_5680521-5680533.png
exit

1
Entering edit mode

What's in the igv.sh? bedToIgv can only generate a batch file for the IGV interface. I don't know how to run this in command line. Thanks!

3
Entering edit mode
10.3 years ago
Gww ★ 2.7k

If you install and setup GBrowse2 you can programmatically generate images for each region (see here). I have used GBrowse2 before to do this with next generation sequencing, it's quite nice because the height of the image is not limited and you can add additional annotations such as coverage plots and gene annotations.

0
Entering edit mode

The images generated by GBrowse2 are pretty and very flexible.Thanks GWW. The only concern is about the difficulties of installing GBrowse2 which, IMO, can be cumbersome for many non-savvy end users.

0
Entering edit mode

@sa9: I agree, it can be challenging depending on your OS. It took a bit effort for me to get it to work on my MacBook.

2
Entering edit mode
10.3 years ago
Ido Tamir 5.2k

most genome browsers have scripting capability [IGB] and [IGV].
GBrowse2 also allows linking and by this scripting.

If you want something really small:

2. and compile them: g++ renderimage.cpp -O3 -o renderimage -lpng -L . -lbam
3. create image with: ./render_image --bam=mybam.bam --options=snps,pairs,arrows,single,faceaway,inversions,linkpairs,colordepth --ref=myref.fa --region="2L:1-200000" --png=2L.a.png

This way you can quickly create small images with the reads from your bam file with many viewing options. like this [example] (click on one button on the right).

0
Entering edit mode

Great! Thanks Ido , I didn't know about lookseq. I will give it a try.

1
Entering edit mode
10.3 years ago

I don't really have a definitive answer, I asked in similar terms a while ago but I think it would be great to have a definitive answer to the question:

Quick Visual Inspection Of Mapping Images For A List Of Regions

cat peaks.txt | while read i; do echo -e "g$i\n." | samtools tview my.bam my.ref.fa.gz; done  The "." at the end is to switch of the dots in the tview configuration. I use Ctrl+C to go to the next region. Is there a better way to do this? Maybe something that generates image files? ADD COMMENT 0 Entering edit mode I must've missed your question when I searched for 'BAM viewer' in BioStar. Thanks avilella. Although this is a practical solution , tview output is not that pretty, don't you think? :) ADD REPLY 0 Entering edit mode @sa9: agreed tview zooms-in at the base-by-base level, which sometimes may not be useful. ADD REPLY 0 Entering edit mode 4.7 years ago @2184687-1231-83: You can use tview for automatized analysis, using the samtools Version: 1.3.1 For instance, with this Bash script while read coor; do samtools tview -d T -p$coor <<Indexed Bam file>>  <<genome>> > \$coor.txt
done < Test_coordinates.txt


The Test_coordinates.txt file could be chr9:86584727 chr15:89155087 chr19:4770712