Question: Calculate numbers of indels <=2 nt within sam
0
gravatar for Rox
3.2 years ago by
Rox1.2k
France / Toulouse / GeT-Plage
Rox1.2k wrote:

Hello everyone,

I am really bad at manipulating sam files... I was looking for a simple way to count, within a sam, the number of indels <=2 nucleotides. I also want it to be position specific, which means for me that, for a single position, if I have a read coverage of 100 and in each read there is a deletion, I only want to count one.

Is this a simple thing to do using samtools or any other scripting method ?

Thanks for your help !

Roxane

script • 995 views
ADD COMMENTlink modified 3.2 years ago by Pierre Lindenbaum129k • written 3.2 years ago by Rox1.2k

I was looking for a simple way to count, within a sam, the number of indels <=2 nucleotides.

This wouldn't be too hard since this information is readily available in the CIGAR strings.

I also want it to be position specific, which means for me that, for a single position, if I have a read coverage of 100 and in each read there is a deletion, I only want to count one.

That makes it a bit more complicated.

I would probably roll my own solution (using python-pysam), but it's perfectly possible that what you are looking for already exists...

ADD REPLYlink written 3.2 years ago by WouterDeCoster44k
0
gravatar for Rahul Sharma
3.2 years ago by
Rahul Sharma600
Germany
Rahul Sharma600 wrote:

Few months ago I did some SNP analyses using mapped sam files. I compared two samples say prep2 and prep3, using samtools and bcftools. Following are the commands I used to investigate all variations including INDELS with a coverage cut-off of 100 and mapping quality cutoff 20. You could grep INDELS from the 'Qualt20DP100_flt.vcf.passed' file and check for position.

 nohup samtools mpileup -t DP -uv -Q 15 -f /media/Storage/Analysis/mm9/genome.fa /media/Storage/Data/bam/prep2.bam /media/Storage/Data/bam/prep3.bam -o prep2-3_mpileup.vcf &       
 bcftools call -mv prep2-3_mpileup.vcf > Only_variations.vcf &
 bcftools filter -s LowQual -e '%QUAL<20 || DP<100' Only_variations.vcf > Qualt20DP100_flt.vcf
 grep -P "\tPASS\t" Qualt20DP100_flt.vcf > Qualt20DP100_flt.vcf.passed

I hope it helps.

Cheers, Rahul

ADD COMMENTlink written 3.2 years ago by Rahul Sharma600

Yes, but I think OP is interested in indels on the read level, also if those aren't reported by a variant caller.

ADD REPLYlink written 3.2 years ago by WouterDeCoster44k
0
gravatar for Brian Bushnell
3.2 years ago by
Walnut Creek, USA
Brian Bushnell17k wrote:

The BBMap package can produce a histogram of such events, like this:

reformat.sh in=mapped.sam indelhist=indelhist.txt

That will show how many insertions and deletions there are of any specific length. You can also produce this directly from BBMap during mapping.

ADD COMMENTlink written 3.2 years ago by Brian Bushnell17k
0
gravatar for Pierre Lindenbaum
3.2 years ago by
France/Nantes/Institut du Thorax - INSERM UMR1087
Pierre Lindenbaum129k wrote:

using bioalcidae: http://lindenb.github.io/jvarkit/BioAlcidae.html with the following script:

while(iter.hasNext()) { 
    var rec=iter.next();
    if(rec.getReadUnmappedFlag()) continue;
    var ref= rec.getAlignmentStart();
    var cigarElements =  rec.getCigar().getCigarElements();
    for(var i=0;i< cigarElements.size();i++)
        {
        var ce = cigarElements.get(i);
        var op = ce.getOperator();
        if( op.isIndelOrSkippedRegion()  && ce.getLength()<2)
            {
            out.println(rec.getContig()+"\t"+ref+"\t"+op.name()+"\t"+ce.getLength());
            }
        if( op.consumesReferenceBases())
            {
            ref += ce.getLength();
            }
        }
    }

and pipe in sort | uniq:

 java -jar dist/bioalcidae.jar -f script.js in.bam | sort | uniq
ADD COMMENTlink written 3.2 years ago by Pierre Lindenbaum129k
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.3.0
Traffic: 678 users visited in the last hour