Question: Automated Plotting Of Seq. Reads In A Bam File
12
gravatar for sa9
6.3 years ago by
sa9750
USA, Cambridge
sa9750 wrote:

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.

viewer bam visualization • 5.8k views
ADD COMMENTlink modified 8 months ago by aiagutierrezdi0 • written 6.3 years ago by sa9750

I think this is a great question, I asked in similar terms a while ago and didn't get a definitive answer: http://biostar.stackexchange.com/questions/5580/quick-visual-inspection-of-mapping-images-for-a-list-of-regions

ADD REPLYlink written 6.3 years ago by 2184687-1231-83-4.8k
7
gravatar for Ryan Dale
6.3 years ago by
Ryan Dale4.6k
Bethesda, MD
Ryan Dale4.6k wrote:

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 COMMENTlink written 6.3 years ago by Ryan Dale4.6k

Simple and effective. Many thanks daler.

ADD REPLYlink written 6.3 years ago by sa9750

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

ADD REPLYlink written 6.3 years ago by Ryan Dale4.6k

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 REPLYlink written 6.3 years ago by J.Rodrigo Flores40

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 REPLYlink written 6.3 years ago by J.Rodrigo Flores40

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 REPLYlink written 6.3 years ago by Ryan Dale4.6k
4
gravatar for Rm
4.9 years ago by
Rm7.5k
Danville, PA
Rm7.5k wrote:

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 COMMENTlink written 4.9 years ago by Rm7.5k

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 REPLYlink written 10 weeks ago by niuyw0
3
gravatar for Gww
6.3 years ago by
Gww2.6k
Canada
Gww2.6k wrote:

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 COMMENTlink written 6.3 years ago by Gww2.6k

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 REPLYlink written 6.3 years ago by sa9750

@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 REPLYlink written 6.3 years ago by Gww2.6k
2
gravatar for Ido Tamir
6.3 years ago by
Ido Tamir4.8k
Austria
Ido Tamir4.8k wrote:

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 COMMENTlink written 6.3 years ago by Ido Tamir4.8k

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

ADD REPLYlink written 6.3 years ago by sa9750
1
gravatar for 2184687-1231-83-
6.3 years ago by
2184687-1231-83-4.8k wrote:

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:

http://biostar.stackexchange.com/questions/5580/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 COMMENTlink written 6.3 years ago by 2184687-1231-83-4.8k

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 REPLYlink written 6.3 years ago by sa9750

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

ADD REPLYlink written 6.3 years ago by 2184687-1231-83-4.8k
0
gravatar for aiagutierrezdi
8 months ago by
aiagutierrezdi0 wrote:

@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 COMMENTlink written 8 months ago by aiagutierrezdi0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 1077 users visited in the last hour