How to quickly determine mutations in a read of a sam file?

Question

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?

Answer 1

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 2

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.

Answer 3

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.