I have two sequenes and I want to get the CIGAR string of a local alignment. I've tried with bwa and bowtie
The first one give me this, all *. When I change a nucleotide in the query, all it does is change the sequence reported in the output (8th column), but nothing else!
$ ./bwa mem ../seq1.fa ../seq2.fa
[M::bwa_idx_load_from_disk] read 0 ALT contigs
@SQ SN:subject LN:1141
@PG ID:bwa PN:bwa VN:0.7.17-r1198-dirty CL:./bwa mem ../seq1.fa ../seq2.fa
[M::process] read 1 sequences (13 bp)...
[M::mem_process_seqs] Processed 1 reads in 0.000 CPU sec, 0.000 real sec
subject 4 * 0 0 * * 0 0 AGCAAAGCAAGTT * AS:i:0 XS:i:0
[main] Version: 0.7.17-r1198-dirty
[main] CMD: ./bwa mem ../seq1.fa ../seq2.fa
[main] Real time: 0.002 sec; CPU: 0.006 sec
With bowtie, if I change a base in the query it gets me the same CIGAR (12M). Doesn't seems to be logic to me.
$ bowtie seq1.fa.fai -f seq2.fa --sam
@HD VN:1.0 SO:unsorted
@SQ SN:subject LN:1141
@PG ID:Bowtie VN:1.2.2 CL:"/Users/juan/Documents/manu/dev/bioinfoch/2/E_TEs/bowtie-1.2.2-macos-x86_64/bowtie-align-s --wrapper basic-0 seq1.fa.fai -f seq2.fa --sam"
subject 0 subject 331 255 12M * 0 0 AGCAAAGCAAGT IIIIIIIIIIII XA:i:0 MD:Z:12 NM:i:0 XM:i:2
# reads processed: 1
# reads with at least one reported alignment: 1 (100.00%)
# reads that failed to align: 0 (0.00%)
Reported 1 alignments
➜ E_TEs
$ bowtie seq1.fa.fai -f seq2.fa --sam
@HD VN:1.0 SO:unsorted
@SQ SN:subject LN:1141
@PG ID:Bowtie VN:1.2.2 CL:"/Users/juan/Documents/manu/dev/bioinfoch/2/E_TEs/bowtie-1.2.2-macos-x86_64/bowtie-align-s --wrapper basic-0 seq1.fa.fai -f seq2.fa --sam"
subject 0 subject 331 255 12M * 0 0 AGCAAAGTAAGT IIIIIIIIIIII XA:i:1 MD:Z:7C4 NM:i:1 XM:i:2
# reads processed: 1
# reads with at least one reported alignment: 1 (100.00%)
# reads that failed to align: 0 (0.00%)
Reported 1 alignments
•
written
by
juan.crescente •
M means that the nucleotide maps to the REF. Maps here means that there are same number of nucleotide in the REF (i.e. No INDELs), but they can be either match or mismatch
Bowtie: Here 12M means the "map" is over 12 nucleotide sequence (i.e. without INDELs), but it can be either match or mismatch.
If you need to see the differences, see the MD tag MD:Z:7C4 (7 Match, one C mismatch and then 4 match)
BWA: BWA is not able to find any match. Note that bwa is looking for full length match of 13bp (vs bowtie of 12bp).
is there a way to get a string represeting the alignment straightforward like MATCH-MATCH-MISMATCH
Parse MD tag. Matches will appear as numbers, mismatches as [ATCG]
Log in to answer this question.
first read has sam flag = 4 your read is unmapped, no cigar string.
I see, with bwa mem. What about the CIGAR in bowtie?
The sequence is to short to work with
bwa mem.Further more a
Min a CIGAR String meansalignment match. This can be a match or a mismatch to the reference. TheMDtag gives you the information whether it is a match or a mismatch to the reference.There are also tools that can convert "normal" CIGAR strings to extended ones. See How to convert from a cigar string to extended cigar string?