If you sort and index your BAM file:
samtools idxstats
. The third column is the number of Illumina reads mapped to each PacBio read. You can then pipe that to awk: samtools idxstats foo.bam | awk '{if($1 != "*" && $3 > 0) tot += 1}END{print tot}
. That will print the number of PacBio reads with at least 1 read mapped to them.
- You'll want
samtools depth
for that
samtools depth -aa foo.bam | awk '{if($1 != lread) {if(lread != "") {print lread, 100*lcov/llen}; lread = $1; lcov=0; llen=0}; llen += 1; if($3>0) {lcov +=1}}END{print lread, 100*lcov/llen}'
That won't give you the number of illumina reads per pacbio read, but you can get that from samtools idxstats
(column 3) and join
the results.
BTW, your usage of "coverage" is non-standard. You want "percentage coverage" rather than "coverage", which actually means "fold-coverage" (i.e., number of reads times their average length divided by the total genome/chromosome/etc. length).
If you instead want "fold-coverage" then the easiest way is to first know the length of your reads. Presuming you have 100 base reads (for paired-end reads, this would be the same, presuming the ends don't overlap):
samtools idxstats foo.bam | awk '{print $1, $3, 100*$3/$2}'
That will give you a table of PacBio read name, number of aligned reads, fold-coverage.