Automated Plotting Of Seq. Reads In A Bam File

Question

[sa9]

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.

[2184687-1231-83-]

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

Answer 1

[Ryan Dale]

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'

    >>> igv.load('http://www.broadinstitute.org/igvdata/1KG/pilot2Bams/NA12878.SLX.bam')
    '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.

[Ryan Dale]

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.

[sa9]

Simple and effective. Many thanks daler.

[Ryan Dale]

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

[J.Rodrigo Flores]

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,

[J.Rodrigo Flores]

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

Answer 2

[Rm]

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
load input.bam
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

[niuyw]

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!

Answer 3

[Gww]

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.

[sa9]

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.

[Gww]

@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.

Answer 4

[Ido Tamir]

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

If you want something really small:

1. download the two files from lookseq [lookseq]
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).

[sa9]

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

Answer 5

[2184687-1231-83-]

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?

[sa9]

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? :)

[2184687-1231-83-]

@sa9: agreed tview zooms-in at the base-by-base level, which sometimes may not be useful.

Answer 6

[aiagutierrezdi]

@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