Identifying Complex Mutations With Samtools

Question

[Jayhesselberth]

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?

[Casey Bergman]

Good question, also relevant to naturally occuring variants: http://www.ncbi.nlm.nih.gov/pubmed/10678838

Answer 1

[Rm]

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

Answer 2

[Jayhesselberth]

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.

[Rm]

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