5
$\begingroup$

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.

$\endgroup$
5
$\begingroup$

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
$\endgroup$
  • $\begingroup$ 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... $\endgroup$ – roblanf Jun 13 '17 at 7:31
  • $\begingroup$ @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. $\endgroup$ – gringer Jun 13 '17 at 10:10
  • $\begingroup$ 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}} $\endgroup$ – terdon Jun 13 '17 at 11:49
  • 1
    $\begingroup$ 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). $\endgroup$ – gringer Jun 13 '17 at 12:20
2
$\begingroup$

Use htsbox:

htsbox samview -p in.bam | less -S
htsbox samview -pS in.sam | less -S

It outputs mapping positions in the PAF 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. 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:

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.