I have 4 bam files and a coordinate file in bed format. I want to plot all 4 bam files in one graph for each coordinate I have, and look for differences in coverage across the samples. How can I do this from Bioconductor packages?
You could use either Gviz or ggbio to do that. An example for using ggbio is here; you want to scroll down to code chunk 28. For Gviz the example is on p24 of the vignette.
But do note that both Gviz and ggbio have steep learning curves. Gviz has a much more comprehensive vignette, so that might be the way to go, but both packages have their strengths. You might also consider that eyeballing a coverage plot might not be the most useful thing to do. Coverage is a function of the number of reads that overlap a position and the total number of reads for a given sample, and plotting just the coverage ignores the total read depth. And eyeballometric measures are...subjective.
You might consider just using rsamtools to extract the read depth at each region defined in your bed file and then making tables of the average read depth, in counts/million. Or you could use simple base R graphics to plot each region. Unless you need something that looks good, fast and dirty might be the way to go.