Mapping Information from SAM/BAM file

Question

I mapped raw illumina reads to longer pacbio reads and I would like to know the following information from my mapping file (SAM/BAM)

• How many PacBio reads are mapped to at least one illumina reads

• A Table showing each pacbio read, the number of illumina reads that mapped, and total coverage (i.e. percentage of mapped regions)

Answer 1

If you sort and index your BAM file:

1. 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.
2. 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.

Answer 2

This infomation can be obtained with

samtools index mybam.bam
samtools idxstats mybam.bam

You'll get a table with one row per PacBio read, the length, the number of mapped reads aligned to it and the number of unmapped reads aligned to it.

The number of PacBio reads mapped to at least one illumina read is the number of rows where the 3rd column is greater than 0

samtools idxstats mybam.bam | awk '$3>0' | wc -l