<p><code>samtools mpileup</code> can do this quickly:</p>

<pre><code>samtools mpileup -f reference.fasta -uv input.sam &gt; variants.vcf
</code></pre>

<p>This will produce a <a href="http://vcftools.sourceforge.net/VCF-poster.pdf" rel="nofollow noreferrer">VCF-formatted</a> file containing information about what variants have been seen in the SAM file, aggregated for all the mapped reads.</p>

CC BY-SA 3.0
2017-05-18T16:23:09.657Z
2017-05-18T16:23:09.657Z
gringer
<p>I wrote a program, <a href="http://asciigenome.readthedocs.io/en/latest/description.html" rel="noreferrer">ASCIIGenome</a>, that I find handy in cases where you want to have a quick look at genomic data. It's a genome browser for the command line. </p>

<p>To view only reads containing mismatches you can use the internal function <code>awk</code>. To filter for reads where the NM tag (number of mismatches) is >0:</p>

<pre><code>ASCIIGenome -fa genome.fa aln.bam
...

[h] for help: awk 'getSamTag("NM") &gt; 0'
</code></pre>

<p>The view on terminal screen may look something like this:<a href="https://i.sstatic.net/qXKkp.png" rel="noreferrer"><img src="https://i.sstatic.net/qXKkp.png" alt="enter image description here"></a></p>

<p>Similarly, to get only reads containing indels you can use <code>awk '$6 ~ "D|I"'</code></p>

<p>Hope this helps and feel free to report bugs &amp; issues.</p>

CC BY-SA 3.0
2017-05-21T19:01:01.607Z
That's really cool. I like it.
ChrisD
CC BY-SA 3.0
2017-05-21T17:47:47.147Z
2017-05-21T17:47:47.147Z
dariober
<p>For qualitative analysis, you're probably better off using something less granular like <a href="http://software.broadinstitute.org/software/igv/" rel="noreferrer">IGV</a> or <a href="http://bioviz.org/igb/overview.html" rel="noreferrer">IGB</a>. However, if you really want to look at a couple of reads:</p>

<p>If you're willing to ignore sequencing errors, you can inspect the CIGAR string or the MD tag, both of which give information on the alignment of a single read.</p>

<p>The CIGAR string gives details on insertions, deletions, clippings, matches and mismatches. From <a href="http://genome.sph.umich.edu/wiki/SAM" rel="noreferrer">Genome Analysis Wiki</a>,</p>

<blockquote>
  <p>The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference. For example:</p>
</blockquote>

<pre><code>RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M
</code></pre>

<p>Most common usage of the CIGAR string uses M (match/mismatch), I (insertion), D (deletion), S (soft clipping), and H (hard clipping). Note that = (match) and X (mismatch) are available as alternatives to the less informative M, but they are less widely used.</p>

<p>The MD tag gives specific details on mismatches and deletions. From the <a href="https://samtools.github.io/hts-specs/SAMtags.pdf" rel="noreferrer">SAMtools tags specification</a>,</p>

<blockquote>
  <p>The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string ‘10A5^AC6’ means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which is different from the aligned read base; the next 5 reference bases are matches followed by a 2bp deletion from the reference; the deleted sequence is AC; the last 6 bases are matches. The MD field ought to match the CIGAR string.</p>
</blockquote>

<p>Note that neither of these will give you any idea of structural variants in short reads, and neither will be particularly readable (or helpful, due to the higher error rate) in long reads.</p>

CC BY-SA 3.0
2017-05-18T03:24:19.127Z
2017-05-18T03:33:12.387Z
2017-05-18T03:33:12.387Z
Scott Gigante
<p>After DNA sequencing, I generated a sam file through alignment of a fastq file. 
Before using well known variant calling programs (eg. Annovar etc.), I want to pick some reads and know what kinds of mutations are there.</p>

<p>Is there any quick way of doing it?</p>

CC BY-SA 3.0
2017-05-18T03:22:39.143Z
You can visualize mismatches with something like [IGB](https://wiki.transvar.org/display/igbman/Visualizing+read+alignments#Visualizingreadalignments-Loadsequencedatatoviewmismatches)
ChrisD
<variant-calling>
How to quickly determine mutations in a read of a sam file?
# How to quickly determine mutations in a read of a sam file?

(Asked by: username_0 on 2017-05-18)

After DNA sequencing, I generated a sam file through alignment of a fastq file.
Before using well known variant calling programs (eg. Annovar etc.), I want to pick some reads and know what kinds of mutations are there.

Is there any quick way of doing it?

**Tags:** variant-calling
> **Comments:**
> **username_1**: You can visualize mismatches with something like [IGB](https://wiki.transvar.org/display/igbman/Visualizing+read+alignments#Visualizingreadalignments-Loadsequencedatatoviewmismatches)  

---

## 3 answers

### Answer by username_4 on 2017-05-18. Score: 8 | Accepted answer

For qualitative analysis, you're probably better off using something less granular like IGV or IGB. However, if you really want to look at a couple of reads:
If you're willing to ignore sequencing errors, you can inspect the CIGAR string or the MD tag, both of which give information on the alignment of a single read.
The CIGAR string gives details on insertions, deletions, clippings, matches and mismatches. From Genome Analysis Wiki,

The sequence being aligned to a reference may have additional bases that are not in the reference or may be missing bases that are in the reference. The CIGAR string is a sequence of of base lengths and the associated operation. They are used to indicate things like which bases align (either a match/mismatch) with the reference, are deleted from the reference, and are insertions that are not in the reference. For example:

RefPos:     1  2  3  4  5  6  7     8  9 10 11 12 13 14 15 16 17 18 19
Reference:  C  C  A  T  A  C  T     G  A  A  C  T  G  A  C  T  A  A  C
Read:                   A  C  T  A  G  A  A     T  G  G  C  T
With the alignment above, you get:
POS: 5
CIGAR: 3M1I3M1D5M

Most common usage of the CIGAR string uses M (match/mismatch), I (insertion), D (deletion), S (soft clipping), and H (hard clipping). Note that = (match) and X (mismatch) are available as alternatives to the less informative M, but they are less widely used.
The MD tag gives specific details on mismatches and deletions. From the SAMtools tags specification,

The MD field aims to achieve SNP/indel calling without looking at the reference. For example, a string ‘10A5^AC6’ means from the leftmost reference base in the alignment, there are 10 matches followed by an A on the reference which is different from the aligned read base; the next 5 reference bases are matches followed by a 2bp deletion from the reference; the deleted sequence is AC; the last 6 bases are matches. The MD field ought to match the CIGAR string.

Note that neither of these will give you any idea of structural variants in short reads, and neither will be particularly readable (or helpful, due to the higher error rate) in long reads.


---

### Answer by username_3 on 2017-05-21. Score: 7

I wrote a program, ASCIIGenome, that I find handy in cases where you want to have a quick look at genomic data. It's a genome browser for the command line. 
To view only reads containing mismatches you can use the internal function awk. To filter for reads where the NM tag (number of mismatches) is >0:
ASCIIGenome -fa genome.fa aln.bam
...

[h] for help: awk 'getSamTag("NM") > 0'

The view on terminal screen may look something like this:
Similarly, to get only reads containing indels you can use awk '$6 ~ "D|I"'
Hope this helps and feel free to report bugs & issues.

> **Comments:**
> **username_1**: That's really cool. I like it.  

---

### Answer by username_2 on 2017-05-18. Score: 4

samtools mpileup can do this quickly:
samtools mpileup -f reference.fasta -uv input.sam > variants.vcf

This will produce a VCF-formatted file containing information about what variants have been seen in the SAM file, aggregated for all the mapped reads.


---

