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
ADD COMMENT
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

ADD REPLY
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'

    >>> 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.
ADD COMMENT
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.

ADD REPLY
0
Entering edit mode

Simple and effective. Many thanks daler.

ADD REPLY
0
Entering edit mode

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

ADD REPLY
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,

ADD REPLY
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

ADD REPLY
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
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
ADD COMMENT
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!

ADD REPLY
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.

ADD COMMENT
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.

ADD REPLY
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.

ADD REPLY
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:

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

ADD COMMENT
0
Entering edit mode

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

ADD REPLY
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

ADD COMMENT

Login before adding your answer.

Traffic: 2457 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6