I am trying to generate Single nucleotide mutations to a sequence, and trying to get SAM files and its related CIGAR sequences.
My test genome is: (ref.fa)
>genome CGGACACACAAAAAGAAAGAAGAATTTTTAGGATCTTTTGTGTGCGAATA GATGGGCCCATCGATTGCCATCCCGGGGAAAATCATCATCATCATCATCA
Then I created a fastq sequence (seq_1.fq) which I inserted a Nucleotide A, to the first line and ran the BWA alignment, then the SAM file
bwa mem ref.fa seq_1.fq > seq_1.sam samtools sort -o seq_1.bam --output-fmt BAM seq_1.sam samtools view seq_1.bam | sam2pairwise
which shows the insertion by giving the CIGAR seq as 29M1I21M:
genome 0 genome 1 60 29M1I21M * 0 0 CGGACACACAAAAAGAAAGAAGAATTTTTAAGGATCTTTTGTGTGCGAATA ||||||||||||||||||||||||||||| ||||||||||||||||||||| CGGACACACAAAAAGAAAGAAGAATTTTT-AGGATCTTTTGTGTGCGAATA
Question: But when I tried to replace a T using G from 'TTTTT' to 'TTGTT' (seq_2.fq) and ran the BWA align and samtools it does not give me the desired CIGAR sequence !!! BUT it just gives it as 50M
bwa mem ref.fa seq_2.fq > seq_2.sam samtools sort -o seq_2.bam--output-fmt BAM seq_2.sam samtools view seq_2.bam | sam2pairwise genome 0 genome 1 60 50M * 0 0 CGGACACACAAAAAGAAAGAAGAATTGTTAGGATCTTTTGTGTGCGAATA |||||||||||||||||||||||||||||||||||||||||||||||||| CGGACACACAAAAAGAAAGAAGAATTGTTAGGATCTTTTGTGTGCGAATA
as it does not recognize the mutation which was introduced to the sequence.
So am missing something here ? How can I generate a CIGAR sequence which states, that there is a mismatch or SNP in the sequence ?
Thank you very much.