Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-01-17T18:48:45.747Z
Answers.LastActivityDate: 2018-01-18T08:01:05.607Z
Answers.LastEditDate: 2018-01-18T08:01:05.607Z
Answers.UserDisplayName: Devon Ryan
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-01-17T18:00:32.357Z
Answers.LastActivityDate: 2018-01-17T18:00:32.357Z
Answers.UserDisplayName: Ian Sudbery
Body: <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>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-01-17T21:25:44.047Z
Comments.Text: By "percentage of mapped regions", do you mean the percentage of the PacBio read that was covered/mapped by illumina reads?
Comments.UserDisplayName: gringer
Tags: <read-mapping><samtools><sam>
Title: Mapping Information from SAM/BAM file
ThreadText: # 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:
  
  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
  
  ```bash
  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
  
  ```
  
  
  ---
  
  