Hi, 

I have paired-end bam files (exomes) for two groups each with multiple individuals (all from one species, non-viral). I am interested to get information as mentioned here http://www.ark-genomics.org/bioinformatics/virome  Would it be possible to do with `Rsamtools`? I am in the process of learning `Rsamtools` and getting stuck with how to process all files together to get final output. My aim is to retrieve genes and SNPs from this species as compared to reference (with same genus). Kindly guide.

Thanks!

The function idxstatsBam() summarized reads per sequence name quickly. I would not call this statistic 'coverage', but it would be sufficient for general QA or visualization purposes.

The rest of your question is too vague -- there is a ton of information at the URL you provide, and there is no magic bullet to produce all of it in a single command.