This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Get A Graphical Representation Of Number Of Reads In Bam File For A Certain Genomic Region

Hi guys and girls,

One of my coworkers is scrambling to redo a figure for a manuscript under review. The reviewers want to see a graphical representation of the number of reads from an RNA seq experiment. My coworkers is a biologists and computers are not his strong suit. We received a BAM file from the company who did the sequeuncing. However, they are not very responsive about getting us the data needed for the figure. I am probably the most computer literate in the lab, but I have never dealt with next-gen sequencing data, so I ma completely at a loss.

From some of the other questions on the forum, I gathered that Samtools could to this. We are interested in plotting the number of reads in "Chr2: 212,000,000 – 213,800,000". I assume, you would do some kind of sliding window in order to bin the reads?

Any pointers would be greatly appreciated!

example of desired output

EDIT: thanks for all the help. Each of the proposed options seems to come with its own little problems. I am sure I'll figure it out though!

bam samtools
• 10,068 views • link
updated by Biostar Biostar • written by Whetting Whetting

4 answers

upload your BAM (or a sub part of the BAM) in the UCSC as a custom track and download the figure(s) as PDF/PS.

• 84 views • link
written by Pierre Lindenbaum Pierre Lindenbaum

I would go with this method, but I would first convert to a coverage plot to reduce the file sizes. I would recommend bedtools genomecov, using the -ibam flag to import the BAM file. With the resulting bedGraph output you can convert it to bigWig (search 'UCSC bigwig'). It is often problematic trying to directly view a BAM file in UCSC, plus a coverage plot is more aesthetically pleasing in a figure.

• 312 views • link
written by Ian Ian

Try IGV, from the Broad institute. You can import indexed .bams to it, and zoom in on the gene, and see the reads pileup up on exons.

• 196 views • link
written by swbarnes2 swbarnes2

Could be a bit tricky given a 1.8Mb region. My IGV stops displaying alignments at regions > 0.1Mb roughly.

• 0 views • link
written by Christof Winter Christof Winter

That setting can be changed under View -> Preferences -> Alignments -> Visibility range threshold. It is set at 30kb by default. With enough memory available you might be able to increase to 2000kb to encompass this region? I tried it on my laptop (for a bam with ~50X coverage) with 2MB of ram allocated for IGV and it was insufficient though. I'll try again later on a bigger memory machine and see if this is possible. Increasing downsampling might help also? But, given the challenges some of the other solutions suggested in this thread might be better.

• 261 views • link
written by Obi Griffith Obi Griffith

Good point, thanks! That works on my workstation (one 5GB RNA-Seq BAM file using ~500MB main memory).

• 0 views • link
written by Christof Winter Christof Winter

How would I get the image out of IGV?

• 0 views • link
written by Whetting Whetting

File -> save image in IGV will save a snapshot of what every you are looking at in IGV. I often just take snapshots using mac's built in screenshot options.

• 0 views • link
written by Obi Griffith Obi Griffith

If you are interested in the number of reads per position to then plot them with a tool of your choice, samtools is indeed your friend:

samtools depth -r chr2:212,000,000-213,800,000 accepted_hits.bam > reads-per-position.txt

will produce 

chr2    212107123    1
chr2    212107124    1
chr2    212107125    1
chr2    212107126    2
chr2    212107127    2
chr2    212107128    2
...

The last column contains the number of reads covering the given position. Note that positions with zero coverage are omitted. You could easily go ahead and turn the output into a BED file and then use bedtools intersect to restrict your positions exons only (using e.g. a GFF or GTF file containing the exon coordinates).

• 134 views • link
written by Christof Winter Christof Winter

is this faster than samtools mpileup? That command gives depth at each position along with some additional information

• 0 views • link
written by Bioinfosm Bioinfosm

For RNA-Seq BAM files that have N's in the CIGAR string to represent introns, samtools mpileup unfortunately also counts the intronic gaps in reads spanning two exons as read depth (which is normally not what you want). Also, by default it ignores reads with an unmapped mate (to get these, add the -A option).

• 0 views • link
written by Christof Winter Christof Winter

what we normally do in our NGS pipeline is to generate a bigwig coverage file for each BAM file, which can be visualized directly on IGV. here is all the code needed:

genomeCoverageBed -bg -ibam input.bam > coverage.bedgraph
bedGraphToBigWig coverage.bedgraph human_hg19.fa.fai coverage.bigwig

you will have to have installed bedtools and the UCSC's bedGraphToBigWig, plus the index of the reference used (human_hg19.fa.fai in our case)

• 86 views • link
written by Jorge Amigo Jorge Amigo

cool, how do you get the image out of IGV?

• 0 views • link
written by Whetting Whetting

File -> save image. Or, take a screenshot.

• 0 views • link
written by Obi Griffith Obi Griffith

true. IGV natively provides a PNG image as output, although it won't be as good as the PDF/PS image that Pierre mentioned.

• 41 views • link
written by Jorge Amigo Jorge Amigo

You can save as SVG as well, so you'll get a vector graphic.

• 0 views • link
written by Christof Winter Christof Winter

Yes. You can choose between .png, .jpeg, and .svg (under 'file format') when you save image.

• 0 views • link
written by Obi Griffith Obi Griffith

Log in to answer this question.

More posts like this