Okay. Is this what you want? Check out Gviz Bioconductor package specifically section 4.9: Alignment track.

For ex:

library(Gviz) #load library
alTrack=AlignmentsTrack("foo_sorted.bam",isPaired=T) #Read bam file
plotTracks(alTrack,from=87072069,to=87073068,chromosome="chr4") #plot the required region; plots coverage, alignment section
plotTracks(alTrack,from=87072069,to=87073068,chromosome="chr4",type="coverage") #plot only coverage, no depth or alignment section

GATK DepthOfCoverage or Bedtools genomecov are probably your best bet

Try this

genomeCoverageBed -ibam test.bam -g my.genome -d > perbase_cov.txt

where my.genome is the genomic file.

Genome files must be tab-delimited and are structured as follows (this is an example for C. elegans):

chrI 15072421
chrII 15279323
...
chrX 17718854
chrM 13794

2nd column is the length of the chromosome.

Refer bedtools

Varun

Since you have "RNA-seq" among your tags I guess this is what you want to plot. In this case I would use bedtools genomeCovergeBed as others have suggested, but with the -split option to correctly represent reads spanning different exons. From the docs:

bedtools genomecov will, by default, screen for overlaps against the entire span of a spliced/split BAM alignment or blocked BED12 feature. When dealing with RNA-seq reads, for example, one typically wants to only screen for overlaps for the portions of the reads that come from exons (and ignore the interstitial intron sequence). The -split command allows for such overlaps to be performed

Thanks for the responses. I actually ended up using:

samtools depth file.bam > coverage.bam

This command gave me the exact numbers I was seeing in IGV on a per-nucleotide basis