<p>If you sort and index your BAM file:</p>

<ol>
<li><code>samtools idxstats</code>. The third column is the number of Illumina reads mapped to each PacBio read. You can then pipe that to awk: <code>samtools idxstats foo.bam | awk '{if($1 != "*" &amp;&amp; $3 &gt; 0) tot += 1}END{print tot}</code>. That will print the number of PacBio reads with at least 1 read mapped to them.</li>
<li>You'll want <code>samtools depth</code> for that</li>
</ol>

<pre class="lang-bash prettyprint-override"><code>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&gt;0) {lcov +=1}}END{print lread, 100*lcov/llen}'
</code></pre>

<p>That won't give you the number of illumina reads per pacbio read, but you can get that from <code>samtools idxstats</code> (column 3) and <code>join</code> the results.</p>

<p>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).</p>

<p>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):</p>

<pre><code>samtools idxstats foo.bam | awk '{print $1, $3, 100*$3/$2}'
</code></pre>

<p>That will give you a table of PacBio read name, number of aligned reads, fold-coverage.</p>

CC BY-SA 3.0
2018-01-17T18:48:45.747Z
2018-01-18T08:01:05.607Z
2018-01-18T08:01:05.607Z
Devon Ryan
<p>This infomation can be obtained with </p>

<pre><code>samtools index mybam.bam
samtools idxstats mybam.bam
</code></pre>

<p>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.</p>

<p>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</p>

<pre><code>samtools idxstats mybam.bam | awk '$3&gt;0' | wc -l
</code></pre>

CC BY-SA 3.0
2018-01-17T18:00:32.357Z
2018-01-17T18:00:32.357Z
Ian Sudbery
<p>I mapped raw illumina reads to longer pacbio reads and I would like to know the following information from my mapping file (SAM/BAM)</p>

<ul>
<li><p>How many PacBio reads are mapped to at least one illumina reads </p></li>
<li><p>A Table showing each pacbio read, the number of illumina reads that
mapped, and total coverage (i.e. percentage of mapped regions)</p></li>
</ul>

CC BY-SA 3.0
2018-01-17T21:25:44.047Z
By "percentage of mapped regions", do you mean the percentage of the PacBio read that was covered/mapped by illumina reads?
gringer
<read-mapping><samtools><sam>
Mapping Information from SAM/BAM file
# Mapping Information from SAM/BAM file

(Asked by: username_0 on 2018-01-17)

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)

**Tags:** read-mapping, samtools, sam
> **Comments:**
> **username_1**: By "percentage of mapped regions", do you mean the percentage of the PacBio read that was covered/mapped by illumina reads?  

---

## 2 answers

### Answer by username_2 on 2018-01-17. Score: 1 | Accepted answer

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.


---

### Answer by username_3 on 2018-01-17. Score: 0

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


---

