Answers.Body: <p><code>samtools stats</code> seems to be able to do most of this, excluding the CIGAR-string parsing stuff (i.e. INDELs):</p>
  
  <pre><code>$ samtools view -h mapped.bam | grep -e '^@' -e 'readName' |
      samtools stats | grep '^SN' | cut -f 2-
  
  raw total sequences:    2
  filtered sequences: 0
  sequences:  2
  is sorted:  1
  1st fragments:  2
  last fragments: 0
  reads mapped:   2
  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
  reads unmapped: 0
  reads properly paired:  0   # proper-pair bit set
  reads paired:   0   # paired-end technology bit set
  reads duplicated:   0   # PCR or optical duplicate bit set
  reads MQ0:  0   # mapped and MQ=0
  reads QC failed:    0
  non-primary alignments: 0
  total length:   13632   # ignores clipping
  bases mapped:   13632   # ignores clipping
  bases mapped (cigar):   13502   # more accurate
  bases trimmed:  0
  bases duplicated:   0
  mismatches: 3670    # from NM fields
  error rate: 2.718116e-01    # mismatches / bases mapped (cigar)
  average length: 6816
  maximum length: 6816
  average quality:    8.9
  insert size average:    0.0
  insert size standard deviation: 0.0
  inward oriented pairs:  0
  outward oriented pairs: 0
  pairs with other orientation:   0
  pairs on different chromosomes: 0
  </code></pre>
  
  <p>For a length filter, an awk length check can be slotted in instead of the grep:</p>
  
  <pre><code>$ samtools view -h mapped.bam |
      awk -F'\t' '{if((/^@/) || (length($10)&gt;500)){print $0}}' |
      samtools stats | grep '^SN' | cut -f 2-
  
  raw total sequences:    131
  filtered sequences: 0
  sequences:  131
  is sorted:  1
  1st fragments:  131
  last fragments: 0
  reads mapped:   131
  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
  reads unmapped: 0
  reads properly paired:  0   # proper-pair bit set
  reads paired:   0   # paired-end technology bit set
  reads duplicated:   0   # PCR or optical duplicate bit set
  reads MQ0:  0   # mapped and MQ=0
  reads QC failed:    0
  non-primary alignments: 0
  total length:   569425  # ignores clipping
  bases mapped:   569425  # ignores clipping
  bases mapped (cigar):   453884  # more accurate
  bases trimmed:  0
  bases duplicated:   0
  mismatches: 113755  # from NM fields
  error rate: 2.506257e-01    # mismatches / bases mapped (cigar)
  average length: 4346
  maximum length: 16909
  average quality:    9.1
  insert size average:    0.0
  insert size standard deviation: 0.0
  inward oriented pairs:  0
  outward oriented pairs: 0
  pairs with other orientation:   0
  pairs on different chromosomes: 0
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T07:31:35.467Z
Answers.Text: all extra points awarded. This is (hopefully) super useful for folks working on long reads, and wanting to summarise the longer bits of their data...
Answers.UserDisplayName: roblanf
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T10:10:00.290Z
Answers.Text: @terdon please don't edit this to make it less readable than it already is. I'm not interested in playing awk golf; I want to show people how it works and make the logic explicit.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T11:49:35.020Z
Answers.Text: Fair enough. I tend to favor brevity in one-liners, but you're right in that it's less clear to the non-expert. Sorry. I would urge you to remove the `\ ` since they are not needed and, IMO, harm readability though. I would also remove the needless parentheses from the awk call for the same reason: `{if(/^@/ || length($10)>500){print}}`
Answers.UserDisplayName: terdon
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T12:20:36.520Z
Answers.Text: I've removed the end of line escapes, but not the brackets. I'm not completely comfortable doing either, because I've had my code do unexpected things for both reasons in the past (although obviously not these exact situations).
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T05:32:15.000Z
Answers.LastActivityDate: 2017-06-13T12:16:07.970Z
Answers.LastEditDate: 2017-06-13T12:16:07.970Z
Answers.UserDisplayName: gringer
Answers.Body: <p>Use <a href="https://github.com/lh3/htsbox" rel="nofollow noreferrer">htsbox</a>:</p>
  
  <pre><code>htsbox samview -p in.bam | less -S
  htsbox samview -pS in.sam | less -S
  </code></pre>
  
  <p>It outputs mapping positions in the <a href="https://github.com/lh3/miniasm/blob/master/PAF.md" rel="nofollow noreferrer">PAF</a> format, which looks something like:</p>
  
  <pre><code>read1  4983  774  4982  +  chr18  80373285  26911072  26915544  3835  4631  60 \
         mm:i:214  io:i:119  in:i:159  do:i:339  dn:i:423
  </code></pre>
  
  <p>This is one long line in terminal. I folded the line for display purposes. For the first 12 fixed fields, see the table in <a href="https://github.com/lh3/miniasm/blob/master/PAF.md" rel="nofollow noreferrer">PAF page</a>. They give mapping positions and identity. The optional fields tell you #mismatches (mm), #insOpens (io), #insertions (in), #delOpens (do) and #deletions (dn). Note that the alignment MUST have an "NM" tag; otherwise the number of matching bases (col. 11) is overestimated and "mm" will always be zero.</p>
  
  <p>You can easily compute summary statistics from PAF if you are familiar with command lines. For example:</p>
  
  <pre class="lang-shell prettyprint-override"><code>htsbox samview -p in.bam | awk '{x+=$10;y+=$11}END{print x/y}'  # identity
  htsbox samview -p in.bam | awk '$2&gt;1000{x+=$10;y+=$11}END{print x/y}'  # identity for &gt;1000bp
  htsbox samview -p in.bam \
    | perl -ane '{$g+=$1 while /[id]n:i:(\d+)/g;$y+=$F[10]}END{print $g/$y,"\n"}' # gap rate
  </code></pre>
  
  <p>If you don't like long one-liners, feel free to write a proper script to collect all statistics in one go.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2017-06-13T14:01:53.477Z
Answers.LastActivityDate: 2017-06-13T14:01:53.477Z
Answers.UserDisplayName: user172818
Body: <p>I have a BAM file, and I have a read ID. What is the simplest way to get mapping statistics of that read in human-readable format?</p>
  
  <p>E.g. I might want: % identity of aligned bases; number of insertions and deletions; number of bases aligned to the reference; etc. </p>
  
  <p>Extra points for one-liners that give detailed summary statistics, and extra-extra points if you can extend the answer to output summaries for all reads longer than some predetermined length.</p>
  
Tags: <sam>
Title: Get the mapping statistics of a single read from a BAM file
ThreadText: # Get the mapping statistics of a single read from a BAM file
  
  (Asked by: username_0 on 2017-06-13)
  
  I have a BAM file, and I have a read ID. What is the simplest way to get mapping statistics of that read in human-readable format?
  
  E.g. I might want: % identity of aligned bases; number of insertions and deletions; number of bases aligned to the reference; etc.
  
  Extra points for one-liners that give detailed summary statistics, and extra-extra points if you can extend the answer to output summaries for all reads longer than some predetermined length.
  
  **Tags:** sam
  
  
  ---
  
  ## 2 answers
  
  ### Answer by username_1 on 2017-06-13. Score: 7
  
  `samtools stats` seems to be able to do most of this, excluding the CIGAR-string parsing stuff (i.e. INDELs):
  
  ```
  $ samtools view -h mapped.bam | grep -e '^@' -e 'readName' |
      samtools stats | grep '^SN' | cut -f 2-
  
  raw total sequences:    2
  filtered sequences: 0
  sequences:  2
  is sorted:  1
  1st fragments:  2
  last fragments: 0
  reads mapped:   2
  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
  reads unmapped: 0
  reads properly paired:  0   # proper-pair bit set
  reads paired:   0   # paired-end technology bit set
  reads duplicated:   0   # PCR or optical duplicate bit set
  reads MQ0:  0   # mapped and MQ=0
  reads QC failed:    0
  non-primary alignments: 0
  total length:   13632   # ignores clipping
  bases mapped:   13632   # ignores clipping
  bases mapped (cigar):   13502   # more accurate
  bases trimmed:  0
  bases duplicated:   0
  mismatches: 3670    # from NM fields
  error rate: 2.718116e-01    # mismatches / bases mapped (cigar)
  average length: 6816
  maximum length: 6816
  average quality:    8.9
  insert size average:    0.0
  insert size standard deviation: 0.0
  inward oriented pairs:  0
  outward oriented pairs: 0
  pairs with other orientation:   0
  pairs on different chromosomes: 0
  
  ```
  
  For a length filter, an awk length check can be slotted in instead of the grep:
  
  ```
  $ samtools view -h mapped.bam |
      awk -F'\t' '{if((/^@/) || (length($10)>500)){print $0}}' |
      samtools stats | grep '^SN' | cut -f 2-
  
  raw total sequences:    131
  filtered sequences: 0
  sequences:  131
  is sorted:  1
  1st fragments:  131
  last fragments: 0
  reads mapped:   131
  reads mapped and paired:    0   # paired-end technology bit set + both mates mapped
  reads unmapped: 0
  reads properly paired:  0   # proper-pair bit set
  reads paired:   0   # paired-end technology bit set
  reads duplicated:   0   # PCR or optical duplicate bit set
  reads MQ0:  0   # mapped and MQ=0
  reads QC failed:    0
  non-primary alignments: 0
  total length:   569425  # ignores clipping
  bases mapped:   569425  # ignores clipping
  bases mapped (cigar):   453884  # more accurate
  bases trimmed:  0
  bases duplicated:   0
  mismatches: 113755  # from NM fields
  error rate: 2.506257e-01    # mismatches / bases mapped (cigar)
  average length: 4346
  maximum length: 16909
  average quality:    9.1
  insert size average:    0.0
  insert size standard deviation: 0.0
  inward oriented pairs:  0
  outward oriented pairs: 0
  pairs with other orientation:   0
  pairs on different chromosomes: 0
  
  ```
  
  > **Comments:**
  > **username_0**: all extra points awarded. This is (hopefully) super useful for folks working on long reads, and wanting to summarise the longer bits of their data...  
  > **username_1**: @username_2 please don't edit this to make it less readable than it already is. I'm not interested in playing awk golf; I want to show people how it works and make the logic explicit.  
  > **username_2**: Fair enough. I tend to favor brevity in one-liners, but you're right in that it's less clear to the non-expert. Sorry. I would urge you to remove the `\ ` since they are not needed and, IMO, harm readability though. I would also remove the needless parentheses from the awk call for the same reason: `{if(/^@/ || length($10)>500){print}}`  
  > **username_1**: I've removed the end of line escapes, but not the brackets. I'm not completely comfortable doing either, because I've had my code do unexpected things for both reasons in the past (although obviously not these exact situations).  
  
  ---
  
  ### Answer by username_3 on 2017-06-13. Score: 2
  
  Use [htsbox](https://github.com/lh3/htsbox):
  
  ```
  htsbox samview -p in.bam | less -S
  htsbox samview -pS in.sam | less -S
  
  ```
  
  It outputs mapping positions in the [PAF](https://github.com/lh3/miniasm/blob/master/PAF.md) format, which looks something like:
  
  ```
  read1  4983  774  4982  +  chr18  80373285  26911072  26915544  3835  4631  60 \
         mm:i:214  io:i:119  in:i:159  do:i:339  dn:i:423
  
  ```
  
  This is one long line in terminal. I folded the line for display purposes. For the first 12 fixed fields, see the table in [PAF page](https://github.com/lh3/miniasm/blob/master/PAF.md). They give mapping positions and identity. The optional fields tell you #mismatches (mm), #insOpens (io), #insertions (in), #delOpens (do) and #deletions (dn). Note that the alignment MUST have an "NM" tag; otherwise the number of matching bases (col. 11) is overestimated and "mm" will always be zero.
  
  You can easily compute summary statistics from PAF if you are familiar with command lines. For example:
  
  ```shell
  htsbox samview -p in.bam | awk '{x+=$10;y+=$11}END{print x/y}'  # identity
  htsbox samview -p in.bam | awk '$2>1000{x+=$10;y+=$11}END{print x/y}'  # identity for >1000bp
  htsbox samview -p in.bam \
    | perl -ane '{$g+=$1 while /[id]n:i:(\d+)/g;$y+=$F[10]}END{print $g/$y,"\n"}' # gap rate
  
  ```
  
  If you don't like long one-liners, feel free to write a proper script to collect all statistics in one go.
  
  
  ---
  
  