Generate pileup for indels only using SAMtools mpileup

Question

[ak221193]

Is there a way to generate a pileup for indels only using the samtools mpileup command? The -I flag can be used to skip indel calling but there is no option to call indels only.

Thanks.

Answer 1

[Devon Ryan]

There's no built-in way to do that, no. You'd could code something to do that in pysam or using the htslib C interface, though this is probably not simple. The general idea would be to either first filter out reads lacking indels or to parse the pileup structure and modify it as you go to remove reads without an indel at a given position.

Answer 2

[dariober]

If you want to see how many reads support an indel at each position, you could count the number of '+' and '-' in the base string (5th column) of the pileup. Not sure if this is what you mean by "pileup for indels" though. This is a crude python implementation of the idea:

samtools mpileup -f myref.fa aln.bam | python -c "
import sys
for line in sys.stdin:
    line= line.strip().split('\t')
    if len(line) < 5:
        print('\t'.join(line))
    else:
        rbase= line[4].replace('^+', '').replace('^-', '') ## Skip +- after ^
        n_ins= rbase.count('+') ##  N. reads supporting insertions
        n_del= rbase.count('-') ##  N. reads supporting deletion
        print('\t'.join(line + [str(n_ins)] + [str(n_del)]))
"

Output is the same as mpileup additional columns for no. of reads supporting insertion and no. reads supporting deletion at each position.