This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Identifying Complex Mutations With Samtools

Is it possible to use samtools to identify complex mutations caused by e.g. UV mutagenesis?

samtools is good at identifying single, isolated mutations (e.g. single transitions and transversions in sequencing reads), but is it capable of identifying neighboring (i.e. consecutive) mutations?

The case I am interested in is UV mutagenesis, which can cause consecutive transitions, going from CT in a wild-type genome to TC at the same positions in in a mutant genome. I have ~50X Illumina coverage of a mutant genome, and can "see" this complex mutation (using samtools tview and zooming to the locus), but can't find arguments that allow me to identify each consecutive mutation via samtools mpileup and bcftools; I'd like to capture these positions in VCF files. I've fiddled with gap opening and extension parameters, but I don't think this is technically a "gapped" problem.

Any insight here?

samtools mutation
• 3,253 views • link
updated by Biostar Biostar • written by Jayhesselberth • Jayhesselberth

2 answers

Quick starting step for you to improve up on:

samtools mpileup -uf Input.bam | bcftools view -cg - | awk -F"\t" '/^$/ { next }; {if ($5 =="."){ ini=0; next }; }{this_coord = $2 ; this_four = $4 ; this_five = $5 }  ini !=1 {diff = this_coord - prev_coord; if (diff==1){ printf "%s\t%s\t%s\t%s\t%s\t%s\t%s\n", this_coord, prev_coord, "diff="diff, prev_four, this_four, prev_five, this_five}} {prev_coord = this_coord ; prev_four = this_four ; prev_five = this_five}'

Potential output:

This_coord Prev_coord difference ref1 var1 ref2 var2
10000   9999    diff=1  C       T       G       A
10353   10352   diff=1  T      A   A        C
• 0 views • link
written by Rm Rm

The problem is that samtools mpileup doesn't call the variants in the first place. All this code does is report them separately if they are called.

• 0 views • link
written by Jayhesselberth • Jayhesselberth

@Jayhesselberth; "if they are called"...you are correct

• 0 views • link
written by Rm Rm

Log in to answer this question.

More posts like this