Generate pileup for indels only using SAMtools mpileup
2
1
Entering edit mode
6.6 years ago
ak221193 ▴ 10

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.

SNP indel SAMtools mpileup • 3.1k views
ADD COMMENT
0
Entering edit mode
6.6 years ago

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.

ADD COMMENT
0
Entering edit mode
6.6 years ago

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.

ADD COMMENT

Login before adding your answer.

Traffic: 1895 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6