Simulating A Read With A Deletion Of 50 Bp Into It.

Question

[Gabrielw]

Im trying to test some indel calling programs with large deletions in sequencing reads.

I ve got this sequence in Chromosome 20 of hg19 :

> CTAGCAAGGG GGCTGTATGG CTTGAGGCCA TAGTCCAGGA CATCATCGGG  50140599
GTATGCGGGT CCGGAGGGTG GGCTGGCGAC CTTATGTGCA TTCGGCTCTT  50140649
ctggaaagag aagggggaag ggggttcttt ttaagcctca aaaacagaac  50140699
TCCCGGTGCA GAGCAATCAC AGGCTGGATT TCGTCTCACG TCGTTTATTT  50140749
TTTTCAAATT CCCACCATGC CAAACCCCAA GCTAGAAGCT CCGTCCCTCA  50140799
CCAGAAGAAA ACTGAGGCCA CTTTGTCCCA GGTCCCTGGA GAAGACAGGT  50140849

thank I created a fastq file with a simulated deletion between positions 50140699 and 50140748:

>@teste.01 ER85B5301EGS5O length=250
CTAGCAAGGGGGCTGTATGGCTTGAGGCCATAGTCCAGGACATCATCGGGGTATGCGGGTCCGGAGGGTGGGCTGGCGACCTTATGTGCATTCGGCTCTTTCCCGGTGCAGAGCAATCACAGGCTGGATTTCGTCTCACGTCGTTTATTTTTTTCAAATTCCCACCATGCCAAACCCCAAGCTAGAAGCTCCGTCCCTCACCAGAAGAAAACTGAGGCCACTTTGTCCCAGGTCCCTGGAGAAGACAGGT
+teste.01 ER85B5301EGS5O length=250
@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@

So I mapped my read using bwasw from BWA, and thats what is mapped in my sam file:

>teste.01    0    chr20    50140700    204    100S150M    *    0    0    CTAGCAAGGGGGCTGTATGGCTTGAGGCCATAGTCCAGGACATCATCGGGGTATGCGGGTCCGGAGGGTGGGCTGGCGACCTTATGTGCATTCGGCTCTTTCCCGGTGCAGAGCAATCACAGGCTGGATTTCGTCTCACGTCGTTTATTTTTTTCAAATTCCCACCATGCCAAACCCCAAGCTAGAAGCTCCGTCCCTCACCAGAAGAAAACTGAGGCCACTTTGTCCCAGGTCCCTGGAGAAGACAGGT    @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@    AS:i:150    XS:i:0    XF:i:3    XE:i:4    XN:i:0
teste.01    0    chr20    50140550    177    100M150S    *    0    0    CTAGCAAGGGGGCTGTATGGCTTGAGGCCATAGTCCAGGACATCATCGGGGTATGCGGGTCCGGAGGGTGGGCTGGCGACCTTATGTGCATTCGGCTCTTTCCCGGTGCAGAGCAATCACAGGCTGGATTTCGTCTCACGTCGTTTATTTTTTTCAAATTCCCACCATGCCAAACCCCAAGCTAGAAGCTCCGTCCCTCACCAGAAGAAAACTGAGGCCACTTTGTCCCAGGTCCCTGGAGAAGACAGGT    @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@    AS:i:100    XS:i:0    XF:i:3    XE:i:3    XN:i:0

So I usde samtools sort, mpileup, and then I used pileup2indel from VARSCAN, but it dindt find any indels, does some one knows what can I do ?

Answer 1

[Michael]

I think there are some problems with your approach. Afaik the bwa - samtools pipeline was not intended to call long but short indels of 1 or maybe 2-3 bases, maybe it can be parameterized differently.

The first problem is that the alignment doesn't work as you expected it to do. The CIGAR string 100S150M tells that there were 150 perfect Matches, but the first 100 bases were skipped. This means that an alignment with the first 100 bases skipped has a better score, than aligning 50 matches, 50 gaps, and then 150 matches. This is likely due to the gap (opening and extension) costs during aligning, such that the score of an alignment with 50 deletions is worse than one having 100 bases skipped.

As there is no alignment to the gapped region, there is no coverage and thus no evidence about a deletion on which indel calling could be based. In order to work, you expect a CIGAR string like 50M50D100M or 50M50I100M (depending on insert or delete occured in the reference).

I think it will be hard to tune bwa (decrease gap extension costs) to call the alignment the way you intend, you need for an alignment tool that allows for larger inserts, you might have more luck with e.g. blat.

The second flaw in that workflow is that there are only two reads (did you really only make two, or was this just an example?). A coverage of 2 seems to me rarely sufficient to call any variant.

You should use a read simulator instead on the reference sequence with the deletion, creating a substantial amount of reads which are then aligned against the reference without deletion.

[Gabrielw]

thanks for the tips ! Can u recommend me a read simulator? thanks again

[Michael]

Here: How To Produce Simulated 'Synthetic' Sequences

either MetaSim (gui) or wgsim (comes with samtools)